Effects of hydrodynamic dispersion on the stability of buoyancy-driven porous media convection in the presence of first order chemical reaction

Carbon dioxide (CO2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {CO}_2$$\end{document}) sequestration is one of several long-term solutions suggested to decrease the amount of greenhouse gases in the atmosphere. Among different methods of carbon dioxide sequestration, the dissolution of CO2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\text {CO}_2$$\end{document} in deep saline aquifers is considered one of the most effective. A significant number of studies are currently being carried out to provide a good understanding of the physical mechanisms involved in this type of storage. The present work focuses on the hydrodynamic part of the problem: setting a model for carbon dioxide-loaded flows in an idealised two-dimensional geometry. It considers the impact of hydrodynamic dispersion in porous media on the development of convective instabilities. Particular attention is paid to the mathematical form of the dispersion tensor widely used in porous media studies, and a new type of bifurcation is investigated. We show that the analysis of bifurcations from the no-flow steady-state solution is a continuous but non-smooth problem, which is a key feature of the analysis. Although the problem is non-smooth, it is also shown that the basic behaviours of linear stability analysis are observed in its solution.


Introduction
Storage of CO 2 in underground geological formations is a potential means of limiting greenhouse gas emissions to the atmosphere whilst continuing the use of fossil fuels [1]. Storage takes place through structural, residual and solubility trapping and through the reaction of CO 2 with minerals in the storage site. The dissolution of CO 2 into the resident brine of the storage site causes an increase in the density of the mixture, which promotes the sinking of the enriched brine-CO 2 mixture whilst reacting with local rock minerals to become a solid carbonate. The resulting convection enhances mixing and reduces reservoir mixing times from thousands to hundreds of years [2], a process commonly referred to as convection enhanced dissolution.
To investigate this problem, researchers have given their attention to the classical problem of Rayleigh-Bénard convection in porous media, greatly improving its understanding. Of significant interest has been the study of the onset time of convection following the introduction of solute from an upper surface [3][4][5][6][7][8]. Associated experiments have exploited the analogy between porous media flow (described by the Darcy equation coupled to a solute transport equation) and flow in a Hele-Shaw cell [2,[9][10][11]. These researchers are interested in understanding a simplified and idealised case, where the problem is reduced to the motion of a fluid through a Hele-Shaw cell and the dispersive transport is based on pure molecular diffusion. The present work is also based on a similar simplified formulation, but instead of using a pure molecular diffusion model, we use a more realistic dispersion relation for flow in porous media. In the literature, a parallel is often drawn between fluid motion in porous media and in Hele-Shaw cells because Darcy's law applies in both cases and dispersive transport is to be observed. However, in this work, it is shown that some fundamental differences exist between these two systems as far as hydrodynamic dispersion is concerned.
In most experimental studies, Hele-Shaw cells are used to analyse the effect of a velocity-dependent dispersion tensor on dispersive transport flows [12]. Theoretical and experimental works [13][14][15][16] have led to a good understanding of dispersive transport in simple geometries, such as pipes and Hele-Shaw cells. In these geometries where a Poiseuille-type flow can be observed, an apparent dispersion of the solute added to the molecular diffusion appears owing to the non-homogeneity of the advection velocity profile. In this type of model, the longitudinal dispersion (along the mean flow direction) varies quadratically with the mean velocity, while the transverse dispersion is zero. In more complex geometries such as porous media, analytical estimation of an effective dispersion tensor cannot be easily obtained. Many numerical and experimental studies have been carried out to determine a proper form of the dispersion tensor in porous media, leading to a linearly dependent longitudinal and transverse dispersion with respect to the Darcy velocity [17,18], commonly referred to as hydrodynamic dispersion.
Hydrodynamic dispersion in porous media has been investigated to understand its impact on the onset and time evolution of the viscous and density fingering instabilities of miscible flow displacements [19,20]. However, limited attention has been given to the effects of hydrodynamic dispersion on the stability and bifurcations from the corresponding no-flow steady-state solution and resulting new flow patterns. A major issue encountered with this type of analysis is that the hydrodynamic dispersion is singular in the limit to the no-flow steady-state solution, in the sense that it is not a sufficiently differentiable function of the fluid velocity when the Darcy velocity tends to zero. Few papers deal with this singularity, but contestable conclusions are drawn regarding the existence of solution in the model framework [21]. One of the main objectives of this work is to investigate the effect of hydrodynamic dispersion on bifurcations from the no-flow steady-state condition of a convective enhanced dissolution process in porous media and their stability. We devote particular attention to the singularity of the dispersion tensor at the bifurcation points, i.e. in the limit where the Darcy velocity tends to zero (no-flow solution).
The present work is motivated by the development of interfacial instabilities arising during the displacement of miscible fluids that can occur in porous media during the process of CO 2 sequestration. However, we do not pretend to analyse the complete chemical interaction between the fluid mixture and the solid matrix. Instead, we present a detailed study of a simplified formulation of the chemical interaction between the fluid mixture and the substrate that is fixed in the underlying solid matrix. The chemical interaction between CO 2 and brine is very complex, with a sequence of reactions taking place leading to acidification and ultimately to a reaction with the host rock (substrate). The reaction timescales span many orders of magnitude [22], requiring significant computational resources for studies that attempt to integrate flow and chemistry. This has motivated fundamental studies where chemistry is represented in a highly simplified manner, either by a first-order reaction, where the dissolved CO 2 reacts with an abundant substrate [7,23], or a second-order reaction, where the substrate may be depleted at the leading edge of the advancing solute field [4,24]. Related studies, addressing the effect of chemical reaction on a Rayleigh-Taylor instability, have employed non-linear reaction terms that allow for chemical front propagation [25,26]. From these works it can be concluded that a chemical reaction with the substrate significantly changes the dynamics of the system through the introduction of a steady state that prevents the solution from penetrating deep into the underlying fluid. This process delays the onset of convection, shortens the evolution time and ultimately alters the spatial patterns of velocity and concentration of solute. In particular, Cardoso and Andres [27] show that, whilst the dissolution of CO 2 in carbonate rocks can persist for a long time, the chemical interactions in silicate-rich rocks may curb this transport drastically and even inhibit it altogether. Although in general chemical reactions delay the onset of convection, in some cases they can have a complete reverse effect depending on the type of reactions taking place in the mixture and surrounding medium. By considering the elementary reaction between three species, A + B → C, in which B, C and the lower phase solvent are insoluble in A but the density of the mixture is linearly proportional to the concentration of each species, Loodts et al. [28] show that a chemical reaction can either enhance or decrease the onset time of convection, depending on the type of density profile building up in time in the reactive solution.
As noted previously, we keep the description of the chemistry deliberately simple in order to understand how a linear reaction term influences the structure of the resulting dynamic system. In this type of model, the solute locally increases the solution density but undergoes a first-order reaction and turns into an inert product with no effect on the solution density. Previous authors have modelled the reaction between dissolved CO 2 and minerals trapped in a fixed substrate as second-order. If the reaction rate is sufficiently rapid, the reaction is confined to a thin zone at the leading edge of the solute field [24,29]. However, if the second-order reaction rate is slow compared to the convective timescale, and the solute is dilute in comparison to the substrate with which it reacts, the problem may instead be described using first-order reaction kinetics, at least over intermediate timescales [23,27,30,31]. In our analysis, we consider cases of deep aquifers, in the sense that the penetration depth of the reaction front is always less deep than the aquifer. Cardoso and Andres [27] conducted a stability analysis of a transient convection diffusive layer undergoing a first-order irreversible chemical reaction in a Hele-Shaw cell; they validated their results with laboratory experiments that mimicked natural geochemical reactions between CO 2 -rich acidic brine and the host rock on the dissolution and spreading of carbon dioxide in the subsurface.
Inevitably, in our first-order chemical reaction formulation, numerous features of real geological formations are neglected, in particular changes in permeability and porosity owing to dissolution or precipitation [32]. Studies by Xu et al. [33] and Zerai et al. [34] suggest that the porosity change caused by these reactions may be assumed to be negligible, as mineral volume change is too small to affect the dynamics of the system. However, for some geological formations, these changes could have some effect on the resulting flow patterns. For simplicity in this work, we neglect any effect due to porosity changes by the chemical reaction between the mixture and the substrate. A similar assumption was made in most of the aforementioned works that deal with the reaction between dissolved CO 2 and minerals trapped in a fixed substrate.

Definition of problem and governing equations
This work considers two-dimensional porous media convection in the presence of a first-order chemical reaction. This is typically described by three parameters (Rayleigh and Damköhler numbers and the domain aspect ratio). The equations governing the system are Darcy's law, the convection-dispersion-reaction equation and the continuity equations [23,30]. The solution density is considered linear in the solute concentration, and the change in density is assumed to be sufficiently small to apply the Boussinesq approximation. Under these considerations, the system of partial differential equations to be solved is where C, p and u = (u x , u z ) are respectively the concentration of solute dissolved in the solution, pressure and Darcy velocity. The parameters K , μ, D, n, α, g and ρ are the medium permeability, solvent viscosity, solute dispersion tensor, medium porosity, reaction rate of the solute, the gravitational acceleration and the density increase per unit concentration of solute. Complex microscopic flow in a porous medium leads to an apparent macroscopic dispersive transport, adding to the molecular diffusion, which is usually modelled by an apparent Fickian dispersion tensor depending non-linearly on the local Darcy velocity of the fluid u [17,18]: where D m is the molecular diffusion coefficient of the solute in the carrying fluid and β L and β T are respectively the longitudinal and transverse dispersion lengths such that β L ≥ β T ≥ 0 ( . always denotes the Euclidean norm on R 2 ). In [23], steady solutions (bifurcation branches from the no-flow solution) and the stability of the preceding system of partial differential equations in the case of pure diffusion, β L = β T = 0 are studied. The obtained structure of primary and secondary solution branches are compared to the reaction-free case in a box of finite depth, i.e. the Rayleigh-Bénard problem [35,36], for which secondary solution branches extend to large values of Ra. Inclusion of the reaction term truncates the solution branches bifurcating from the base state as Ra increases, and under certain flow conditions some modes states can be suppressed, which would otherwise exist in the reaction-free case. Kim and Choi [37] extended Ward et al.'s [23] stability analysis at the steady-state limit to the instability phenomenon for the time-evolving transient period. Their analysis shows that the onset of instability motion can occur during the transient period even if the system is stable at the steady state. In the analysis presented in [23] and [37], as well as in most works dealing with the study of stability and bifurcation maps for this type of problem, only molecular diffusion is considered instead of the hydrodynamic dispersion in porous media. One of the main objectives of this work is to investigate the effect of hydrodynamic dispersion on the stability of steady-state solutions and the corresponding bifurcation map of the enhanced dissolution problem considered in [23].
The problem domain is defined by a two-dimensional rectangle with horizontal length L and vertical depth 2H (x ∈ [0, L] and z ∈ [−H, H ]), where gravity acts in the vertical direction (Fig. 1). The aquifer is considered to be surrounded by impermeable rocks, and the solute (CO 2 ) is dissolved through its upper surface. On the upper boundary of the domain we enforce whilst at the three other boundaries, no-flux boundary conditions are applied:

Non-dimensional stream function formulation
By defining the velocity field in terms of its stream function formulation, , that automatically satisfies the continuity equation (3), i.e.
the pressure field can be eliminated from the above system of partial differential equations, resulting in a new set of governing equations for the unknowns ( , C), subject to the boundary conditions = 0 on every boundary.
Equations (10) and (11) and the corresponding boundary conditions are written in dimensionless form using the following dimensionless variables: The three dimensionless parameters that characterise the problem (without dispersion) are the Rayleigh and Damköhler numbers and the aspect ratio of the domain, which are defined here as and the dimensionless parameters characterising dispersion are Three timescales can be considered in this problem; they are defined by the preceding set of governing equations: the time needed by the solute to reach the bottom of the domain by simply sinking (buoyancy effects) τ b , the characteristic molecular diffusion timescale τ d and the reaction timescale τ r , where the characteristic buoyancy velocity is given by Darcy's law as K gC ρ/μ. Taking into account that the effective molecular diffusion coefficient in a porous medium must include the porosity, these timescales can be represented as from which it can be seen that the Rayleigh and Damköhler numbers can be defined in terms of the ratios between these timescales: Dropping the primes on the resulting set of dimensionless governing equations, we have the closed problem where the non-linear dispersive flux J = D∇C is given by subject to the boundary conditions = 0 on every boundary.
Integrating Eq. (16) over the problem domain and imposing the boundary conditions yields where we define the global flux F S accounting for the solute flux across the upper boundary and the total mass of dissolved solute M (both non-dimensional) as We see that the total mass of solute in the domain increases through solute dissolution across the upper boundary and decreases through the first-order chemical reaction. The efficiency of the convective enhanced dissolution can thus be evaluated by the computation of F S /Ra. dispersion on the bifurcations from the no-flow steady-state condition (trivial solution) of (15) and (16). We will focus on the study of neutral stability curves and bifurcations. The corresponding eigenvalue problems encountered in linear stability analysis are solved using MATLAB routines with finite-difference schemes. The steady-states bifurcation diagram is investigated numerically with a finite-element code based on Newton's method and arclength continuation.
Let φ = C be the state vector and λ = (Ra, Da, L) be the vector of parameters. For a given value of β T and β L , Eqs. (15) and (16) and the boundary conditions define a non-linear time-dependent problem of the form where and The study of steady states is therefore reduced to the non-linear problem We can numerically compute paths and branches of steady-state solutions of (29). For each steady state along the paths, a linear analysis enables us to determine the linear stability of those solutions. Provided that the Jacobian of F exists, the linearisation about a solution φ of (29) gives the following equation for a perturbation φ : subject to the homogeneous boundary conditions = 0 on every boundary. Setting we have Therefore, the linear stability can be deduced from the computation of the eigenvalues σ of the general problem given by G ≡ 0 0 0 1 and the Jacobian ∂F/∂φ. The corresponding eigenvectorsφ associated to eigenvalues with non-negative real part gives the shape of the instabilities that are likely to appear and the path to new steady states.
For steady state, (23) gives Therefore, the efficiency of the sequestration for different values of Ra (with Da fixed) and different steady-state solutions can be evaluated with the evaluation of either F S /Ra or M, with the numerical computation of the latter easier and more accurate.
The no-flow steady state φ 0 = 0 C 0 can be obtained analytically as which is a solution of (29) for every value of the parameter λ. This solution corresponds to a pure reaction-diffusion of the solute in the domain. For this solution, we have where is an operator related to hydrodynamic dispersion, provided that ∇ · J is -differentiable. The analytical evaluation of the dissolved mass of the pure reaction-diffusion process is straightforward and given by which increases with the aspect ratio L and decreases with Ra. M 0 is the lower bound of the capture efficiency since the uptake of the solute increases as soon as convection occurs in the domain.

The case of pure molecular diffusion
As a reference, we consider here the simple case where hydrodynamic dispersion is neglected, i.e. when β L = β T = 0, which implies and Here, the Fickian dispersion tensor is reduced to the molecular diffusion term. Setting with the wavenumber we can carry out a linear stability analysis by taking into account the homogeneous boundary conditions. In this case, the Jacobian can be rewritten as so that the general eigenvalue problem becomes where the eigenvectorsφ are only functions of z. We numerically compute the eigenvalues and eigenvectors of the preceding problem using a first-order finite-difference scheme to discretise the z-derivative operators appearing in the Jacobian and at the bottom boundary condition. Boundary conditions at the lower and upper boundaries are directly included in the discretised matrices. The problem is solved for a set of Rayleigh numbers and wavenumbers, and we look for the appearance of eigenvalues with non-negative real part. For L = π and Da = 0.1, increasing the Rayleigh number has a destabilising effect on the no-flow state with successive bifurcation of modes n = 2, 3, 4, . . . from the no-flow steady state. Figure 2 shows the neutral stability curve, where n ∈ N * corresponds to the different bifurcation modes, i.e. mode n corresponds to the onset of n convection cells in the domain. Figure 3 shows (a) streamlines and (b) concentration contours of mode 4, computed at the critical Rayleigh number of the perturbed unstable no-flow steady state (for more details about this case see [23]).
The Jacobian, and therefore the linear analysis, is invariant to the transformation showing in particular that even modes 2n in a domain of size L bifurcate at the same critical value of the parameter as modes n in a domain of size L/2. For instance, although mode 1 does not exist in a domain of size L = π , it does exist in a domain of size L/2 and bifurcates at the same critical value of mode 2 for L = π . It is also observed that the no-flux condition is satisfied between convection cells; therefore, a horizontal translation of the modes by one cell must also be a solution of the problem. Owing to the symmetric behaviour across the centreline of each of the recirculation cells (Fig. 3), the resulting solution obtained from the translation of one cell corresponds to a mirror image of the previous solution. By doing this translation, we find the solution on the other bifurcating branch.

Non-linear analysis
For a given set of parameters, the steady non-linear problem is solved using Newton's method, which ensures convergence towards both stable and unstable steady states [38]. To compute the paths of the steady states, arclength continuation was used, which enabled the solution to follow a path, even when turning points were encountered [38].
As observed in the linear stability analysis for even modes, the concentration profile and velocity field are a symmetric mirror image with respect to the line x = L/2, whilst odd modes introduce asymmetry in the concentration profile and velocity field. Two linear functionals will therefore be used to plot bifurcation diagrams, according to the type of mode, as   Fig. 3. The corresponding dissolved mass M given in Fig. 4b shows the capture efficiency due to the convective enhanced dissolution of this mode. In comparison to the no-flow steady state, the recirculation enhances mixing and significantly increases the uptake of the solute, which is of key importance in practical applications.

Effect of dispersion
In this case, we have If the operator O(φ 0 ) is well defined, it follows that O(φ 0 ) = 0 since the flux J is symmetric in . Therefore, in such a case the hydrodynamic dispersion would have no local effect on the bifurcations from the no-flow steady state (unchanged critical values of the parameters and the same unstable eigenvectors). However, the dispersion tensor considered in porous media is not a sufficiently differentiable function of the fluid velocity in φ 0 . We can indeed easily check that ∇ · J is not Frechet differentiable in φ 0 if β L = 0 or β T = 0. Although the model is singular at φ 0 , it is frequently used in studies of transport in porous media because it remains a good macroscopic empirical law and both theoretical and numerical studies have proved its consistency. Therefore, it is important to analyse the properties of this model near the no-flow steady state where the singularity leads to drastic changes in the bifurcations diagram and particularly in the stability regions of the different steady states.

The case of Hele-Shaw cells
Experimental studies of dispersive transport are often carried out in Hele-Shaw cells, where the fluid is confined between two close parallel plates. The dispersion tensor generally used in Hele-Shaw cells takes into account the Taylor-Aris effective dispersion effect (adding a longitudinal dispersion to the molecular diffusion due to the Poiseuille flow between the two plates of the cell). If b is the gap width of the Hele-Shaw cell and D m the molecular diffusion, the dispersion coefficient in cases of small but finite D m can be written in dimensionless form as [15,16] where R is a rotation matrix and A the Taylor-Aris dispersion: This leads to the frame-independent form In contrast to the dispersion tensor used in porous media, this tensor is quadratic in the velocity and therefore regular in the no-flow state. In this case, since D HS does not change when is changed to − , we have immediately O(φ 0 ) = 0, and the problem reduces to the same problem considered in Sect. 3.1. Dispersion, therefore, has no effect on the critical Rayleigh number at bifurcations from the no-flow solution in Hele-Shaw cells, although the bifurcation diagram changes. In the remainder of the article, we will consider the dispersion tensor for porous media.

Linear stability analysis
In the case of flow in porous media, the operator O(φ) is not defined at φ = φ 0 since the divergence of the dispersive flux J is not Frechet differentiable there, having Gateaux derivatives. Because of this singularity, the problem of stability of the no-flow solution can no longer be written as a general linear eigenvalue problem since the Jacobian does not exist at the no-flow solution.
Following the basic idea of linear stability analysis, we consider small perturbations about the no-flow solution, i.e.
with the vector φ as previously defined at the beginning of Sect. 3; then we have the set of equations for the perturbation φ = C : where and subject to the same homogeneous boundary conditions as previously mentioned. As expected, the right-hand side of the advection-dispersion-reaction equation is non-linear in the perturbation. This particularity is directly linked to the fact that the divergence of the dispersive flux J is Gateaux differentiable (homogeneous of degree one) but not Frechet differentiable in φ 0 , i.e. the operator O is not defined at φ 0 . However, this non-linear term is still positive homogeneous of O(1), so that it is still possible to look for perturbations of the form The preceding representation of the perturbed field shows that exponential growth of instabilities is likely to occur under unstable conditions, as will be seen subsequently in the transient simulations. Then we have The main aim in this type of analysis is to find values of Ra for which there exist non-trivial solutions to the preceding set of equations with a non-negative growth rate σ . We do this by analysing the solution of the original full non-linear problem and the transient evolution of the perturbation problem.
Contrary to the case without dispersion, the preceding system is not invariant to the horizontal transformation since in this case the symmetry at a cross section of the recirculation cells is lost (Fig. 6) and the resulting new solution from a cell translation is no longer a mirror image of the previous solution, suggesting the loss of symmetry in the bifurcation diagram and the disappearance of pitchfork bifurcations, as observed in the following non-linear analysis.

Non-linear analysis (non-smooth bifurcations)
In this case, bifurcations from the no-flow steady-state solution also exist but for different critical Rayleigh numbers, and they are no longer classic pitchfork bifurcations encountered in the regular case of a well-defined Jacobian. Figure 5 shows the non-smooth bifurcation of mode 4 with β L = 10, β T = 7, and L = π . In this figure, the dotted line below the solid inclined line with a positive angle corresponds to the symmetric image of a solid inclined line with a negative angle. The two solid inclined lines are the bifurcation branches from the no-flow solution in the case of dispersion, whilst the continuous dashed curve corresponds to the pitchfork bifurcation of the pure diffusion case. As can be seen from the figure, the bifurcation diagram is no longer symmetric with respect to the horizontal axis at the bifurcation point. Furthermore, we observe that the convective velocity is reduced by hydrodynamic dispersion and that the critical Rayleigh number is shifted to a higher value with respect to the case of pure diffusion.
Dispersion also affects the shape of the various recirculation modes. The cells seem to organise in pairs of contra-rotating convection rolls, as shown in Fig. 6, corresponding to different values of β L and β T . The shapes of the bifurcation diagram and solution branches reveal that two non-collinear solutions of the perturbation problem with a non-negative growth rate are to be found at the new critical Rayleigh number (one for each bifurcating branch, instead of only one unstable eigenvector in the case without dispersion). Indeed, just as in the case without dispersion, lateral no-flux conditions are observed between cells, so that any horizontal translation of the modes by one cell must be a solution of the problem. By doing this translation, we find the solution on the other bifurcating branch; however, in this case, the translation is not equivalent to a mirror image of the other solution. The efficiency of the mass dissolution M of both bifurcation branches is the same owing to the invariant property of the translation, just like the case without dispersion.
Similar to the cell translation approach, it is possible to construct new recirculation flow patterns in wider and narrower domains by copying or deleting cells at the edges of a given solution. For example, by copying at the end of the domain the first cell in the mode 4 diagram for L = π (Fig. 6), we obtain a mode 5 for L = (5/4)π , or by copying the first two cells we can construct a mode 6 for L = (3/2)π . Equivalently, by deleting at the end of the domain the last one or two cells, we obtain mode 3 for L = (3/4)π and mode 2 for L = (1/2)π . In particular, by adding or subtracting two cells, a bifurcation and flow behaviour identical to the original cell configuration is obtained. In our numerical solutions, this reproduction approach was used to verify the consistency of the results. It is important to observe that in the limit to the no-flow steady-state solution, i.e. when the velocity field tends to zero, the hydrodynamic dispersion tends to the molecular diffusion. However, the corresponding divergence of the dispersive flux J is singular in such a limit, in the sense that it is not Frechet differentiable at zero flow velocity, leading to continuous bifurcation curves with discontinuous first derivatives at the corresponding bifurcation points (Fig. 5).
Computing numerically the branches of the bifurcation diagram enables us to find the new critical Rayleigh numbers. The shift δ(β L , β T ) = Ra c − Ra β L =β T =0 c , with respect to the case of pure diffusion, was computed for different modes, taking L = π and Da = 0.1. Table 1 summarises the results for variable values of β ≡ β L and two different fixed ratios of β L /β T , equal to 1 and 10, where the ratio of β L /β T = 10 is considered a rule of thumb in hydrogeological models (see [20]). From this table it can be inferred that the effects of dispersion on the shift δ seem to depend quadratically on β, and tend to stabilise the bifurcation modes, i.e. δ > 0; however, it has an opposite effect on mode 2, making it bifurcate earlier when β L /β T = 10. Figure 7 shows the streamlines of mode 2 for β T = 1 and β L = 10, (a) on the upper branch and (b) on the lower branch.  As expected, when the aspect ratio L is multiplied or divided by a factor 2, modes show exactly the same dependency on dispersion as modes with twice as much or fewer convection cells, respectively. For instance, the destabilising effect on mode 2 observed when L = π and β L /β T = 10 (Fig. 11a) is also found on mode 1 for L = π/2 (Fig. 8) as well as on mode 4 for L = 2π . Equally, for L = π/2, mode 2 is found to show the same behaviour as mode 4 with L = π (Fig. 5), having dispersion in this case a stabilising effect on mode 2 for L = π/2 instead of a destabilising effect for L = π . We used this repetitive behaviour of different modes in different domains to confirm the unexpected destabilising effect of dispersion on mode 2 for L = π and β L /β T = 10, from which it can be inferred that dispersion has a destabilising effect on modes 2 ( j−1) for L = 2 ( j−2) π , with j as an integer number. All other modes found in this work on different domain sizes are stabilised by dispersion.
Because the system is Gateaux but not Frechet differentiable, a first-order highly non-linear term appears in the perturbation problem, causing the observed degeneracy of the bifurcation pattern. The non-linearity of the perturbation problem and the singular behaviour of the dispersive flux at the bifurcation points make it difficult to infer from the governing equations the observed behaviour of the bifurcation map and corresponding bifurcation points. For this reason, the reported bifurcation points and branches were found numerically using three different approaches: solving the transient form of the perturbation problem (Fig. 13), predicting them by the arclength continuation of Newton's solution of the steady-state non-linear problem at the bifurcation branches as the flow velocity tends to zero (Fig. 11a), as well as, where possible, when the solution of those repetitive modes has exactly the same bifurcation and flow behaviour at different domain sizes. Increasing dispersion changes the critical values of the parameters and affects the whole bifurcation diagram. Figure 8 shows paths of steady state after bifurcation towards mode 1 with increasing dispersion, for L = π/2, Da = 0.1, β L = β and β T = β/10. The crossings of bifurcation branches on the right of the horizontal axis do not account for any bifurcation from the no-flow steady state but are a projection effect in the figure due to the linear form used; if the mode 1 component is indeed null when crossing the axis, the steady state contains a non-zero mode 2 component (secondary bifurcations from mode 2). Figure 9 shows steady-state branches of bifurcations of modes 1, 2 and 3 at low Ra for L = π/2, Da = 0.1 and β L = 10, β T = 1. The symmetric component of the path of steady states after bifurcation towards mode 1 stems from the interaction with mode 2. As can be observed, in this case, dispersion leads to a new type of bifurcation.
Dispersion also tends to diminish the efficiency of convective enhanced dissolution. In fact, although the global flux F S presents an additional lateral component due to dispersion when β T = 0, the concentration gradient at the top of the domain is weaker owing to the added lateral dispersion process, thereby diminishing F S . Therefore, the intensity of the convection process is diminished by hydrodynamic dispersion, slightly diminishing the efficiency of mixing. Figure 10 shows the effect of dispersion on the total dissolved mass M after bifurcation towards mode 4, for L = π , Da = 0.1 and β L = 10, β T = 1, with a loss of efficiency up to 5 %. The obtained value of M in the case of a first-order chemical reaction without dispersion (dashed line in figure) is in agreement with previous simulations by Ward et al. [23] and Andres and Cardoso [31], which are also identical (at leading order) to the total dissolved mass of the classical one-sided Rayleigh-Bénard-Darcy convection in the case of a deep domain, i.e. without the effect of the lower boundary [11]. On the other hand, the obtained value of M with dispersion (thin solid line in figure) is consistent with the values reported by Hidalgo and Carrera [20], where the effect of dispersion is considered but chemical reaction is neglected. As noted by Ward et al. [30], the dissolution flux is determined by the boundary layer immediately beneath the distributed solute source and is independent of the chemical reaction between the mixture and the substrate since the length scale of the corresponding chemical reaction is larger than the boundary layer thickness, but, as evidenced by our results, it is affected by dispersion. The plumes beneath the boundary layer accommodate the flux provided by the boundary layer and distribute it according to the corresponding recirculation regions, which are affected by dispersion.

Regularised dispersion tensor
As mentioned earlier, a smooth and -symmetric dispersion tensor would have no effect on the bifurcations from the no-flow solution (such as the Hele-Shaw cell case). We tried to regularise the dispersion tensor near the no-flow state by making it locally quadratic in the velocity: where We found that the critical values of the Rayleigh numbers were indeed identical to the case without dispersion, but since the present regularisation does not affect the tensor far from the no-flow solution, the resulting bifurcation map is a mixture of the regular case (near the no-flow state) and of the dispersive case (far from the no-flow state). We observe an intermediary zone where the two maps smoothly match. Figure 11 shows the bifurcation of mode 2: (a) without regularisation (η = 0) and (b) when the dispersion tensor is regularised with η = 10 −5 , in both cases β T = 1, β L = 10 and L = π . In the singular case, i.e. η = 0, the critical Rayleigh number is smaller than the corresponding value for the non-dispersive case (destabiliser effect), as previously predicted. On the other hand, in the regularised case, the obtained critical Rayleigh number is, as expected, the same as in the case without dispersion.
The bifurcation diagram after regularisation shows a subcritical pitchfork bifurcation and two saddle-node bifurcations. The singular case can be seen as a limit case where the saddle-node bifurcations appear infinitely close to the no-flow steady state. However, as soon as η > 0, the no-flow solution remains stable until the subcritical pitchfork bifurcation happens. Therefore, the critical Rayleigh number of the singular case is not likely to be found using a limit process η → 0.
It is interesting to observe that the regularised solution in Fig. 11b also confirms the destabilising effect previously reported of the non-regularised mode 2 for L = π . For a far-field solution with respect to the bifurcation point, over a mode 2 branch, where the regularisation has no significant effect, the bifurcation point looks to be shifted with respect to the non-dispersive case, with δ < 0 (destabilising effect). However, close to the bifurcation point, near the no-flow solution, no shift (δ = 0) with respect to the non-dispersive case is observed, as expected, with a matching region between them defining the two saddle-node bifurcations. By regularising the dispersion tensor, we showed that the shift in the values of the critical Rayleigh number observed from the numerical solution of the non-linear problem and the concurrent change in the stability was not captured.

Transient behaviour
The preceding non-linear analysis suggests that dispersion is likely to change the value of the critical Rayleigh numbers at bifurcations from the no-flow steady state. In this section, we study the transient evolution of a perturbation near the no-flow solution to confirm whether exponential growth of instability occurs beyond the new critical values of the Rayleigh number. Furthermore, it is shown that the growth rate behaves smoothly near the critical Rayleigh number. The corresponding time-dependent problem is solved numerically using a version of Gear's method.

Perturbation problem
We first analyse the evolution of small perturbations φ to the no-flow steady-state solution that satisfies the governing equations (53) and (54)  coefficients of β L = 10, β T = 1 (Fig. 11a). For the dispersive case, this value of Ra is beyond the critical value, and perturbations are unstable with an exponential growth rate σ ≈ 0.13 (Fig. 12). This approach makes it possible to compute σ for different values of the Rayleigh number. In Fig. 13, for the dispersive case with β L = 10 and β T = 1, corresponding to the bifurcation of mode 2, it can be observed that near the critical Ra number, σ can be approximated by a straight line with a slope of ∂σ ∂Ra ≈ 0.17. This behaviour is encountered in smooth systems where the eigenvalues of the linear stability problem smoothly cross the imaginary axis at the bifurcation point, implying that the non-linear system of equations governing the time evolution of the perturbations is a smooth system. Consequently, in the dispersive case, the existence of the bifurcation can be characterised, even though the bifurcation itself is not smooth at the critical point, as previously shown.

Non-linear evolution
Here, we consider the full non-linear system (15), (16), taking into account the higher-order terms. We set the initial condition as As expected, for the same parameters as previously, the higher-order terms stabilise the solution at the bifurcating upper branch after an initial exponential transient growth with the same growth rate σ ≈ 0.13 as encountered in the perturbation problem (Fig. 14).

Conclusions
Although the stability of convection dissolution problems without dispersion are described by classical bifurcation theory and linear stability analysis, we show here that singularities occur at the no-flow solution when hydrodynamic dispersion is taken into account, owing to the lack of differentiability of the dispersion tensor at the no-flow condition. This singularity is typical of porous media and does not occur in simpler geometries such as Hele-Shaw cells, where the dispersive tensor has a smooth quadratic form. Despite the lack of differentiability of the dispersion tensor, we show that typical stability and bifurcation results of regular problems are also encountered in the corresponding stability analysis of the no-flow solution in porous media. However, the stability problem can no longer be written as a general linear eigenvalue problem since the corresponding Jacobian does not exist at the no-flow solution. Because the stability problem is no longer defined by a linear eigenvalue problem, degeneracy of the bifurcation instabilities occurs, causing non-smooth bifurcations. In this case, bifurcations from the no-flow steady state still exist, but for different critical Rayleigh numbers, and bifurcations are no longer the classic pitchfork bifurcations encountered in the regular case.
In the case of hydrodynamic dispersion, the evolution of a perturbation to the no-flow solution is governed by a non-linear system of partial differential equations. Numerical simulations of the transient evolution near the noflow solution show exponential growth of instability for values of Ra beyond the new critical values. Furthermore, we show that the growth rate behaves smoothly near the critical Rayleigh number, and consequently the analysis of the evolution of a perturbation can be used to characterise the existence of a bifurcation, even though the bifurcation itself is not smooth at the critical point. Hydrodynamic dispersion also changes the shape of the whole bifurcation diagram, which leads to new convective steady states having slightly less solute dissolution flux than the corresponding diffusion-only case. Consequently, hydrodynamic dispersion tends to diminish the efficiency of convective enhanced dissolution. The intensity of the whole convection process is weakened slightly, thereby diminishing the efficiency of mixing.
Hydrodynamic dispersion tends to stabilise bifurcation modes; however, it has a destabilising effect on modes 2 ( j−1) for L = 2 ( j−2) π , with j as an integer number. Owing to the non-linearity of the perturbation problem and the singular behaviour of the dispersive flux at the bifurcation points, it is difficult to infer from the governing equations the observed behaviour of the bifurcation map and corresponding bifurcation points. For this reason, the reported bifurcation points were found by solving the transient form of the perturbation problem, and at the same time their existence was verified by the arclength continuation of Newton's solution of the steady-state non-linear problem at the bifurcation branches as the flow velocity tends to zero. As expected, when the aspect ratio L is multiplied or divided by a factor of 2, modes show exactly the same dependency on dispersion as modes with twice as many or fewer convection cells respectively. We also used this repetitive behaviour of different modes in different domains to confirm the destabilising effect of dispersion on modes 2 ( j−1) for L = 2 ( j−2) π .
There is significant evidence that chemical reaction with the substrate changes the dynamics of the system through the introduction of a steady state that prevents the solute penetrating deep into the underlying fluid, delays the onset of convection, shortens the evolution time and alters the spatial patterns of velocity and concentration of solute. However, chemical reaction does not affect the critical value of the Ra number at bifurcation points, nor does it modify the solute dissolution flux, which is determined by a very small boundary layer immediately beneath the distributed solute source. Chemical reaction does, however, completely change the bifurcation map. On the other hand, the inclusion of hydrodynamic dispersion shifts the critical value of Ra and reduces the dissolution flux by redistributing the dissolved solute along the boundary layer. From these considerations it can be concluded that the present analysis of the effect of hydrodynamic dispersion on the stability of convection-enhanced dissolution in the case of non-reactive systems will result in the same critical values of the Ra number and dissolution flux as those found here for the reacting case, but with a completely different bifurcation map.