key: cord-0010153-wrl6026n
authors: Freeman, Brubeck Lee; Cleall, Peter John; Jefferson, Anthony Duncan
title: An indicator‐based problem reduction scheme for coupled reactive transport models
date: 2019-08-27
journal: Int J Numer Methods Eng
DOI: 10.1002/nme.6186
sha: b2bedd9406f6529cc4b36b446f8913d3df4c5e91
doc_id: 10153
cord_uid: wrl6026n

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.

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-diffusion 1 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-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 ions [15] [16] [17] [18] or calcium leaching 19, 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 al 23 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 cost 25 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 Knabner 26 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.

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 C ) and chemical concentration (c i ) 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.

Transport in the porous medium is described using mass balance equations for liquid moisture, moisture vapor, and chemical species. Defining the domain ∈ R d with boundary and noting that the time interval considered is given as t ∈ [t 0 , T], 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 n denotes the porosity, S w is the degree of liquid moisture saturation, w is the liquid density, S g is the degree of gas saturation, and v is the moisture vapor density. t denotes a time derivative, v w is the liquid moisture velocity, t (m v ) is the rate of moisture transfer from liquid to vapor, and J v is the moisture vapor diffusion described here by Fick's law. 32 That is,

where g is the gas density and D mv is the moisture vapor diffusivity tensor given by 22

where A v = 1, B v = 1.667, and f v = 0.01. D v0 = 2.47 × 10 −5 m 2 /s is the moisture vapor diffusion coefficient in air, T r = 273 K is the reference temperature, and P g 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, w is the viscosity of the fluid, and K rw is the relative permeability of the moisture phase given by 8 

where A w = 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).

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 p is the degree of saturation of precipitated or sorbed material, p is the mass of the precipitated or sorbed material, and J i d 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, D mol i 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. D d i 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 condition 10, [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 Bachmat 1 is adopted here, where the mechanical dispersion tensor is given as

where L and T are the longitudinal and transverse dispersivities, respectively, and kj is the Kronecker delta ( kj = 1 if k = j or kj = 0 if k ≠ j).

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 follows 35 :

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.

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 w , q v , and q c are the prescribed fluxes of the moisture, moisture vapor, and chemical species, respectively, and c and c are the convective transfer coefficients of the moisture and chemical species, respectively. n 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

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 ∈ U, such that

is the vector of primary variables, and the space for the trial functions is defined as

In this study, the continuous Galerkin weighted-residual method 32 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 c ) is interpolated from the vector of nodal variables in the standard manner (ie, P c = ∑ N i P c i = NP c ). 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 scheme 33 for time discretization leads to

This set of nonlinear equations can then be solved using a standard Newton-Raphson procedure 33 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 ( t+1 k+1 ):

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 al 12 :

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.

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.

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 D d i is a function of the liquid moisture velocity and the dispersivities (see Equation (14)). Noting that both D d i and v w 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, Δc ind , based on the relative difference between the known diffusion coefficients of the current non-indicator species D i and of the indicator species D ind and D jnd , as follows:

in which g(D i , Δc ind ) 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 D i and D ind 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) 

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 c ), 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 t + Δt, as follows:

where c i is the concentration of a chemical species, c gr 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 Fetter 42 

, where R 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) ( f c ), 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 f c 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 n, w , and D mol i 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 c , as follows:

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 c err i is the concentration correction vector for species i and M i is an approximate tangent coefficient matrix. Two forms of M i were considered as follows:

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.

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

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

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.

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 al 35 is considered. This problem is based on experiments carried out by Francy 43 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 .

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. 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.

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 .

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.

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. 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 (D ind,m ) 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 D ind,m that produced the most accurate solutions.

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 RE = |(c prs i − c full i )∕c b i |. 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 0 1 2 Diffusion coefficient 0.8D ind <D i <1.5D ind D ind,u <6D ind,l D ind,u <16D ind,l a range Abbreviation: PRS, problem reduction scheme. a 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%. 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.

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. 

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. b This recommendation is based on the fact that PRS0 was found to have an upper limit to its applicable diffusion coefficient range of 1.5D ind (see Table 5 ).

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 L 2 -error norm was employed in these convergence studies:

where e is the error vector and c r 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.

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 al 44 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.

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.

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. 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.

Diffusion experiments, as well as accompanying numerical simulations, presented by Song et al 36 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 − , SO 4 2− , and Ca 2+ ). The reaction rates for the chemical species were considered as given by

where k a and k d 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 Ca 2+ is obtained from t (S p ) Ca = r 1 +r 6 ). Song et al 36 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 Δt = 36 s, 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 (D ind,u < 6D ind,l ); 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 ).

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. 

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. 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 + , Al 3+ , and K + were chosen as indicator species.

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 .

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. 

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 w = 0.09. The problem has been extended beyond that considered by Baroghel-Bouny et al, 35 with the inclusion of species Ca 2+ and SO 4 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 − , Ca 2+ , and Cl − as indicator species.

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 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.

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 . 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 (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 c ) to simulate moisture flow, whereas Gawin et al 32 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.

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 al 35 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. 

Introduction to Modelling of Transport Phenomena in Porous Media

A model for simulating transport of reactive multispecies components: model development and demonstration

A geochemical transport model for redox-controlled movement of mineral fronts in groundwater flow systems: a case of nitrate removal by oxidation of pyrite

Modeling of multicomponent reactive transport in groundwater: 1. Model development and evaluation

PHT3D: a reactive multicomponent transport model for saturated porous media

PhreeqcRM: a reaction module for transport simulators based on the geochemical model PHREEQC

Analysis of the Thmc behaviour of compacted swelling clay for radioactive waste isolation

Inclusion of some aspects of chemical behavior of unsaturated soil in thermo/hydro/chemical/mechanical models. I: model development

Inclusion of some aspects of chemical behavior of unsaturated soil in thermo/hydro/chemical/mechanical models. II: application and transport of soluble salts in compacted bentonite

Diffusive reactive transport of multicomponent chemicals under coupled thermal, hydraulic, chemical and mechanical conditions

Reactive transport of chemicals in unsaturated soils: numerical model development and verification

Multi-component reactive transport modeling of natural attenuation of an acid plume at a uranium mill tailings site

Hydrochemical modelling for preliminary assessment of minewater pollution

Red mud and fly ash for remediation of mine sites contaminated with As, Cd, Cu, Pb and Zn

Modelling chloride diffusion in concrete effect of fly ash and slag

Factors influencing chloride transport in concrete structures exposed to marine environments

Assessment and prediction of RC structure service life by means of durability indicators and physical/chemical models

Modelling of salt crystallization in building materials with microstructure -poromechanical approach

Coupled chemo-mechanical deterioration of cementitious materials. Part I: modeling

Modeling deterioration of cementitious materials exposed to calcium leaching in non-isothermal conditions

Chemical-diffusive modelling of the self-healing behavior in concrete

A coupled thermo-hygro-chemical model for characterising autogenous healing in ordinary cementitious materials

Use of parallel computing and visualisation techniques in the simulation of large scale geoenvironmental engineering problems

Accuracy of operator splitting for advection-dispersion-reaction problems

A critical evaluation of recent developments in hydrogeochemical transport models of reactive mutichemical components

A reduction scheme for coupled multicomponent transport-reaction problems in porous media: generalization to problems with heterogeneous equilibrium reactions

A general reduction scheme for reactive transport in porous media

On the behavior of approaches to simulate reactive transport

Solute transport with multiple equilibrium controlled or kinetically controlled chemical reactions

A formulation for decoupling components in reactive transport problems

A reaction-based paradigm to model reactive chemical transport in groundwater with general kinetic and equilibrium reactions

Hygro-thermo-chemo-mechanical modelling of concrete at early ages and beyond. Part I: hydration and hygro-thermal phenomena

Porosity development in a thermo-hygral finite element model for cementitious materials

The finite element method in the static and dynamic deformation and consolidation of porous media

Modelling of isothermal coupled moisture-ion transport in cementitious materials

Modeling of chloride diffusion in concrete immersed in CaCl 2 and NaCl solutions with account of multi-phase reactions and ionic interactions

The treatment of multi-component diffusion and ion pairs in diagenetic fluxes

Transport of reacting solutes in porous media: relation between mathematical nature of problem formulation and chemical nature of reactions

Analytical and variational numerical methods for unstable miscible displacement flows in porous media

Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-stokes equations

Finite element methods for modeling water flow in variably saturated porous media: numerical oscillation and mass-distributed schemes

Modelling of Chloride Ions Ingress in Partially Water Saturated Mortars

Stabilized reduced order models for the advection-diffusion-reaction equation using operator splitting

A reduced finite element formulation based on POD method for two-dimensional solute transport problems

A closed-form equation for predicting the hydraulic conductivity of unsaturated soils

The authors gratefully acknowledge primary financial support from the Engineering and Physical Sciences Research Council under grant EP/P02081X/1 through the Resilient Materials for Life (RM4L) project. The authors would also like to acknowledge support from the Natural Environment Research Council under grant NE/L013908/1 through the "In situ recovery of resources from waste repositories" project. Information on the data underpinning the results presented here, including how to access them, can be found in the Cardiff University data catalog at http://doi.org/10.17035/d.2019. 0079832549.