Abstract
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 leadingorder 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 preindustrial 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 CO_{2} emissions is carbon capture and storage (CCS). There are a number of options for longterm 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). Longterm 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 CO_{2}, 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 CO_{2} in water and the subsequent chemical reactions, and to offer systematic approximations that may be exploited in multiscale 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: 1.1a 1.1b 1.1c 1.1d 1.1e 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 fifthorder 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 muchsimplified sets of equations. Our analysis identifies six different time scales within the system, each exhibiting a different dominant balance of terms representing locally ratelimiting 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.1a–e) 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 CO_{2}(aq) remains below carbon dioxide’s solubility limit (without the need to impose a maximum value on its concentration), and that all the CaCO_{3} 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), [CaCO_{3}]=[CaCO_{3}(aq)]+[CaCO_{3}(s)], k_{1} to k_{5} and k_{−1} to k_{−5} are, respectively, the forward and reverse rate constants for the five reactions (1.1a–e) and H(.) is the Heaviside step function (H(x)=1 if x>0, H(x)=0 if x<0): 2.1a 2.1b 2.1c 2.1d 2.1e 2.1f 2.1g 2.1h 2.1i 2.1j
Equations (2.1) can be simplified by removing the explicit distinction between dissolved and precipitated CaCO_{3} and by imposing massconservation constraints consistent with the initial conditions. From equations (2.1i,j), 2.2 equations (2.1g,h) and (2.2) can alternatively be written as 2.3aand 2.3b Therefore, equivalently to equations (2.1), the system comprises equations (2.1a–f) and (2.3), and whenever a more detailed split between aqueous and solid CaCO_{3} is required, this is given in terms of the solubility limit S by 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.1h–j) sum to zero, it follows that (where subscript zero indicates an initial value), denoting conservation of calcium. Similarly, conservation of carbon, hydrogen and oxygen are, respectively, given by and
By constructing an alternative formulation of the four constraints, such that [CO_{2}(g)],[H^{+}],[CaCO_{3}] (denoted X, Y and Z, respectively) and [H_{2}O] are expressed as functions of the remaining five species and the initial concentrations, the system (2.1a–f) and (2.3) simplifies to the following system of five coupled equations: 2.5a 2.5b 2.5c 2.5d 2.5e where 2.6a 2.6b 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 k_{1} and k_{−1} are estimates based on CO_{2} 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 wellmixed scenario (Rochelle et al. 2004) and similar in speed to reaction (1.1c). Also, k_{5} and k_{−5} are estimates based on reaction (1.1e) being relatively slow, with K_{5}≈1660 at 25^{°}C (Drever 2002). As CaCO_{3} dissociates to one molecule of Ca^{2+} and one molecule of CO, S (in M) is given by the square root of the solubility product constant, which, for CaCO_{3} in pure water at 25^{°}C, is 4.96×10^{−9} (Tro 2008).
The behaviour of the system will be illustrated with initial concentrations [Ca^{2+}]_{0}=10 mg l^{−1}=2.50×10^{−4} M, [HCO 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 ), 2.7 and so 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 [CO_{2}(g)] to 0.065 M. This is sufficient to increase [CO_{2}(aq)] to carbon dioxide’s solubility limit, which at 25^{°}C and 1 atm is about 0.033 M. Any more injected CO_{2}(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.
(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 . In table 4, we define 18 dimensionless parameters, denoted with Greek symbols, with the initial concentration of water and chosen to be the arbitrary reference concentration and time scale, respectively. These parameters are evaluated using data in tables 1 and 2 and [CO_{2}(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×10^{8}). This is a precursor to taking the distinguished asymptotic limit , with all hatted parameters remaining O(1) in the limit.
As the initial conditions (prior to CO_{2} injection) have been chosen to satisfy the equilibrium equations (as in equation (2.7)), reactions (1.1b–e) yield the following identities: 2.8 Substituting the dimensionless variables, parameters and identities (2.8) into equations (2.5) yields the following dimensionless system, containing 14 independent parameters: 2.9a 2.9b 2.9c 2.9d 2.9e The five variables all equal zero at .
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 [CO_{2}(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 CO_{2}(g) leads to increases in [CO_{2}(aq)], [H_{2}CO_{3}], [H^{+}] and [HCO]. The relative increase in [H^{+}] is much larger than the increase in [HCO], as reaction (1.1d) converts some of the generated HCO into CO and further H^{+}. The ratio of a reaction’s concentrations is constant for any equilibrium, so from reaction (1.1d), Therefore, a large increase in [H^{+}] (from [H^{+}]_{0} to [H^{+}]) and a small increase in [HCO] must be balanced by a (large) decrease in [CO] (as shown in figure 1a), as expected (Orr et al. 2005). Similarly, from reaction (1.1e), Therefore, a drop in [CO] must be balanced by an increase in [Ca^{2+}] (as shown in figure 1a) and a decrease in [CaCO_{3}]. Therefore, CO_{2} cannot be sequestered as CaCO_{3} through these reactions alone.
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 exponent 3.1 provides a useful measure of the rate of change of each variable (as 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 t_{1}–t_{6} in figure 1, and discussed in more detail below): (i) CO_{2}(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) H_{2}CO_{3} and HCO reach an approximately linear evolution in time (Ψ_{D} and Ψ_{F} approach 1), with CO (G) varying quadratically, and Ca^{2+} (H) cubically in time, although the changes are all still negligible compared with the initial concentrations, (iii) and (iv) the rate of change of CO declines rapidly (Ψ_{G} falls to zero), (v) H_{2}CO_{3}, HCO and CO (D, F and G) reach equilibrium while Ca^{2+} evolves linearly in time (Ψ_{H}≈1) and (vi) Ca^{2+} 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 CO_{2}(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 leadingorder behaviour of the system at each time scale in the distinguished limit (hatted parameters in table 4 remaining O(1)).
The different leadingorder 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 largetime limits of the leadingorder approximations, given by equations (A3a), (A23a–c) and (A26), give estimates of the equilibrium concentrations, which in dimensional terms are 5.1a 5.1b 5.1c 5.1d 5.1e These nontrivial 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 CO_{2}(g) that is sequestered (i.e. dissolves into the water) is controlled only by the dissolution constant K_{1}=k_{1}/k_{−1}; the further process of dissolved CO_{2} being removed by reacting with water is less important. In the present example, with K_{1}=1, half of the injected CO_{2}(g) dissolves in the water. The proportion of injected and dissolved CO_{2} that then reacts with the water to form carbonic acid is given simply by the ratio of this reaction’s rate constants (K_{2}; 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 leadingorder 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), as . Therefore, a decrease in must be accompanied by a decrease in [CaCO_{3}]. Therefore, the amount of CO_{2} 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 and [CaCO_{3}]_{0} are all negligible compared with (which is true for these parameter values with error O(α^{−1/4})), we find 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 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 wellmixed scenario, CO_{2}(aq) reaches equilibrium within O(10^{−11} s), H_{2}CO_{3}, HCO and CO reach equilibrium within O(10^{−2} s) and Ca^{2+} reaches equilibrium in O(10^{3} 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.1a–d). 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 wellmixed scenario is given by equation (A25),
The results obtained above remain valid while the six time scales remain distinct, i.e. 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 CO_{2}(g) into a dilute system, based on the mass action law. We have also used the method of matched asymptotic expansions to derive leadingorder 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 nontrivial 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 CO_{2} 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 , 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 CO_{2}(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 ratelimiting, particularly for the extremely rapid dissolution of CO_{2} (time scale 1). In this wellmixed system, however, the reaction with calcium is ultimately ratelimiting. The present model provides a useful foundation on which to build more sophisticated multiscale models, incorporating, for example, additional reactions that buffer pH, interactions with dissolved salt, catalysts or other metals (such as olivine (Mg_{2}SiO_{4}) and serpentine (Mg_{3}Si_{2}O_{5}(OH)_{4})), as well as spatial effects relevant to aquifers or other applications.
Acknowledgements
This study was supported by EPSRC grant EP/F012098/1.
Appendix A
A.1. Time scale 1:
Writing , , , , and gives the following leadingorder behaviour of equations (2.9) on this time scale as : A 1a A 1b A 1c A 1d A 1e with at . equilibrates rapidly, with changes in the remaining variables slaved sequentially to this. Equations (A 1) have the solutions A 2a A 2b A 2c A 2d A 2e where . 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 reaches O(α^{1/2}), i.e. , at which time , , , and . The leadingorder largetime limits of these variables, valid for (i.e. between time scales t_{1} and t_{2}), form the initial conditions for the next time scale. These limits are, in dimensionless unscaled variables, A 3a A 3b A 3c A 3d A 3e
A.2. Time scale 2: = O(α^{−1})
We now define O(1) variables using , , , , and to yield the following leadingorder reduction of equations (2.9): A 4a A 4b A 4c A 4d A 4e Therefore, . Also, d()/d, implying that (the constant of integration is zero, by matching with equations (A 3)). Substituting this expression for into equation (A 4c), solving for and matching the smalltime limit of this expression with equation (A 3c), gives where . Likewise, solving equations (A 4d,e) and matching with equations (A 3d,e) give and Again, looking at the largest neglected terms in equations (A 4), we find that these approximations are valid until reaches O(α^{1/4}), at which time , , , and . The leadingorder largetime limits of these variables, valid for (i.e. between time scales t_{2} and t_{3}), form the initial conditions for the next time scale. These limits are (A 3a) plus A 5a A 5b A 5c A 5d broadly consistent with the exponents Ψ_{B}≈0, Ψ_{D}≈1, Ψ_{F}≈1, Ψ_{G}≈2 and Ψ_{H}≈3 for in figure 1c.
A.3. Time scale 3:
On this time scale, we define O(1) variables using , , , , and , to yield the following reduction of equations (2.9): A 6a A 6b A 6c A 6d with as in equation (A 4a). As the lefthand sides of equations (A 6a,b) are O(1), the O(α^{1/4}) terms on the righthand sides must sum to zero, i.e. A 7
In the limit , the leadingorder behaviour of equations (A 6) on this time scale is given by equations (A 6c,d) plus A 8aand A 8b
Therefore, d()/d, implying (matching with equations (A 5a,b)). Substituting this into equation (A 7) gives A 9aand A 9b Substituting equation (A 9b) into equation (A 6c), and solving for gives where and . is given by The constants of integration for and can be deduced by matching with equations (A 5c,d). The neglected terms in equations (A 6) indicate that these approximations are valid until reaches O(α^{1/4}), at which point, is O(1), , and are O(α^{1/4}), and is O(α^{1/2}). The leadingorder largetime limits of these variables, valid for (i.e. between time scales t_{3} and t_{4}), are given by equations (A 3a), (A 5a,b), A 10 consistent with Ψ_{G}≈1, Ψ_{H}≈2 in figure 1c (although these exponents do not exhibit clear plateaux for ).
A.4. Time scale 4:
On the next time scale, we define O(1) variables using , with , , , , and . 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 A 11aand A 11b
Then, in the limit , the leadingorder behaviour is given by A 12a A 12b A 12c A 12d Therefore, d()/d, and so . This, together with equation (A 7), gives equations (A 9). Substituting equation (A 9b) into equation (A 11b) gives A 13a and substituting equation (A 13a) into equation (A 12d) gives A 13b which match with equations (A 10) for . These approximations are valid until reaches O(α^{1/4}), at which point and are O(1), and , and are O(α^{1/4}). The largetime limits of these variables, valid for (i.e. between time scales t_{4} and t_{5}), and consistent with Ψ_{G}≈0, Ψ_{H}≈1, are given by equations (A 3a), (A 5a,b), plus A 14
A.5. Time scale 5:
On this time scale, we define O(1) variables using , with , , , , and . The O(α^{3/4}) terms in the ODEs for D, F and G give A 15aand A 15b the O(α^{1/2}) terms give A 16aand A 16b and the O(α^{1/4}) terms give A 17aand A 17b
Then, the leadingorder behaviour on this time scale is given by A 18a A 18b A 18c A 18d
Therefore, d()/d, and from equation (A 15a), the evolution of and on this time scale is given by A 19aand A 19b and so has the implicit solution A 20 where . From equations (A 15b), (A 16b) and (A 18d), respectively, the evolution of and is given by A 21a A 21b A 21c
Thus, , and rapidly equilibrate for , and evolves linearly in time, consistent with figure 1c. As as , the smalltime limits of equations (A 19) and (A 21b) are A 22a A 22b A 22c ensuring matching with equations (A 9) and (A 14). These approximations are valid until reaches O(α^{1/2}), at which point, , , and are O(1), and is O(α^{1/2}). The largetime limits for D and F, valid for (i.e. between time scales t_{5} and t_{6}), are A 23aand A 23b The corresponding largetime limits of G and H are A 23c and A 23d (with F satisfying equation (A 23b)), consistent with Ψ_{B}=Ψ_{D}=Ψ_{F}=Ψ_{G}=0 and Ψ_{H}≈1.
A.6. Time scale 6:
On this final time scale, we define O(1) variables , with , , , , and . Substituting these expansions into equations (2.9) gives identities (A15–A17). Also, O(1) balances are given by Therefore, the solutions to D, F and G on this time scale are given by equations (A 23a–c). Also, A 24 and the solution to equation (A 24), which includes the first correction term for G and matches with equation (A 23d), is A 25 Thus, for , the largetime limits are equations (A 3a), (A 23a–c) and 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 largetime limits, to give composite asymptotic formulae. Otherwise, the expressions have been patched together to provide a good fit to the full simulations as follows: and where is given implicitly by Ω≡γ+ν, Δ≡ψ+ψδ/μ and .
Footnotes
 Received July 5, 2009.
 Accepted November 11, 2009.

This is an openaccess 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.
 © 2009 The Royal Society