Numerical Simulation for a Fractal MIM Model for Solute Transport in Porous Media

A fractal mobile-immobile (MIM in short) model for solute transport in heterogeneous porous media is investigated from numerics. An implicit finite difference scheme is set forth for solving the coupled system, and stability and convergence of the scheme are proved based on the estimate of the spectral radius of the coefficient matrix. Numerical simulations with different parameters are presented to reveal the solute transport behaviors in the fractal case.


Introduction
Solute transport in porous media is a complicated process involving in physical/chemical and biological reactions with fluid mechanics, and the traditional models are the advection-dispersion equations and the MIM solute transport models. The MIM model describes the hydrodynamic behavior in the mobile zone and the mass transfer process between the mobile zone and the immobile zone, which can characterize the physical/chemical non-equilibrium of solute transport in heterogeneous porous media. Although the physical and chemical non-equilibrium models are based on different concepts, they can be described by the same mathematical equation in dimensionless form, see Pang and Close (1999), Toride, Leij and Van Genuchten (1995) for instance. A MIM solute transport undergoing linear sorption without degradations and source/sink reactions in the studied region is given as follows Meerschaert (2009), Li, Wen, Zhu andJakada (2020), Lu, Wang, Zhao and Rathore et al (2018), Schumer and Benson (2003), Van Genuchten and Wagenet (1989)): (1.1) where C 1 , C 2 are the dimensionless solute concentrations in the equilibrium and non-equilibrium sites respectively, P > 0 is the Pelect number, and R ≥ 1 is the retardation factor due to the sorption, and β ∈ (0, 1) is a partitioning coefficient between the equilibrium and non-equilibrium phases, and ω > 0 is the first-order mass transfer rate.
The system (1.1) is a classical integer-order MIM model for solute transport in porous media which has been studied and applied widely by hydrogeologists not only in laboratory but also in field tests. However, there were some researches in the last decades constantly indicated that fractional differential equations could be more suitable than those of classical models to describe non-Darcian flow or anomalous diffusion in some special environment, especially in low-permeability porous media, see Caputo and Plastino (2004), Hansbo (2001), Miller and Low (1963), Obembe, Hossain and Abu-Khamsin (2017), Raghavan (2011), Zhou, Yang andZhang (2018) for instance. The solute mass transfer or the chemical reaction in a heavy heterogeneous porous media is not an instantaneous process but a longtime dynamical behavior due to the memory effect, in which case fractional diffusion equations incorporating with the memory effect are expected to describe the anomalous diffusion processes in the porous media, see Baeumer and Meerschaert (2007), Benson, Wheatcraft and Meerschaert (2000), Gerolymatou, Vardoulakis and Hilfer (2006), Kelly and Meeschaert (2019), Zhang, Benson andReeves (2009), Zhou, Yang andZhang (2019), for instance.
The initial value condition is assumed to be which means the solute concentration in the mobile and immobile zones are both zero at the initial stage. The boundary condition at x = 0 is given as which implies that the left-hand side of the region in the mobile is a source, and keeps a constant distribution. The condition at x = 1 is impermeable, which is given by As a result we get a determined system composed by (1.2) with (1.4)-(1.6), which is a novel fractal MIM model for describing reactive solute transport processes with dynamical sub-diffusion behaviors. Its solution can be expressed by the inverse Laplace transform, but the computational cost is very expensive. The finite difference methods are usually utilized for solving fractional differential equations, see Du, Gunzburger, Lehoucq and Zhou (2012), Li, Sun, Jia and Du (2016), Liu, Zhuang, Anh, Turner and Burrage (2007), Liu, Zhuang and Burrage (2014), Meerschaert and Tadjeran (2004), for instance. As for the fractal MIM solute transport system (1.2), there are no studies not only for the solution's properties in mathematics but also for numerical methods. The motivation of this paper is to set forth a numerical method and give numerical simulations for the coupled system. By discretizing the fractional derivatives and integer-order derivatives in (1.2), an implicit finite difference scheme is established, and the unconditional stability and convergence of the scheme are proved by estimating the spectral radius of the coefficient matrix under natural conditions, and numerical simulations are performed with different parameters to reveal the subdiffusion behaviors of the solute in fractal cases.
The rest of the paper is organized as follows.
In section 2, an implicit finite difference scheme for solving the system (1.2) is put forward, and the stability and convergence of the scheme are proved in Section 3, and numerical simulations are presented in Section 4. Concluding remarks are given in section 5.

The Finite Difference Scheme
Let m, n be positive integers, and h = 1/m, τ = T/n be grid steps to discretize the domain. Denote as the approximations. By the general finite difference method as used to fractional diffusion equations (Li, Sun, Jia and Du (2016), Liu, Zhuang, Anh, Turner and Burrage (2007), Meerschaert and Tadjeran (2004)), we have (2.2) (2.3) By omitting the truncated terms we get an implicit difference equations given as and the initial boundary value conditions are discretized as By rearranging (2.4) we get the difference scheme in the matrix form: where M 11 , M 12 , M 21 and M 22 are all m − 1-order matrices given by And the matrices N and N 0 in (2.5) are all 2(m − 1)-order defined by where I is the m − 1-order identity matrix, O denotes the m − 1-order zero matrix, and and And the matrix Ψ k j is defined by for j = 1, · · · , k − 1 and k = 2, · · · , n − 1.

Stability and Convergence
According to the background of solute transport in porous media and the physical/chemical laws, the model parameters satisfy the following condition: Under the natural condition (3.1), we will prove the coefficient matrix M given by (2.6) is strictly diagonal dominant, and the finite difference scheme (2.5) is uniquely solvable, and the stability and convergence of the scheme are obtained by estimating the spectral radius of the coefficient matrix.
Proposition 3.1 For the coefficient matrix M given by (2.6), there holds with ∑ 2m−2 j=1, j i M i j < 0 for i = 1, 2, · · · , 2m − 2. Proof Obviously, the coefficients A, B, C, D, E, and F defined in (2.3) are all positive under the natural condition (3.1), and the assertion is valid by the following verifications.
we also get the assertion (3.7).
Corollary 3.1 By this proposition, the coefficient matrix M of the difference scheme (2.5) is strictly diagonal dominant, and the difference scheme is uniquely solvable.
By this proposition we also get an estimates of the spectral radius of the matrix M by direct computations.
Lemma 3.4 For the coefficient matrix M given by (2.6) there hold here ∥M∥ ∞ = max 1<i<2m−2 {M ii }, and ρ(M) is the spectral radius of M, and there exists a matrix norm ∥ · ∥ * such that where M −1 denotes the inverse matrix of M.

Stability and Convergence
In this subsection, we discuss the stability and convergence of the difference scheme (2.5) utilizing the norm ∥ · ∥ * given in Lemma 3.4, and we denote it as ∥ · ∥ for convenience of writing.
Theorem 3.1 Under the natural condition (3.1), the implicit finite difference scheme (2.5) is of unconditional stability.
Proof By the linear difference scheme (2.5), we get where U 0 denotes the initial function with noises, E k = U k − U k denotes the solutions difference at the k-th level, and k = 0, 1, · · · . By (3.11) and Lemma 3.4, there holds Suppose that there are E j < E 0 for j = 1, 2, · · · , k. Then we get also by (3.11) and Lemma 3.4 (3.13) On the other hand, for given α, γ ∈ (0, 1), by Lemmas 3.2-3.3, the sum of norms ∥N∥ ∞ + ∥N 0 ∥ ∞ + ∑ k−1 j=1 ∥Ψ k j ∥ ∞ is bounded independent of k ∈ N, and thanks to the equivalence of the norm in finite-dimensional space, there exists a positive constant C > 0 such that E k+1 ≤ C E 0 , (3.14) which implies that the assertion of this theorem is valid by the inductive principle. The proof is over.
Denote the errors in the solutions by and u i,k j , v i,k j j = 1, 2, are the exact and numerical solutions at (x i , t k ) of the problem respectively.
Theorem 3.2 The difference solution of (2.5) is convergent to the exact solution as h, τ → 0 for any finite time T < ∞, and there holds where C also denotes positive constant independent of k ∈ N. Proof Also by the difference scheme (2.5), and noting to (2.1) and (2.2), we have where R k , k = 1, 2, · · · , denote the truncated terms in the solutions' approximation satisfying , k = 1, 2, · · · , (3.18) here C > 0 also denotes a positive constant independent of k ∈ N.

Numerical Tests
In this section we give numerical simulations using the finite difference scheme (2.5) for the fractal MIM solute transport model, and we will investigate on the transport behaviors by taking different model parameters. It is noted that the computational cost is O(m 3 ) for solving a m-order matrix equation, and it becomes O(m 3 n 2 ) for the scheme (2.5) for solving the fractal MIM model due to the time-fractional derivatives, where m, n are the numbers of the discretized space and time domains, respectively.

Example 1
The model parameters are chosen as P = 1, β = 0.5, ω = 0.5, λ = 0.05, µ = 0.5, R 1 = R 2 = 2 and α = γ = 0.5. Taking the space and time steps as h = 1/20 and τ = 1/20, the finite difference scheme (2.5) is applied to solve the system, and numerical solutions are plotted in Figure 1. In order to observe the solute transport behaviors in the fractal cases, we choose h = 1/20, τ = 1/100 and α = γ = 0.5, and the system is solved numerically with different partition parameters and mass transfer rates, and the simulation results are plotted in Figures 2-3.
From Figures 2-3 it can be seen that the partition parameter and the mass transfer rate both have important influences on the solute transport process. The variation of v increases largely and u decreases largely when the partition parameter goes to large. In addition, the mass transfer rate has a great impact on the immobile zone. As the mass transfer rate decreases, the peak value of the variation of v gradually goes down. Next, we investigate on the dynamics  Figure 4, as well as the case of taking α = 1 where the solute transport in the mobile zone is assumed to be governed by the integer-order diffusion. From Figure 4 it can be seen that the solute transport becomes slow as the fractional order α goes to small, especially when the transport time t > 5. It is worth noting that the peak value also decreases as the differential order α decreases. Finally, the solute concentration at x = 0.5 during the time interval of [0, 50] are plotted in Figure 5, where we fix the fractional order in the mobile zone be α = 0.5, and take the immobile fractional orders be γ = 0.2, 0.4, 0.6 and γ = 0.8 respectively. As we can see, the immobile fractional order also plays an important role in the solute transport, and the peak value of the solute concentration in the immobile also decreases as the order γ becomes small.

Example 2
In this example the model parameters are chosen as P = 1, β = 0.5, ω = 1.5, λ = 0.05, µ = 0.1, R 1 = R 2 = 2 and α = 0.8, γ = 0.5. By using the same time-space steps as done in Ex.1, the numerical solutions are plotted in Figure 6, and the solute transport behaviors are plotted in Figures 7-8 with different partition parameters and mass transfer rates, and the long time behaviors with different fractional orders at x = 0.5 for t ∈ [0, 50] are plotted in Figures 9-10, respectively.

Conclusion
A fractal MIM model for reactive solute transport in heterogeneous porous media is investigated from numerics. A stable and convergent finite difference scheme for solving the forward problem is established, and numerical simulations under different conditions are presented. The numerical computations show that the fractal MIM model can be utilized to characterize the subdiffusion behaviors of a solute transport in a long time interval, and the fractional order, the partition parameter and the mass transfer rate all have important influences on the solute transportation. In other words, the fractal transport model is more suitable for describing a long time subdiffusion process, and the solute transport in the immobile zone can not be ignored when coping with heavily heterogeneous porous media. We are to deal with 2D fractal MIM model and focus on inverse problems related with the fractal MIM solute transport models in the near future.