Solution (source code)

= Solution

Use the <Dupuit approximation>: the saturated region is shallow enough that <hydrostatic pressure> is $p=\rho g(h-z)$ and flow is predominantly horizontal. <Darcy's law> then gives the horizontal <Darcy velocity> $u_x=-u_b(1+\beta z)h_x$, where $u_b=k_0\rho g/\mu$. Integrating through the saturated depth gives the <volume flux per unit width>
$$
q=\int_0^h u_x\,dz=-u_b\left(h+\frac\beta2h^2\right)h_x.
$$
The stored water volume per unit area is $\phi h$, so <mass conservation> yields the <unconfined aquifer with depth-dependent permeability> equation
$$
\boxed{\phi h_t=u_b\partial_x\left[\left(h+\frac\beta2h^2\right)h_x\right]+R,\qquad h(0,t)=0,\qquad q(L,t)=0.}
$$
Define the positive river discharge by $Q(t)=-q(0,t)$. A useful check on all subsequent results is the integrated <mass conservation> law
$$
\phi\frac d{dt}\int_0^Lh\,dx=RL-Q.
$$