# A model of carbon dioxide dissolution and mineral carbonation kinetics

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

## 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 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: 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 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): 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 CaCO3 and by imposing mass-conservation 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.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 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.1hj) 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 [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: 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 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 CO, 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, [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 [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 . 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 [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 , 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: 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 [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 [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 [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 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 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 HCO reach an approximately linear evolution in time (ΨD and ΨF approach 1), with CO (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 CO declines rapidly (ΨG falls to zero), (v) H2CO3, HCO and CO (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 (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 5.1a 5.1b 5.1c 5.1d 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), as . Therefore, a decrease in 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 and [CaCO3]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 well-mixed scenario, CO2(aq) reaches equilibrium within O(10−11 s), H2CO3, HCO and CO 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),

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

## Acknowledgements

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

## Appendix A

### A.1. Time scale 1:

Writing , , , , and gives the following leading-order 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 leading-order large-time limits of these variables, valid for (i.e. between time scales t1 and t2), 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 leading-order 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 small-time 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 leading-order large-time limits of these variables, valid for (i.e. between time scales t2 and t3), 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 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. A 7

In the limit , the leading-order 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 leading-order large-time limits of these variables, valid for (i.e. between time scales t3 and t4), 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 leading-order 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 large-time limits of these variables, valid for (i.e. between time scales t4 and t5), 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 leading-order 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 small-time 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 large-time limits for D and F, valid for (i.e. between time scales t5 and t6), are A 23aand A 23b The corresponding large-time 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 23ac). 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 large-time limits are equations (A 3a), (A 23ac) 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 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: and where is given implicitly by Ωγ+ν, Δ≡ψ+ψδ/μ and .