An indicator-based problem reduction scheme for coupled reactive transport models.

A number of effective models have been developed for simulating chemical transport in porous media; however, when a reactive chemical problem comprises multiple species within a substantial domain for a long period of time, the computational cost can become prohibitively expensive. This issue is addressed here by proposing a new numerical procedure to reduce the number of transport equations to be solved. This new problem reduction scheme (PRS) uses a predictor-corrector approach, which "predicts" the transport of a set of non-indicator species using results from a set of indicator species before "correcting" the non-indicator concentrations using a mass balance error measure. The full chemical transport model is described along with experimental validation. The PRS is then presented together with an investigation, based on a 16-species reaction-advection-diffusion problem, which determines the range of applicability of different orders of the PRS. The results of a further study are presented, in which a set of PRS simulations is compared with those from full model predictions. The application of the scheme to the intermediate-sized problems considered in the present study showed reductions of up to 82% in CPU time, with good levels of accuracy maintained.

INTRODUCTION

The prediction of chemical solute transport behavior in porous media is of great importance in a wide range of engineering applications. To this end, a significant number of numerical transport models have been developed. Typically, these are coupled models that consider advection and hydrodynamic dispersion, the latter of which comprises both mechanical dispersion and self‐diffusion1 of the chemical species coupled with heat flow and, often, the mechanical behavior of the medium. In addition to these flow and deformation processes, the solute can also be considered to be either reactive or nonreactive depending on the nature of the problem considered.

The application of these models has varied considerably, with many studies concentrating on geochemical problems, such as modeling groundwater systems,2, 3, 4, 5, 6 assessing the performance of engineered barriers,7, 8, 9, 10, 11 or attenuation of mine water tailings.12, 13, 14 The application of these models to cementitious materials has most often investigated the ingress of chloride ions15, 16, 17, 18 or calcium leaching19, 20; however, recently, these models have also been used for the investigation of self‐healing concrete.21, 22

One of the major disadvantages of these models is that they can become computationally expensive when the chemical system becomes relatively large. Cleall et al23 suggest that the magnitude of a problem is governed by three main aspects, namely, the domain size, timescale, and the complexity of the analysis, which includes the number of variables, the degree of coupling between them, the number of processes considered, and the degree of nonlinearity of the system. The degree of nonlinearity depends on the underlying physical processes; for example, the system will be highly nonlinear if the reactions are sensitive to small changes in concentration or the permeability of the porous medium is sensitive to the degree of saturation. To deal with this issue, several techniques have been proposed to improve the computational efficiency of the associated solution process, including “operator splitting” techniques and reformulation of the coupled system. Operator splitting divides or “splits” a time step into a transport calculation and a reaction calculation, to be solved sequentially, with many models iterating between the two. This is an approximation of the time integral of the governing equation,24 which therefore changes the numerical framework from a global implicit approach to a sequential iterative or a sequential non‐iterative approach. Such calculation splitting methods have been found to reduce the computational cost25 and have been used by many authors.2, 3, 4, 8, 9, 10, 11 A disadvantage of the sequential non‐iterative approach is its propensity to introduce splitting errors (for example, mass balance errors), whereas the sequential iterative approach tends to require prohibitively small time steps and a relatively large number of iterations to achieve convergence.26, 27 Consequentially, these operator splitting approaches have been found to increase computational demand for certain “chemically difficult” cases (for example, cases with high kinetic reaction rates).27, 28

A number of authors have taken an alternative approach, concentrating instead on improving the efficiency of global implicit approaches. Most studies have focused on reformulating the system by decoupling a number of the partial differential equations (PDEs) and eliminating certain local (spatially invariant) equations, including both ordinary differential equations and algebraic equations.26, 27, 29, 30, 31 A common assumption used in such formulations is that all chemical species have the same diffusion coefficient.26, 27 This has been justified by the fact that predicted responses are relatively insensitive to differences in diffusion constants because mechanical dispersion is normally the dominant transport mechanism.27 An alternative view was expressed by Thomas et al,10 who suggested that this assumption is not valid for some chemical systems. The presence of minerals, whose concentrations can reach zero, may also cause issues when such schemes are used,26 although Kräutle and Knabner26 have introduced a moving boundary condition (BC) and a complimentary function to alleviate these difficulties.

The reduction in computational demand that can be achieved by the operator splitting and reformulation approaches described above is often limited by the chemistry of a problem; for example, only limited reductions are possible when a large number of kinetic reactions with high reaction rates are considered. To address this issue, an indicator‐based multi‐order problem reduction scheme (PRS) is proposed, which allows for greater reduction in problem size that is independent of the nature of reactions and which explicitly considers species‐dependent diffusion coefficients. The scheme works by decoupling a number of the nonlinear PDEs from the global system of equations. In this new PRS, a small number of indicator species (≤ 3) are selected for full computation, and then, on a stepwise basis, the response of other species is inferred from these indicators using a predictor‐corrector approach. The reduction in computational demand is achieved through the reduction in size of the nonlinear global system of equations and by the replacement of nonlinear PDEs with stepwise predictor‐corrector computations for non‐indicator species. Three different orders of the PRS are considered, each applicable to different chemical systems. The new PRS is presented along with a study to establish the range of applicability of each order considered. This is followed by a set of validation examples in which the reactive transport of chemical species through mortar and concrete specimens is simulated.

THEORETICAL FORMULATION

Mathematical model

The theoretical model is based on the approach of Gawin et al.32 The porous medium is assumed to be composed of three main phases, namely, solid skeleton, liquid moisture, and gas, the latter of which consists of dry air and moisture vapor. It is also assumed that dissolved chemical species are present in the liquid moisture phase and that precipitation can occur. Following the approach of Chitez and Jefferson,33 the combined gas pressure of dry air and moisture vapor is assumed to remain constant at atmospheric pressure but that the moisture vapor may diffuse through the gas phase.

The primary variables considered here are the capillary pressure (P ) and chemical concentration (c ) of each species. The advantages of using capillary pressure as the primary moisture variable have been discussed by Gawin et al.32 These include the fact that it is convenient for coupling moisture flow to mechanical behavior, due to the direct relationship that exists between the pressure and the components of stress, and that it provides a physically meaningful driving force for moisture flow.

Macroscopic balance equations, as derived from the volume‐averaging theorem and the hybrid mixture theory,34 are presented in the following sections. Heat transport is included in the model, governed by the enthalpy balance equation; however, the problems considered here are isothermal, and so, this part of the model has not been used in the present study.

Governing equations

Moisture transport

Transport in the porous medium is described using mass balance equations for liquid moisture, moisture vapor, and chemical species. Defining the domain with boundary and noting that the time interval considered is given as , the mass balance of the liquid moisture and moisture vapor can be written as Summing these up gives the total mass balance for moisture, as follows: where denotes the porosity, is the degree of liquid moisture saturation, is the liquid density, is the degree of gas saturation, and is the moisture vapor density. denotes a time derivative, is the liquid moisture velocity, is the rate of moisture transfer from liquid to vapor, and J is the moisture vapor diffusion described here by Fick's law.32 That is, where is the gas density and is the moisture vapor diffusivity tensor given by22 where A  = 1, B  = 1.667, and f  = 0.01. D  = 2.47 × 10−5 m2/s is the moisture vapor diffusion coefficient in air, T  = 273 K is the reference temperature, and P and P atm denote the gas pressure and atmospheric pressure, respectively.

The flow of the liquid moisture through the medium can be described by Darcy's law.32 That is, where K in is the intrinsic permeability tensor of the medium, μ is the viscosity of the fluid, and K rw is the relative permeability of the moisture phase given by8 where A  = 4 in the present work.

The remaining constitutive relationships used in the moisture transport component of the model are given in the Appendix (see Table A1; see also the Nomenclature table and the notation on the Greek letters and subscripts/superscripts used in this paper).

Table A1

Constitutive relationships used in the moisture transport component

Constitutive relation	Expression	Values
Kelvin's law22	P_c=−ρ_wRTM_wlnP_vP_vs	Mw = 18 kg/kmol
		R = 8314.5 J/K kmol
Dalton's law22	P_g  = Pda+Pv	–
Antoine's law22	ρ_vs=b_1·10^b_2−b_3b_4+T−b_5	b1 = 133.32, b2 = 8.07,
		b3 = 1730.63, b4 = 233.43
		b5 = 273
Moisture retention22, 46	S_w=P_Ca^11−m+1^−m	a = 1.867 × 107 Pa
		m = 0.44

Chemical transport

The mass balance of a dissolved chemical species, i, can be written as where the final term on the left‐hand side is the source/sink (SS) term due to chemical reactions, S is the degree of saturation of precipitated or sorbed material, is the mass of the precipitated or sorbed material, and is the dispersive flux given by the Poisson‐Nernst‐Planck equations, which can be written as Equation (9),35 where the first term on the right‐hand side represents mechanical dispersion and the remaining right‐hand side terms account for diffusion.

In the above equations, is the coefficient of molecular diffusion, z is the charge of an ion, ε is the dielectric permittivity, is the charge density, F is Faraday's constant, and is the electrical potential. is the mechanical dispersion tensor, defined in Equation (14). The second term in the bracket of Equation (9) represents dispersion due to the local electric field , which can be calculated from Equation (10). It can be noted that the first term in Equation (10) is negligible,35, 36 as is the charge density ; thus, Equation (10) reduces to the following charge neutrality condition10, 35, 36, 37: and since the pore solution is initially charge neutral, this condition is ensured through the “no electrical current” condition,10, 36, 37 as follows: where is the diffusive flux of species i.

Equation (12) can be used to eliminate the electrical potential gradient, giving the dispersive flux of species i as

Hydrodynamic dispersion is the sum of molecular diffusion and mechanical dispersion. The definition of hydrodynamic dispersion of Bear and Bachmat1 is adopted here, where the mechanical dispersion tensor is given as where and are the longitudinal and transverse dispersivities, respectively, and is the Kronecker delta ( if k = j or ).

Chemical reactions

The chemical reactions considered here are all assumed to be nonequilibrium reactions, ie, transport of the ions may be faster than the rate of reactions such that binding or precipitation is not instantaneous.38 In addition, it is assumed that the rates of the reaction can be computed using a Freundlich‐type isotherm, as follows35: in which μ and λ are rate parameters and τ is a characteristic time that accounts for nonequilibrium behavior. It can be noted that molar concentrations are used to calculate the reaction rates throughout.

Boundary conditions

In order to solve the system, both the initial conditions and BCs are required. The initial conditions considered here define the values of all variables at time t = t 0 throughout the domain and on the boundary, as follows:

The BCs considered here are of the Cauchy type and the Dirichlet type. Those of the former type are a combination of imposed fluxes (which, in isolation, describe Neumann BCs)32 and convective fluxes, which are a function of the difference in variables between the sample and the environment. Applied to the governing mass balance equations of moisture and of dissolved chemical species, the Cauchy‐type BCs describe the rate of mass transfer from the environment to the sample and are given as where q , q , and q are the prescribed fluxes of the moisture, moisture vapor, and chemical species, respectively, and β and γ are the convective transfer coefficients of the moisture and chemical species, respectively. is the unit normal vector to the boundary.

The Dirichlet‐type BCs fix the value of the variable on the boundary and are given by

NUMERICAL SOLUTION

In the present study, governing Equations (3) and (8) are discretized using the finite element method. The resulting variational problem, after the application of the Gauss‐Green divergence theorem, may be written as follows: find , such that where is the vector of primary variables, and the space for the trial functions is defined as .

In this study, the continuous Galerkin weighted‐residual method32 is employed such that the weight functions (W) are chosen to be equal to the shape functions (N), with the primary variables being the capillary pressure and the species' concentrations. The resulting system of discretized equations is as follows: where the superior dot denotes a time derivative and each variable (eg, P ) is interpolated from the vector of nodal variables in the standard manner (ie, ).

It is recognized that the continuous Galerkin method does not guarantee local mass conservation when applied to Equations (21) and (22), and the system described by Equation (23) may be subject to spurious oscillations,39 particularly for an advection‐dominant case. To address this issue, a number of stabilization techniques may be used, including mass lumping, the streamline‐upwind Petrov‐Galerkin method, and enrichment of the Galerkin method.39, 40, 41 Alternatively, a different method could be employed for spatial discretization, such as the discontinuous Galerkin method.39 In the present work, however, no spurious oscillations were observed in the problems considered, and therefore, it did not prove necessary to employ any of the above stabilization techniques.

The global matrices are given by where superscript “ne” is the number of elements and the element matrices are given as follows:

The global system can be written in compact form as

Applying an implicit Euler backward difference scheme33 for time discretization leads to

This set of nonlinear equations can then be solved using a standard Newton‐Raphson procedure33 based on a first‐order Taylor series expansion of the mass/energy balance error, which leads to the following incremental‐iterative update of the primary variable vector ( ): where is the approximation error, given here as

Without loss of generality, bilinear quadrilateral elements were used throughout this study.

Since convergence is not always guaranteed with the Newton‐Raphson procedure,39 the stability of a solution was checked using the following Courant‐Friedrichs‐Lewy condition, as suggested by Zhu et al12:

PROBLEM REDUCTION SCHEME

The balance equations governing reactive transport in porous media (Equation (27)) are often highly coupled and nonlinear, and, as such, the computational demand associated with their solution can become prohibitively expensive, particularly when solutions are required for large domains and relatively long time periods. The present approach addresses this issue by decoupling a number of the nonlinear PDEs from the global system of equations.

The proposed PRS is a predictor‐corrector approach that employs one or more indicator species and then, in the “predictor” step, computes the transport of the other (non‐indicator) species by interpolation. The “corrector” step then refines the predicted concentrations using error approximation. Reactions involving both indicator and non‐indicator species are dealt with on a point‐wise basis (at nodal points) at the end of each time step. Before describing the interpolations and detailed processes used in the PRS, the overall solution algorithm is presented to show where the PRS fits into the transient solution procedure.

Algorithm

Formulation

The aim of the PRS is to greatly reduce the computational cost of solving multispecies chemical transport problems by reducing the number of primary variables solved in the coupled system. In this scheme, a full solution is undertaken for the capillary pressure and a selected number of indicator species, and then, a predictor‐corrector approach is used to compute the concentration of non‐indicator species. The predictor first computes the concentration of non‐indicator species from the concentration of the indicators, and then, a correction is applied to these values using error approximation derived from an appropriate balance equation.

Predictor step

The concentrations of the indicator species at each time step are computed from the solution of the coupled equation system (Equation (27) and Algorithm 1). Once these indicator concentrations are known, the predictor step is applied to compute the transport of non‐indicator species. To derive this predictor step, the governing mass balance equation for a conservative chemical species i is considered, which, neglecting the charge neutrality condition, is given as where is a function of the liquid moisture velocity and the dispersivities (see Equation (14)).

Noting that both and are the same for all species, it can be seen from Equation (32) that, for a conservative chemical species, the difference between the rates of transport of chemical species depends on their diffusion coefficients, concentrations, and concentration gradients.

In the present study, Lagrangian polynomial interpolation is used for diffusion, with the interpolated diffusion being used to weight the change in concentration of an indicator species over a time step, Δ, based on the relative difference between the known diffusion coefficients of the current non‐indicator species and of the indicator species and , as follows: in which ( is a diffusion interpolation function, which is incorporated into the final interpolation function below (Equation (34)), n ind is equal to the number of indicator species, subscripts “ind” and “jnd” denote specific indicator species, subscript i denotes non‐indicator species, and and represent the diffusion coefficient of any species and a specific indicator, respectively. It can be noted that, in the case of n ind = 1, Equation (34) simplifies to (.

Equation (33) does not account for the effect of the difference in concentration gradients between indicator and interpolated non‐indicator species. This is addressed by introducing a concentration gradient function ( f ), and adding terms to account for chemical reactions and the charge neutrality condition leads to Equation (34) for predicting the concentration of a non‐indicator species at time , as follows: where c is the concentration of a chemical species, refers to a concentration gradient, and superscript t denotes time.

The penultimate term in Equation (34) represents the SS term due to chemical reactions, and the final term ( f) represents the diffusion due to the charge neutrality condition, which is elaborated in Section 4.3. In the scenario in which a reaction is “sufficiently fast” in comparison with the rates of the transport processes, this may be treated under the local equilibrium assumption.38 As such, the change in concentration over a time step should be multiplied by a retardation factor, as detailed in the work of Fetter42 (ie, Equation (34) becomes , where denotes the retardation factor for an equilibrium reaction and the SS term has been retained to account for any kinetic reactions).

A number of options were considered for the concentration gradient function in Equation (34) (  ), including using the ratio of non‐indicator‐to‐indicator gradients from the previous time step and using a weighted measure of concentrations from the sources. Both these options proved to be flawed and to have problem‐dependent levels of accuracy. A more successful approach was to derive from a “source influence solution” (SIS), which is a linear steady‐state solution of the diffusion problem for each chemical species that is carried out at the beginning of the solution procedure (thereby having a negligible effect on the overall solution time) and is given as noting that , , and are constants, for a uniformly saturated medium Equation (35) simplifies to the following Poisson equation:

The outputs from the SISs are concentrations of each chemical species throughout the domain, which are subsequently employed in f , as follows:

Corrector step

The corrector refines the concentrations computed from the predictor step using the mass balance approximation error (Equation (30)) for each non‐indicator species in turn. The concentration correction for each species i is then based on a first‐order Taylor series expansion of a single species form of Equation (29), as follows: where is the concentration correction vector for species i and is an approximate tangent coefficient matrix.

Two forms of were considered as follows:

Option 1: ;

Option 2: , in which is the lumped diagonal storage matrix and is a diagonal form of the flux matrix.

The use of diagonal lumped matrices implies that the concentrations can be updated on a point‐wise basis. The concentrations are refined as

It was found that the two options gave similar results in terms of the accuracy of solutions, and therefore, the simpler option 2 was adopted for all subsequent computations. In this study, the corrector step was treated as non‐iterative since it was found that sufficient accuracy could be obtained using a single corrector step, as illustrated in Sections 6 and 7 of this paper.

Summary of the predictor‐corrector scheme

In the present work, three different orders of the generalized reduction scheme are investigated, denoted PRS0, PRS1, and PRS2, which use one, two, and three indicator species, respectively. Recalling Equations (34) and (39), the predictor‐corrector scheme is summarized as

BCs and charge neutrality condition

The initial conditions for the PRS can be defined in the same way as in the full model (ie, using Equation (16) for all species). The BCs for the PRSs are also defined in the same way as for the full model and should be of the same type for both “indicator” and “non‐indicator” species, with the latter being calculated on a point‐wise basis.

Another key consideration is how the charge neutrality condition is satisfied in the PRSs, since the transport of non‐indicator species is not calculated in the reduced system of governing equations. In the present approach, the diffusive flux due to the electric field is considered explicitly by using concentrations from the previous time step and moving these to the right‐hand side of the governing equations in a similar manner to the way moisture flow under gravity is included in liquid transport computations. For non‐indicator species, this can then be calculated on a point‐wise basis and subtracted at the end of a PRS predictor step (Equation (40)). The diffusive flux due to the electric field is therefore given by

Selection of indicator species

The choice of indicator species is an important aspect of the PRS. For PRS1 and PRS2, species with the highest and lowest diffusion coefficients are always defined as indicators. This ensures that the computed responses of all non‐indicator species are bounded by those of fully computed species. For PRS2, which requires a third indicator species, it was found, following an error analysis and sensitivity study (see Section 6.2), that the greatest accuracy was achieved by using a third indicator species with a mean diffusion coefficient. In cases where such a real species is not available, artificial species can be used. Artificial indicator species may also be used for other reasons; for example, as mentioned above, PRS1 and PRS2 use indicator species with the highest and lowest diffusion coefficients. However, the presence of any reactions associated with one of the bounding species, with reaction rates of the same order as (or higher than) those of the transport processes, may alter its rate of transport, such that it may no longer be the highest/lowest. In such situations, a nonreactive artificial indicator should be used in its place, thereby maintaining the solution bound. The criterion used to select the single indicator species required for PRS0 is given in Section 6.2

Another aspect of reactive transport problems, which is relevant to the choice of the indicator species, is that the transport of different species can take place over very different timescales. When this occurs due to a dominant (or “sufficiently fast”38) reaction linked to a particular non‐indicator species, this is dealt with by modeling the process as an equilibrium reaction. If this is due to other factors, then indicator species could be chosen to represent each different timescale.

It is recognized that other factors may also affect the choice of indicators; however, the factors considered in the present work encompass a wide range of real problems.

The number of indicator species, and associated order of the scheme, is also a key consideration when employing the PRS. The choice of the number of indicator species to use depends upon the problem under consideration and specifically upon the range of diffusion coefficients of the chemical species in the system. In general, the smaller the range, the lower the order of the scheme required, with PRS0 being applicable to problems in which the diffusion coefficients lie in a narrow range, as quantified in Section 6.2. An exception to this rule is problems with a significant degree of advection, for which PRS0 is not appropriate (unless the diffusion coefficients are equal for all species), since solutions for the non‐indicator species transport are not bounded, which can lead to an overestimation of the advection.

The specific selection criterion used to determine the indicators for each order of the PRS is given in Section 6.2, along with the expected accuracy for each degree of the scheme.

VALIDATION OF THE FULL MODEL

Before investigating the behavior of the PRS, the validity of the full model is demonstrated. To this end, a non–steady‐state diffusion problem reported by Baroghel‐Bouny et al35 is considered. This problem is based on experiments carried out by Francy43 on cement discs of dimension 120(d) mm × 20(h) mm. In these tests, the left‐hand side of the sample was exposed to a salt solution, whereas the remaining sides were sealed. The transport of Na+, OH−, K+, and Cl− ions was considered, and nonequilibrium chloride binding was also taken into account, using Equation (15), in addition to the instantaneous formation of Friedel's salt. It is assumed that, since chloride ions sorb onto the solid mass, hydroxide ions are released to preserve charge neutrality. At the end of the test, measurements were taken of the free chloride and total chloride content (Tcc) of the sample at different locations.35, 43 The problem setup can be seen in Figure 1.

Figure 1

Finite element mesh and problem geometry (not to scale) [Colour figure can be viewed at http://wileyonlinelibrary.com]

The time period considered was 7 days, and, following a mesh and time‐step convergence study, a time step of Δt = 3.6 s was selected, giving a total number of time steps of 168 000, along with a mesh of 20 bilinear quadrilateral elements with a maximum element size of 1 mm. To reflect experimental conditions, the sample was assumed to be initially saturated. The model parameters, BCs, and diffusion coefficients of the chemical species are given in Table 1. It should be noted that, in this example, it was found that a tortuosity factor, D , was needed to correctly predict the chemical transport; this factor takes into account the tortuous pathways of the medium and is simply multiplied by the species diffusion coefficients.

Table 1

Model parameters

Parameter	Value	Species	Initial conc (kg/kg)	Boundary conc	D mol
				(kg/kg)	(10−11 m2/s)
n	0.13	Na+	0.000299	0.01352	1.33
γ c (kg/m2s)	6.5 × 10−3	OH−	0.001105	0.00188	5.3
D τ	0.5	K+	0.002028	0.00319	1.96
μ	2.61	Cl−	0.0	0.01954	2.1
λ	0.61
τ (s)	36 000
p C 0 (kN/m2)	1.34 × 103

The concentration profile and the Tcc, as predicted by the full model, are compared with the corresponding experimental results in Figure 2. The numerical results are considered sufficiently close to the experimental results to validate the full model. The CPU time for the simulation was 1382 seconds.

Figure 2

Cl− concentrations and total chloride content (Tcc) profiles as predicted by the full model at t = 7 days (where Tcc is measured in kg/m3 of mortar) [Colour figure can be viewed at http://wileyonlinelibrary.com]

INVESTIGATION OF THE RANGE OF APPLICABILITY

The accuracy and range of applicability of the PRS with three different orders of the scheme (namely, PRS0, PRS1, and PRS2) are studied by considering a wick action test on a mortar sample in which the transport of 16 chemical species is simulated. The analysis undertaken in this study considered chemical reactions between the ions and the cement matrix, as well as advective and dispersive transport. It was decided to use an artificial set of chemical species for this study in order to have control over the range and spread of the diffusion coefficients considered.

The reactions concerned the adsorption of the chemical ions onto the cement matrix, described by the nonequilibrium Freundlich isotherm (Equation (15)). The time period considered was 24 hours, and the initial concentration of each ion, as well as the sorbed chemical mass for each species in the sample, was assumed to be zero. The mortar sample was assumed to be initially saturated, prior to the left‐hand side of the specimen being exposed to the chemical solution and the right‐hand side being exposed to an environmental humidity of 60%, with all remaining sides being sealed, thereby ensuring one‐dimensional transport. A nonuniform mesh of 25 bilinear quadrilateral elements was used along the length of the specimen, with a maximum element size of 4 mm. A time step of Δt = 36 s was used, giving a total number of time steps of 2400. As with the previous validation case, the mesh and time step were adopted following a convergence study.

The diffusion coefficients for all species can be seen in Table 2, whereas model parameters, including the boundary mass transfer coefficients for all chemical species, are given in Table 3.

Table 2

Diffusion coefficients

Species	D mol	Species	D mol
	(10−10 m2/s)	(cont'd)	(10−10 m2/s)
1	0.25	9	6
2	0.5	10	7
3	1	11	8
4	1.5	12	9
5	2	13	10
6	3	14	12
7	4	15	14
8	5	16	16

Table 3

Model parameters

Parameter	Value	Parameter	Value
μ	16	K in (m2)	35 × 10−21
λ	2.0	c b a (kg/kg)	1 × 10−3
τ (s)	2000	z a (−)	1
n	0.13	c 0 a (kg/kg)	0.0
β c (m/s)	2.5 × 10−3	p C 0 (kN/m2)	1.34 × 103
γ c (kg/m2s)	1 × 10−4

For all species.

Full model results

To determine the accuracy of the reduction schemes, an analysis of the full problem was undertaken. The calculated profiles of chemical concentration, moisture content, and sorbed masses from this full analysis at time t = 24 h are given in Figure 3. The responses are plotted in groups of species in order of increasing diffusion coefficient (D mol). The relative difference between the four responses in each of the groups reflects the relative spread of diffusion coefficients across the group, with group 1 (species 1‐4) representing a 6‐fold increase in diffusion coefficient and group 4 having a relative increase of only 1.6. The moisture profile shows that some drying of the specimen has occurred over the 24 hour analysis period, but much of the sample remains saturated or near saturated. The sorbed mass profiles are very similar to the concentration profiles but with less penetration into the sample. The CPU time for the simulation was 312 seconds.

Figure 3

Normalized concentration, sorbed mass profiles, and saturation profile as predicted by the full model at t = 24 h

PRS results

In this section, the results from a set of analyses of the problem considered in Section 6.1, undertaken with the three orders of the PRS, are compared with the results from the full analysis described above. The indicator species chosen for each of the solutions are presented in Table 4. An artificial species labeled “A” has been chosen for PRS0 and PRS2 in order to allow the use of an indicator with a diffusion coefficient corresponding to the mean value of the species diffusion coefficients.

Table 4

Indicator species chosen and the corresponding diffusion coefficients (10−10 m2/s)

PRS	Indicator 1	Indicator 2	Indicator 3
0	A (6.2)	–	–
1	3 (1)	9 (6)	–
2	3 (1)	A (7.03)	16 (16)

Abbreviation: PRS, problem reduction scheme.

This problem is used to determine the range of applicability of each of the reduction schemes based on a maximum allowable error in any one chemical species at any one time. The tolerance measure considered is the mean absolute percentage error (MAPE), defined as

An acceptable MAPE of 3% was selected for PRS0, and 1% was selected for both PRS1 and PRS2 to enable the definition of the range of applicability of each scheme. It is acknowledged that these tolerance values are a matter of judgement and that different problems will require different levels of accuracy. Here, 1% is considered to be a reasonable tolerance value for most practical purposes, and 3% is thought to be an acceptable error for the type of relatively coarse approximation associated with PRS0.

For PRS0, it is possible to consider species with progressively larger or smaller diffusion coefficients until a profile exceeds the tolerance; however, this is not possible for PRS1 and PRS2 as they use the extremes of the diffusion coefficient range as indicator species. To overcome this issue, the range was successively varied on a trial‐and‐error basis until the maximum range that meets the tolerance was found.

The resulting ranges, over which each PRS achieves the selected tolerance, are given in Table 5 (where superscripts u and l indicate upper and lower indicators, respectively). The intermediate indicator () for PRS2 was chosen to be the mean diffusion coefficient, which was based on the results of a series of analyses aimed at finding the value of that produced the most accurate solutions.

Table 5

Diffusion coefficient ranges over which each scheme is applicable

PRS	0	1	2
Diffusion coefficient	0.8Dind<Di<1.5Dind	Dind,u<6Dind,l	Dind,u<16Dind,l a
range

Abbreviation: PRS, problem reduction scheme.

The maximum mean absolute percentage error was 1.13%, but this was deemed close enough to the 1% tolerance; the average was only 0.63%.

To further explore the influence of changing the diffusion coefficient of an indicator species, a sensitivity analysis was undertaken on the value of the diffusion constant used for the intermediate indicator for PRS2. Using the problem considered in Section 6.1, it was found that changing this coefficient by 20% resulted in an increase in the MAPE of 0.13% and the mass balance error of 0.17%, which suggests that the scheme is relatively insensitive to the choice of the intermediate indicator.

Profiles showing the two species with the largest departure from the results of the full analysis, within the range of applicability given in Table 5, are presented for each PRS in Figure 4. The relative percentage error plots for each PRS can be seen in Figure 5, where the relative percentage error is given as . It can be seen from the profiles that PRS results are very close to those of the full model, particularly for PRS1 and PRS2. The profiles corresponding to the PRS without the corrector step show a larger departure from the results of the full analysis; however, the profiles are still relatively close to one another. This shows that, for this particular problem, the corrector has a limited effect, although, as is shown in Section 7.1, this is not the case for many other problems.

Figure 4

Worst case profiles in the applicable range for (A) PRS0 (species 12 and 13), (B) PRS1 (species 5 and 6), and (C) PRS2 (species 4 and 5) at t = 24 h [Colour figure can be viewed at http://wileyonlinelibrary.com]

Figure 5

Relative error profiles for (A) PRS0, (B) PRS1, and (C) PRS2

In addition to the relative error in concentrations, another key consideration is the mass balance error. To assess this, the total accumulated mass balance error was calculated by comparing the mass of each chemical in the system obtained from the PRS analysis with the results from the corresponding full analysis, with the latter being taken as the correct solution for this purpose. The average relative mass balance errors were 1.14% for PRS0, 4.24% for PRS1, and 1.98% for PRS2. In this case, the mass balance error for PRS0 is the smallest, even though PRS0 predicts the concentrations less accurately. This is because PRS0 has areas of overprediction and underprediction, leading to a small overall error in mass balance.

The results of this investigation led to the PRS selection criteria for the PRS given in Tables 6 and 7. These recommendations are based on the range of diffusion coefficients for each scheme, which was found to give an acceptable MAPE, and the corresponding indicator species selected. The selection of the order of the PRS should also take into account the fact that PRS0 is an inappropriate choice for problems with a significant degree of advection, as discussed in Section 4.4. Clearly, it would be possible to use higher‐order PRSs, but the range of orders examined here (ie, 0‐2) is considered to cover the majority of practical examples.

Table 6

Selection of order of the problem reduction scheme (PRS)

Diffusion coefficient rangea	PRS
Du<1.9Dl	0b
Du<6.0Dl	1
Du<16Dl	2

Where and indicate the upper and lower diffusion coefficients, respectively.

Where advection is significant, PRS1 is the appropriate choice.

Table 7

Selection of indicator species for each order of the problem reduction scheme (PRS)

PRS	Indicator numbera
	1	2	3
0	0.67Du b	–	–
1	Dl	Du	–
2	Dl	Dmean	Du

Where the selected indicator species should have diffusion coefficients corresponding to the tabulated values; in the case that no species exists with such a diffusion coefficient or the species is highly reactive, an artificial species with the corresponding diffusion coefficient should be used.

This recommendation is based on the fact that PRS0 was found to have an upper limit to its applicable diffusion coefficient range of 1.5 (see Table 5).

Convergence

The above problem (described in Section 6.1) was also used to examine the convergence properties of the model with respect to mesh refinement (ie, h‐refinement) and to the order of the PRS (or PRS order). The following ‐error norm was employed in these convergence studies: where is the error vector and is the vector of reference concentrations.

In the absence of an analytical solution to this problem, the reference solution for the h‐convergence study was obtained by using a very fine mesh (h = 3.125 × 10−4 mm); for example, the reference solution for PRS0 is the PRS0 solution obtained using this very fine mesh.

The reference solution for PRS order convergence was the solution from the full model. The results of the investigation are given in Figure 6. It is noted that uniform meshes were used for all of the analyses in this convergence study and that the mesh size for PRS order convergence was h = 6.25 × 10−4 mm.

Figure 6

Mesh size and problem reduction scheme (PRS) order convergence study results

The results presented in Figure 6 show that the full model and all PRSs exhibit satisfactory h‐convergence. It noticeable that the order of h‐convergence is greatest for PRS0 but least for the full solution; however, this is believed to relate to the fact that the reference solutions for each set of results are those obtained from a fine mesh with the order of the scheme being considered. The “error” measures are therefore relative and not absolute.

The results also show the convergence with respect to PRS order. It is interesting to note that the relative errors of the PRS solutions are similar to those of the reduced‐order model of the reaction‐advection‐diffusion equation found in the work of McLaughlin et al44 and of the solute transport problem found in the work of Luo et al.45 In addition to this, the convergence with respect to PRS order shows similarities to the convergence of the aforementioned reduced‐order models with respect to the number of basis functions employed.

VALIDATION OF THE PRS

Having established a range of applicability for each PRS order, it is considered desirable to validate the final procedure using data reported by other authors. To this end, a problem for each reduction order has been considered in addition to an advection‐dominant problem for PRS1 and PRS2. The first two problems consider the diffusion and reactions of chemical species in a mortar sample and are based on alternative numerical solutions presented by previous authors.35, 36 The third example is a hypothetical scenario and considers two‐dimensional diffusive‐reactive transport in a mortar specimen based on the example in the work of Zhu et al.12 The advection‐dominant problem is an extension of the advective example reported in the work of Baroghel‐Bouny et al.35

As before, the meshes and time steps used for all of the analyses reported in this section resulted from convergence studies.

PRS0

The diffusion case reported by Baroghel‐Bouny et al,35 which was used to validate the full model, is considered here over a time period of 12 hours, using the same mesh and time step as reported in Section 5. The total number of time steps was 12 000. The sample was assumed to be initially saturated.

The narrow range of the diffusion coefficients of the dominant species in this problem suggests that PRS0 is applicable. The chosen indicator species was Na+, and the diffusion coefficients for non‐indicator species K+ and Cl− lie within the established range of validity of PRS0; however, the other species, OH−, lies outside the nominal applicability range, but since there is little transport of OH−, this was deemed acceptable in this case. The model parameters, BCs, and diffusion coefficients of the chemical species can be seen in Table 1, with the one difference being that a value of τ = 360 s was used for chloride binding.

The concentration profiles as predicted by the full model and PRS0 can be seen in Figure 7, along with the Tcc.

Figure 7

Chemical concentrations and total chloride content (Tcc) profiles as predicted by the full model and PRS0 at t = 12 h

It can be seen from the profiles that the PRS results show a very close match to those of the full model, with the exception of the OH− profile, which involves a relatively small overall change in concentration. The reason for the greater difference observed in the predicted OH− profiles is that the diffusion coefficient lies outside the applicable range of the scheme; as such, the transport is being overpredicted by PRS0, and the local peak of OH− ions dissipates more quickly than in the full model simulations. The average relative mass balance error in this analysis is 2.05%.

The Cl− concentration profiles predicted by the full model and PRS0 with and without the corrector step are shown in Figure 8. It may be seen from the profiles that, in this example, the corrector step has a significant effect on the accuracy of the solution, with the profile predicted by PRS0 without the corrector showing a much greater departure from that predicted by the full model. In addition to this, neglecting the corrector step increases the average relative mass balance error in this analysis to 8.31%. The CPU times using the full model and PRS0 were 101 and 62 seconds, respectively.

Figure 8

Cl− concentration profiles as predicted by the full model, PRS0, and PRS0 without corrector at t = 12 h [Colour figure can be viewed at http://wileyonlinelibrary.com]

PRS1

Diffusion experiments, as well as accompanying numerical simulations, presented by Song et al36 are considered here. In these tests, a concrete slab was cured for 90 days, after which a series of 100(d) mm × 50(h) mm cylindrical cores were taken from the slab. The sides and bottom of the specimens were then sealed, and the remaining surface was exposed to a salt solution for 6 months. The simulation considers the transport of six chemical species (OH−, K+, Na+, Cl−, SO4 2−, and Ca2+). The reaction rates for the chemical species were considered as given by where k and k represent the adsorption and desorption rates, respectively. The Freundlich‐type isotherm was used, and the adsorption was calculated based on a nonequilibrium assumption. The SS term for the mass balance equation (Equation (8)) for each species is given as the sum of the SS terms due to each relevant reaction (eg, the SS for Ca2+ is obtained from ). Song et al36 also considered the reaction of the cement matrix minerals; however, the associated changes were negligible, and these reactions were therefore not included in the present computations.

The model parameters are given in Table 8. The BCs are given in Table 9 along with the diffusion coefficients and the D es factors that account for electrostatic double‐layer effects.36 The sample was assumed to be initially saturated, and the time period considered was 2 months. The time‐step size was Δ, giving a total number of time steps of 144 000, and a mesh of 52 bilinear quadrilateral elements was used with an element size of 2 mm. PRS1 was applicable since the diffusion coefficient range satisfies the condition (); the OH− and K+ species were chosen as the indicator species in this case as they are the species with the highest and lowest effective diffusion coefficients, respectively (where the effective diffusion coefficient is given as the product of D mol and D es).

Table 8

Model parameters

Parameter	Value
n	0.13
γ c (kg/m2s)	1 × 10−4
p C 0 (kN/m2)	1.34 × 103

Table 9

Boundary conditions, initial conditions, and diffusion coefficients of chemical species

Species	Initial	Bound	D mol	D es	Eq	k a (10−8)	k d (10−8)	λ
	conc	conc	(10−9 m2/s)
	(kg/kg)	(kg/kg)
Na+	0.001978	0.0	1.33	0.25	r1	18.0	12.0	0.35
OH−	0.004573	0.0	5.3	0.25	r2	24.0	144.0	0.35
K+	0.007215	0.0	1.96	0.0875	r3	330.0	1375.0	0.35
Cl−	0.0	0.01775	2.1	0.25	r4	1.275	6.42	0.2
SO4 2−	0.000192	0.0	1.07	1.0	r5	0.75	5.4	0.2
Ca2+	0.00004	0.01	0.79	0.4	r6	1000.00	1200.0	0.35

The concentration and sorbed mass profiles as predicted by the full model and PRS1 can be seen in Figure 9. The results of the PRS solution are a close match to those of the full model. The largest difference can be seen in the Cl− profile, which is due to the fact that Cl− ions are very reactive in this case, being involved in four of the six reactions. The average relative mass balance error in this analysis is 2.23%. The CPU times using the full model and PRS1 were 409 and 198 seconds, respectively.

Figure 9

Chemical concentrations and sorbed mass profiles as predicted by the full model and PRS1 at t = 2 months (where sorbed mass concentrations are measured in wt% of concrete, and sorbed masses containing Cl− are measured in Cl− concentration of the sorbed mass in wt% of concrete)

PRS2

The third simulation is based on a study by Zhu et al.12 This involves the transport of 10 chemical species and the precipitation of six minerals and considers one immobile solid species. This is a two‐dimensional problem with a point and a line source for the chemical species (marked A and B, respectively, in Figure 10; point A corresponds to the coordinates (0, 0), whereas the two ends of line B correspond to the coordinates (15, 4) and (15, 8), respectively). The geometry of the problem is given in Figure 10. The BCs impose a zero flux on all sides. The time period considered was 24 hours, Δt = 36 s, giving a total number of time steps of 2400, and a mesh of 1000 bilinear quadrilateral elements was used with an element size of 0.5 mm (as shown in Figure 10). It was assumed that the sample was initially saturated.

Figure 10

Problem geometry (not to scale) [Colour figure can be viewed at http://wileyonlinelibrary.com]

The Freundlich‐type isotherm (Equation (15)) was used, and the reactions were calculated based on a nonequilibrium assumption. The reaction rates for the six solid minerals considered are given as

The model parameters are given in Table 10. The BCs, initial conditions, reaction rates, and diffusion coefficients are given in Table 11. PRS2 was chosen to model this example since the diffusion coefficient range is greater than the range of applicability of both PRS0 and PRS1. H+, Al3+, and K+ were chosen as indicator species.

Table 10

Model parameters

Parameter	Value
n	0.3
γ c (kg/m2s)	1 × 10−4
p C 0 (kN/m2)	1.34 × 103

Table 11

Boundary conditions, initial conditions, and diffusion coefficients of chemical species

Species	Initial	Boundary	Boundary	D mol	Eq	k a	k d	λ
	conc	conc A	conc B	(10−9 m2/s)		(10−7)	(10−8)
	(kg/kg)	(kg/kg)	(kg/kg)
H+	0.000028	0.00005	0.00010	9.311	r 1	2.6	29.6	0.61
Ca2+	0.0003164	0.00048	0.00096	0.792	r 2	0.6	8.67	0.2
Mg2+	0.0010230	0.001944	0.003888	0.706	r 3	1.4	9.0	0.07
HCO3 −	0.00061	0.000030	0.000060	1.185	r 4	2.8	2.0	0.11
Al3+	0.000837	0.00135	0.00270	0.541	r 5	2.1	6.0	0.43
SO4 2	0.016896	0.024	0.048	1.065	r 6	0.425	74.0	0.35
Fe3+	0.0019920	0.00279	0.00558	0.604
K+	0.0000612	0.000078	0.000156	1.957
Cl−	0.0010295	0.001775	0.00355	2.032
Na+	0.0018515	0.000345	0.000690	1.334
SiO2	0.015	–	–	–

The dissolved concentration profiles as predicted by the SIS can be seen in Figure 11 for two example species. The dissolved and precipitated mass concentration profiles as predicted by the full model are given in Figure 12 for selected species. An example of the transient behavior is provided in Figure 13.

Figure 11

Dissolved concentrations predicted by the source influence solution for (A) Al3+ and (B) HCO3 2−

Figure 12

Dissolved and precipitated mass concentration profiles as predicted by the full model for (A) Ca+, (B) SO4 2−, (C) Al(OH)3, and (D) CaCO3 at t = 24 h

Figure 13

HCO3 − concentration after (A) 2.5 hours, (B) 5 hours, (C) 7.5 hours, and (D) 24 hours

It can be seen from the profiles presented in Figure 14 that, for this example, both the dissolved concentration profiles and the precipitated solids were accurately predicted by the PRS. Figure 14 shows the profiles of a group of species that includes species whose diffusion coefficient is furthest away from an indicator. The results shown are typical of the results of all considered species. The average relative mass balance error in this analysis is 2.77%. The CPU times using the full model and PRS2 were 8710 and 1543 seconds, respectively.

Figure 14

Longitudinal profiles taken along the x‐axis as predicted by the full model and PRS2 at t = 24 h

Advection dominant—PRS1 and PRS2

The final example is based on an advection‐dominant problem reported by Baroghel‐Bouny et al.35 The problem setup is the same as the diffusion case; however, in this test, the sample was first dried at a relative humidity of 4%, giving a uniform initial moisture content of S  = 0.09. The problem has been extended beyond that considered by Baroghel‐Bouny et al,35 with the inclusion of species Ca2+ and SO4 2−, and the chemical reaction considered has been changed to the formation of NaCl given by the nonequilibrium Freundlich isotherm, ie,

The time period considered was 3 hours, and the time step chosen was Δt = 0.9 s, giving a total number of time steps of 12 000. The finite element mesh and problem geometry were kept the same as for the diffusion case. The model parameters, BCs, and diffusion coefficients of the chemical species are given in Table 12. PRS1 and PRS2 were chosen to model this example, with OH−, Ca2+, and Cl− as indicator species.

Table 12

Model parameters

Parameter	Value	Species	Initial conc	Boundary conc	D mol
			(kg/kg)	(kg/kg)	(10−9 m2/s)
n	0.13	Na+	0.003312	0.013524	1.33
μ	0.3	OH−	0.012274	0.000000	5.3
λ	0.61	K+	0.024980	0.000000	1.96
τ (s)	72 000	Cl−	0.000000	0.038624	2.1
K in (m2)	10 × 10−21	Ca2+	0.000000	0.010000	0.79
p C 0 (kN/m2)	233.6 × 103	SO4 2−	0.003000	0.000000	1.07
p C b (kN/m2)	1.3 4 × 103

The concentration, sorbed mass, and saturation profiles as predicted by the full model and by both PRS1 and PRS2 can be seen in Figure 15, whereas transient profiles corresponding to a point 1 mm from the exposed face are given in Figure 16. The results of the PRS simulations are a good match to those of the full model, with the biggest differences being seen in the Cl− and K+ profiles for PRS1. The average relative mass balance error in this analysis is 2.24% for PRS1 and 1.23% for PRS2. The CPU times using the full model, PRS1, and PRS2 were 146, 83, and 103 seconds, respectively.

Figure 15

Concentration, sorbed mass, and saturation profiles as predicted by the full model, PRS1, and PRS2 at t = 3 h (where the black line indicates the position of the wetting front)

Figure 16

Transient concentration, sorbed mass, and saturation profiles as predicted by the full model, PRS1, and PRS2 at x = 1 mm

COMPUTATIONAL COST

The purpose of the reduction schemes is to reduce the computational cost of solving reactive chemical transport problems. In the above examples, the computational time for solving the problems with PRSs was substantially less than that of the associated full solutions. In order to quantify the reduction in computational cost, CPU times from each of the example problems presented in Sections 6 and 7 using the PRS are compared with those from a full solution. This analysis was performed on a PC with an Intel Core i7‐7700HQ at 2.80 GHz and 15.9 GB of usable RAM. The CPU time of the time‐step loop was measured for a number of runs, and an average was taken. The results for each of the problems considered are given in Tables 13 and 14.

Table 13

Normalized CPU times and percentage reduction for example problems

Example problem	Full model time (−)	PRS time (−)			Reduction (%)
		0	1	2	0	1	2
4 ion	1	0.616	–	–	38.37	–	–
6 ion	1	–	0.485	–	–	51.55	–
6 ion adv	1	–	0.571	0.709	–	42.93	29.14
10 ion	1	–	–	0.177	–	–	82.28
16 iona	1	0.226	0.280	0.461	77.39	71.97	53.92

Times given relate to the problem reduction scheme (PRS) applied over the applicable ranges and the corresponding full model time.

Table 14

Average Newton iterations for example problems

Example problem	Full model average	PRS average iterations (−)
	iterations (−)	0	1	2
4 ion	2	2	–	–
6 ion	2	–	2	–
6 ion adv	2	–	2	2
10 ion	2	–	–	2
16 ion	2	2	2	2

Abbreviation: PRS, problem reduction scheme.

It can be seen from the table that the PRS achieves significant reductions in computational times for the example problems, with reductions being up to 82% of CPU time. The smallest reduction is 29% for the advection‐dominant case (where the highly nonlinear moisture transport is responsible for the bulk of the computational cost). It can be seen from Table 14 that the implementation of the reduction scheme does not increase the number of Newton iterations required per time step. Table 14 also shows that the reduction in computational cost achieved by the scheme is a result of reducing the size of the nonlinear system through the reduction of the number of coupled nonlinear PDEs to be solved, as opposed to reducing the nonlinearity of the PDEs or improving the convergence of the Newton‐Raphson procedure. The nature of the predictor‐corrector approach (Equations (40) and (41)) makes this part of the model readily parallelizable, and, in addition to this, the model is compatible with domain decomposition methods, both of which would reduce the computational cost further.

Finally, it is noted that the reduction in total solution times gained from using the proposed PRS may depend on the moisture transport model employed and on whether or not temperature is included as a primary variable. In the present case, the authors considered isothermal problems and have used a single degree of freedom (P ) to simulate moisture flow, whereas Gawin et al32 employed two degrees of freedom in their moisture transport model. The authors expect that the reduction gained, relative to a full model solution, from reducing the number of chemical species that are primary variables would be different if a temperature‐dependent two‐variable moisture model had been used. However, we believe that the use of the PRS will always result in a substantial reduction in CPU time (for any practical engineering problem) when the number of primary chemical species is reduced significantly.

CONCLUSIONS

The coupled computational model described in this paper is capable of simulating moisture and chemical transport in porous media, as demonstrated in a validation exercise using experimental data from the works of Baroghel‐Bouny et al35 and Francy.43

The computational cost of large‐scale multispecies chemical transport problems can become prohibitively expensive, and therefore, efficient computational procedures and/or PRSs are required to solve such problems within reasonable times.

The approach of undertaking full solutions for a limited number of chemical (indicator) species and interpolating the response of the other (non‐indicator) species provides a direct and effective method for reducing the computational size of multispecies chemical transport problems.

The new predictor‐corrector PRS, based on Lagrangian interpolation and mass balance correction, is able to solve multispecies coupled chemical transport problems with good accuracy and efficiency, with errors in chemical concentration and mass being within 1% and 4%, respectively, of those from a full solution. The accuracy of the PRS depends on the order of the scheme, with higher‐order schemes being applicable to problems that have larger ranges of diffusion coefficients. The ranges of applicability of zero‐, first‐, and second‐order schemes (with one, two, and three indicator species, respectively), in terms of a reference (lower) diffusion constant D, are from D to 1.9D, from D to 6D, and from D to 16D, respectively, with an accuracy level of 1% (or 3% for PRS0).

The reduction in computational cost of an example 10‐ion transport problem using the PRS is 82% relative to the corresponding cost of a full solution.

Overall, the new numerical procedure described in this paper has the potential to offer significant reductions in computational demand for highly coupled multispecies reactive transport models. Furthermore, the particular nature of the PRS algorithm makes it suitable for parallelization, and the scheme is compatible with domain decomposition methods. This implies that the proposed PRS is capable of overcoming existing barriers to the exploitation of supercomputing facilities in the solution of such coupled problems.