A model of carbon dioxide dissolution and mineral carbonation kinetics

Mark J. Mitchell, Oliver E. Jensen, K. Andrew Cliffe, M. Mercedes Maroto-Valer


The kinetics of the dissolution of carbon dioxide in water and subsequent chemical reactions through to the formation of calcium carbonate, a system of reactions integral to carbon sequestration and anthropogenic ocean acidification, is mathematically modelled using the mass action law. This group of reactions is expressed as a system of five coupled nonlinear ordinary differential equations, with 14 independent parameters. The evolution of this system to equilibrium at 25°C and 1 atm, following an instantaneous injection of gaseous carbon dioxide, is simulated. An asymptotic analysis captures the leading-order behaviour of the system over six disparate time scales, yielding expressions for all species in each time scale. These approximations show excellent agreement with simulations of the full system, and give remarkably simple formulae for the equilibrium concentrations.

1. Introduction

The quantity of carbon dioxide in the atmosphere has increased by about a third since the start of the Industrial Revolution in the late eighteenth century, from around 280 parts per million by volume (ppmv) in 1850 to around 380 ppmv to date, and is currently rising by about 1.7 ppmv per year (Lal 2008). When carbon dioxide reacts with water, carbonic acid is formed, from which hydrogen ions dissociate, increasing the acidity of the system. Therefore, in addition to any greenhouse effect, anthropogenic carbon dioxide emissions to the atmosphere can increase the acidity of the atmosphere and precipitation. Around 30–40% of anthropogenic carbon dioxide emitted to the atmosphere dissolves into the oceans (Millero 1995), where the reaction with sea water has increased ocean acidity by 0.1 pH units since pre-industrial times (Intergovernmental Panel on Climate Change (IPCC) 2007). Ocean ecosystems are affected both through acidification and by associated reductions in carbonate ion concentrations (Orr et al. 2005). Around 20 per cent of anthropogenic carbon dioxide emissions to the atmosphere are absorbed by the terrestrial biosphere (Feely et al. 2004), and the reaction between the absorbed carbon dioxide and soil moisture can alter soil acidity. Carbon dioxide emissions to the atmosphere can therefore increase the acidity of land, sea and air.

One promising solution to reduce CO2 emissions is carbon capture and storage (CCS). There are a number of options for long-term storage, of which saline aquifers have the largest capacity (IPCC 2005) and may therefore have the most potential. Saline aquifers are very large deep porous geological formations saturated with brine, and are often rich in different metals. Although the initial principal trapping mechanism is often an impermeable cap rock above the underground reservoir (‘structural trapping’), solubility trapping and mineral trapping are also important in the longer term (IPCC 2005). Solubility trapping occurs when carbon dioxide dissolves into the brine solution, and mineral trapping occurs when the dissolved carbon dioxide reacts with the water to eventually form stable carbonate compounds such as calcium carbonate and magnesium carbonate (Lagneau et al. 2005; Druckenmiller et al. 2006). Long-term storage may also include the generation of stable carbonate compounds in an industrial process above ground, or the injection of carbon dioxide into the deep oceans.

Among the numerous factors to be taken into account in developing computational models of CCS, the chemical reactions that take place between injected CO2, water and the various metals which may be present are fundamental. These interactions are strongly coupled with the associated transport processes, for example, via precipitates building up in the pore spaces and affecting permeability, or dissolved carbon dioxide changing the density of groundwater leading to buoyant fluxes. A variety of sophisticated modelling and simulation packages have been developed for fluid flows with chemical reactions, including TOUGH2/TOUGHREACT (Xu et al. 2006), PHREEQC (Parkhurst & Appelo 1999), FEHM (Zyvoloski et al. 1991), Geochemist's Workbench (Bethke 2002) and CHILLER (Reed 1982; Reed & Spycher 1998). In order to provide improved theoretical understanding of the kinetics of the dissolution of CO2 in water and the subsequent chemical reactions, and to offer systematic approximations that may be exploited in multi-scale modelling tools, this paper seeks to simplify the complex kinetic modelling of these key reactions, by using the method of matched asymptotic expansions to identify the distinct time scales over which the reactions take place and to provide simple expressions for the resulting dynamic and equilibrium concentrations. We consider the following system of reversible chemical reactions (IPCC 2005), where the system is initially at equilibrium with all species at low concentrations compared with water, and is then subject to an instantaneous injection of gaseous carbon dioxide: Embedded Image 1.1a Embedded Image 1.1b Embedded Image 1.1c Embedded Image 1.1d Embedded Image 1.1e Embedded Image 1.1f

In reaction (1.1a), gaseous carbon dioxide is dissolved in water, reacting to form carbonic acid (1.1b). Hydrogen ions dissociate from the carbonic acid, to give bicarbonate (1.1c), and then a carbonate ion (1.1d ), which then reacts with a calcium cation to form calcium carbonate (1.1e). Some of this calcium carbonate precipitates out of the solution (1.1f ). We do not include any further reactions with buffers that remove some of the hydrogen ions from the system; this removal will only drive the system of reactions further forward, and so result in more hydrogen being produced. We therefore obtain the minimum quantity of hydrogen ions that is generated.

We construct a mathematical model for the kinetics of this system of reactions using the mass action law and known rate and equilibrium constants. Although the model involves coupled fifth-order nonlinear ordinary differential equations (ODEs), our asymptotic approach yields remarkably simple and accurate approximations for both transient and equilibrium concentrations. This is possible because the rate constants differ by many orders of magnitude, allowing the essence of the system’s behaviour to be captured by much-simplified sets of equations. Our analysis identifies six different time scales within the system, each exhibiting a different dominant balance of terms representing locally rate-limiting reactions. Approximations are matched between time scales and are validated against numerical simulations. We anticipate that the predictions of our study will be of value in simulations relating to CCS or ocean acidification involving processes that overlap with the slower time scales under consideration here; our results provide useful approximations of the more rapid reaction processes and will help reduce the computational cost of resolving reactions with rapid rate constants that lead to numerical stiffness.

We present the governing equations and parameters in §2, and demonstrate the accuracy of asymptotic approximations against simulations in §3. The asymptotic analysis is outlined in §4; technical details appear in appendix A. Section 5 contains a summary of the key predictions and an explanation of their range of validity.

2. Mathematical modelling

(a) Full dimensional model

Reactions (1.1ae) are assumed to follow the mass action law, with the additional assumptions that water is abundant, that the components are always well mixed and that the temperature, pressure, salinity and ionic strength are constant and uniform, and any impact upon the rates of reaction caused by varying one or more of these factors can be captured by varying the rate constants appropriately. It is also assumed that CO2(aq) remains below carbon dioxide’s solubility limit (without the need to impose a maximum value on its concentration), and that all the CaCO3 that is formed in solution remains dissolved until its concentration reaches calcium carbonate’s molar solubility limit S; thereafter, it precipitates instantaneously. We also assume that the reactions take place in a homogeneous environment (so any nucleation that occurs will also be homogeneous).

The kinetics of equations (1.1) is therefore modelled by the following rate equations, where t is time (in seconds), [ ] denotes concentration (measured in moles per litre, or M), [CaCO3]=[CaCO3(aq)]+[CaCO3(s)], k1 to k5 and k−1 to k−5 are, respectively, the forward and reverse rate constants for the five reactions (1.1ae) and H(.) is the Heaviside step function (H(x)=1 if x>0, H(x)=0 if x<0): Embedded Image 2.1a Embedded Image 2.1b Embedded Image 2.1c Embedded Image 2.1d Embedded Image 2.1e Embedded Image 2.1f Embedded Image 2.1g Embedded Image 2.1h Embedded Image 2.1i Embedded Image 2.1j

Equations (2.1) can be simplified by removing the explicit distinction between dissolved and precipitated CaCO3 and by imposing mass-conservation constraints consistent with the initial conditions. From equations (2.1i,j),Embedded Image 2.2 equations (2.1g,h) and (2.2) can alternatively be written as Embedded Image 2.3aand Embedded Image 2.3b Therefore, equivalently to equations (2.1), the system comprises equations (2.1af) and (2.3), and whenever a more detailed split between aqueous and solid CaCO3 is required, this is given in terms of the solubility limit S byEmbedded Image 2.4

The four different elements within this system (C, O, H and Ca) must be conserved throughout the system’s evolution, and so the system of nine equations has four constraints and five degrees of freedom. Noting that equations (2.1hj) sum to zero, it follows thatEmbedded Image (where subscript zero indicates an initial value), denoting conservation of calcium. Similarly, conservation of carbon, hydrogen and oxygen are, respectively, given byEmbedded Image and Embedded Image

By constructing an alternative formulation of the four constraints, such that [CO2(g)],[H+],[CaCO3] (denoted X, Y and Z, respectively) and [H2O] are expressed as functions of the remaining five species and the initial concentrations, the system (2.1af) and (2.3) simplifies to the following system of five coupled equations: Embedded Image 2.5a Embedded Image 2.5b Embedded Image 2.5c Embedded Image 2.5d Embedded Image 2.5e where Embedded Image 2.6a Embedded Image 2.6b Embedded Image 2.6c This simplification reduces the number of parameters from 19 (the 10 rate constants and nine initial concentrations) to 18.

(b) Parameter estimates and initial conditions

The system to be analysed is based on some of the shallower wells of the Botucatu aquifer in Brazil, which has been well characterized (Meng & Maynard 2001; e.g. Wells 56, 89 and 149). These shallower wells have temperatures close to 25°C, and so estimates of the magnitudes of the rate constants at 25°C and 1 atm are given in table 1. Note that k1 and k−1 are estimates based on CO2 distributing itself ‘approximately equally between its gas phase and aqueous phase at ordinary temperatures’ (Stumm & Morgan 1996), and reaction (1.1a) being fast in a well-mixed scenario (Rochelle et al. 2004) and similar in speed to reaction (1.1c). Also, k5 and k−5 are estimates based on reaction (1.1e) being relatively slow, with K5≈1660 at 25°C (Drever 2002). As CaCO3 dissociates to one molecule of Ca2+ and one molecule of COEmbedded Image, S (in M) is given by the square root of the solubility product constant, which, for CaCO3 in pure water at 25°C, is 4.96×10−9 (Tro 2008).

View this table:
Table 1.

Parameter estimates at 25°C and 1 atm.

The behaviour of the system will be illustrated with initial concentrations [Ca2+]0=10 mg l−1=2.50×10−4 M, [HCOEmbedded Image M and pH=6, i.e. [H+]0=10−6 M. From these three initial concentrations, the remaining initial concentrations of the species are calculated using the fact that a reaction’s equilibrium constant equals both its ratio of rate constants and its ratio of equilibrium concentrations. For example, for reaction (1.1d ),Embedded Image 2.7 and so Embedded Image M. Similar calculations for the other four reactions give the remaining initial equilibrium concentrations (table 2), with water’s concentration not deviating from its large pure value. Although aquifer conditions can be far from 25°C, 1 atm and these chosen initial concentrations, the results obtained here are extendable to a wider range of scenarios, as explained in §5.

We assume that the system is perturbed by an instantaneous injection of gaseous carbon dioxide that increases [CO2(g)] to 0.065 M. This is sufficient to increase [CO2(aq)] to carbon dioxide’s solubility limit, which at 25°C and 1 atm is about 0.033 M. Any more injected CO2(g) will therefore remain undissolved, and will not affect any of the other species. We seek the evolution of the system to a new equilibrium.

View this table:
Table 2.

Initial concentrations.

(c) Dimensionless model

In table 3, we define dimensionless variables B, D, F, G and H, representing deviations of concentrations from initial values, each a function of a rescaled time variable Embedded Image. In table 4, we define 18 dimensionless parameters, denoted with Greek symbols, with the initial concentration of water and Embedded Image chosen to be the arbitrary reference concentration and time scale, respectively. These parameters are evaluated using data in tables 1 and 2 and [CO2(g)]0=0.065 M, and the resulting values, spanning 23 orders of magnitude, are also given in table 4. The parameters are also expressed as the product of a corresponding O(1) parameter (indicated by a hat) and a power of the parameter α(=1.95×108). This is a precursor to taking the distinguished asymptotic limit Embedded Image, with all hatted parameters remaining O(1) in the limit.

View this table:
Table 3.

Dimensionless variables.

View this table:
Table 4.

Dimensionless parameters.

As the initial conditions (prior to CO2 injection) have been chosen to satisfy the equilibrium equations (as in equation (2.7)), reactions (1.1be) yield the following identities:Embedded Image 2.8 Substituting the dimensionless variables, parameters and identities (2.8) into equations (2.5) yields the following dimensionless system, containing 14 independent parameters: Embedded Image 2.9a Embedded Image 2.9b Embedded Image 2.9c Embedded Image 2.9d Embedded Image 2.9e The five variables all equal zero at Embedded Image.

3. The evolution of the system to equilibrium

Numerical integration of equations (2.5) using the parameter values given in tables 1 and 2 (with [CO2(g)]0=0.065 M) yields the results shown in figure 1a. (The evolution of calcium is given as the increase from its initial value.) As expected, the injection of CO2(g) leads to increases in [CO2(aq)], [H2CO3], [H+] and [HCOEmbedded Image]. The relative increase in [H+] is much larger than the increase in [HCOEmbedded Image], as reaction (1.1d) converts some of the generated HCOEmbedded Image into COEmbedded Image and further H+. The ratio of a reaction’s concentrations is constant for any equilibrium, so from reaction (1.1d),Embedded Image Therefore, a large increase in [H+] (from [H+]0 to [H+]Embedded Image) and a small increase in [HCOEmbedded Image] must be balanced by a (large) decrease in [COEmbedded Image] (as shown in figure 1a), as expected (Orr et al. 2005). Similarly, from reaction (1.1e),Embedded Image Therefore, a drop in [COEmbedded Image] must be balanced by an increase in [Ca2+] (as shown in figure 1a) and a decrease in [CaCO3]. Therefore, CO2 cannot be sequestered as CaCO3 through these reactions alone.

Figure 1.

(a) Simulation results expressed in dimensional variables, satisfying equations (2.5); (b) simulation of the dimensionless system (2.9) (continuous lines), with leading-order approximations from appendix B (dashed lines, where visible); and (c) the time exponent ΨJ of each species’ evolution, from simulation of equations (2.9). In each graph, the vertical lines t1t6 indicate the six time scales: t1=α−3/2, t2=α−1, t3=α−3/4, t4=α−1/2, t5=α−1/4, t6=α1/4.

The evolution of the system is plotted using dimensionless variables (2.9) in figure 1b with continuous lines, along with asymptotic approximations (dashed lines) that will be derived below. There is close agreement between the simulations and the approximations. For each species J (J=B,D,F,G,H), the exponentEmbedded Image 3.1 provides a useful measure of the rate of change of each variable (as Embedded Image when ΨJ is constant). Figure 1c shows how the ΨJ exhibit stepwise changes with time.

The plateaux in figure 1c separate six different time scales in this system (indicated by t1t6 in figure 1, and discussed in more detail below): (i) CO2(aq) evolves linearly in time (ΨB≈1), and then rapidly reaches equilibrium (ΨB=0), the remaining species initially evolving with successively higher powers of time, (ii) H2CO3 and HCOEmbedded Image reach an approximately linear evolution in time (ΨD and ΨF approach 1), with COEmbedded Image (G) varying quadratically, and Ca2+ (H) cubically in time, although the changes are all still negligible compared with the initial concentrations, (iii) and (iv) the rate of change of COEmbedded Image declines rapidly (ΨG falls to zero), (v) H2CO3, HCOEmbedded Image and COEmbedded Image (D, F and G) reach equilibrium while Ca2+ evolves linearly in time (ΨH≈1) and (vi) Ca2+ reaches equilibrium.

4. Asymptotic analysis

We seek asymptotic approximations of the solutions of equations (2.9) in each of the six time scales identified above; the approximations are matched with each other where the time scales overlap. Complete details are given in appendix A. Briefly, in each time scale, the five variables are expressed as the product of a corresponding O(1) variable (indicted by a hat), and an appropriate scaling factor, given as a power of the parameter α. This parameter is chosen as it is a measure of the quantity of injected CO2(g); as this is a nonlinear system, different initial conditions can produce qualitatively different outcomes, and our predictions are subject to the constraint (5.3) below. The hatted variables and parameters (from table 4) and all their corresponding scaling factors of α are substituted into equations (2.9), and the terms in successive powers of α−1/4 are balanced, giving the leading-order behaviour of the system at each time scale in the distinguished limit Embedded Image (hatted parameters in table 4 remaining O(1)).

The different leading-order approximations for each species are summarized in appendix B, and are combined into a single plot as the dashed lines that are graphically indistinguishable from the continuous lines in figure 1b. There is uniformly close agreement between the simulations and these approximations for the large, but finite, value of α chosen.

5. Results

(a) Equilibrium concentrations

The large-time limits of the leading-order approximations, given by equations (A3a), (A23ac) and (A26), give estimates of the equilibrium concentrations, which in dimensional terms are Embedded Image 5.1a Embedded Image 5.1b Embedded Image 5.1c Embedded Image 5.1d Embedded Image 5.1e These non-trivial algebraic relationships emerge from our systematic asymptotic analysis; their simplicity reflects the wide numerical range of dimensionless rate constants in table 4. Equation (5.1a) shows that the proportion of injected CO2(g) that is sequestered (i.e. dissolves into the water) is controlled only by the dissolution constant K1=k1/k−1; the further process of dissolved CO2 being removed by reacting with water is less important. In the present example, with K1=1, half of the injected CO2(g) dissolves in the water. The proportion of injected and dissolved CO2 that then reacts with the water to form carbonic acid is given simply by the ratio of this reaction’s rate constants (K2; equation (5.1b)). From equation (5.1c), the proportion of extra carbonic acid that ends up as bicarbonate is governed (through a quadratic relationship) by the carbonic acid dissociation constant. For the present parameter regime, the leading-order approximation for the final quantity of carbonate ions (5.1d) turns out to be independent of its initial value.

The concentrations of all the remaining species can be calculated using equations (5.1) and the conservation of mass constraints. From equations (5.1e) and (2.6c),Embedded Image as Embedded Image. Therefore, a decrease in Embedded Image must be accompanied by a decrease in [CaCO3]. Therefore, the amount of CO2 bound to calcium decreases as a result of these reactions. In this example, the concentration of calcium carbonate decreases to 2.2×10−10 M.

Also, substituting equations (5.1d,e) into equation (2.6b), and assuming Embedded Image and [CaCO3]0 are all negligible compared with Embedded Image (which is true for these parameter values with error O(α−1/4)), we findEmbedded Image 5.2 which equals 5.1×10−5 M. Therefore, with no buffers, the pH has dropped from 6 to 4.3.

(b) The six time scales

So far, the six time scales used to construct the asymptotic approximations have been defined as powers of α. Greater insight comes from examining the rate constants that appear in the solutions, such as the time scale Embedded Image implicit in (A2a). The six time scales, expressed in terms of the original dimensional quantities, are given in table 5, together with the corresponding values at 25°C and 1 atm. In this well-mixed scenario, CO2(aq) reaches equilibrium within O(10−11 s), H2CO3, HCOEmbedded Image and COEmbedded Image reach equilibrium within O(10−2 s) and Ca2+ reaches equilibrium in O(103 s). Therefore, in practice, the rates of the first four reactions will be determined by the neglected transport and mixing processes, and may be assumed to be in local equilibrium as described by equations (5.1ad). The reaction with calcium takes longer, and in practice will therefore be determined by transport, mixing and reaction processes. In dimensional terms, its evolution in the well-mixed scenario is given by equation (A25),Embedded Image

View this table:
Table 5.

Summary of the six time scales.

The results obtained above remain valid while the six time scales remain distinct, i.e.Embedded Image 5.3 The validity of the model can therefore be assessed under different initial concentrations, and under different conditions (temperature, pressure, salinity, surface area, etc.), via the impact on the rate constants, provided this ordering is preserved.

6. Discussion

We have formulated a nonlinear mathematical model of reactions (1.1) following an instantaneous injection of CO2(g) into a dilute system, based on the mass action law. We have also used the method of matched asymptotic expansions to derive leading-order approximations to the solutions of equations (1.1), which agree closely with numerical simulations (figure 1). We identified six distinct time scales within this system (table 5), over which there are different dominant balances of terms. Both the transient evolution (appendix B) and the ultimate equilibria (5.1) can be captured by a remarkably simple but non-trivial set of functions in terms of the initial concentrations and equilibrium constants only. For example, we are able to predict explicitly the level of acidification arising from gaseous CO2 injection (5.2).

The full nonlinear model, when expressed in dimensionless variables, comprises a set of five ODEs with 14 independent parameters. We chose one of these (1/α) as the small parameter around which we constructed an asymptotic expansion. We assigned to all the remaining parameters a power of α (table 4), and formally took the distinguished limit in which the parameters varied proportionally to the appropriate power of α (even though these relationships may not have physical justification). In the limit Embedded Image, the six resulting time scales become increasingly widely separated in magnitude, and the asymptotic approximation should converge to the exact solution of the ODEs. In practice, we make this comparison at a finite, but large, value of α (using physically realistic parameter values), still demonstrating good agreement. To extend the present approximation to other parameter regimes (such as different temperatures and pressures, or smaller initial injections of CO2(g), for which α is smaller in magnitude), it is essential that the six time scales (stated in physical variables in table 5) remain reasonably widely separated. In the present example, the closeness of the expressions for time scales 4 and 5 indicates where the present approximation may first break down.

At 25°C and 1 atm, the first five time scales are rapid (less than 1 s), and mixing processes not captured in the present model are likely to be rate-limiting, particularly for the extremely rapid dissolution of CO2 (time scale 1). In this well-mixed system, however, the reaction with calcium is ultimately rate-limiting. The present model provides a useful foundation on which to build more sophisticated multi-scale models, incorporating, for example, additional reactions that buffer pH, interactions with dissolved salt, catalysts or other metals (such as olivine (Mg2SiO4) and serpentine (Mg3Si2O5(OH)4)), as well as spatial effects relevant to aquifers or other applications.


This study was supported by EPSRC grant EP/F012098/1.

Appendix A

A.1. Time scale 1: Embedded Image

Writing Embedded Image, Embedded Image, Embedded Image, Embedded Image, Embedded Image and Embedded Image gives the following leading-order behaviour of equations (2.9) on this time scale as Embedded Image: Embedded Image A 1a Embedded Image A 1b Embedded Image A 1c Embedded Image A 1d Embedded Image A 1e with Embedded Image at Embedded Image. Embedded Image equilibrates rapidly, with changes in the remaining variables slaved sequentially to this. Equations (A 1) have the solutions Embedded Image A 2a Embedded Image A 2b Embedded Image A 2c Embedded Image A 2d Embedded Image A 2e where Embedded Image. The species evolve in successively higher powers of time, as expected from figure 1c. Examination of the neglected terms in equations (A 1) reveals that these approximations are valid until Embedded Image reaches O(α1/2), i.e. Embedded Image, at which time Embedded Image, Embedded Image, Embedded Image, Embedded Image and Embedded Image. The leading-order large-time limits of these variables, valid for Embedded Image (i.e. between time scales t1 and t2), form the initial conditions for the next time scale. These limits are, in dimensionless unscaled variables, Embedded Image A 3a Embedded Image A 3b Embedded Image A 3c Embedded Image A 3d Embedded Image A 3e

A.2. Time scale 2: Embedded Image = O(α−1)

We now define O(1) variables using Embedded Image, Embedded Image, Embedded Image, Embedded Image, Embedded Image and Embedded Image to yield the following leading-order reduction of equations (2.9): Embedded Image A 4a Embedded Image A 4b Embedded Image A 4c Embedded Image A 4d Embedded Image A 4e Therefore, Embedded Image. Also, d(Embedded Image)/dEmbedded Image, implying that Embedded Image (the constant of integration is zero, by matching with equations (A 3)). Substituting this expression for Embedded Image into equation (A 4c), solving for Embedded Image and matching the small-time limit of this expression with equation (A 3c), givesEmbedded Image where Embedded Image. Likewise, solving equations (A 4d,e) and matching with equations (A 3d,e) giveEmbedded Image and Embedded Image Again, looking at the largest neglected terms in equations (A 4), we find that these approximations are valid until Embedded Image reaches O(α1/4), at which time Embedded Image, Embedded Image, Embedded Image, Embedded Image and Embedded Image. The leading-order large-time limits of these variables, valid for Embedded Image (i.e. between time scales t2 and t3), form the initial conditions for the next time scale. These limits are (A 3a) plus Embedded Image A 5a Embedded Image A 5b Embedded Image A 5c Embedded Image A 5d broadly consistent with the exponents ΨB≈0, ΨD≈1, ΨF≈1, ΨG≈2 and ΨH≈3 for Embedded Image in figure 1c.

A.3. Time scale 3: Embedded Image

On this time scale, we define O(1) variables using Embedded Image, Embedded Image, Embedded Image, Embedded Image, Embedded Image and Embedded Image, to yield the following reduction of equations (2.9): Embedded Image A 6a Embedded Image A 6b Embedded Image A 6c Embedded Image A 6d with Embedded Image as in equation (A 4a). As the left-hand sides of equations (A 6a,b) are O(1), the O(α1/4) terms on the right-hand sides must sum to zero, i.e.Embedded Image A 7

In the limit Embedded Image, the leading-order behaviour of equations (A 6) on this time scale is given by equations (A 6c,d) plus Embedded Image A 8aand Embedded Image A 8b

Therefore, d(Embedded Image)/dEmbedded Image, implying Embedded Image (matching with equations (A 5a,b)). Substituting this into equation (A 7) gives Embedded Image A 9aand Embedded Image A 9b Substituting equation (A 9b) into equation (A 6c), and solving for Embedded Image givesEmbedded Image where Embedded Image and Embedded Image. Embedded Image is given byEmbedded Image The constants of integration for Embedded Image and Embedded Image can be deduced by matching with equations (A 5c,d). The neglected terms in equations (A 6) indicate that these approximations are valid until Embedded Image reaches O(α1/4), at which point, Embedded Image is O(1), Embedded Image, Embedded Image and Embedded Image are O(α1/4), and Embedded Image is O(α1/2). The leading-order large-time limits of these variables, valid for Embedded Image (i.e. between time scales t3 and t4), are given by equations (A 3a), (A 5a,b),Embedded Image A 10 consistent with ΨG≈1, ΨH≈2 in figure 1c (although these exponents do not exhibit clear plateaux for Embedded Image).

A.4. Time scale 4: Embedded Image

On the next time scale, we define O(1) variables using Embedded Image, with Embedded Image, Embedded Image, Embedded Image, Embedded Image, Embedded Image and Embedded Image. Similarly to time scale 3, the resulting O(α1/2) terms in the ODEs for D and F must sum to zero, which gives equation (A 7), as must the O(α1/4) terms in the ODEs for D, F and G, which give Embedded Image A 11aand Embedded Image A 11b

Then, in the limit Embedded Image, the leading-order behaviour is given by Embedded Image A 12a Embedded Image A 12b Embedded Image A 12c Embedded Image A 12d Therefore, d(Embedded Image)/dEmbedded Image, and so Embedded Image. This, together with equation (A 7), gives equations (A 9). Substituting equation (A 9b) into equation (A 11b) gives Embedded Image A 13a and substituting equation (A 13a) into equation (A 12d) gives Embedded Image A 13b which match with equations (A 10) for Embedded Image. These approximations are valid until Embedded Image reaches O(α1/4), at which point Embedded Image and Embedded Image are O(1), and Embedded Image, Embedded Image and Embedded Image are O(α1/4). The large-time limits of these variables, valid for Embedded Image (i.e. between time scales t4 and t5), and consistent with ΨG≈0, ΨH≈1, are given by equations (A 3a), (A 5a,b), plusEmbedded Image A 14

A.5. Time scale 5: Embedded Image

On this time scale, we define O(1) variables using Embedded Image, with Embedded Image, Embedded Image, Embedded Image, Embedded Image, Embedded Image and Embedded Image. The O(α3/4) terms in the ODEs for D, F and G give Embedded Image A 15aand Embedded Image A 15b the O(α1/2) terms give Embedded Image A 16aand Embedded Image A 16b and the O(α1/4) terms give Embedded Image A 17aand Embedded Image A 17b

Then, the leading-order behaviour on this time scale is given by Embedded Image A 18a Embedded Image A 18b Embedded Image A 18c Embedded Image A 18d

Therefore, d(Embedded Image)/dEmbedded Image, and from equation (A 15a), the evolution of Embedded Image and Embedded Image on this time scale is given by Embedded Image A 19aand Embedded Image A 19b and so Embedded Image has the implicit solutionEmbedded Image A 20 where Embedded Image. From equations (A 15b), (A 16b) and (A 18d), respectively, the evolution of Embedded Image and Embedded Image is given by Embedded Image A 21a Embedded Image A 21b Embedded Image A 21c

Thus, Embedded Image, Embedded Image and Embedded Image rapidly equilibrate for Embedded Image, and Embedded Image evolves linearly in time, consistent with figure 1c. As Embedded Image as Embedded Image, the small-time limits of equations (A 19) and (A 21b) are Embedded Image A 22a Embedded Image A 22b Embedded Image A 22c ensuring matching with equations (A 9) and (A 14). These approximations are valid until Embedded Image reaches O(α1/2), at which point, Embedded Image, Embedded Image, Embedded Image and Embedded Image are O(1), and Embedded Image is O(α1/2). The large-time limits for D and F, valid for Embedded Image (i.e. between time scales t5 and t6), are Embedded Image A 23aand Embedded Image A 23b The corresponding large-time limits of G and H are Embedded Image A 23c and Embedded Image A 23d (with F satisfying equation (A 23b)), consistent with ΨB=ΨD=ΨF=ΨG=0 and ΨH≈1.

A.6. Time scale 6: Embedded Image

On this final time scale, we define O(1) variables Embedded Image, with Embedded Image, Embedded Image, Embedded Image, Embedded Image, Embedded Image and Embedded Image. Substituting these expansions into equations (2.9) gives identities (A15–A17). Also, O(1) balances are given byEmbedded Image Therefore, the solutions to D, F and G on this time scale are given by equations (A 23ac). Also,Embedded Image A 24 and the solution to equation (A 24), which includes the first correction term for G and matches with equation (A 23d), isEmbedded Image A 25 Thus, for Embedded Image, the large-time limits are equations (A 3a), (A 23ac) andEmbedded Image A 26 (with F satisfying equation (A 23b)).

Appendix B

The combined expressions for the five species, given below, are plotted in figure 1b. Wherever possible, they are constructed by summing the different time scales’ approximations, and subtracting the intervening large-time limits, to give composite asymptotic formulae. Otherwise, the expressions have been patched together to provide a good fit to the full simulations as follows:Embedded Image and Embedded Image where Embedded Image is given implicitly by Embedded Image Ωγ+ν, Δ≡ψ+ψδ/μ and Embedded Image.


    • Received July 5, 2009.
    • Accepted November 11, 2009.
  • This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


View Abstract