Iterative schemes for surfactant transport in porous media

TL;DR

This paper develops and analyzes iterative numerical schemes for modeling surfactant transport in variably saturated porous media, focusing on convergence, efficiency, and performance through various linearization techniques.

Contribution

It introduces and compares monolithic and splitting schemes based on Newton, Picard, and L-scheme linearizations for coupled surfactant and water flow equations.

Findings

Schemes converge under certain conditions

Splitting schemes are computationally efficient

Linear systems have manageable condition numbers

Abstract

In this work we consider the transport of a surfactant in a variably saturated porous media. The water flow is modelled by the Richards equations and it is fully coupled with the transport equation for the surfactant. Three linearization techniques are discussed: the Newton method, the modified Picard and the L-scheme. Based on these, monolithic and splitting schemes are proposed and their convergence is analyzed. The performance of these schemes is illustrated on four numerical examples. For these examples, the number of iterations and the condition numbers of the linear systems emerging in each iteration are presented.

Tables

Table 1. Table 1 : Parameters involved in all the examples

$T_{max}$	10 h
$\Omega$	$[0,1]\times [0,1]$
$L_{\Psi }$	.1
$L_c$	.005
Van Genuchten parameters
${\theta }_s$	.42
${\theta }_r$	.0026
$n$	2.9
$\alpha$	.95
$a$	.044
$b$	.4745
Surface tension parameters
$\zeta$	2.4901
$\sigma 0$	73 mN/m
$K_s$	.12 cm/min
$D_0$	6.0e-04
Accuracy requirement
$ϵ$	e-07

Table 2. Table 2 : Example 1A: unsaturated medium, fixed τ = 1 / 10 𝜏 1 10 \tau=1/10

	Monolithic			NonLinS				AltLinS
	Newton			Newton				Newton
					cond. #				cond. #
dx	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/10	20	511.0045	40	333.4035		5.9916	20	333.4019		5.9916
1/20	20	2.2933e+03	40	1.5040e+03		6.2079	20	1.5040e+03		6.2079
1/40	20	9.4458e+03	40	6.1312e+03		6.3234	20	6.1312e+03		6.3234
1/80	20	3.8371e+04	40	2.4774e+04		6.3816	20	2.4774e+04		6.3817
	L Scheme			L Scheme				L Scheme
					cond. #				cond. #
dx	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/10	277	183.4223	540	177.4742		2.1356	264	177.4725		2.1356
1/20	300	812.5650	650	796.5765		2.1839	316	796.5755		2.1839
1/40	363	3.3450e+03	750	3.2584e+03		2.2092	368	3.2584e+03		2.2092
1/80	510	1.3585e+04	850	1.3191e+04		2.2220	421	1.3191e+04		2.2220
	Picard			Picard				Picard
					cond. #				cond. #
dx	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/10	100	326.8280	40	177.4610		5.9916	20	177.4601		5.9916
1/20	110	1.4667e+03	40	796.5170		6.2079	20	796.5129		6.2079
1/40	120	6.0380e+03	40	3.2581e+03		6.3234	20	3.2581e+03		6.3234
1/80	130	2.4522e+04	40	1.3190e+04		6.3816	20	1.3190e+04		6.3817

Table 3. Table 3 : Example 1A: unsaturated medium, fixed dx=1/40

	Monolithic			NonLinS				AltLinS
	Newton			Newton				Newton
					cond. #				cond. #
$\tau$	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/10	20	9.4458e+03	40	6.1312e+03		6.3234	20	6.1312e+03		6.3234
1/20	40	4.7275e+03	80	3.2581e+03		6.3234	40	3.2580e+03		6.3234
1/40	80	2.3677e+03	160	1.7024e+03		6.3234	80	1.7024e+03		6.3234
1/80	160	1.1876e+03	320	870.8016		6.3234	160	870.8010		6.3234
	L Scheme			L Scheme				L Scheme
					cond. #				cond. #
$\tau$	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/10	363	3.3450e+03	750	3.2584e+03		2.2092	368	3.2584e+03		2.2092
1/20	570	1.7540e+03	1300	1.7026e+03		2.2092	633	1.7026e+03		2.2092
1/40	1048	898.9759	2160	870.2808		2.2092	1050	870.8979		2.2092
1/80	1914	455.3332	3520	440.6573		2.2092	1700	440.8161		2.2092
	Picard			Picard				Picard
					cond. #				cond. #
$\tau$	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/10	120	6.0380e+03	40	3.2581e+03		6.3234	20	3.2581e+03		6.3234
1/20	220	3.0216e+03	80	1.7025e+03		6.3234	40	1.7025e+03		6.3234
1/40	400	1.5132e+03	160	870.8263		6.3234	80	870.8251		6.3234
1/80	640	758.9936	320	440.8018		6.3234	160	440.8015		6.3234

Table 4. Table 4 : Example 1B: variably saturated medium, fixed τ = 1 / 10 𝜏 1 10 \tau=1/10

	Monolithic			NonLinS				AltLinS
	Newton			Newton				Newton
					cond. #				cond. #
dx	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/10	28	2.2753e+11	50	5.7734e+09		1.1251	-	5.5783e+09		1.0117
1/20	-	1.2345e+12	-	4.6521e+09		1.0126	-	.6521e+09		1.0126
1/40	-	4.5159e+12	-	5.2321e+09		1.0124	-	5.2321e+09		1.0124
1/80	-	.7232e+13	-	5.5219e+09		1.0123	-	5.5219e+09		1.0123
	L Scheme			L Scheme				L Scheme
					cond. #				cond. #
dx	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/10	175	247.8672	440	239.2940		1.3314	264	239.8408		1.3314
1/20	314	1.0576e+03	650	1.0432e+03		1.3338	316	1.0432e+03		1.3338
1/40	352	4.2256e+03	750	4.1291e+03		1.3328	368	4.1291e+03		1.3328
1/80	408	1.6902e+04	910	1.6437e+04		1.3323	421	1.6437e+04		1.3323
	Picard			Picard				Picard
					cond. #				cond. #
dx	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/10	-	4.5478e+11	50	5.7735e+09		1.1251	-	5.5783e+09		.0117
1/20	-	2.4690e+12	-	4.6521e+09		1.0126	-	4.6521e+09		1.0126
1/40	-	9.0318e+12	-	5.2321e+09		1.0124	-	5.2321e+09		1.0124
1/80	-	3.4465e+13	-	5.5219e+09		1.0123	-	5.5219e+09		1.0123

Table 5. Table 5 : Example 1B: variably saturated medium, fixed dx=1/40

	Monolithic			NonLinS				AltLinS
	Newton			Newton				Newton
					cond. #				cond. #
$\tau$	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/10	-	4.5159e+12	-	5.2321e+09		1.0124	-	5.5783e+09		1.0117
1/20	-	4.5194e+12	-	2.1747e+10		1.0062	-	4.6521e+09		1.0126
1/40	80	4.5265e+12	200	1.0442 e+10		1.1325	-	5.2321e+09		1.0124
1/80	160	4.5406e+12	400	4.3494e+10		1.1325	-	5.5219e+09		1.0123
	L Scheme			L Scheme				L Scheme
					cond. #				cond. #
$\tau$	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/10	352	4.2256e+03	750	4.1291e+03		1.3328	368	4.1291e+03		1.3328
1/20	627	2.2518e+03	1300	2.1862e+03		1.3328	633	2.1890e+03		1.3328
1/40	1100	1.1624e+03	2160	1.1258e+03		1.3328	1050	1.1266e+03		1.3328
1/80	1900	589.7690	3520	570.6119		1.3328	1700	571.0523		1.3328
	Picard			Picard				Picard
					cond. #				cond. #
$\tau$	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/10	-	2.4690e+12	-	5.2321e+09		1.0124	-	5.2321e+09		1.0124
1/20	-	9.0388e+12	-	1.0442e+10		1.0062	-	1.0442e+10		1.0062
1/40	-	9.0529e+12	200	2.1748e+10		1.1325	-	2.0860e+10		1.0031
1/80	-	9.0811e+12	400	4.3494e+10		1.1325	-	4.1698e+10		1.0015

Table 6. Table 6 : Example 2A: unsaturated medium, fixed τ = 1 / 50 𝜏 1 50 \tau=1/50

	Monolithic			NonLinS				AltLinS
	Newton			Newton				Newton
					cond. #				cond. #
dx	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/50	100	7.6243e+09	350	22.2624		3.3828e+09	200	22.1804		3.1779e+09
1/100	100	4.2369e+10	350	25.9419		5.2695e+10	200	31.8223		5.5839e+10
1/150	100	1.0394e+11	350	42.8667		2.8324e+11	200	33.0562		2.8324e+11
1/200	100	1.9614e+11	350	56.4999		9.1099e+11	200	42.8581		9.1662e+11
	L Scheme			L Scheme				L Scheme
					cond. #				cond. #
dx	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/50	2100	2.7971e+09	3700	4.2830		1.6426e+09	2400	4.2489		1.4735e+09
1/100	2100	2.5556e+10	3800	5.0706		2.6387e+10	2700	5.0005		2.6331e+10
1/150	2100	7.6270e+10	3900	6.3840		1.4791e+11	3100	6.3731		1.4165e+11
1/200	2100	1.3652e+11	4100	8.1437		4.5556e+11	3200	8.2452		4.5159e+11
	Picard			Picard				Picard
					cond. #				cond. #
$\tau$	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/50	465	4.3612e+09	400	22.3247		1.7215e+09	200	22.1612		7.6372e+08
1/100	470	3.2796e+10	450	25.9427		2.6395e+10	200	26.0764		1.3243e+10
1/150	480	9.0401e+10	500	33.4899		1.4173e+11	200	33.0714		7.0961e+10
1/200	490	1.7576e+11	500	42.7698		4.5570e+11	200	42.737		2.2773e+11

Table 7. Table 7 : Example 2A: unsaturated medium, fixed dx=1/50

	Monolithic			NonLinS				AltLinS
	Newton			Newton				Newton
					cond. #				cond. #
$\tau$	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/50	100	7.6243e+09	300	22.2624		3.3828e+09	200	22.1804		3.1779e+09
1/100	200	2.7797e+09	350	23.0671		8.0658e+08	400	22.9854		3.8663e+08
1/150	300	1.4883e+09	400	21.3871		3.8252e+08	600	21.9674		2.0727e+08
1/200	400	9.2045e+08	500	21.2462		1.9030e+08	800	21.3715		1.8499e+08
	L Scheme			L Scheme				L Scheme
					cond. #				cond. #
$\tau$	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/50	2100	2.7971e+09	3700	4.2830		1.6426e+09	2400	4.2489		1.4735e+09
1/100	4100	8.0703e+08	7700	4.1966		3.6673e+08	4200	4.1408		1.8802e+08
1/150	5900	3.8356e+08	11250	4.1439		1.6219e+08	5800	4.1052		8.4601e+07
1/200	7700	2.2641e+08	14600	4.0893		9.1518e+07	7200	4.0870		4.7818e+07
	Picard			Picard				Picard
					cond. #				cond. #
$\tau$	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/50	465	4.3612e+09	400	22.3247		1.7215e+09	200	22.1612		7.6372e+08
1/100	880	1.3673e+09	600	23.0000		1.9895e+08	400	21.5095		1.9878e+08
1/150	1300	6.6937e+08	900	22.4168		9.6242e+07	600	21.2971		8.8586e+07
1/200	1700	4.0868e+08	1200	21.2490		5.0376e+07	750	21.2490		5.0376e+07

Table 8. Table 8 : Example 2B: variably saturated medium, fixed τ 𝜏 \tau = dx/100

	Monolithic			NonLinS				AltLinS
	Newton			Newton				Newton
					cond. #				cond. #
dx	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/10	-	2.2435e+10	-	20.9641		3.2267e+08	-	20.9795		3.2267e+08
1/20	-	1.2309e+11	-	20.9642		5.3778e+08	-	20.9800		5.3778e+08
1/40	-	4.5128e+11	-	20.9644		1.2100e+09	-	20.9826		1.2100e+09
	L Scheme			L Scheme				L Scheme
					cond. #				cond. #
dx	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/10	10000	1.1768e+04	4000	5.000		6.2544e+03	3000	5.000		6.2544e+03
1/20	26000	2.5237e+04	16000	5.000		1.4067e+04	6000	5.000		1.4067e+04
1/40	-	5.1895e+04	44000	5.000		2.8441e+04	12000	5.000		2.8440e+04
	Picard			Picard				Picard
					cond. #				cond. #
dx	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/10	-	4.4870e+10	-	20.9641		3.2267e+08	-	20.9726		3.2267e+08
1/20	-	2.4617e+11	-	20.9642		5.3778e+08	-	20.9732		5.3778e+08
1/40	-	9.0255e+11	-	20.9644		1.2100e+09	-	20.9752		1.2100e+09

Table 9. Table 9 : Example 2B: variably saturated medium, fixed dx=1/10

	Monolithic			NonLinS				AltLinS
	Newton			Newton				Newton
					cond. #				cond. #
$\tau$	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/1000	-	2.2435e+10	-	20.9641		3.2267e+08	-	20.9795		3.2267e+08
1/2000	-	2.2444e+10	-	20.9641		6.4533e+08	-	20.9799		6.4533e+08
1/4000	-	2.2461e+10	20000	20.9640		1.2907e+09	-	20.9807		1.2907e+09
	L Scheme			L Scheme				L Scheme
					cond. #				cond. #
$\tau$	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/1000	10000	1.1768e+04	4000	5.000		6.2544e+03	3000	5.000		6.2544e+03
1/2000	14000	5.9591e+03	6000	5.000		3.2036e+03	6000	5.000		3.2036e+03
1/4000	20000	3.0483e+03	12000	5.000		1.6697e+03	12000	5.000		1.6697e+03
	Picard			Picard				Picard
					cond. #				cond. #
$\tau$	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/1000	-	4.4870e+10	-	20.9641		3.2267e+08	-	20.9726		3.2267e+08
1/2000	-	2.4617e+11	-	20.9642		6.4533e+08	-	20.9731		6.4533e+08
1/4000	-	9.0255e+11	-	20.9644		1.2907e+09	-	20.9739		1.2907e+09

Table 10. Table 10 : Example 3: heterogeneous variably saturated medium, fixed τ 𝜏 \tau = T/10

	Monolithic			NonLinS				AltLinS
	Newton			Newton				Newton
					cond. #				cond. #
dx	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/10	-	1.7934e+09	-	240.1589		6.6555e+07	-	240.1589		6.6555e+07
1/20	-	3.1287e+09	-	1.2146e+03		2.6668e+08	-	1.2146e+03		2.6668e+08
1/40	-	4.7011e+10	-	3.3058e+03		1.0721e+09	-	3.3058e+03		1.0721e+09
	L Scheme			L Scheme				L Scheme
					cond. #				cond. #
dx	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/10	602	3.0502e+07	223	13.3371		1.8192e+07	135	13.3371		1.8192e+07
1/20	849	1.0614e+08	222	47.4216		6.3261e+07	152	47.4216		6.3261e+07
1/40	949	4.2918e+08	230	224.9096		2.9895e+08	153	224.9096		2.9895e+08
	Picard			Picard				Picard
					cond. #				cond. #
dx	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
1/10	-	1.5365e+09	-	308.0140		2.1527e+07	-	308.0140		2.1527e+07
1/20	-	5.4220e+09	-	946.9821		9.3559e+07	-	946.9821		9.3559e+07
1/40	-	4.0185e+10	-	4.3966e+03		3.5695e+08	-	4.3966e+03		3.5695e+08

Table 11. Table 11 : Example 3: heterogeneous variably saturated medium, fixed dx= 1/10

	Monolithic			NonLinS				AltLinS
	Newton			Newton				Newton
					cond. #				cond. #
$\tau$	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
T/$10$	-	1.7934e+09	-	240.1589		6.6555e+07	-	240.1589		6.6555e+07
T/${10}^2$	-	6.7021e+08	-	483.2030		6.6082e+06	-	483.2030		6.6082e+06
T/${10}^3$	-	4.6611e+08	-	3.4913e+03		6.6008e+05	-	3.4913e+03		6.6008e+05
	L Scheme			L Scheme				L Scheme
					cond. #				cond. #
$\tau$	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
T/$10$	602	3.0502e+07	223	13.3371		1.8192e+07	135	13.3371		1.8192e+07
T/${10}^2$	6880	3.0503e+06	2350	2.3236		1.6678e+06	1431	2.3236		1.6678e+06
T/${10}^3$	23188	3.0528e+05	26890	1.5385		1.5358e+05	7950	1.5385		1.5358e+05
	Picard			Picard				Picard
					cond. #				cond. #
$\tau$	# iterations	condition #	# iterations	Richards		Transport	# iterations	Richards		Transport
T/$10$	-	1.5365e+09	-	308.0140		2.1527e+07	-	308.0140		2.1527e+07
T/${10}^2$	-	8.1784e+09	-	341.0311		2.5362e+06	-	341.0311		2.5362e+06
T/${10}^3$	-	1.0572e+10	-	366.8343		2.8550e+05	-	366.8343		2.8550e+05

Equations

$$\frac{\partial \theta (\Psi ,c)}{\partial t}-\nabla \cdot (K(\theta (\Psi ,c))\nabla (\Psi +z))=H_1$$

$$\frac{\partial \theta (\Psi ,c)c}{\partial t}-\nabla \cdot (D\nabla c-{\mathbf{u}}_{\mathbf{w}}c)+R(c)=H_2,$$

$$\Psi ={\Psi }_0\text{ and }c=c_0\text{ at }t=0.$$

$$\theta (\Psi )=\{\begin{array}{ll}{\theta }_r+({\theta }_s-{\theta }_r){\left(\frac{1}{1+(-\alpha \Psi {)}^n}\right)}^{\frac{n-1}{n}},& \Psi \le 0\\ {\theta }_s,& \Psi >0,\end{array}$$

$$K(\theta (\Psi ))=\{\begin{array}{ll}{K}_s{\theta }_e^{\frac{1}{2}}{\left[1-{\left(1-{\theta }_e^{\frac{n}{n-1}}\right)}^{\frac{n-1}{n}}\right]}^2,& \Psi \le 0\\ K_s,& \Psi >0,\end{array}$$

$$\theta(\Psi,c):=\theta\Big{(}\frac{\gamma(c)}{\gamma_{0}(c_{0})}\Psi\Big{)},\qquad\text{ with }\qquad\frac{\gamma(c)}{\gamma_{0}(c_{0})}=\frac{1}{1-b\log(c/a+1)}.$$

$$\theta(\Psi,c)=\begin{cases}\theta_{r}+(\theta_{s}-\theta_{r})\left[1/\Big{(}1+\big{(}-\alpha(\frac{1}{1-b\log(c/a+1)})\Psi\big{)}^{n}\Big{)}\right]^{\frac{n-1}{n}},&\Psi\leq 0\\ \theta_{s},&\Psi>0,\end{cases}$$

$$K(\theta (\Psi ,c))=\{\begin{array}{ll}{K}_s{\theta }_e^{\frac{1}{2}}{\left[1-{\left(1-{\theta }_e^{\frac{n}{n-1}}\right)}^{\frac{n-1}{n}}\right]}^2,& \Psi \le 0\\ K_s,& \Psi >0.\end{array}$$

$$<{\partial }_t\theta (\Psi ,c),v_1>+<K(\theta (\Phi ,c))\nabla (\Psi +z),\nabla v_1>=<H_1,v_1>$$

$$<{\partial }_t(\theta (\Psi ,c)c),v_2>+<D\nabla c+{\mathbf{u}}_{\mathbf{w}}c,\nabla v_2>=<H_2,v_2>$$

$$V_h:=\{v_h\in H_0^1(\Omega )s.t.v_{h|T}\in {\mathbb{P}}_1(T),\text{ for any }T\in T_h\},$$

$$\begin{array}{rl}& <\theta ({\Psi }_h^n,c_h^n)-\theta ({\Psi }_h^{n-1},c_h^{n-1}),v_h>\\ +\tau <K(& \theta ({\Psi }_h^n,c_h^n))(\nabla ({\Psi }_h^n)+{\mathbf{e}}_{\mathbf{z}}),\nabla v_h>=\tau <H_1,v_h>\end{array}$$

$$\begin{array}{rl}& <\theta ({\Psi }_h^n,c_h^n)c_h^n-\theta ({\Psi }_h^{n-1},c_h^{n-1})c_h^{n-1},w_h>\\ +& \tau <D\nabla c_h^n+{\mathbf{u}}_{\mathbf{w}}^{\mathbf{n}-1}{c}_h^n,\nabla w_h>=\tau <H_2,w_h>.\end{array}$$

$$\{\begin{array}{l}{F}_1^{lin}({\Psi }^{n,j+1},c^{n,j+1})=0,\\ F_2^{lin}({\Psi }^{n,j+1},c^{n,j+1})=0.\end{array}$$

$$\{\begin{array}{l}{F}_1({\Psi }^{n,j+1},c^{n,j})=0,\text{ followed by }\\ F_2({\Psi }^{n,j+1},c^{n,j+1})=0.\end{array}$$

$$\{\begin{array}{l}{F}_1^{lin}({\Psi }^{n,j+1},c^{n,j})=0,\text{ followed by }\\ F_2^{lin}({\Psi }^{n,j+1},c^{n,j+1})=0.\end{array}$$

$$\begin{array}{r}<\theta ({\Psi }_h^{n,j},c_h^{n,j})-\theta ({\Psi }_h^{n-1},c_h^{n-1}),v_h>+<\frac{\partial \theta }{\partial \Psi }({\Psi }_h^{n,j},c_h^{n,j})({\Psi }_h^{n,j+1}-{\Psi }_h^{n,j}),v_h>\\ +\tau <K(\theta ({\Psi }_h^{n,j},c_h^{n,j}))(\nabla ({\Psi }_h^{n,j+1})+{\mathbf{e}}_{\mathbf{z}}),\nabla v_h>\\ +\tau <\frac{\partial K}{\partial \Psi }(\theta ({\Psi }_h^{n,j},c_h^{n,j}))(\nabla ({\Psi }_h^{n,j+1})+{\mathbf{e}}_{\mathbf{z}})({\Psi }_h^{n,j+1}-{\Psi }_h^{n,j}),\nabla v_h>=\tau <H,v_h>\end{array}$$

$$\begin{array}{rl}& <\theta ({\Psi }_h^{n,j},c_h^{n,j})c_h^{n,j+1}-\theta ({\Psi }_h^{n-1},c_h^{n-1})c_h^{n-1},w_h>\\ & +<\frac{\partial \theta }{\partial c}({\Psi }_h^{n,j},c_h^{n,j})(c_h^{n,j+1}-c_h^{n,j}),v_h>\\ & +\tau <D\nabla c_h^{n,j+1}+{\mathbf{u}}_{\mathbf{w}}^{\mathbf{n}-1}{c}_h^{n,j+1},\nabla w_h>=\tau <H_c,w_h>.\end{array}$$

$$\begin{array}{rl}& <\theta ({\Psi }_h^{n,j},c_h^{n,j})-\theta ({\Psi }_h^{n-1},c_h^{n-1}),v_h>\\ & +<\frac{\partial \theta }{\partial \Psi }({\Psi }_h^{n,j},c_h^{n,j})({\Psi }_h^{n,j+1}-{\Psi }_h^{n,j}),v_h>\\ +\tau <K(\theta (& {\Psi }_h^{n,j},c_h^{n,j}))(\nabla ({\Psi }_h^{n,j+1})+{\mathbf{e}}_{\mathbf{z}}),\nabla v_h>=\tau <H,v_h>\end{array}$$

$$\begin{array}{rl}& <\theta ({\Psi }_h^{n,j},c_h^{n,j})c_h^{n,j+1}-\theta ({\Psi }_h^{n-1},c_h^{n-1})c_h^{n-1},w_h>\\ & +<\frac{\partial \theta }{\partial c}({\Psi }_h^{n,j},c_h^{n,j})(c_h^{n,j+1}-c_h^{n,j}),w_h>\\ +\tau <D\nabla & c_h^{n,j+1}+{\mathbf{u}}_{\mathbf{w}}^{\mathbf{n}-1}{c}_h^{n,j+1},\nabla w_h>=\tau <H_c,w_h>.\end{array}$$

$$\begin{split}\theta(\Psi^{n,j+1}_{h},c^{n,j+1}_{h})\rightarrow\theta(\Psi^{n,j}_{h},c^{n,j}_{h})&+\Big{(}\frac{\partial\theta}{\partial\Psi}\Big{)}(\Psi^{n,j}_{h},c^{n,j}_{h})(\Psi^{n,j+1}_{h}-\Psi^{n,j}_{h})\\ &+\Big{(}\frac{\partial\theta}{\partial c}\Big{)}(\Psi^{n,j}_{h},c^{n,j}_{h})(c^{n,j+1}_{h}-c^{n,j}_{h}).\end{split}$$

$$\begin{array}{rl}& <\theta ({\Psi }_h^{n,j},c_h^{n,j})-\theta ({\Psi }_h^{n-1},c_h^{n-1}),v_h>+L_1<{\Psi }_h^{n,j+1}-{\Psi }_h^{n,j},v_h>\\ & \tau <K(\theta ({\Psi }_h^{n,j},c_h^{n,j}))(\nabla ({\Psi }_h^{n,j+1})+{\mathbf{e}}_{\mathbf{z}}),\nabla v_h>=\tau <H,v_h>,\end{array}$$

$$\begin{array}{rl}<\theta ({\Psi }_h^{n,j},c_h^{n,j})c_h^{n,j+1}-& \theta ({\Psi }_h^{n-1},c_h^{n-1})c_h^{n-1},w_h>\\ +L_2<c_h^{n,j+1}-c_h^{n,j},& w_h>+\tau <D\nabla c_h^{n,j+1}+{\mathbf{u}}_{\mathbf{w}}^{\mathbf{n}-1}{c}_h^{n,j+1},\nabla w_h>\\ & =\tau <H_c,w_h>.\end{array}$$

$$\begin{array}{rl}<\theta ({\Psi }_h^{n,j+1},& c_h^{n,j})-\theta ({\Psi }_h^{n-1},c_h^{n-1}),v_h>\\ +\tau <K(\theta ({\Psi }_h^{n,j},c_h^{n,j})& )(\nabla ({\Psi }_h^{n,j+1})+{\mathbf{e}}_{\mathbf{z}}),\nabla v_h>=\tau <H,v_h>\end{array}$$

$$\begin{array}{rl}<\theta ({\Psi }_h^{n,j+1},& c_h^{n,j+1})c_h^{n,j+1}-\theta ({\Psi }_h^{n-1},c_h^{n-1})c_h^{n-1},w_h>+\tau <D\nabla c_h^{n,j+1}\\ & +{\mathbf{u}}_{\mathbf{w}}^{\mathbf{n}-1}{c}_h^{n,j+1},\nabla w_h>=\tau <H_c,w_h>.\end{array}$$

$$\begin{array}{rl}& <\theta ({\Psi }_h^{n,j},c_h^{n,j})-\theta ({\Psi }_h^{n-1},c_h^{n-1}),v_h>\\ +<{\theta }^{\prime }({\Psi }_h^{n,j},c_h^{n,j})(& {\Psi }_h^{n,j+1}-{\Psi }_h^{n,j}),v_h>+\tau <K(\theta ({\Psi }_h^{n,j},c_h^{n,j}))(\nabla ({\Psi }_h^{n,j+1})\\ & +{\mathbf{e}}_{\mathbf{z}}),\nabla v_h>+\tau <\frac{\partial K}{\partial \Psi }(\theta ({\Psi }_h^{n,j},c_h^{n,j}))(\nabla ({\Psi }_h^{n,j+1})\\ & +{\mathbf{e}}_{\mathbf{z}})({\Psi }_h^{n,j+1}-{\Psi }_h^{n,j}),\nabla v_h>=\tau <H,v_h>\end{array}$$

$$\begin{array}{rl}<\theta ({\Psi }_h^{n,j+1},& c_h^{n,j})c_h^{n,j+1}-\theta ({\Psi }_h^{n-1},c_h^{n-1})c_h^{n-1},w_h>\\ +<\frac{\partial \theta }{\partial c}({\Psi }_h^{n,j+1},& c_h^{n,j})(c_h^{n,j+1}-c_h^{n,j}),v_h>+\tau <D\nabla c_h^{n,j+1}\\ & +{\mathbf{u}}_{\mathbf{w}}^{\mathbf{n}-1}{c}_h^{n,j+1},\nabla w_h>=\tau <H_c,w_h>.\end{array}$$

$$\begin{array}{rl}& <\theta ({\Psi }_h^{n,j},c_h^{n,j})-\theta ({\Psi }_h^{n-1},c_h^{n-1}),v_h>\\ & +<\frac{\partial \theta }{\partial \Psi }({\Psi }_h^{n,j},c_h^{n,j})({\Psi }_h^{n,j+1}-{\Psi }_h^{n,j}),v_h>\\ +\tau <K(\theta ({\Psi }_h^{n,j},& c_h^{n,j}))(\nabla ({\Psi }_h^{n,j+1})+{\mathbf{e}}_{\mathbf{z}}),\nabla v_h>=\tau <H,v_h>\end{array}$$

$$\begin{array}{rl}<\theta ({\Psi }_h^{n,j+1},& c_h^{n,j})c_h^{n,j+1}-\theta ({\Psi }_h^{n-1},c_h^{n-1})c_h^{n-1},w_h>\\ & +<\frac{\partial \theta }{\partial c}({\Psi }_h^{n,j+1},c_h^{n,j})(c_h^{n,j+1}-c_h^{n,j}),w_h>\\ +\tau <D\nabla & c_h^{n,j+1}+{\mathbf{u}}_{\mathbf{w}}^{\mathbf{n}-1}{c}_h^{n,j+1},\nabla w_h>=\tau <H_c,w_h>.\end{array}$$

$$\begin{array}{rl}& <\theta ({\Psi }_h^{n,j},c_h^{n,j})-\theta ({\Psi }_h^{n-1},c_h^{n-1}),v_h>+L_1<{\Psi }_h^{n,j+1}-{\Psi }_h^{n,j},v_h>\\ & \tau <K(\theta ({\Psi }_h^{n,j},c_h^{n,j}))(\nabla ({\Psi }_h^{n,j+1})+{\mathbf{e}}_{\mathbf{z}}),\nabla v_h>=\tau <H,v_h>\end{array}$$

Full text

Iterative schemes for surfactant transport in porous media††thanks: Project founded by VISTA, a collaboration between the Norwegian Academy of Science and Letters and Equinor.

Davide Illiano davide.illiano@uib.no Department of Mathematics, University of Bergen

Iuliu Sorin Pop sorin.pop@uhasselt.be Department of Mathematics, University of Bergen

Faculty of Sciences, University of Hasselt

Florin Adrian Radu florin.radu@uib.no@university.edu Department of Mathematics, University of Bergen

Abstract In this work we consider the transport of a surfactant in variably saturated porous media. The water flow is modelled by the Richards equations and it is fully coupled with the transport equation for the surfactant. Three linearization techniques are discussed: the Newton method, the modified Picard and the L-scheme. Based on these, monolithic and splitting schemes are proposed and their convergence is analyzed. The performance of these schemes is illustrated on four numerical examples. For these examples, the number of iterations and the condition numbers of the linear systems emerging in each iteration are presented.

Keywords Richards equation, reactive transport, linearization schemes, L-scheme, modified Picard, Newton method, splitting solvers.

1 Introduction

Many societal relevant problems are involving multiphase flow and multicomponent, reactive transport in porous media. Examples in this sense appear in the enhanced oil recovery, geological $CO_{2}$ storage, diffusion of medical agents into the human body, or water or soil pollution. In all these situations, experimental results are difficult and expensive to obtain, therefore numerical simulations become a key technology. Together with lab experiments and field data, they can help us comprehend these complex phenomena. The mathematical models for problems as mentioned above are (fully or partially) coupled, non-linear, possible degenerate partial differential equations. In most cases, deriving explicit solutions is not possible, whereas developing appropriate algorithms for finding numerical solutions is a challenge in itself. Here we investigate robust and efficient methods for solving the nonlinear problems obtained after performing an implicit time discretization, the focus being on iterative, splitting or monolithic schemes for fully coupled flow and transport.

Of particular interest here is a special case of multiphase, reactive flow in porous media, namely the surfactant transport in soil [2, 19, 31, 22, 24, 25]. Surfactants, which are usually organic compounds, are commonly used for actively combating soil and water pollution [16, 36, 11, 41, 12]. They contain both hydrophobic and hydrophilic groups and are dissolved in the water phase, being transported by diffusion and convection. Typically, the surfactants are employed in soil regions near the surface (vadose zone), where water and air are present in the pores. Consequently, the outcoming mathematical model accounts the transport of at least one species (the surfactant, but often also the contaminant) in a variably saturated porous medium. Whereas the dependence of the species transport on the flow is obvious, one can encounter the reverse dependence as well when surfactants are affecting the interfacial tension between water and air, leading to a dependency of the water flow on the concentration of surfactant. In other words, one has to cope with a fully coupled flow and transport problem, and not only with a one-way coupling, i.e. when only the transport depends on the flow, as mostly considered in reactive transport [33].

Whereas the surfactant transport is described by a reaction-diffusion-convection equation, water flow in variably saturated porous media is modelled by the Richards equation [7, 18]. The main assumption in this case is that the air remains in contact with the atmosphere, having a constant pressure (the atmospheric pressure, here assumed zero). This allows reducing the flow model to one equation, the Richards’ equation. In mathematical terms, this equation is degenerate parabolic, whose solution has typically low regularity [3].

From the above, and adopting the pressure head as the main unknown in the Richards’ equation, we study here different linearization schemes for the model

[TABLE]

and

[TABLE]

holding for $\vec{x}\in\Omega$ ($z$ being the vertical coordinate of $\vec{x}$, pointing against gravity) and $t\in(0,T]$. Here $\Omega$ is a bounded, open domain in $\mathbb{R}^{d}$ ($d=1,2$ or $3$) having a Lipschitz continuous boundary $\partial\Omega$ and $T>0$ is the final time. Further, $\theta(\cdot,\cdot)$ denotes the water content, and is a given function depending on the pressure head $\Psi$ and of the surfactant concentration $c$. Also, $K(\cdot)$ is the hydraulic conductivity, $D>0$ the diffusion/dispersion coefficient. Finally, ${{\bf{u_{w}}}}:=-K(\theta,c)\nabla(\Psi+z)$ is the water flux, $R(\cdot)$ the reaction term expressed as a function of the concentration $c$, and $H_{1},H_{2}$ are the external sinks/sources. Initial and boundary conditions, which are specified below, complete the system.

We point out that the water content and the hydraulic conductivity, $\theta(\cdot,\cdot)$ and $K(\cdot)$ are given non-linear functions. They are medium- and surfactant-dependent and are determined experimentally (see [18]). Specific choices are provided in Section 2.

To solve numerically the system (1) – (2) one needs to discretize in time and space. We refer to [15] for a practical review of numerical methods for the Richards equation. Due to the low regularity of the solution and the need of relatively large time steps, the backward Euler method is the best candidate for the time discretization. Multiple spatial discretization techniques are available, such as the Galerkin Finite Element Method (FEM) [30, 5, 37], the Mixed Finite Element Method (MFEM) [4, 34, 42, 44], the Multi-Point Flux Approximation (MPFA) [23, 6, 1] and the Finite Volume Method (FVM) [9, 13, 14].

Since the time discretization is not explicit, the outcome is a sequence of non-linear problems, for which a linearization step has to be performed. Widely used linearization schemes are the quadratic, locally convergent Newton method and the modified Picard method [10]. For both, the convergence is guaranteed if the starting point is close to the solution. Since for evolution equations the initial guess is typically the solution at the previous time, this may induces severe restrictions on the time step size (see [35]). Among alternative approaches we mention the L-scheme (see [45, 38, 32, 28]) and the modified L-scheme [29], both being robust w.r.t. the mesh size, but converging linearly. In particular, the L-scheme converges for any starting point, and the restriction on the time step, if any, is very mild. The modified L-scheme makes explicit use of the choice of the starting point as the solution obtained at the previous time, and has an improved convergence behaviour if the changes in the solutions at two successive times are controlled by the time step. Neverhteless, the modified L-scheme involves computation of derivatives while the L-scheme does not. Finally, the robustness of the Newton method is significantly increased if one considers combinations of the Picard and the Newton methods [8], and in particular of the L-scheme and the Newton scheme [28].

We conclude this discussion by mentioning that in this paper we adopt the FEM and the MPFA, but the iterative schemes presented here can be applied in combination with any other spatial discretization. The focus is on effectively solving the flow and transport system (1) – (2), and in particular on the adequate treating of the coupling between the two model components (the flow and the reactive transport). The schemes are divided in three main categories: monolithic (Mon), non-linear splitting (NonLinS) and alternate splitting (AltS). Subsequently, we denote e.g. by Mon-NE, the monolithic scheme obtained by applying the Newton method as linearization. The nonlinear splitting schemes (NonLinS) should be understood as solving at each time step first the flow equation until convergence, by using the surfactant concentration from the last iteration, and then with the obtained flow solving the transport equation until convergence. The procedure can be continued iteratively, this being the usual or classical splitting method for transport problems. The convergence of NonLinS does not depend on the linearization approach used for each model component (Newton, Picard or L-scheme), because we assume that the nonlinear subproblems are solved exactly, i.e. until convergence. Finally, the alternate splitting methods (AltS) have a different philosophy. Instead solving each subproblem until convergence within each iteration, one performs only one step of the chosen linearization. For example, AltS-NE will perform one Newton step for each model component, and iterate. These schemes are illustrated in Figures 1, 2.

All the schemes can be analysed theoretically, and we do this exemplary for Mon-LS, i.e. for the monolithic approach combined with the L-scheme. Based on comparative numerical tests performed for academic and benchmark problems, we see that the alternate methods can save substantial computational time, while maintaining the robustness of the L-scheme.

The remaining of the paper is organized as follows. In Section 2 we establish the mathematical model and the notation used and present the iterative monolithic and splitting schemes. In Section 3 we prove the convergence of the $Mon-LS$ scheme and briefly discuss the convergence of the other schemes. Section 4 presents four different numerical examples. They are inspired by the cases already studied in the literature [28, 24]. Section 5 concludes this work.

2 Problem formulation, discretization and iterative schemes

We solve the fully coupled system (1)–(2), completed by homogeneous Dirichlet boundary conditions for both $\Psi$ and $c$ and the initial conditions:

[TABLE]

We use the van Genuchten-Mualem parameterization [17]

[TABLE]

[TABLE]

where $\theta_{r}$ and $\theta_{s}$ represent the values of the residual and saturated water content, $\theta_{e}=(\theta(Psi)-\theta_{r})/(\theta_{s}-\theta_{r})$ is the effective water content, $K_{s}$ is the conductivity and $\alpha$ and $n$ are model parameters depending on the soil.

Observe that in the expression above for $\theta$ the influence of the surfactant on the water flow is neglected. As reported in [20, 24, 40], the surface tension between water and air does depend on the surfactant concentration $c$, implying the same for the function $\theta$ above. The following parametrization is proposed in [24]

[TABLE]

Here $\gamma$ and $\gamma_{0}$ are the surface tensions at concentrations $c$ and $c_{0}$, the second being a reference concentration. The parameters $a$ and $b$ depend on the fluid and the medium. We refer to [39, 40] for details about the scaling factor in (5).

This gives the following expressions for $\theta$ and $K$

[TABLE]

[TABLE]

This shows that the flow component also depends on the reactive transport, implying that the model is coupled in both directions.

In the following we proceed by discretizing the equations (1) and (2) in time and space. We will use common notations in functional analysis. We denote by $L^{2}(\Omega)$ the space of real valued, squared integrable function defined on $\Omega$ and $H^{1}(\Omega)$ its subspace, containing the functions having also the first order derivatives in $L^{2}(\Omega)$. $H_{0}^{1}(\Omega)$ is the space of functions belonging to $H^{1}(\Omega)$ and vanishing on $\partial\Omega$. Further, we denote by $<\cdot,\cdot>$ the $L^{2}(\Omega)$ scalar product (and by $\left\lVert\cdot\right\rVert$ the associated norm) or the pairing between $H1_{0}$ and its dual $H^{-1}$. Finally, by $L^{2}(0,T;X)$ we mean the Bochner space of functions taking values in the Banach-space $X$, the extension to $H^{1}(0,T;X)$ being straightforward.

With this we state the weak formulation of the problem related to (1) – (2):

Problem P: Find $\Psi,c\in L^{2}(0,T;H_{0}^{1}(\Omega))\cap H^{1}(0,T;H^{-1}(\Omega))$ such that

[TABLE]

and

[TABLE]

hold for all $v_{1},v_{2}\in H_{0}^{1}(\Omega)$ and almost every $t\in(0,T]$.

We now combine the backward Euler method with linear Galerkin finite elements for the discretization of Problem P. We let $N\in\mathbb{N}$ be a strictly positive natural number and the time step $\tau\ :=\ T/N$. Correspondingly, the discrete times are $t_{n}\ :=\ n\tau\ (n\in\{0,1,\dots,N\})$. Further, we let $T_{h}$ be a regular decomposition of $\Omega$, $\overline{\Omega}=\underset{T\in T_{h}}{\cup}T$ into $d$-dimensional simplices, with $h$ denoting the mesh diameter. The finite element space $V_{h}\subset H_{0}^{1}(\Omega)$ is defined by

[TABLE]

where $\mathbb{P}_{1}(T)$ denotes the space of linear polynomials on $T$ and $v_{h|T}$ the restriction of $v_{h}$ to $T$.

For the fully discrete counterpart of Problem P we let $n\geq 1$ be fixed and assume that $\Psi^{n-1}_{h},c^{n-1}_{h}\in V_{h}$ are given. The solution pair at time $t_{n}$ solves

Problem Pn: Find $\Psi^{n}_{h},c^{n}_{h}\in V_{h}$ such that for all $v_{h},w_{h}\in V_{h}$ there holds

[TABLE]

and

[TABLE]

${{\bf{e_{z}}}}$ denotes the unit vector in the direction opposite to gravity.

Remark 1

Observe that ${\bf{u^{n-1}_{w}}}$ appears in the convective term in (12). This choice is made for the ease of presentation. Nevertheless, all calculations carried out in this paper were doubled by ones where $\bf{u^{n}_{w}}$ has replaced ${\bf{u^{n-1}_{w}}}$. The differences in the results were marginal.

Observe that Problem Pn is a coupling system of two elliptic, nonlinear equations. In the following we discuss different iterative schemes for solving this system.

2.1 Iterative linearization schemes

We discuss monolithic and splitting approaches for solving Problem Pn, combined with either the Newton-method, or the modified Picard [10] or the L-scheme [32, 28]. In the following the index $n$ always refers to the time step, whereas $j$ denotes the iteration index. As a rule, the iterations start with the solution at the last time, $t_{n-1}$.

In the monolithic approach one solves the two equations of the system (11)-(12) at once, and combined with a linearization method. Formally, this becomes

Problem PMonn,j+1: Find $\Psi^{n,j+1}$ and $c^{n,j+1}$ such that

[TABLE]

$F_{k}^{Lin}$ is a linearization of the expression $F_{k}$ ($k=1,2$) appearing in the system (11)-(12). Depending on the used linearization technique, one speaks about a monolithic-Newton scheme (Mon-Newton), or monolithic-Picard (Mon-Picard) or monolithic-L-scheme (Mon-LS). These three schemes will be presented in detail below.

In the iterative splitting approach one solves each equation separately and then iterates between these using the results obtained previously. We distinguish between two main splitting ways: the nonlinear slitting and the alternate splitting. This is schematized in Figure 1 and Figure 2 respectively. The former becomes

Problem PNonLinSn,j+1: Find $\Psi^{n,j+1}$ and $c^{n,j+1}$ such that

[TABLE]

For the linearization of $F_{1}$ and $F_{2}$ one can use one of the three linearization techniques mentioned before. In contrast, in the alternate splitting one performs only one linearization step per iteration, see also Figure 2. The alternate splitting scheme becomes

Problem PAltSn,j+1: Find $\Psi^{n,j+1}$ and $c^{n,j+1}$ such that

[TABLE]

Depending on which linearization is used, one speaks about alternate splitting Newton (AltS-NE) or alternate splitting L-scheme (AltS-LS). Both schemes are presented in detail below.

2.1.1 The monolithic Newton method (Mon-Newton)

We recall that the Newton scheme is quadratically, but only locally convergent. The monolithic Newton method applied to (11)-(12) gives

Problem PMon-Newtonn,j+1: Let $\Psi^{n-1}_{h},c^{n-1},\Psi^{n,j}_{h},c^{n,j}_{h}\in V_{h}$ be given,

find $\Psi^{n,j+1}_{h},c^{n,j+1}_{h}\in V_{h}$ such that for all $v_{h},w_{h}\in V_{h}$ one has

[TABLE]

and

[TABLE]

2.1.2 The monolithic Picard method (Mon-Picard)

The modified Picard method was initially proposed by Celia [10] for the Richards equation. It is similar to the Newton method in dealing with the nonlinearity in the saturation, but not in the permeability. Being a modification of the Newton method, modified Picard method is only linearly convergent [35]. The monolithic Picard method applied to (11)-(12) becomes

Problem PMon-Picardn,j+1: Let $\Psi^{n-1}_{h},c^{n-1}_{h},\Psi^{n,j}_{h},c^{n,j}_{h}\in V_{h}$ be given,

find $\Psi^{n,j+1}_{h},c^{n,j+1}_{h}\in V_{h}$ such that for all $v_{h},w_{h}\in V_{h}$ one has

[TABLE]

and

[TABLE]

The equations (18) and (19) have been expressed as functions of only the unknown pressure $\Psi_{h}^{n,j+1}$ and concentration $c^{n,j+1}_{h}$, respectively. To achieve this, in the former we used only the derivative of $\theta$ with respect to $\Psi$ and only the derivative of $\theta$ with respect to $c$ in the latter.

Alternatively, both of the partial derivatives can be involved in the formulation,

[TABLE]

2.1.3 The monolithic L-scheme (Mon-LS)

The monolithic L-scheme for solving (11)–(12) becomes

Problem PMon-LSn,j+1: Let $\Psi^{n-1}_{h},\Psi^{n,j}_{h},c^{n-1}_{h},c^{n,j}_{h}\in V_{h}$ be given and

with $L_{1},L_{2}>0$ large enough (as specified below), find $\Psi^{n,j+1}_{h},c^{n,j+1}_{h}\in V_{h}$ s.t. for all $v_{h},w_{h}\in V_{h}$

[TABLE]

[TABLE]

$L_{1}$ and $L_{2}$ are free to be chosen parameters but should be large enough to ensure the convergence of the scheme, see Sec. 3. In practice, the values of $L_{1},L_{2}$ are connected to $\displaystyle\max_{\Psi}\left\lVert\frac{\partial\theta}{\partial\Psi}\right\rVert$, $\displaystyle\max_{c}\left\lVert\frac{\partial\theta}{\partial c}\right\rVert$.

The L-scheme does not involve the computations of derivatives, and the linear systems to be solved within each iteration are better conditioned compared to the ones given by Newton or Picard method (see [28]). Moreover, this scheme is (linearly) convergent for any initial guess for the iteration.

2.1.4 The non-linear splitting approach (NonLinS)

The non-linear splitting approach for solving (11)–(12) becomes

Problem PNonLinSn,j+1: Let $\Psi^{n-1}_{h},c^{n-1},\Psi^{n,j}_{h},c^{n,j}_{h}\in V_{h}$ be given, find $\Psi^{n,j+1}_{h}\in V_{h}$ s.t.

[TABLE]

holds true for all $v_{h}\in V_{h}$. Then, with $\Psi^{n,j+1}_{h}$ obtained, find $c^{n,j+1}_{h}\in V_{h}$ such that for all $w_{h}\in V_{h}$ it holds

[TABLE]

As for the monolithic schemes, one can apply the different linear iterative schemes to obtain fully linear versions of the splitting approach. This is done first to solve (23) and, once a solution to (23) is available, this is employed in the linearization of (24).

2.1.5 The alternate Newton method (AltS-Newton)

In the alternate Newton method applied to (11)-(12) one solves

Problem PAltS-Newtonn,j+1: Let $\Psi^{n-1}_{h},c^{n-1},\Psi^{n,j}_{h},c^{n,j}_{h}\in V_{h}$ be given,

find $\Psi^{n,j+1}_{h}\in V_{h}$ s.t.

[TABLE]

holds true for all $v_{h}\in V_{h}$. Then, with $\Psi^{n,j+1}_{h}$ obtained above, find $c^{n,j+1}_{h}\in V_{h}$ such that for all $w_{h}\in V_{h}$ one has

[TABLE]

2.1.6 The alternate Picard method (AltS-Picard)

The alternate Picard method applied to (11)-(12) becomes

Problem PAltS-Picardn,j+1: Let $\Psi^{n-1}_{h},c^{n-1},\Psi^{n,j}_{h},c^{n,j}_{h}\in V_{h}$ be given,

find $\Psi^{n,j+1}_{h}\in V_{h}$ s.t.

[TABLE]

hold true for all $v_{h}\in V_{h}$. Then, with $\Psi^{n,j+1}_{h}$ obtained above, find $c^{n,j+1}_{h}\in V_{h}$ such that for all $w_{h}\in V_{h}$ one has

[TABLE]

2.1.7 The alternate L-scheme (AltS-LS)

The alternate L-scheme for solving (8–9) becomes

Problem PAltS-LSn,j+1: Let $\Psi^{n-1}_{h},c^{n-1},\Psi^{n,j}_{h},c^{n,j}_{h}\in V_{h}$ be given, find $\Psi^{n,j+1}_{h}\in V_{h}$ s.t.

[TABLE]

hold true for all $v_{h}\in V_{h}$. Then, with $\Psi^{n,j+1}_{h}$ obtained above, find $c^{n,j+1}_{h}\in V_{h}$ such that for all $w_{h}\in V_{h}$ one has

[TABLE]

Remark 2

(Stopping criterion) For both monolithic and splitting schemes, one stops the iteration process whenever

[TABLE]

where $\epsilon_{1},\epsilon_{2}$ are small numbers. Here we took $10^{-07}$ or $10^{-08}$.

3 Convergence analysis

In this section we analyse the convergence of the monolithic L-scheme introduced through Problem PMon-LSn,j+1. We restrict the analysis to this iteration, but mention that the convergence analysis for the other (monolithic or splitting) schemes introduced above, can be done in a similar fashion. We start by defining the errors

[TABLE]

obtained at iteration $j+1$. The scheme is convergent if both errors vanish when $j\rightarrow\infty$.

The convergence is obtained under the following assumptions:

• (A1)

There exist $\alpha_{\Psi}>0$ and $\alpha_{c}\geq 0$ such that for any $\Psi_{1},\Psi_{2}\in\mathbb{R}$ and $c_{1},c_{2}\in\mathbb{R}_{+}$

[TABLE]

Furthermore, there exist two constants $\theta_{m}\geq 0$ and $\theta_{M}<\infty$ such that $\theta_{m}\leq\theta\leq\theta_{M}$.

• (A2)

The function $K(\theta(\cdot,\cdot))$ is Lipschitz continuous, with respect to both variables, and there exist two constants $K_{m}$ and $K_{M}$ such that $0\leq K_{m}\leq K\leq K_{M}<\infty$.

• (A3)

There exist $M_{u},M_{\Psi},M_{c}\geq 0$ such that

$\left\lVert\bf{u^{n}_{w}}\right\rVert_{L^{\infty}}\leq M_{u}$, $\left\lVert\nabla\Psi^{n}\right\rVert_{L^{\infty}}\leq M_{\Psi}$ and $\left\lVert c^{n}\right\rVert_{L^{\infty}}\leq M_{c}$ for all $n\in\mathbb{N}$.

Remark 3

(A2) is satisfied in most realistic situations. (A3) is a pure technical one, being satisfied when data is sufficiently regular, which is assumed to be the case for the present analysis. The inequality (32) in (A1) is a coercivity assumption. It is in particular satisfied if $\theta$ only depends on $\Psi$, and for common relationships $\theta--\Psi$ encountered in the engineering literature.

Theorem 1

Let $n\ \in\ \{1,2,\dots N\}$ be given and assume (A1)-(A3) be satisfied. If the time step is small enough (see (42) below), the monolithic L-scheme in (29) and (30) is linearly convergent for any $L_{1}$ and $L_{2}$ satisfying (41).

Proof 1

We follow the ideas in [32, 28] and start by subtracting (11) from (29) to obtain the error equation

[TABLE]

Testing now the above equation with $v_{h}=e_{\Psi}^{j+1}$, one obtains

[TABLE]

By (A2) and after some algebraic manipulations we further get

[TABLE]

Using now (A1), (A3), the Lipschitz continuity of $K$ and twice the Young and Cauchy-Schwarz inequalities, for any $\delta_{0}>0$ and $\delta_{1}>0$, from (35) one obtains

[TABLE]

Similarly, subtracting (12) from (30) and choosing $w_{h}=e_{c}^{j+1}$ in the resulting one gets

[TABLE]

This can be rewritten as

[TABLE]

Using again (A1), (A3) and the Cauchy-Schwarz and Young inequalities, from (38) it follows that for any $\delta_{2},\delta_{3},\delta_{4}>0$ one has

[TABLE]

Summing adding (36) to (39) and using (A1) one gets

[TABLE]

Choosing $\delta_{0}=\delta_{2}=\dfrac{\alpha_{\Psi}}{2}$, $\delta_{1}=\dfrac{K_{m}}{2}$, $\delta_{3}=\theta_{m}$ and $\delta_{4}=\dfrac{D}{2}$ in (40), and assuming that

[TABLE]

and the time step $\tau$ satisfies the mild conditions

[TABLE]

where $C_{\Omega}$ denotes the Poincare constant, then we obtain

[TABLE]

Finally, by using the Poincare inequality two times we get from (43)

[TABLE]

From (42), (44) implies that the errors are contracting and therefore the monolithic L-scheme (29) - (30) is convergent.

Remark 4

The convergence rate resulting from (44) does not depend on the spatial mesh size. Also, observe that this convergence is obtained for any initial guess. Based on this, the method is globallly convergent, which is in contrast to the Newton or (modified) Picard schemes, converging only locally. It can be observed that, the larger the time step, the smaller the constants $L_{1}$ and $L_{2}$ are, resulting in a faster convergence. For small steps instead the convergence rate can approach 1. On the other hand, if the time step is small enough, one may reach the regime where the Newton scheme becomes convergent (see [35]). Alternatively, one may first perform a number of L-scheme iterations, and use the resulting as an initial guess for the Newton scheme (see [28]), or consider the modified L-scheme in [29]. In either situations, the convergence behaviour was much improved.

Remark 5

The convergence of the modified Picard and Newton method applied to the Richards equation has been already proved in [35]. Such results can be extended to the coupled problems considered here.

4 Numerical examples

In this section we consider four test cases for the proposed linearization schemes, inspired by the literature [28, 24]. The schemes have been implemented in the open source software package MRST [27], an open source toolbox based on Matlab, in which multiple solvers and models regarding flows in porous media are incorporated.

Example 1A: flow and transport in strictly unsaturated media

We start our numerical studies with a manufactured problem admitting an analytical solution [28]. The unit square $\Omega$ is divided into two sub-domains: $\Omega_{up}$ and $\Omega_{down}$. The two regions are defined as: $\Omega_{up}=[0,1]\times[1/4,1]$ and $\Omega_{down}=[0,1]\times[0,1/4]$. Dirichlet boundary conditions, $\Psi=-3$, and no-flow Neumann boundary conditions are imposed on $\Gamma_{D}=[0,1]\times{1}$ and $\Gamma_{N}=\partial\Omega/\Gamma_{D}$, respectively. A constant initial pressure $p^{0}_{up}=-2$, and a non-constant $p^{0}_{down}=-y-1/4$ are defined in the upper and in the lower part of the domain, $\Omega_{up}$ and $\Omega_{down}$. The van Genuchten parameters are presented in Table 1.

Further, for both Richards and transport equations, we have a source term, $f(x,y)=.006\cos(4/3\pi y)\sin(x)$, on $\Omega_{up}$. No external forces or sources, are defined in the lower region, i.e. $f=0$ on $\Omega_{down}$. Finally, the initial condition for the concentration is given by $c^{0}=1$ and the boundary conditions by $c_{|{\Gamma_{D}}}=4$.

We performed simulations using regular meshes, consisting of squares, whose sides were of length $dx=[1/10,\ 1/20,\ 1/40,\ 1/80]$. We consider also varying time steps of sizes $\tau=[1/10,\ 1/20,\ 1/40,\ 1/80]$. In Fig. 3(a) we are plotting the pressure and concentration profiles at the final time $T=1$. We point out that in this first example we are always in the strictly unsaturated regime, implying that Richards equation is a regular. All the proposed iterative schemes were converging for this example. In Fig. 4 is given the total number of iterations for the different schemes.

More details regarding the total number of iterations and the condition number of the linear systems are presented in Tables 2, 3. The condition number is computed at the first iteration of each algorithm and with respect to the Euclidean norm. In Table 2, we fixed a time step $\tau=1/10$ and we investigated different mesh sizes, precisely $dx=[1/10,1/20,1/40,1/80]$. In Table 3 we use a constant $dx=1/40$ and varying the time step sizes $\tau=[1/10,1/20,1/40,1/80]$. We point out that the alternate schemes are converging much faster than the classical ones. We also remark the high differences in the condition numbers, the L-scheme based algorithms being much better conditioned.

Example 1B: flow and transport in variably saturated porous media

For the second example we use the same domain, mesh sizes, boundary conditions and parameters, but we allow a saturated/unsaturated regime by changing the initial condition for the pressure. We consider a subdivision of $p^{0}$ between upper and lower regions, precisely: $p^{0}_{up}=-2$ and $p^{0}_{down}=-y+1/4$. This new expression for $p^{0}_{down}$ gives a positive pressure in the lower part of the domain (saturated region). For this example the Richards equation is now degenerate parabolic, therefore more challenging for the numerical schemes. Furthermore, we introduce this time also a reaction term $R(c)$ in the transport equation, given by $R(c)\ :=\ c/(1+c)$.

At the iteration $j+1$, the term $R(c)$ is linearized in the following way:

[TABLE]

In Fig. 10(a) we show again the pressure and concentration profiles at the final step $T=1$. The main differences to the previous example, i.e. Figure (3(a)) are in the values of the pressure. We can observe again a discontinuity in the pressure profile but, more importantly, it is evident a jump from negative to positive values. Such results were expected considering the initial pressure imposed on the domain.

In Fig. 6 are presented the total number of iterations. We remark that in this case only the L-scheme based algorithms are converging. It is also interesting to observe that the difference in the number of iterations between the more commonly used non-linear splitting approach (NonLinS) and alternate splitting (AltLinS) approach. The alternate method appears to be a valid alternative to the common formulation. It produces equally accurate results, requiring fewer iterations.

As for the previous example we present in the Tables 4, 5 the precise numbers of iterations and condition numbers for each algorithm implemented for different mesh diameters and time steps. Each segment ($-$) in the tables below, implies that the method failed to converge for such particular combination of time step and space mesh. As already observed in Fig. 6 the L-scheme based solvers are the only ones converging in all cases. Moreover, the linear systems associated with the L-scheme are better conditioned than the ones for Picard or Newton methods. We finally remark that, as expected, for smaller time steps the Newton and Picard schemes converge, see Table 5. Anyhow, the condition numbers of the systems associated to the linearized problems, remain ill-conditioned. This can cause further numerical problems.

Example 2A: well in unsaturated porous media

Our next example is inspired from [24]. We consider same domain (e.g. the unit square), boundary and initial condition and parameters as in the first numerical example (1A). The medium is again strictly unsaturated. We include now, in the upper part of the domain, a well and inject water with a specific concentration of the external component. No analytical solution is available for this case. Due to the higher complexity of the problem we use more refined meshes, precisely $dx=[1/50,1/100,1/150,1/200]$. The pressure at the well is set to $p_{W}=-10$ and the concentration of the surfactant to $c_{W}=10$.

In Fig. 7, we present the different profiles of pressure and concentration at the initial time $t^{0}$ and at final time $T=1$ day.

Once more, in Fig. 8 we compare the different solving algorithms. We study the numbers of iterations and the conditions numbers of the linearized systems. As for the first example, the media being unsaturated, the Richards equation does not degenerate and all the schemes converge. We can observe, in the Tables 6, 7, that the monolithic Newton method is the fastest, in term of numbers of iterations. We remark that the alternate splitting approach (AltLinS), once more, requires fewer iterations than the classical splitting algorithm (NonLinS) for all the linearization schemes. The linear systems resulting by applying the L-scheme based solvers are better conditioned compared with the other solvers.

Example 2B: well in variably saturated porous media

This fourth numerical example is obtained by changing the initial condition for pressure in the example 2A. We use the same $p^{0}$ as in example 1B. The profiles of pressure and concentration at the beginning and end of the simulation, i.e. at $T=1$ hour, are presented in Fig. 11. We can observe smaller changes, compared to the previous example, due to a smaller time interval (1 hour versus 1 day).

In Fig. 11 we present the total number of iterations for the different schemes applied to example 2B. Similar to the example 1B, due to the degeneracy of the Richards equation, many of the considered schemes show convergence problems. In the Tables 8,9 we study the convergence of the schemes and the condition number of the associated linear systems. The results are very similar with the previous examples, with the L-scheme based solvers being the most robust one for all the cases and with the alternate method being faster than the classical splitting schemes.

4.1 Example 3: Highly heterogeneous porous medium

Physical porous media, such underground reservoirs, are characterized by highly heterogeneous properties. Here we investigate a domain $\Omega\subset\mathbb{R}^{2}$ with highly heterogeneous porosity and permeability as presented in Fig. 12.

We consider the problem already studied in Example 1B, using the same initial conditions and parameters. We set the external forces to be equal to zero so that the flow and transport are governed only by the initial and boundary conditions. The problem investigated is highly heterogeneous and, thanks to the initial pressure $p^{0}$, also variably saturated, presenting a discontinuity in the water content. Due to the low conductivity, Fig. 12(b), a large time interval is considered, $T=10^{5}$.

We can observe, in Tables 10 and 11 and in Fig. 13, the total numbers of iterations required by each solving algorithm and the condition numbers of the associated linearized systems. For this particular problem, it is interesting to notice that, even tough we are investigating a large time domain, the L-schemes converges for a time step $\tau=T/10$. Considering the results presented in the previous examples, where for smaller time steps all the schemes converged, we tested also $\tau=[T/10^{2},T/10^{3},T/10^{4}]$. Anyhow, due to the highly complex domain, both Newton method and modified Picard did not converge. We did not test smaller time steps because the resulting number of iterations would have been much larger than the one obtained with the original $\tau$ and the L-scheme.

5 Conclusions

In this paper we considered surfactant transport in variably saturated porous media. The water flow and the transport are in this case fully coupled. Three linearization techniques were considered: the Newton method, the modified Picard and the L-scheme. Based on these, monolithic and splitting schemes were designed, analyzed and numerically tested. We conclude that the only quadratic convergent scheme is the monolithic Newton, that the L-scheme based solvers are the most robust ones and produce well-conditioned linear systems and that the alternative schemes are faster than the classical splitting approaches.

Although we recognized the existence of improved Newton solvers, e.g. [21, 26, 43, 46], we believe that, also due to its extreme simplicity, the L-scheme is a powerful tool and can be particularly useful in degenerative cases here investigated.

Acknowledgments 1

*The research of D. Illiano was funded by VISTA, a collaboration between the Norwegian Academy of Science and Letters and Equinor, project number 6367, project name: adaptive model and solver simulation of enhanced oil recovery. The research of I.S. Pop was supported by the Research Foundation-Flanders (FWO), Belgium through the Odysseus programme (project G0G1316N) and Equinor through the Akademia grant.

We thank our colleagues from the Sintef research group. In particular Olav Moyner PhD, for assistance with the implementation of the numerical examples in MRST, the toolbox based on Matlab developed at Sintef itself.*