5 Approximate Numerical Solution

[page 2333, §1]

5.1 Initial conditions

[2333.1.1] The preceding theoretical considerations are confirmed by numerical solutions of eq. (23) with a simple jump-type hysteresis, as formulated in eq. (29). [2333.1.2] The initial conditions are monotone and of the form

$s(z,0)=S^{{\rm out}}+(S^{{\rm P}}-S^{{\rm out}})(1-\mathrm{\Theta}(z-z^{*}))$

(41)

where $S^{{\rm out}}$ is given in Table 1. [2333.1.3] Note that all initial conditions are monotone and do not have an overshoot. [2333.1.4] Two cases are, labelled A and B, will be illustrated in the figures, namely

	$\displaystyle S^{{\rm P}}=0.6554=S^{}\text{ with width }z^{}=0.25,$		(A)
	$\displaystyle S^{{\rm P}}=0.9540\neq S^{}\text{ with width }z^{}=0.015.$		(B)

[2333.1.5] Case A is chosen to illustrate travelling wave solutions with $c_{{\mathrm{im}}}=c_{{\mathrm{dr}}}$ . [2333.1.6] The second initial condition illustrates a general overshoot solution with two wave speeds $c_{{\mathrm{im}}}\neq c_{{\mathrm{dr}}}$ . [2333.1.7] Moreover it illustrates one possibility for the limit $z^{*}\to 0$ and $S^{{\rm P}}\to 1-{S_{{\mathbb{O}\rm r}}}^{{\mathrm{im}}}$ in which a very thin saturated layer (that may arise from sprinkling water on top of the column) initiates an overshoot profile. [2333.1.8] This could provide a partial answer to the question of how saturation overshoot profiles can be initiated.

5.2 Numerical Methods

[2333.2.1] The numerical solution of eq. (29) was performed using the open source toolkit for computational fluid mechanics Open $\nabla$ FOAM [82]. [2333.2.2] This requires to develop an appropriate solver routine. The solver routine employed here was derived from the solver [page 2334, §0] scalarTransportFoam of the toolkit. [2334.0.1] Two different discretizations of eq. (29) have been implemented. One is based on a direct discretization of $\partial f_{i}(S)/\partial z$ using the fvc::div-operators, the other on $f_{i}^{\prime}(S)\times$ $\partial S/\partial z$ using fvc::grad-operators. [2334.0.2] For the discretization of the time derivatives an implicit Euler scheme was chosen, for the discretization of $\partial S/\partial z$ an explicit least square scheme, and for the discretization of the second order term an implicit Gauss linear corrected scheme has been selected. [2334.0.3] The second order term had to be regularized by replacing the function $D_{i}(S)$ with $\max _{{S\in[0,1]}}\{ D_{i}(S),10^{{-3}}\}$ . [2334.0.4] This avoids oscillations at the imbibition front. [2334.0.5] The discretized system was solved with an incomplete Cholesky conjugate gradient solver in Open $\nabla$ FOAM . [2334.0.6] The numerical solutions were found to differ only in the numerical diffusion of the algorithms. [2334.0.7] The divergence formulation seems to be numerically more stable and accurate.

5.3 Results

[2334.1.1] Equation (38) can have one solution, several solutions or no solution. [2334.1.2] This can be seen numerically, but also analytically. [2334.1.3] Depending on the values of the parameters the velocity of the imbibition front may be larger or smaller than that of the drainage front. [2334.1.4] With the parameters from Table 1 the overshoot saturation for a travelling wave as computed from (38) is found to be $S^{*}=0.6554$ and its velocity is $c_{{\mathrm{dr}}}(S^{*})=c_{{\mathrm{im}}}(S^{*})=4.2$ .

[2334.2.1] The travelling wave solution expected from the graphical construction in Figure 1a) and b) for the initial condition A with $S^{{\rm P}}=S^{*}=0.6554$ and width $z^{*}=0.25$ is displayed in Figure 1c) at the initial time $t=0$ and after dimensionless time $t=0.053222$ corresponding to $4450$ s. [2334.2.2] It confirms a travelling wave of the form (36) whose velocity $c^{*}=4.2$ agrees perfectly with that predicted from eq. (38). [2334.2.3] We have also checked that the solution does not change with grid refinement.

[2334.3.1] For different initial conditions the numerical solutions also agree with the theoretical considerations. [2334.3.2] Saturation overshoot is found also when initially a thin saturated layer with steplike or linear profile is present. [2334.3.3] These solutions, however, do not form travelling waves. [2334.3.4] Instead they have $c_{{\mathrm{im}}}\neq c_{{\mathrm{dr}}}$ for their imbibition and drainage front as predicted analytically. [2334.3.5] For $S^{{\rm P}}>S^{*}$ the overshoot region grows linearly with time at the rate $c_{{\mathrm{im}}}-c_{{\mathrm{dr}}}$ . [2334.3.6] Such an overshoot solution is generated by the second initial condition B and compared with the travelling wave solution in Figure 2. [2334.3.7] The non-travelling overshoot for initial condition B with $z^{*}=0.015$ and $S^{{\rm P}}=0.954$ is shown as the solid profiles at times $t=0,0.02392,0.053222$ corresponding to $0,2000,4450$ s. [2334.3.8] The profile quickly approaches an overshoot solution with plateau value $S\approx 0.682$ close to the saturation of the Welge tangent construction. [2334.3.9] This is higher than the plateau value $S^{*}=0.6554$ . [2334.3.10] The difference will be difficult to observe experimentally because of the large uncertainties in the measurements [54, 55]. [2334.3.11] The velocities of the numerical imbibition and drainage fronts in the case of initial conditions B with $z^{*}=0.015$ and $S^{{\rm P}}=0.954$ (solid curves in Fig. 2) again agree perfectly with $c_{{\mathrm{im}}}(0.682)=4.27$ and $c_{{\mathrm{dr}}}(0.682)=2.813$ from the theoretical analysis.

[2334.3.12] Equation (29) with initial and boundary conditions for $S$ and an initial condition for $\partial S/\partial t$ is conjectured to be a well defined semigroup of bounded operators on $L^{1}([0,L])$ on a finite interval $[0,T]$ of time. [2334.3.13] The conjecture is supported by the fact that each of the equations (23) individually defines such a semigroup, and because multiplication by $\Xi$ or $(1-\Xi)$ are projection operators.

[2334.3.14] Equation (29) with initial and boundary conditions for $S$ and an initial condition for $\partial S/\partial t$ is conjectured to be a well defined semigroup of bounded operators on $L^{1}([0,L])$ on a finite interval $[0,T]$ of time. [2334.3.15] The conjecture is supported by the fact that each of the equations (23) individually defines such a semigroup, and because multiplication by $\Xi$ or $(1-\Xi)$ are projection operators.

[2334.4.1] Other values of the initial saturation $S^{{\rm P}}$ have also been investigated. [2334.4.2] A rich phenomenology of possibilities indicates that the existence or nonexistence of overshoot solutions depends very sensitively on the parameters of the problem and the initial conditions.

[2334.5.1] Note also the asymmetric shape of the overshoot solutions (travelling or not). [2334.5.2] This asymmetry resembles the asymmetry seen in experiment [54, 55]. [2334.5.3] While the leading imbibition front is steep, the trailing drainage front is more smeared out. [2334.5.4] This results from the capillary flux term $D_{i}(S)$ whose values at $S^{{\rm out}}$ are smaller than at $S^{{\rm in}}$ .

Figure 2: Numerical solutions $S(z,0)$ , $S(z,0.02392),S(z,0.053222)$ of eq. (29) at dimensionless times $t=0,0.02392,0.053222$ with parameters from Table 1 for initial condition A (dashed) and B (solid). For case A (shown as a short dashed line) the profile for $t=0.053222$ is the same as in Fig.1c). For the solid profiles (case B) the imbibition front moves with speed $c_{{\mathrm{im}}}=4.27$ while the trailing drainage front moves with $c_{{\mathrm{dr}}}=2.813$ .

Authors:	R. Hilfer and R. Steinle
originally published in:	European Physical Journal: Special Topics, Vol. 223, pages 2323-2338 (2014)
originally submitted:	August 18th, 2014