## Appendix C Numerical Model  
To solve the boundary value problem, an explicit centered finite difference scheme has been developed. The time derivative of the water content is expressed as  
\[\frac{\partial w}{\partial t}=\frac{\partial w}{\partial s}\frac{\partial s} {\partial t}=\frac{\partial\ (\theta\rho_{\rm w}+(\theta_{sat}-\theta)\rho_{\rm v})}{\partial s}\ \frac{\partial s}{\partial t}=\left[(\rho_{\rm w}-\rho_{\rm v})\frac{\partial \theta}{\partial s}+(\theta_{\rm sat}-\theta)\frac{\partial\rho_{\rm v}}{ \partial s}\right]\!\frac{\partial s}{\partial t}\] (A6)  
Approximating the derivatives by increments,  
\[\frac{\partial s}{\partial t}\approx\frac{\Delta s}{\Delta t}=\frac{s_{\rm i} ^{\rm j}-\dot{s}_{\rm i}^{\rm j-1}}{\Delta t}\] (A7)  
where \(\dot{s}_{\rm i}^{\rm j}\) refers to the suction in the ith grid point at jth time step. Consequently, the increase of the moisture content in each time step is calculated by  
\[\frac{\Delta w}{\Delta t}=\left[(\rho_{\rm w}-\rho_{\rm v})\frac{\partial\theta }{\partial s}+(\theta_{\rm sat}-\theta)\frac{\partial\rho_{\rm v}}{\partial s }\right]_{\rm i}^{\rm j-1}\!\frac{\Delta s}{\Delta t}\] (A8)  
On the other hand, the water flow rate (a scalar magnitude in this 1D model) can be written as  
\[g=g_{\rm w}+g_{\rm v}=\left[\frac{\kappa\ K}{\mu}\ \rho_{\rm w}-\tau_{\rm v}\ ( \theta_{\rm sat}-\theta)\ D_{\rm v}\frac{\partial\rho_{\rm v}}{\partial s} \right]\!\frac{\partial s}{\partial z}\] (A9)  
The water flow rate is calculated at the center of two consecutive grid points; i.e., in \(i+\frac{1}{2}\) and \(i-\frac{1}{2}\) (see Figure 11). Therefore, it can be assumed that  
\[\dot{g}_{\rm i+\frac{1}{2}}^{\rm j-1}=\left[\frac{\kappa\ K}{\mu}\ \rho_{\rm w}-\tau_{\rm v}\ ( \theta_{\rm sat}-\theta)\ D_{\rm v}\frac{\partial\rho_{\rm v}}{\partial s} \right]_{\rm i+\frac{1}{2}}^{\rm j-1}\frac{s_{\rm i+1}^{\rm j-1}-s_{\rm i}^{ \rm j-1}}{\Delta z}\] (A10)  
\[\dot{g}_{\rm i-\frac{1}{2}}^{\rm j-1}=\left[\frac{\kappa\ K}{\mu}\rho_{\rm w} -\tau_{\rm v}\ (\theta_{\rm sat}-\theta)\ D_{\rm v}\frac{\partial\rho_{\rm v}}{\partial s} \right]_{\rm i-\frac{1}{2}}^{\rm j-1}\frac{s_{\rm i}^{\rm j-1}-s_{\rm i-1}^{ \rm j-1}}{\Delta z}\] (A11)  
Therefore, the divergence of the water flow density can be approximated by  
\[\frac{\partial g}{\partial z}=\frac{\Delta g}{\Delta z}\dot{s}_{\rm i}^{\rm j -1}\approx\frac{\dot{g}_{\rm i+\frac{1}{2}}^{\rm j-1}-\dot{g}_{\rm i-\frac{1}{ 2}}^{\rm j-1}}{\Delta z}\] (A12)Considering this approximation, the suction of the ith grid point at time step \(\mathrm{j}-1\left(s_{\mathrm{i}}^{\mathrm{j}-1}\right)\) is updated to \(s_{\mathrm{i}}^{\mathrm{j}}\) as  
\[s_{\mathrm{i}}^{\mathrm{j}}=s_{\mathrm{i}}^{\mathrm{j}-1}-\frac{\frac{\Delta g }{\Delta z}\left|\frac{\mathrm{j}-1}{\left[\left(\rho_{\mathrm{w}}-\rho_{ \mathrm{v}}\right)\frac{\partial\rho}{\partial s}+\left(\theta_{\mathrm{sat}}- \theta\right)\frac{\partial\rho_{\mathrm{v}}}{\partial s}\right]_{\mathrm{i}}^ {\mathrm{j}-1}}\ \Delta t \tag{13}\]  
The described formulation constitutes a forward-time central-space (FTCS) scheme (with an evaluation of the spatial differences at the midpoints of the grid for a better evaluation of the non-linear storage term).  
The computational time interval is defined by the initial and final time; i.e., \"Tini\" and \"Tffin\". The time step of the interval, \(\Delta t\), was defined in the code, which was an arbitrarily chosen value and remained constant throughout the calculation process. The time was updated after each time step until the final value of the interval (\"Tfin\") was reached, when the calculation process ended.