## 4 Numerical Solution to the Viscous Free Interface Model
### Discretization by an Optimality Condition  
The second numerical method for solving conditions (3), (4) and (5) with positive viscosity \(\
u>0\) uses a formulation as an optimal problem set on the whole interval \((0,h)\), valid under the assumption \(\partial_{Z}S(t,Z)\leq 0\),  
\[\left\{\begin{aligned} &\min(\partial_{t}U+S-\partial_{Z}(\
u \partial_{Z}U),U)=0&&\text{for all }Z\in(0,h),\\ &\partial_{Z}U=0&&\text{at }Z=h,\\ & U=0&&\text{at }Z=0,\end{aligned}\right. \tag{26}\]  
where we consider a no-slip boundary condition at the bottom. Note that this condition becomes relevant whenever the static/flowing interface reaches the bottom, a situation encountered in our simulations as shown in Figure 10c. In (26), the static/flowing interface position \(b(t)\) no longer appears explicitly, but has to be deduced from the velocity profile as  
\[b(t)=\inf\left\{Z\in[0,h]\text{ such that }U(t,Z)>0\right\}. \tag{27}\]  
We split the space domain \((0,h)\) into \(n_{Z}\) cells of length \(\Delta Z\) with \(n_{Z}\Delta Z=h\), and denote the center of the cells by \(Z_{j}=(j-1/2)\Delta Z\), for all \(j=1\ldots n_{Z}\). The discrete times \(t^{n}\), for all \(n\geq 0\), are related by \(t^{n+1}=t^{n}+\Delta t^{n}\), where \(\Delta t^{n}\) is the time step (chosen according to the CFL condition (30) below) and \(t^{0}=0\). We discretize the problem (26) using a finite difference scheme by writing for the discrete unknowns \(U_{j}^{n}\simeq U(t^{n},Z_{j})\)  
\[\min\left(\frac{U_{j}^{n+1}-U_{j}^{n}}{\Delta t^{n}}+S(t^{n},Z_{j})-\frac{v_{ j+1/2}^{n}(U_{j+1}^{n}-U_{j}^{n})-v_{j-1/2}^{n}(U_{j}^{n}-U_{j-1}^{n})}{ \Delta Z^{2}},U_{j}^{n+1}\right)=0, \tag{28}\]  
for all \(j=1\ldots n_{Z}\). The boundary conditions, which are discretized as \(U_{n_{Z}+1}^{n}=U_{n_{Z}}^{n}\) at the free surface (\(Z=h\)) and as \(U_{0}^{n}=-U_{1}^{n}\) at the bottom (\(Z=0\)), are used to provide the ghost values involved in the discretization of the diffusive term. The problem (28) is solved in two steps as  
\[\begin{split}\frac{U_{j}^{n+1/2}-U_{j}^{n}}{\Delta t^{n}}+S(t^{n},Z_{j})-\frac{v_{j+1/2}^{n}(U_{j+1}^{n}-U_{j}^{n})-v_{j-1/2}^{n}(U_{j}^{n}-U_{ j-1}^{n})}{\Delta Z^{2}}=0,\\ U_{j}^{n+1}=\max(U_{j}^{n+1/2},0).\end{split} \tag{29}\]  
Owing to the explicit discretization of the diffusive term, we use the CFL condition  
\[2\sup_{j}v_{j+1/2}^{n}\frac{\Delta t^{n}}{\Delta Z^{2}}\leq 1. \tag{30}\]  
We could also consider an implicit discretization to avoid any CFL condition, but each time step would be more computationally demanding. Finally, the thickness of the static layer can be evaluated as  
\[b^{n}=\left(\min\left\{j\in\left\{1\ldots n_{Z}\right\}\text{ such that }U_{j}^{n}>0\right\}-1\right)\times\Delta Z. \tag{31}\]  
However, to prevent \(b^{n}\) from being influenced by small values of \(U\) and for better accuracy, we prefer to use  
\[b^{n}=\max\!\left(Z_{k^{n}}-\sqrt{2U_{k^{n}}^{n}/C_{0}},0\right),\qquad k^{n} =\min\left\{j\in\left\{1,\ldots,n_{Z}\right\}\text{ such that }U_{j}^{n}\geq C_{0}\Delta Z^{2}/2\right\}, \tag{32}\]  
Figure 8: Velocity \(U_{1}^{n+1}\) versus \(b^{n+1}\). The chosen value for \(b^{n+1}\) is determined by the intersection of the curve with the horizontal axis.  
where \(C_{0}\) is an appropriate constant of the order of \(S/\
u\), in relation to Equations (4a), (4b) and (7).  
In our computations, we use 200 space cells. With the value \(h=0.02\) m, this leads to \(\Delta Z=10^{-4}\) m. Since \(S/\
u\simeq 1.35\) ms\({}^{-2}/(5\times 10^{-5}\) m\({}^{2}\)s\({}^{-1})=2.7\times 10^{4}\) m\({}^{-1}\)s\({}^{-1}\), this gives a velocity threshold of \((S/\
u)\Delta Z^{2}/2\simeq 1.4\times 10^{-4}\) ms\({}^{-1}\). This value is much smaller than the accuracy of experimental velocity measurements which is around 1 cm/s.