Solution (source code)

= Solution

Define the <reduced gravity>
$$
g'=g\frac{\rho_1-\rho_2}{\rho_2},
$$
and use the Boussinesq limit $\rho_1\simeq\rho_2$ except in buoyancy. The balances become
$$
g'_t+ug'_x=-g'\frac{w_e}{h},
$$
$$
h_t+uh_x+hu_x=w_e-w_d,
$$
$$
u_t+uu_x+g'h_x+\frac h2g'_x=-u\frac{w_e}{h}.
$$
Thus, for $\mathbf q=(g',h,u)^T$,
$$
\boxed{
\mathbf q_t+
\begin{pmatrix}
u&0&0\\
0&u&h\\
h/2&g'&u
\end{pmatrix}\mathbf q_x
=
\begin{pmatrix}
-g'w_e/h\\
w_e-w_d\\
-uw_e/h
\end{pmatrix}},
$$
which is the required <entraining shallow-water layer> system.