2.1 Streamfunction formulation and non-dimensionalisation

Under the assumption that the rate of change of porosity is much slower than the transport of dissolved solute, (2.1) may be approximated as

This is a standard simplification (see Ritchie & Pritchard (Reference Ritchie and Pritchard2011) and references therein), which allows us to define a streamfunction $\hat{\unicode[STIX]{x1D713}}(\hat{x},\hat{z},\hat{t})$ so that $\hat{u} =-\unicode[STIX]{x2202}\hat{\unicode[STIX]{x1D713}}/\unicode[STIX]{x2202}\hat{z}$ and ${\hat{w}}=\unicode[STIX]{x2202}\hat{\unicode[STIX]{x1D713}}/\unicode[STIX]{x2202}\hat{x}$ .

We now define dimensionless variables as

where $T^{\prime }$ and $C^{\prime }$ are perturbations to the base-state temperature and concentration fields respectively. We can then eliminate $\hat{p}$ to obtain the dimensionless governing equations

subject to the boundary conditions

The dimensionless parameters ${\mathcal{R}}_{C}$ , ${\mathcal{R}}_{T}$ , $Le$ , $\unicode[STIX]{x1D706}$ , $k_{0}$ and $\unicode[STIX]{x1D6FF}$ are defined as

The parameters ${\mathcal{R}}_{C}$ and ${\mathcal{R}}_{T}$ can be either positive or negative, with positive values encouraging instability. We can define

where ${\mathcal{R}}_{C}^{solute}$ and ${\mathcal{R}}_{T}^{temp}$ are the Rayleigh numbers that emerge naturally from the pure solutal and pure thermal problems respectively. The Lewis number $Le$ is known to be greater than unity for liquids. The specific heat ratio $\unicode[STIX]{x1D706}$ satisfies $0<\unicode[STIX]{x1D706}\leqslant 1$ , and for the purposes of this study we will take it to be constant. The dimensionless dissolution/precipitation rate (Damköhler number) $k_{0}\geqslant 0$ ; estimates for the rates of geochemical processes may vary by many orders of magnitude (Phillips Reference Phillips2009, § 2.8). Finally, the magnitude of the density ratio $\unicode[STIX]{x1D6FF}$ , which determines the rate of matrix evolution relative to the solutal and thermal dynamics, may be assumed to be much less than unity since mass concentrations of solute are typically much smaller than the density of the solid mineral (Phillips Reference Phillips2009, § 2.8). In contrast to Ritchie & Pritchard (Reference Ritchie and Pritchard2011), $\unicode[STIX]{x1D6FF}$ is not necessarily positive since we may take ${\hat{C}}_{1}\lessgtr {\hat{C}}_{0}$ .

3.1 Dissolution and precipitation but no porosity evolution

Before considering the full linear stability problem, it is useful to recapitulate some general features of the stability diagram in the $({\mathcal{R}}_{C},{\mathcal{R}}_{T})$ plane when $\unicode[STIX]{x1D6FF}=0$ , as presented in detail in Pritchard & Richardson (Reference Pritchard and Richardson2007).

With $\unicode[STIX]{x1D6FF}=0$ , the terms in (3.13) that are linear in $q$ disappear, and it becomes a cubic in $q^{2}$ with solutions $q=\pm \text{i}n\unicode[STIX]{x03C0}$ for $n\in \mathbb{N}$ . The solvability condition for this simpler system is given by

where we define $\unicode[STIX]{x1D70E}=\unicode[STIX]{x1D70E}_{0}$ as the growth rate for $\unicode[STIX]{x1D6FF}=0$ and

This is identical to the solvability condition given by equations (3.7) and (3.8) of Pritchard & Richardson (Reference Pritchard and Richardson2007). Instability will occur whenever $\text{Re}(\unicode[STIX]{x1D70E}_{0+})>0$ with $n\in \mathbb{N}$ , i.e. whenever (i) $b<0$ or (ii) $b>0$ and $c<0$ (so that $\sqrt{b^{2}-4ac}>b$ ). In all cases, the system is most unstable for $n=1$ , so only these results are presented below.

For a given wavenumber $\unicode[STIX]{x1D6FC}$ , we define the critical values ${\mathcal{R}}_{T}={\mathcal{R}}_{T,0}^{b}({\mathcal{R}}_{C},\unicode[STIX]{x1D6FC})$ on which $b=0$ , and ${\mathcal{R}}_{T}={\mathcal{R}}_{T,0}^{c}({\mathcal{R}}_{C},\unicode[STIX]{x1D6FC})$ on which $c=0$ . There are then two boundaries of interest: (i) the boundary ${\mathcal{R}}_{T}={\mathcal{R}}_{T,0}^{b\,min}=\min _{\unicode[STIX]{x1D6FC}}{\mathcal{R}}_{T,0}^{b}({\mathcal{R}}_{C},\unicode[STIX]{x1D6FC})$ above which there exist $\unicode[STIX]{x1D6FC}$ such that $b<0$ and (ii) the boundary ${\mathcal{R}}_{T}={\mathcal{R}}_{T,0}^{c\,min}=\min _{\unicode[STIX]{x1D6FC}}{\mathcal{R}}_{T,0}^{c}({\mathcal{R}}_{C},\unicode[STIX]{x1D6FC})$ above which there exist $\unicode[STIX]{x1D6FC}$ such that $c<0$ . (The minima occur at $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x03C0}$ for $k_{0}=0$ , but this is not true in general; see Pritchard & Richardson (Reference Pritchard and Richardson2007).) The unstable region is then the union of the regions ${\mathcal{R}}_{T}>{\mathcal{R}}_{T,0}^{b\,min}$ and ${\mathcal{R}}_{T}>{\mathcal{R}}_{T,0}^{c\,min}$ . Figure 2 illustrates the evolution of these stability boundaries in the $({\mathcal{R}}_{C},{\mathcal{R}}_{T})$ plane as $k_{0}$ is increased from $0$ .

Figure 2.

Linear stability boundaries for $\unicode[STIX]{x1D706}=0.25$ , $Le=10$ and (a) $k_{0}=0$ , (b) $k_{0}=5$ , (c) $k_{0}=20$ . The solid lines represent stability boundaries on which $\unicode[STIX]{x1D70E}_{0+}=0$ , the dashed lines represent stability boundaries on which $\text{Re}(\unicode[STIX]{x1D70E}_{0+})=0$ but $\text{Im}(\unicode[STIX]{x1D70E}_{0+})\neq 0$ , and the dotted lines represent the continuation of ${\mathcal{R}}_{T,0}^{b\,min}$ and ${\mathcal{R}}_{T,0}^{c\,min}$ above the boundary. The numbered arrows in (b) and (c) represent the transects taken when examining the stability of the full problem. (Replotted from Pritchard & Richardson (Reference Pritchard and Richardson2007).)

When $k_{0}=0$ (figure 2 a), the stability boundary is given by ${\mathcal{R}}_{T,0}^{b\,min}$ in the solutally stable (SS) regime ${\mathcal{R}}_{C}<0$ and by ${\mathcal{R}}_{T,0}^{c\,min}$ in the solutally unstable (SU) regime ${\mathcal{R}}_{C}>0$ . In the SS regime, $\text{Im}(\unicode[STIX]{x1D70E}_{0+})\neq 0$ on the stability boundary, and so the onset of convection is by an oscillatory mode, whereas in the SU regime, $\text{Im}(\unicode[STIX]{x1D70E}_{0+})=0$ on the stability boundary, and the transition to instability is by a monotonically growing mode. Furthermore, in the SS regime, ${\mathcal{R}}_{T,0}^{b\,min}$ increases as ${\mathcal{R}}_{C}$ decreases: a greater destabilising temperature gradient is required to overcome a greater stabilising concentration gradient and allow the system to become unstable.

As the Damköhler number increases, the boundary ${\mathcal{R}}_{T,0}^{b\,min}$ moves up everywhere; meanwhile, ${\mathcal{R}}_{T,0}^{c\,min}$ moves upwards and rightwards in the SU regime ( ${\mathcal{R}}_{C}>0$ ) and downwards in the SS regime ( ${\mathcal{R}}_{C}<0$ ). The consequence is that when the Damköhler number is increased to $k_{0}=5$ (figure 2 b), the system becomes more stable in both regimes. However, when the Damköhler number is increased further to $k_{0}=20$ (figure 2 c), the system becomes more stable only in the SU regime. In the SS regime, the boundary ${\mathcal{R}}_{T,0}^{c\,min}$ has now passed entirely below the boundary ${\mathcal{R}}_{T,0}^{b\,min}$ , and so the transition to instability is by a monotonic rather than an oscillatory mode in both regimes. In particular, the stability boundary in the SS regime has now moved down from its original position, and the ‘reaction’ term now destabilises the system by weakening the stabilising solutal perturbations. Hence, a lower value of ${\mathcal{R}}_{T}$ is required for instability in the SS regime when $k_{0}=20$ than when $k_{0}=0$ .

The critical Damköhler number $k_{0}^{crit}$ at which the boundary ${\mathcal{R}}_{T,0}^{c\,min}$ passes entirely below the boundary ${\mathcal{R}}_{T,0}^{b\,min}$ is given by

3.2 Linear stability analysis of the full problem

When $\unicode[STIX]{x1D6FF}\neq 0$ , an analytical solution of (3.17) in terms of $\unicode[STIX]{x1D70E}$ is not available. Therefore, we followed the approach of Ritchie & Pritchard (Reference Ritchie and Pritchard2011) for solutal convection and investigated the stability problem numerically using a continuation method, tracking $\unicode[STIX]{x1D70E}$ as ${\mathcal{R}}_{T}$ was changed while keeping ${\mathcal{R}}_{C}$ and all other parameters fixed. A numerical solution of (3.17) with $\unicode[STIX]{x1D6FF}=0$ and for a large value of ${\mathcal{R}}_{T}$ was used as an initial guess.

The task of tracking each solution branch is numerically very laborious, largely because the analytical solutions for $q_{i}$ are prohibitively complicated and so nested numerical solutions of (3.13) and (3.21) are required. As with the solutal convection problem (Ritchie & Pritchard Reference Ritchie and Pritchard2011), we found that to track a single branch it was necessary to reduce ${\mathcal{R}}_{T}$ in steps no larger than $10^{-2}$ , with smaller steps required around the bifurcation points discussed below. Additionally, to track the bifurcation point accurately, we needed to further simplify the determinant (3.17). Using column operations, the determinant can be simplified to

where $\det (\unicode[STIX]{x1D63E})$ is omitted here for brevity (see Corson (Reference Corson2012, Appendix A.3)). When $\unicode[STIX]{x1D6FF}\neq 0$ ,

Therefore, to track each solution branch, we only require to solve $\det (\unicode[STIX]{x1D63E})=0$ . Nevertheless, each set of plots within figures 3 and 4 took around 12 h to produce on a desktop computer. This effectively precluded the thorough investigation of the full problem defined by (3.21); instead, we employed a less formal Galerkin approach in order to both search parameter space and provide insight into the solution structure.

Figure 3.

Behaviour of $\unicode[STIX]{x1D70E}$ for $k_{0}=5$ , $Le=10$ , $\unicode[STIX]{x1D706}=0.25$ , $\unicode[STIX]{x1D6FC}=3$ , and (a,b) ${\mathcal{R}}_{C}=-15$ , $\unicode[STIX]{x1D6FF}=-0.1$ close to $\unicode[STIX]{x1D6E5}^{b}=0$ , (c,d) ${\mathcal{R}}_{C}=30$ , $\unicode[STIX]{x1D6FF}=0.1$ close to $\unicode[STIX]{x1D6E5}^{c}=0$ . Panels (a,c) show $\text{Re}(\unicode[STIX]{x1D70E})$ , while (b,d) show $\text{Im}(\unicode[STIX]{x1D70E})$ . The solid lines are full numerical solutions to the linear problem, the dashed lines are Galerkin solutions and the dotted lines in (a,c) show the results for $\unicode[STIX]{x1D6FF}=0$ . The Roman numerals in (a,b) are to aid in identifying the corresponding real and imaginary parts of each $\unicode[STIX]{x1D70E}$ branch.

Figure 4.

Behaviour of $\unicode[STIX]{x1D70E}$ close to $\unicode[STIX]{x1D6E5}^{c}=0$ , for $k_{0}=20$ , $Le=10$ , $\unicode[STIX]{x1D706}=0.25$ , $\unicode[STIX]{x1D6FC}=3$ , and (a,b) ${\mathcal{R}}_{C}=-50$ , $\unicode[STIX]{x1D6FF}=-0.1$ , (c,d) ${\mathcal{R}}_{C}=10$ , $\unicode[STIX]{x1D6FF}=0.1$ , (e,f) ${\mathcal{R}}_{C}=50$ , $\unicode[STIX]{x1D6FF}=0.1$ . Panels (a,c,e) show $\text{Re}(\unicode[STIX]{x1D70E})$ , while (b,d,f) show $\text{Im}(\unicode[STIX]{x1D70E})$ . The solid lines are full numerical solutions to the linear problem, the dashed lines are Galerkin solutions and the dotted lines in (a,c,e) show the results for $\unicode[STIX]{x1D6FF}=0$ . The Roman numerals in (a,b) are to aid in identifying the corresponding real and imaginary parts of each $\unicode[STIX]{x1D70E}$ branch.

3.3 Galerkin approach to the linear stability problem

We start with the system of equations (3.7)–(3.10), subject to the boundary conditions

and to the usual non-triviality condition that the solution is not identically zero. Rearranging (3.10), we see that

and we can eliminate $\unicode[STIX]{x1D6F7}(z)$ from (3.9) to obtain the system

subject to the same boundary conditions.

We seek approximate solutions of the form

It should be noted that the boundary conditions are automatically satisfied and that the non-triviality condition has been imposed by normalising the $\sin (\unicode[STIX]{x03C0}z)$ component of $\unicode[STIX]{x1D6F9}(z)$ . In the limit $\unicode[STIX]{x1D6FF}=0$ , this approximation will become exact with $\unicode[STIX]{x1D6F9}_{2}=\unicode[STIX]{x1D6E9}_{2}=\unicode[STIX]{x1D712}_{2}=0$ .

With six unknowns, we can choose to satisfy six integral conditions. The natural ones are obtained by extracting the first and second Fourier sine components of (3.25)–(3.27),

for $n=1$ and for $n=2$ . This yields a system of six algebraic equations,

We can then eliminate $\unicode[STIX]{x1D6F9}_{2}$ , $\unicode[STIX]{x1D6E9}_{1}$ , $\unicode[STIX]{x1D6E9}_{2}$ , $\unicode[STIX]{x1D712}_{1}$ and $\unicode[STIX]{x1D712}_{2}$ , and find that $\unicode[STIX]{x1D70E}$ satisfies the sextic

where the coefficients $a_{i}$ are omitted here for brevity. It is straightforward to locate all of the complex roots of (3.40) numerically, for example using the fsolve command in Maple, and thus to track all of the solutions through parameter space.

3.4 Eigenfunction structure and instability mechanism

As well as providing the complex growth rate $\unicode[STIX]{x1D70E}$ , the Galerkin stability analysis also provides the quantities $W_{2}$ , $\unicode[STIX]{x1D6E9}_{1}$ , $\unicode[STIX]{x1D6E9}_{2}$ , $\unicode[STIX]{x1D712}_{1}$ and $\unicode[STIX]{x1D712}_{2}$ which define the corresponding vertical eigenfunction. To elucidate the mechanism involved in the instability, it is useful to examine these quantities.

We will consider three different scenarios, as described in the previous section: first, the combination of a destabilising solutal gradient with either a stabilising or a destabilising thermal gradient; second, a destabilising thermal gradient with a stabilising solutal gradient and a subcritical Damköhler number, so that the onset of instability is oscillatory; third, a destabilising thermal gradient with a stabilising solutal gradient and a supercritical Damköhler number.

In our formulation, when the thermal gradient is destabilising, $\hat{T}_{1}-\hat{T}_{0}<0$ , so $T^{\prime }<0$ in regions that become warmer than the base state and $T^{\prime }>0$ in regions that become cooler. Similarly, when the solutal gradient is stabilising, ${\hat{C}}_{1}-{\hat{C}}_{0}<0$ , so $C^{\prime }<0$ in regions of higher concentration and $C^{\prime }>0$ in regions of lower concentration. Therefore, in order to avoid confusion when considering the spatial structure of the fastest-growing eigenfunction, we define $\tilde{T}=-\text{sgn}({\mathcal{R}}_{T})T^{\prime }$ and $\tilde{C}=\text{sgn}({\mathcal{R}}_{C})C^{\prime }$ . With this reformulation, $\tilde{T}>0$ and $\tilde{C}>0$ will always correspond to regions of higher temperatures and concentrations respectively, while $\tilde{T}<0$ and $\tilde{C}<0$ will correspond to regions of lower temperatures and concentrations respectively. Using (3.24) and the definitions of $\tilde{T}$ and $\tilde{C}$ , we also define

so that $\tilde{\unicode[STIX]{x1D719}}>0$ always corresponds to regions in which the porosity is greater than $\unicode[STIX]{x1D719}_{0}$ , while $\tilde{\unicode[STIX]{x1D719}}<0$ always corresponds to regions in which the porosity is lower than $\unicode[STIX]{x1D719}_{0}$ .

4.1 Stabilising solutal gradient, destabilising thermal gradient, subcritical Damköhler number

In this and the subsequent subsection, we consider the combination of a stabilising solutal gradient with a destabilising thermal gradient. Here, we consider this combination together with a subcritical Damköhler number; we recall that in this case, the transition to instability is through a growing oscillatory mode.

Figures 11 and 12 show the evolution of the concentration perturbation $\tilde{C}$ , the temperature perturbation $\tilde{T}$ and the porosity $\unicode[STIX]{x1D719}$ for the left half of the domain, for ${\mathcal{R}}_{C}=-15$ , ${\mathcal{R}}_{T}=299$ and $k_{0}=5$ . These figures should be read in conjunction with figure 13(a), which shows the evolution of the maximum absolute velocity and the maximum and minimum porosity within the domain. In figure 11, one should note the different magnitudes of $\tilde{C}$ , $\tilde{T}$ and $\unicode[STIX]{x1D719}$ in each subset of plots. The solid lines in each plot are the streamlines, which correspond to evenly spaced values of the streamfunction in steps of (a,c) 0.04, (b) 0.003 and (d) 1.6. The flow direction is indicated by the arrow at $x=0.5$ in the plot of $\unicode[STIX]{x1D719}$ .

At $t=2$ (figure 11 a), all of the fields are in phase; that is, regions of $\tilde{C}>0$ and $\tilde{T}>0$ coincide with regions of lower porosity ( $\unicode[STIX]{x1D719}<\unicode[STIX]{x1D719}_{0}$ ), and regions of $\tilde{C}<0$ and $\tilde{T}<0$ coincide with regions of higher porosity ( $\unicode[STIX]{x1D719}>\unicode[STIX]{x1D719}_{0}$ ). At this time, we have upward flow in $\tilde{C}>0$ and $\tilde{T}>0$ regions, and downward flow in $\tilde{C}<0$ and $\tilde{T}<0$ regions. However, by $t=2.2$ (figure 11 b), the flow direction has switched, but the concentration, temperature and porosity fields have remained in phase. The concentration and temperature fields move a half-cell by $t=2.3$ (figure 11 c), so again we have upward flow in $\tilde{C}>0$ and $\tilde{T}>0$ regions, but now the concentration and temperature fields are out of phase with the porosity field. The fields remain out of phase, and by $t=3$ (figure 11 d), nonlinear effects have become significant, as seen by the plumes of high- and low-concentration perturbations.

The oscillatory onset can also be seen in the inset of figure 13(a). We can see that the maximum velocity grows as convection tries to establish itself, but then drops as the flow direction changes. This occurs twice before the phase shift at $t=2.3$ , after which the velocity grows rapidly as convection develops.

Since nonlinear effects enter the system so rapidly, there is not a chance to see the upward movement of the circulation and temperature fields and downward movement of the concentration and porosity fields predicted by the two-term Galerkin approximation to the linear stability analysis. However, the concentration and temperature fields are out of phase with the porosity field, as predicted.

Figure 13.

Maximum absolute velocity $|\boldsymbol{u}|$ (solid) and maximum and minimum porosity $\unicode[STIX]{x1D719}$ (dashed) throughout the whole domain for (a) $k_{0}=5$ , ${\mathcal{R}}_{C}=-15$ , ${\mathcal{R}}_{T}=299$ ; (b) $k_{0}=20$ , ${\mathcal{R}}_{C}=-15$ , ${\mathcal{R}}_{T}=245$ .

Figure 12 shows the evolution of the concentration perturbation $\tilde{C}$ , temperature perturbation $\tilde{T}$ and porosity $\unicode[STIX]{x1D719}$ for the left half of the domain. The solid lines in each plot are streamlines, which correspond to evenly spaced values of the streamfunction in steps of 2 at each time. To aid with the visualisation of the flow, additional streamlines are plotted with dashed lines. (These correspond to further evenly spaced values of the streamfunction, but the spacing is not the same at each time.) The flow direction is indicated by the arrow at $x=0.5$ in the plot of $\unicode[STIX]{x1D719}$ . Figure 12 should be read in conjunction with figure 13(a), which shows the evolution of the maximum absolute velocity, and the maximum and minimum porosity.

At $t=100$ (figure 12 a), the concentration and temperature fields are very similar to those at the onset of convection ( $t=3$ , figure 11 d), but the porosity field has begun to evolve, with porosity increasing where $\tilde{C}<0$ and $\tilde{T}<0$ and decreasing where $\tilde{C}>0$ and $\tilde{T}>0$ . Therefore, the porosity is reduced slightly near the top of the layer, and correspondingly the centres of the circulation cells have moved down towards the bottom of the layer. This region of reduced porosity is pierced by higher-porosity channels which correspond with the downflow regions of the circulation. Here, downflow is faster than upflow, as seen by the more tightly packed streamlines.

The porosity continues to evolve and by $t=500$ (figure 12 b), the region of lower porosity at the top of the domain has extended downwards, flattening the circulation cells even more. High-porosity channels still pierce the low-porosity region in areas of downward flow, but now upflow under the low-porosity barrier is faster than downflow.

Figure 13(a) shows the evolution of the maximum absolute velocity, and we can see that convection is established and maintained until $t\approx 750$ , when the porosity evolution has forced the circulation into too small a region to be maintained (figure 12 c). We then enter a regime of slow porosity evolution, which strives to reach a state where convection can be maintained again.

At $t=1000$ (figure 12 d), the circulation is still contained in the lower half of the domain, but additional counter-rotating cells have formed in the high-porosity region at the bottom of the channels. By $t=5000$ (figure 12 e), these new circulation cells have expanded upwards and now extend the full height of the domain, confining the original circulation cells even more strongly to the bottom of the domain. It should be noted that further pairs of counter-rotating circulation cells have formed at the top of the domain at the top of the high-porosity channel, and by $t=10\,000$ (figure 12 f), these circulation cells have started to dominate the flow.

From the maximum-porosity curve in figure 13(a), we can split the slow-evolution regime into two periods: one in which the maximum porosity decreases ( $t\approx 750$ –10 000) and one in which it increases ( $t\approx 10\,000$ –18 000). These periods correspond to times when the system is dominated by faster upflow and faster downflow respectively.

By $t=15\,000$ (figure 12 g), the new circulation cells have extended to cover most of the domain and the high-porosity channel is now more evident. This channel continues to be eroded ( $t=17\,500$ , figure 12 h) until vigorous convection can be re-established ( $t=17\,600$ , figure 12 i). The downward movement of the circulation cells continues, however, and by the end of the simulation ( $t=20\,000$ , figure 12 j), the very vigorous convection is localised in high-porosity regions near the bottom of the domain (note the different scale on the porosity plot).

The maximum porosities of $\unicode[STIX]{x1D719}=0.8$ observed at the end of this simulation may not be physically realistic in the context of naturally occurring porous media, where $\unicode[STIX]{x1D719}$ rarely exceeds 0.6 (Nield & Bejan Reference Nield and Bejan2006, § 1.2). At such high porosities, extensions of Darcy’s law (such as the Brinkman equation) would be required (Nield & Bejan Reference Nield and Bejan2006, § 1.5). This high porosity is a result of our assumption that the porous medium is entirely soluble. We would expect a more complete model to predict lower maximum porosities; however, our results strongly suggest that even in such a model, complete depletion of the soluble component from parts of the rock matrix may be expected.

Figure 14.

Evolution for ${\mathcal{R}}_{C}=-15$ , ${\mathcal{R}}_{T}=245$ and $k_{0}=20$ . The shading indicates the solutal concentration perturbation $\tilde{C}$ (left), the temperature perturbation $\tilde{T}$ (centre) and the porosity $\unicode[STIX]{x1D719}$ with its horizontal average $\overline{\unicode[STIX]{x1D719}}$ (right). The solid lines in each plot are streamlines, which correspond to evenly spaced values of the streamfunction in steps of 1 at each time. Additional streamlines are plotted with dashed lines. (These correspond to further evenly spaced values of the streamfunction, but the spacing is not the same at each time.) The flow direction is indicated by the arrow at $x=0.5$ in the plot of $\unicode[STIX]{x1D719}$ .

4.2 Stabilising solutal gradient, destabilising thermal gradient, supercritical Damköhler number

We now consider the combination of a stabilising solutal gradient with a destabilising thermal gradient together with a supercritical Damköhler number, so that the transition to instability is direct. Figure 14 shows the evolution of the concentration perturbation $\tilde{C}$ , temperature perturbation $\tilde{T}$ and porosity $\unicode[STIX]{x1D719}$ for the left half of the domain, for ${\mathcal{R}}_{C}=-15$ , ${\mathcal{R}}_{T}=245$ and $k_{0}=20$ . The solid lines in each plot are streamlines, which correspond to evenly spaced values of the streamfunction in steps of 1 at each time. To aid with the visualisation of the flow, additional streamlines are plotted with dashed lines. (These correspond to further evenly spaced values of the streamfunction, but the spacing is not the same at each time.) The flow direction is indicated by the arrow at $x=0.5$ in the plot of $\unicode[STIX]{x1D719}$ . Figure 14 should be read in conjunction with figure 13(b), which shows the evolution of the maximum absolute velocity, and the maximum and minimum porosity.

Convection is apparent by $t\approx 2$ (not shown). At this time, the fields are out of phase, with regions of $\tilde{C}>0$ and $\tilde{T}>0$ coinciding with regions of higher porosity ( $\unicode[STIX]{x1D719}>\unicode[STIX]{x1D719}_{0}$ ), and regions of $\tilde{C}<0$ and $\tilde{T}<0$ coinciding with regions of lower porosity ( $\unicode[STIX]{x1D719}<\unicode[STIX]{x1D719}_{0}$ ). Once convection has been established, the flow very quickly becomes dominated by nonlinear effects. Therefore, as in the subcritical case, there is no chance to see the vertical displacement predicted by the Galerkin analysis. However, the concentration and temperature fields are out of phase with the porosity field, as predicted.

Since the diffusion of solute is much slower than the diffusion of heat, the concentration perturbations are advected with the flow. This is the picture that has emerged by $t=100$ (figure 14 a), when we see plumes of low concentration being drawn towards the bottom of the domain, and those of high concentration being drawn towards the top of the domain, by the circulation. The porosity field continues to evolve, with porosity increasing in regions that coincide with $\tilde{C}<0$ and $\tilde{T}<0$ , and decreasing in regions that coincide with $\tilde{C}>0$ and $\tilde{T}>0$ . Therefore, the porosity is reduced near the top of the layer, and correspondingly the centres of the circulation cells have moved downwards. However, this region of low porosity is pierced by higher-porosity channels in which fluid is carried towards the bottom of the layer. Furthermore, downflow is faster than upflow, as can be seen by the more tightly packed streamlines in these regions.

From figure 13(b), we can see that the maximum velocity increases as convection is established, and then suddenly convection is halted at $t\approx 130$ . During this time, the porosity contrast has increased substantially (figure 13 b), but once convection is halted, we enter a regime of slow porosity evolution. This is identical behaviour to that seen for a subcritical Damköhler number in § 4.1, although in this case the period of vigorous convection is much shorter.

Figure 14(b) shows the fields at $t=130$ just as the convection is halted. By this time, both the concentration and the temperature perturbations are greatly reduced, and, although the porosity field looks similar to that at $t=100$ , the porosity contrast between the regions of high and low porosity has increased. As a result of the region of low porosity extending across the top half of the domain, the circulation has been confined closer to the bottom of the layer. Consequently, the vigorous convection cannot continue, resulting in this much slower circulation.

By $t=500$ (figure 14 c), the circulation has been confined even closer to the bottom of the domain, and there is now no circulation in the top half of the domain. Small additional circulation cells have appeared at the bottom of the domain in the high-porosity regions, and these grow in size ( $t=1000$ , figure 14 d), until by $t=2000$ (figure 14 e), they extend the full height of the domain. It should be noted that throughout this evolution, the region of high porosity at the bottom of the layer has become thinner, so now most of the layer is dominated by lower porosity. However, channels of slightly higher porosity still pierce the domain, but these now contain fluid moving upwards towards the top of the domain. Furthermore, the top of each channel now contains a pair of counter-rotating circulation cells.

By $t=5000$ (figure 14 f), the circulation cells that appeared at the top of the higher-porosity channels have increased in size, and, in the region of downflow, the porosity has increased. These circulation cells continue to expand downwards while the lower cells shrink (figure 14 g), until by $t=7000$ (figure 14 h), the upper cells almost fill the domain.

From the maximum-porosity curve in figure 13(b), we can see that the slow evolution can be split into three regions: from $t\approx 130$ –2000 the maximum porosity remains constant, then from $t\approx 2000$ –6000 it decreases, before increasing again from $t\approx 6000$ . These regions correspond to periods during which the system is dominated by faster downflow in high-porosity regions, then faster upflow, and finally faster downflow respectively. A similar pattern was seen in § 4.1 for a subcritical Damköhler number.

By $t=9000$ (figure 14 i), the centres of the circulation cells have moved downwards, and the porosity has increased further in the downflow channel. The porosity continues to increase until, by $t=10\,000$ (figure 14 j), vigorous convection is able to return.