I recently came across this problem in one of my differential equation texts:
Explain why the characteristic curves generate the solution surface of $P(x,y) \frac{\partial z}{\partial x} + Q(x,y) \frac{\partial z}{\partial y} = R(x,y,z)$ assuming $P,Q$ are Lipschitz in $x,y$ and $R$ is Lipschitz in $z$.
Now I know the proof of this fact, but I think this is not the desired answer. Is there perhaps a (geometrical) explanation of why this might be the case?
I think it might be possible to reason as follows: One can show that $(P,Q,R)$ lies in the tangent plain of any point on the solution surface. $(P,Q,R)$ is by definition also tangent to any point on the characteristic curves...
Let's look at the proof of the fact and try to extract the geometric meaning in its steps. We are trying to solve: $$ P(x, y)u_x + Q(x, y) u_y - R(x, y, u) = 0 \label{1}\tag{1}. $$ The characteristic equations are $$ \begin{cases} (a) \ \ (\dot x(s), \dot y (s)) = (P(x, y), Q(x, y)) \\ (b) \ \ \dot z(s) = R(x, y, z). \end{cases} $$ We can write this as just a single equation for the triple $(x(s), y(s), z(s)) = \gamma(s)$: $$ \dot \gamma(s) = (P(x, y), Q(x, y), R(x, y, z)) $$ (note that all arguments on the right hand side depend on the parameter $s$). These are our characteristic curves.
Now we can argue as you say. The point is exactly the tangency you mentioned. The solution surface is the graph of a solution $z = u(x, y)$ in $(x, y, z)$-space. The first step is to recast the PDE \eqref{1} into a geometric interpretation. For this we use that a normal vector to a graph of a function $z = f(x, y)$ is given by $(f_x, f_y, -1)$. So the normal vector to the graph of $z = u(x, y)$ is $$ (u_x, u_y, -1). $$ If we dot product this with $\vec{V} = (P, Q, R)$, as you say, we see that $$ \vec{V}\cdot (u_x, u_y, -1) = Pu_x + Qu_y - R = 0 $$ by the PDE \eqref{1}. So, saying that $u$ solves \eqref{1} is equivalent to saying that $\vec{V}$ is tangent to the graph of $u$. Now let's look at the characteristic curves. The next point is that the characteristic curves are integral curves of $\vec{V}$. This we have actually already shown above: $\gamma$ satisfies $$ \dot \gamma(s) = (P, Q, R) = \vec{V}(\gamma(s)). $$ The last equality just follows by rewriting things in terms of our new $\vec{V}$ notation.
Thus, the curves $\gamma$ must lie entirely within the solution surface, since they are integral curves of the vector field $\vec{V}$ which, we computed above, is tangent to the solution surface. As for the question of why they must cover the whole solution surface: given any $(x_0, y_0, z_0 = u(x_0, y_0))$ on the solution surface, we have the vector field $\vec{V}$ here and can solve the characteristic ODE $\dot\gamma = \vec{V}(\gamma)$ with the initial condition $\gamma(0) = (x_0, y_0, z_0)$.