Nodal solutions for Neumann systems with gradient dependence

We consider the following convective Neumann systems:\begin{equation*}\left(\mathrm{S}\right)\qquad\left\{\begin{array}{ll}-\Delta_{p_1}u_1+\frac{|\nabla u_1|^{p_1}}{u_1+\delta_1}=f_1(x,u_1,u_2,\nabla u_1,\nabla u_2) & \text{in}\;\Omega,\\ -\Delta _{p_2}u_2+\frac{|\nabla u_2|^{p_2}}{u_2+\delta_2}=f_2(x,u_1,u_2,\nabla u_1,\nabla u_2)&\text{in}\;\Omega, \\ |\nabla u_1|^{p_1-2}\frac{\partial u_1}{\partial \eta }=0=|\nabla u_2|^{p_2-2}\frac{\partial u_2}{\partial \eta}&\text{on}\;\partial\Omega,\end{array}\right.\end{equation*}where $\Omega$ is a bounded domain in $\mathbb{R}^{N}$ ($N\geq 2$) with a smooth boundary $\partial\Omega$,$\delta_1,\,\delta_2 >0$ are small parameters, $\eta$ is the outward unit vector normal to $\partial \Omega,$ $f_1,\,f_2:\Omega\times\mathbb{R}^2\times\mathbb{R}^{2N}\rightarrow \mathbb{R}$ are Carath\'{e}odory functions that satisfy certain growth conditions, and $\Delta _{p_i}$ ($1<p_i<N,$ for $i=1,2$) are the $p$-Laplace operators $\Delta _{p_i}u_i=\mathrm{div}(|\nabla u_i|^{p_i-2}\nabla u_i)$,for every $\,u_i\in W^{1,p_i}(\Omega).$ In order to prove the existence of solutions to such systems, we use a sub-supersolution method. We also obtain nodal solutions by constructing appropriate sub-solution and super-solution pairs. To the best of our knowledge, such systems have not been studied yet.