By analogy with nonlinear elastodynamics (or with the $p$-system of fluid dynamics), we introduce the velocity variable $v = u_t$ and the strain variable $\theta = u_x$. Using the equality of mixed derivatives and the second-order PDE in OP, we rewrite the problem as a nonlinear system of conservation laws of the form ${\bf u}_t + {\bf f}({\bf u})_x = {\bf 0}$,
where
$$
{\bf u} = \begin{pmatrix}
\theta\\
v
\end{pmatrix},
\qquad
{\bf f}({\bf u}) = -\! \begin{pmatrix}
v\\
\tfrac{1}{3}(1+\theta)^3
\end{pmatrix} .
$$
As suggested in OP, we assume that $v = F(\theta)$, where $F$ is a smooth function to be determined. Injecting this Ansatz in the PDE system leads to
$$
\begin{pmatrix}
1 & -F'(\theta)\\
F'(\theta) & -(1+\theta)^2
\end{pmatrix}
\begin{pmatrix}
\theta_t \\
\theta_x
\end{pmatrix} =
\begin{pmatrix}
0 \\
0
\end{pmatrix} .
$$
Nontrivial solutions can be obtained provided that the determinant of the matrix above vanishes, i.e. $F'(\theta)^2 = (1+\theta)^2$. Therefore, the (smooth) function $F$ satisfies $F'(\theta) = \pm(1+\theta)$, which solutions are $F(\theta) = c \pm (\theta+\frac12 \theta^2)$ for some constant $c$. Note that this constant equals $v \mp \int_0^\theta (1+\vartheta)\,\text d \vartheta$ which is a Riemann invariant of the system.
With this expression of $v=F(\theta)$, the system of conservation laws amounts to the single PDE
$$
\theta_t \mp (1 + \theta)\, \theta_x = 0 \, .
$$
Using the boundary condition $\theta(x_0,0) = h'(x_0)$, one deduces that the characteristic curves are the lines $x = x_0 \mp (1+\theta) t$ along which $\theta = h'(x_0)$ is constant, which can be written in implicit form as
$$
\theta = h'\big(x \pm (1+\theta) t\big)\, .
$$
Since $\theta$ is constant along the characteristic lines, the velocity $v = F(\theta)$ is constant along those lines too. Therefore, the variable $u$ is increasing along those lines with a constant rate. More precisely, we have
$$
\frac{\text d}{\text d t} u\big(x_0 \mp (1+\theta) t, t\big) = \mp (1+\theta) \theta + F(\theta)
= \mp\tfrac12 \theta^2 + c \, ,
$$
where $\theta = h'(x_0)$ is constant.
Finally, one obtains
$$
u = h\big(x \pm (1+\theta) t\big) \mp\tfrac12 \theta^2 t + c t \, ,
$$
where $\theta$ satisfies the implicit equation $\theta = h'(x \pm (1+\theta) t)$.
Note that this expression is different from the one in OP (which obviously doesn't match the boundary condition), and that it is only valid until the breaking time $t_b = \pm 1/\inf h''$.
The concept of simple wave solution refers to an Ansatz of the form ${\bf u}(x,t) = {\bf v}(\xi)$, where $\xi$ is a smooth function of $(x,t)$ and ${\bf v}(\xi)$ follows an integral curve. Along the integral curve, the Riemann invariant $v \mp \int_0^\theta (1+\vartheta)\,\text d \vartheta$ is constant, and ${\bf v}'(\xi)$ is the eigenvector of the Jacobian matrix ${\bf f'}({\bf v}(\xi))$ corresponding to the eigenvalue $\mp (1+\theta(\xi))$. The quasi-linear PDE thus yields the scalar conservation law $\xi_t \mp (1+\theta(\xi))\, \xi_x = {0}$ for $\xi$, which is constant along its characteristic lines deduced from the initial condition $\theta(\xi(x,0)) = h'(x)$. Note that $v$ and $\theta$ are constant along these lines too, which is expressed by the PDE on $\theta$ above.
Note: This example was used in the introduction of the following paper: F. John, "Delayed singularity formation in solutions of nonlinear wave equations in higher dimensions", Comm. Pure Appl. Math. 29, 1976, 649-682. doi:10.1002/cpa.3160290608 In: Moser, J. (eds) Fritz John. Collected Papers Volume 1, Birkhäuser, 1985. doi:10.1007/978-1-4612-5406-5_37
