Assessment of a non-traditional operator split algorithm for simulation of reactive transport
Introduction
The operator split (OS) technique provides a convenient means of solving the advection-diffusion-reaction (ADR) equation:where c is the solute concentration, D the diffusion coefficient, the advective velocity, x the spatial variable and t is time. The function represents kinetic processes within the system; for some applications, non-linearities in the system are confined to . For the purposes of this manuscript, we will interpret as a reaction.
The standard two-step OS scheme used to simulate (1) involves splitting the reaction part from the equation, yielding a simpler linear PDE and a linear or non-linear system of ordinary differential equations (ODE) [5], [9], [12], [20], [24], [25], [28], [30]:Traditionally, the linear advection-diffusion equation (2) is solved using a standard finite difference or finite element approach [20], [30], while the reaction equation (3) is solved using an ODE solution technique that is appropriate for the particular form of [20], [30]. We refer to this scheme as the DA-R scheme.
The use of OS schemes have become widespread because of the flexibility offered by choosing different solution techniques to suit the various terms in the conservation equation [24]. OS schemes are also advantageous because the decoupling of reaction and transport discretisations permits different reaction types to be simulated with minimal recoding [13], [14]. Obtaining numerical solutions using an OS technique can increase the computational efficiency compared to solving the full system and also provides opportunities for parallel computation [16], [19]. OS schemes are particularly attractive for coupled multi-species problems or problems with non-linear kinetics [11], [18]. Possible extensions of the traditional OS technique, such as the alternating OS scheme [20], [25], [30], the iterative OS scheme [18] and specialised schemes that combine OS with LU factorisation [9] have been proposed. This work, however, focuses on a sequential OS scheme as this is the most widely used approach [8], [14], [32], [35]. Previous investigations have shown that solving the ADR equation with the traditional two-step OS scheme can introduce an error, which is minimised by choosing sufficiently small [20], [25], [30].
Traditional OS schemes can provide several advantages over solving the full conservation equation [30]. However, if the transport part of the equation is solved with standard Eulerian techniques, then the overall accuracy is limited by the well-known grid Péclet number restriction. The grid Péclet number restriction for mass-lumped linear finite element or finite difference solutions of the advection-diffusion equation requires in order to prevent spurious oscillations and excessive numerical diffusion [8]. For advection-dominated problems, the grid Péclet number restriction can be quite severe, requiring excessively fine grids to prevent numerical errors.
In addition to the traditional DA-R scheme (2), (3), previous analyses have used a three-step splitting scheme, whereby the diffusion, advection and reaction steps are all separated (D-A-R) [13], [14], [32], [35]:The D-A-R procedure has been successfully used to simulate the ADR equation, where the advection term is typically simulated using particle tracking or front tracking techniques [13], [14], [32], [35]. These methods have been invoked to help overcome the grid Péclet number restriction of traditional Eulerian approaches. The advantage of the D-A-R method is that the original PDE is split into two linear PDEs and a linear or non-linear ODE system. Each of these steps can be solved by a solution technique appropriate to the nature of each term.
Overcoming the traditional difficulties associated with solving the ADR equation under advection-dominated conditions has been an active area of research. Investigations of alternative Eulerian algorithms that minimise numerical error for advection-dominated conditions have yielded a wide range of improved schemes. A survey of the developments of alternative difference schemes can be found in Yeh [6], Leonard [34] and Christon et al. [33]; similar developments based on weighted residual strategies are surveyed by Yeh [6] and Donea and Huerta [2].
Recent advances in solving hyperbolic conservation laws with central schemes and flux limiters provide high accuracy alternatives to standard Eulerian solutions to the advective transport equation [33]. In this work, the recent central scheme proposed by Kurganov and Tadmor (KT) [21] is used to solve the advection-reaction part of the ADR equation. The KT discretisation permits a semi-discrete form, which is explicitly integrated with a Runge–Kutta method. It is therefore convenient to simply incorporate any reaction terms into this integration process thereby temporally integrating both the advection and reaction terms simultaneously. Following this argument, we propose a “non-traditional” OS scheme, whereby the ADR equation is split:This non-traditional OS scheme (D-AR) at first appears more complicated than the traditional OS schemes as the single linear or non-linear PDE has been split into two PDEs: the first is a linear diffusion equation and the second is a linear or non-linear advection-reaction equation. The reason for splitting a single PDE into two other PDEs is that the KT algorithm requires explicit temporal integration [21]. In general, temporal integration of the diffusion term cannot be completed explicitly because of the well-known stability criteria due to the stiffness of the discretised form of the diffusion operator [1], [8]. Therefore, to maintain generality, the diffusion equation is to be solved implicitly using a standard finite difference or finite element approach.
The impetus for splitting the ADR equation according to (7), (8) originated through simulating cell migration problems. For models of cell migration by haptotaxis [23], [26] and chemotaxis [22], the advection term is non-linear. In these cases, the advective velocity is proportional to the spatial gradient of a signalling (chemotaxis) or adhesive (haptotaxis) substance. The advective velocity can be temporally and spatially variable and is determined as a part of the solution. This complication results in a non-linear advection term [22]. Under these conditions, traditional splitting (2), (3) does not convert the non-linear PDE into a linear transport PDE and linear or non-linear ODE pair. This is an important point that distinguishes the application of OS strategies to non-linear advection problems from classical linear advection problems such as those encountered in contaminant transport analyses within aquifer systems [6], [8].
Previous work by Marchant et al. [23] investigated the influence of adding diffusion to the pure haptotactic tumour invasion problem described by Perumpanani et al. [26]. Marchant et al.’s numerical simulations were limited by stability constraints owing to the explicit temporal integration of the diffusion equation within the KT algorithm. Therefore, only small amounts of diffusion were added [23]. Recent work by Landman et al. [22] successfully used an OS scheme analogous to the proposed non-traditional splitting (7), (8) to perform numerical simulations of a combined diffusion and chemotaxis cell invasion problem. Landman et al. simulated any combination of diffusion or chemotaxis without any stability constraints associated with the diffusive migration [22]. This previous work successfully investigated the interactions of diffusive and chemotactic migration, and all simulations were performed with a non-traditional splitting scheme similar to that which is to be analysed here.
The objective of this work is to explore the robustness of D-AR splitting scheme for solving the ADR equation. The explicit KT central scheme is used to solve the advection-reaction part of the conservation equation. The diffusion term is split and solved implicitly. Because the D-AR OS scheme is different from all standard OS approaches, the convergence behaviour of the proposed OS scheme is investigated. In addition, the mass balance error (MBE) incurred by the proposed D-AR scheme is investigated analytically and compared to the MBE for standard DA-R and D-A-R approaches. This analytical error analysis is particularly important as previous analytical analyses of OS errors have been restricted to standard sequential and alternating forms of traditional two-step OS schemes [20], [25], [30].
Section snippets
Implicit finite difference solution of the diffusion equation
The first step in the D-AR OS scheme is to solve (7). The diffusion equation is solved using a standard implicit finite difference approach [1], [3]. The domain is uniformly discretised with nodes spaced apart. The nodes are indexed . Temporal integration can be performed with a time-weighted scheme, which yields the discretised form of the diffusion equation at a central node:where the subscript i
Discussion and conclusion
This work has extended the applicability of the KT central scheme to be used in conjunction with an OS technique to simulate ADR systems. The algorithm presented here is likely to be of interest to researchers analysing reactive transport problems, which are ubiquitous throughout all applied science disciplines. The KT central scheme is typically implemented with an explicit method for the temporal integration of the discretised equations. Therefore, to develop a general algorithm for the
Acknowledgements
Matthew Simpson and Kerry Landman are supported by a National Health and Medical Research Council project grant ID237144. Kerry Landman also appreciates the support of the Particulate Fluids Processing Center, an Australian Research Council Special Research Center. In addition, we thank D. Andrew Barry for his thorough review comments and Anjani Kumar and for his initial comments.
References (35)
- et al.
Alternative split-operator approach for solving chemical reaction/groundwater transport models
Adv. Water Resour.
(1996) - et al.
Temporal discretisation errors in non-iterative split-operator approaches to solving chemical reaction/groundwater transport models
J. Contam. Hydrol.
(1996) - et al.
Analysis of split-operator methods for nonlinear and multispecies groundwater chemical transport models
Math. Comput. Simul.
(1997) - et al.
Comparison of split-operator methods for solving coupled chemical non-equilibrium reaction/groundwater transport models
Math. Comput. Simul.
(2000) - et al.
Operator-splitting procedures for reactive transport and comparison of mass balance errors
J. Contam. Hydrol.
(2004) - et al.
Coupling of transport and chemical processes in numerical transport models
Geoderma
(1989) - et al.
A note on splitting errors for advection reaction equations
Appl. Numer. Math.
(1995) - et al.
Critical assessment of the operator-splitting technique in solving the advection-dispersion-reaction equation 1: first-order reaction
Adv. Water Resour.
(1995) - et al.
New high-resolution central schemes for non-linear conservation laws and convection-diffusion equations
J. Comput. Phys.
(2000) - et al.
Critical assessment of the operator-splitting technique in solving the advection-dispersion-reaction equation 2: Monod-kinetics and coupled transport
Adv. Water Resour.
(1995)
A two parameter family of travelling waves with a singular barrier arising from the modelling of extracellular matrix mediated cellular invasion
Physica D
Generalized solution to multi-dimensional multi-species transport equations coupled with a first-order reaction network involving distinct retardation factors
Adv. Water Resour.
Efficient implimentation of essentially non-oscillatory shock-capturing schemes
J. Comput. Phys.
Analytical solutions for chemical transport with simultaneous adsorption, zero-order production and first-order decay
J. Hydrol.
The Mathematics of Diffusion
Finite Element Methods for Flow Problems
Applied Finite Element Analysis
Cited by (38)
Survival, extinction, and interface stability in a two-phase moving boundary model of biological invasion
2023, Physica D: Nonlinear PhenomenaPattern formation and front stability for a moving-boundary model of biological invasion and recession
2023, Physica D: Nonlinear PhenomenaThe effect of geometry on survival and extinction in a moving-boundary problem motivated by the Fisher–KPP equation
2022, Physica D: Nonlinear PhenomenaLocal meshless method for PDEs arising from models of wound healing
2017, Applied Mathematical ModellingCitation Excerpt :In general, the discretization of convection term in convection-dominated IBVPs must be tackled with care [20]. Conventional discretization approaches are frequently unsuitable for such types of problems [19] due to excessive numerical smearing and nonphysical oscillations near sharp gradients if not complemented with some stabilization procedure. Additional difficulties include strong couplings and nonlinearities of advection-dominated systems of PDEs [21].
Are in vitro estimates of cell diffusivity and cell proliferation rate sensitive to assay geometry?
2014, Journal of Theoretical BiologyModels of collective cell motion for cell populations with different aspect ratio: Diffusion, proliferation and travelling waves
2012, Physica A: Statistical Mechanics and its Applications