## 4 Numerical Solution to the Viscous Free Interface Model
### Variable Viscosity  
The viscoplastic description of granular materials involves fundamentally a Drucker-Prager yield stress proportional to the pressure, the coefficient being the static friction parameter \(\mu_{s}=\tan\delta\). A full rheological law including this yield stress (i.e., defined for all values of the strain rate and not only close to zero) has been proposed in [28], the so called \(\mu(I)\) rheology. As described in [32], this rheology can be interpreted as a decomposition (i.e., (2)) of the deviatoric stress tensor in a rate-independent pure plastic part proportional to \(\mu_{s}\) and a viscous part with pressure- and rate-dependent dynamic viscosity \(\eta\) given by  
\[\eta=(\mu(I)-\mu_{s})\frac{p_{dyn}}{2\|D\|}, \tag{37}\]  
where \(p_{dyn}\) is the dynamic pressure, \(D\) is the strain rate tensor, \(\|D\|^{2}=\frac{1}{2}\sum_{ij}D_{ij}^{2}\), and \(I\) is the inertial number defined by  
\[I=\frac{2\|D\|d}{\sqrt{p_{dyn}/\rho_{s}}}, \tag{38}\]  
with \(d\) the grain diameter as before and \(\rho_{s}\) the grain density. The kinematic pressure \(p\) is related to \(p_{dyn}\) by \(p=p_{dyn}/\rho\), with \(\rho=\phi\rho_{s}\) the density of the granular material, \(\phi\) being the volume fraction. If we consider a slope aligned velocity field depending only on the normal coordinate \(Z\) (simple shear flow), the arguments of Section 2.3 show that the solution to the viscoplastic model (1) and (2) is described by the system (3)-(5) with the constant and uniform source term \(S\) given by (14) and \(\
u=\eta/\rho\). Such a flow has hydrostatic pressure \(p=g\cos\theta(h-Z)\) and shear rate \(\|D\|=\partial_{Z}U/2\). Thus with (37), the kinematic viscosity becomes  
\[\
u=(\mu(I)-\mu_{s})\frac{g\cos\theta(h-Z)}{\partial_{Z}U}, \tag{39}\]  
with  
\[I=\frac{\partial_{Z}Ud}{\sqrt{\phi g\cos\theta(h-Z)}}. \tag{40}\]  
Then the term appearing in (3) is  
\[\
u\partial_{Z}U=(\mu(I)-\mu_{s})g\cos\theta(h-Z), \tag{41}\]  
and the part of this term proportional to \(\mu_{s}\) indeed balances in (3) the corresponding part in \(S\) from (14). The nonlinearity is given according to [28] as  
\[\mu(I)=\mu_{s}+(\mu_{2}-\mu_{s})\frac{I}{I+I_{0}}, \tag{42}\]with \(\mu_{2}>\mu_{s}\) the friction at large strain rate and \(I_{0}\) is a constant that depends on the material used and other device specificities, see ([27], Appendix A). Here we take the value \(I_{0}=0.279\) proposed by [27]. Note that the boundary condition (4c) is automatically satisfied and can therefore be skipped. For the discretization, we use the method of Section 4.2 with the discrete unknowns \(U_{j}^{n}\simeq U(t^{n},Z_{j})\), \(Z_{j}=(j-1/2)\Delta Z\) for \(j=1\ldots n_{Z}\), \(n_{Z}\Delta Z=h\). The Equations (3) and (41) with (42) and (40) are discretized with finite differences under the conservative form  
\[\frac{U_{j}^{n+1/2}-U_{j}^{n}}{\Delta t^{n}}+S(t^{n},Z_{j})-\frac{\Phi_{j+1/2} ^{n}-\Phi_{j-1/2}^{n}}{\Delta Z}=0, \tag{43}\]  
for all \(j=1\ldots n_{Z}\), with  
\[\Phi_{j+1/2}^{n}=(\mu(I_{j+1/2})-\mu_{s})g\cos\theta(h-Z_{j+1/2}),\qquad I_{j+ 1/2}=\frac{d\left(U_{j+1}^{n}-U_{j}^{n}\right)/\Delta Z}{\sqrt{\phi g\cos \theta(h-Z_{j+1/2})}} \tag{44}\]  
for \(j=0\ldots n_{Z}\), with \(Z_{j+1/2}=Z_{j}+\Delta Z/2=j\Delta Z\), and once \(U_{j}^{n+1/2}\) has been computed, we set  
\[U_{j}^{n+1}=\max(U_{j}^{n+1/2},0). \tag{45}\]  
Since \(\Phi_{n_{Z}+1/2}^{n}=0\), there is no need to define \(U_{n_{Z}+1}^{n}\), in accordance with the loss of the free surface boundary condition (4c). On the left boundary (bottom), we use as before the no-slip condition, \(U_{0}^{n}=-U_{1}^{n}\). We use the CFL condition \(2\
u_{\max}\Delta t^{n}\leq\Delta Z^{2}\), with \(\
u_{\max}\) defined as the maximum value of \(\
u\) computed with (39). The interface position \(b^{n}\) is computed according to (32). As previously, we take \(\mu_{s}=\tan(26^{\circ})\), \(h=0.02\) m, \(b^{0}=0.005\) m, and \(d=7\times 10^{-4}\) m. As in [12], we take \(\phi=0.62\).