Solution of eikonal equation is locally the distance from a hypersurface, up to a constant

Consider the Eikonal equation (with right handside 1)

$$\sum_{i=1}^{n}(\frac{\partial u}{\partial x_i})^2=1$$

I want to see why any solution to this is locally the sum of a distance function from a hypersurface plus a constant. I think I came up with a way but I wonder if there is a simpler way to argue this. Here is what I did:

Assume $u(x)$ is a solution which is smooth in some domain. Take a level set $N= u^{-1}(c)$ (restricted to this domain). Then we want to solve a first order non-linear partial differetial equation with the initial condition $u|_{N}=c$. The characteristic equation for this pde is

$$\dot{x} = 2p$$

$$\dot{p} = 0$$

$$\dot{u} = 2p^2 = 2$$

We know that the full solution to this equation is a function $u(x)$ whose graph $(x,u(x),\nabla u(x))$ in fiber bundle of the jet of functions over the domain is the set of characteristics passing through $N$ (since $\nabla u$ is not the tanget space of the initial manifold). This means that locally in coordinates $p$ does not change since if $g(w,t)$ are the characteristics curves:

$$p(g(w,t)) = p(w) + \int_{0}^{t} \dot{p}(g(w,s))ds = p(w)$$

for all $t$. This also means that if project a characteristic curve starting at $(x_0,u_0,p_0)$ down to the base manifold it is a straight line of the form

$$x(t) = x_0 + 2tp_0 = x_0 + 2t\nabla u(x_0)$$

This proves that locally $\nabla u$ are aligned in straight lines and is of norm 1 that is level sets of $u$ are parallel translates of each other. Now it is easy to see that $u(x)$ can be locally written as distance from $N$ + c. Infact again by the characteristic equations

$$u(t) = c + 2t$$

where $d(x(t),x_0)=2t$ by the equation above and the fact that characteristic lines projected down are straight lines.

Your proof makes sense, although I feel that talking about fiber bundles is somehow unnecessary. First, I would observe that $u$ is $1$-Lipschitz, since its gradient is $1$. This implies that for any $b<c$, the level surfaces satisfy $$\operatorname{dist}(u^{-1}(b),u^{-1}(c))\ge c-b \tag1$$ On the other hand, starting from a point of $u^{-1}(b)$ and following the trajectory of steepest ascent $\dot x=\nabla u(x)$, we have $\frac{d}{dt} u(x(t))=1$. Therefore, $x(t)\in u^{-1}(c)$ when $t=c-b$. The length of trajectory $x([0,c-b])$ is $c-b$. Contrasting with (1), we see that
• every point of $u^{-1}(b)$ is at distance $c-b$ from $u^{-1}(c)$.
The same applies with $b$ and $c$ interchanged, except we follow the path of steepest descent then.
Therefore, for every $a$, the function $u$ agrees with the signed distance function of $u^{-1}(a)$, plus $a$.