The Elder Problem with Reactive Infiltration Effects

The occurrence of buoyancy-driven flow and reactive solute transport in a fluid-saturated porous medium can be induced by either natural processes or human activities. Typical examples include the groundwater salinization in carbonate-rock aquifers and the acid treatment of oil wells in petroleum drilling industry. In this paper, the classical Elder problem of buoyancy-driven convection in two-dimensional porous media is extended to include the local chemical interactions between the solute in the pore liquid (e.g. salt such as NaCl or acids such as HCl and HCl/HF mixtures) and the solid space of the porous medium (e.g. minerals such as calcite and dolomite). Effects of the geochemical processes on the flow and mass transport are investigated. For the reactive, strongly solute-driven case in the regime dominated by the diffusive mass transport, a decrease in the net solute concentration is found as compared to the non-reactive case. This decrease is pronounced at higher values of the Damköhler number when the solute reaction rate becomes larger than the solute diffusion rate. Furthermore, the flow structure is affected by products generated by the chemical reaction when the Rayleigh number for the products is sufficiently high. In that case, numerical simulations show the formation of diluted fluid tongues exhibiting damped periodic oscillations about a nonzero mean. The results are obtained using the pseudospectral numerical method, verified against the analytical solution for the non-reactive, purely diffusive case.


Introduction
A fluid flow driven by buoyancy and gravitational forces in fluid-saturated porous media is a physical phenomenon commonly observed under natural or laboratory conditions.Buoyancy-driven flow arises due to significant fluid density changes induced by thermal or concentration gradients.When a denser fluid layer is placed onto a less dense one in the porous medium with sufficient permeability, the denser fluid can flow downwards to displace the lighter fluid, which in turn flows upwards to replace the denser fluid, and the convective circulation ensues.
Using mathematical modelling techniques, buoyancy-driven flow in porous media was studied in many publications (for reviews, see e.g.Holzbecher (1998); Phillips (2009) ).One of the first laboratory and numerical experiments describing this type of fluid flow was performed by Elder (1967a, b).In Elder (1967a), he investigated steady convection in a horizontal slab of a homogeneous porous material saturated with a viscous fluid where buoyant instabilities were created by sudden heating a portion of the lower boundary of the slab.Some attention was also paid to the role of end effects and mass discharge in the system.Extending his original model, Elder (1967b) numerically computed a time-dependent solution featuring a single plume.The other solutions to the Elder problem were later found using advanced numerical procedures (Frolkovič and De Schepper 2000;Johannsen 2003;van Reeuwijk et al. 2009).
To eliminate truncation errors associated with spatial discretization (Arakawa 1966), van Reeuwijk et al. (2009) implemented the pseudospectral method to solve a saline analogue of the Elder problem.Unlike the original thermal Elder problem, the convective circulation pattern in the saline Elder problem is initiated by a saline concentration gradient, where the salt is dissolved by freshwater at a portion of the upper boundary of the system.In van Reeuwijk et al. (2009), all the known stable steady-state solutions to the saline Elder problem were reproduced, and a bifurcation diagram for the steady states up to the saline Rayleigh number of 400 was determined.The bifurcation diagram confirmed the existence of three stable steady-state solutions reported previously using numerical continuation methods (Johannsen 2003).Moreover, van Reeuwijk et al. (2009) recommended the Elder problem at the Rayleigh number of 60, where only one steady-state solution exists, as a benchmark test case for density-dependent groundwater flow and mass transport models.For the purely diffusive case, corresponding to zero Rayleigh number, the analytical solution was derived and tested against the numerical pseudospectral solution with excellent agreement.
Buoyancy-driven flow plays a crucial role in a wide range of hydrogeological processes such as saline contamination (Bear et al. 1999;Gilman and Bear 1996;Zimmermann et al. 2006;Lin et al. 2009), carbon dioxide sequestration (Riaz et al. 2006;Ranganathan et al. 2012;Huppert and Neufeld 2014;Celia et al. 2015) and heat transfer in geothermal systems (Horne and O'Sullivan 1974;Cheng 1979;Oldenburg and Pruess 1998).Such flows through porous media are often accompanied by local chemical reactions between the liquid phase and a solid structure of the porous matrix.Adopting a system configuration similar to that of the saline Elder problem, Post and Prommer (2007) considered a reactive multicomponent transport with complex chemical reactions based on ion exchange and calcite dissolution and precipitation.They investigated the feedback mechanisms between geochemical reactions and buoyancy-driven flow in the system.The main result of their work was that the flow patterns developed during saltwater intrusion into freshwater were significantly different in comparison to the non-reactive saline Elder case for relatively small density contrasts between saltwater and freshwater.However, the importance of this feedback was shown to decrease with the increasing density contrast.
Although Post and Prommer (2007) did not take into account any porosity and permeability changes of the porous medium due to calcite dissolution or precipitation, these changes may affect the development of the flow and mass transport in a system.This was shown for calcite dissolution during seawater intrusion in a coastal carbonate-rock aquifer using the Henry problem geometry by Laabidi and Bouhlila (2016).Specifically, a reactive solute flowing through a porous medium tends to follow relatively narrow and high-permeability zones and dissolve them further, thereby increasing the permeability even more.This local reactive-infiltration instability is generally caused by an unstable dissolution front.The local stability of dissolution fronts moving through porous media has been addressed in a number of studies (e.g.Chadam et al. 1986Chadam et al. , 2001;;Sherwood 1987;Hinch and Bhatt 1990;Szymczak and Ladd 2013)).
More recently, Ritchie and Pritchard (2011) carried out a theoretical study of geochemical convection in a fluid-saturated porous medium with evolving porosity and permeability.Kinetics of the dissolution reactions was very similar to that introduced by Chadam et al. (1986Chadam et al. ( , 2001)).The convective motion was evoked by imposing a vertically varying equilibrium solubility.Linear stability analysis of the full problem was employed to describe the onset of convection, and numerical simulations were used to investigate the long-term behaviour of the system.In the strongly buoyantly-unstable regime at high Rayleigh numbers, it was shown that the reaction process stabilised the system by removing the destabilising reactive solute.Over very long-time scales, horizontally-layered porosity structures, penetrated by narrow vertical channels of enhanced porosity, were developed.
In this paper, we present a model of buoyancy-driven flow and reactive solute transport in a fluid-saturated porous medium (Sect.2).Adopting the Elder problem geometry, we consider the simplest model in which local chemical interactions between a solute and a mineral structure of the solid can occur, with the reactive solute and dissoluble minerals consumed and reaction products generated in the liquid.A second-order reaction rate for solid mineral dissolution is assumed following Lund and Fogler (1976) and Hinch and Bhatt (1990), and any dissolution effects on porosity and permeability are neglected as in Post and Prommer (2007).In Sect.3, we describe a pseudospectral numerical approach to tackle the considered problem.In Sect.4, for a weakly reactive case without buoyancydriven flow and for a reactive case with solute-driven flow, we analyse and explain effects of the geochemical processes on the flow and mass transport.Furthermore, we investigate the effects occurring when the flow is driven not only by a concentration gradient of the solute but also by a concentration gradient of the products.

Mathematical Formulation
In a two-dimensional Cartesian coordinate system x * = (x * , z * ) T , we consider a homo- geneous and isotropic porous medium, originally saturated with freshwater, in the region We consider a solute (e.g.salt) in the pore liquid with the concentration c * (x * , t * ) [mol m −3 ] , which reacts with a solid structure of the porous medium and in turn dissolves it.The porous matrix consists of dissoluble minerals (e.g.calcite) with the concentration m * (x * , t * ) [mol m −3 ] .For simplicity, there is a bimolecular reaction between a solute molecule in the liquid phase and a mineral molecule in the solid phase, which generates reaction products in the liquid phase represented by a single concentration quantity n * (x * , t * ) [mol m −3 ] .The solute and the products are perfectly mixed in the liquid phase and do not interact themselves.The liquid motion through the porous medium is governed by Darcy's law for the volume flux u * (x * , t * ) [m s −1 ] .The pressure field in the liquid phase is described by p * (x * , t * ) [kg m −1 s −2 ] .All dimensional variables and operators are denoted by the asterisk.
Boundary conditions for the concentrations of the solute and products are as follows.At the top, we impose that the solute concentration has a fixed profile: for |x * | > H * , we have c * (x * , H * , t * ) = c * 0 and for |x * | ≤ H * , we have c * (x * , H * , t * ) = c * 1 > c * 0 (see Fig. 1a).Here, c * 0 and c * 1 represent the reference concentrations of the solute in the liquid.The above boundary conditions, applied to the system with the initial solute concentration c * 0 , yield a vertical concentration gradient.In case of the concentration-dependent fluid density, it can result in a sudden emergence of buoyant instabilities.Along the whole upper boundary, we consider the concentration of the products n * 0 .At the bottom, the solute concentration c * 0 and the concentration of the products n * 0 are prescribed.While these conditions entail removing the solute and products from the bottom of the system, there is no mass loss across the vertical boundaries, which is guaranteed by zero Neumann conditions for the solute and products.In addition, all the boundaries are impermeable to the liquid, so there is no advective flux across any boundary.
Taking into account advection and diffusion of the solute, and the chemical interactions between the solute and the dissoluble minerals, the local mass conservation law for the solute is governed by where is the porosity (liquid fraction) of the porous medium, D * c [m 2 s −1 ] is the molecu- lar diffusion coefficient for the solute of the liquid phase and k * [m 3 mol −1 s −1 ] is the rate constant for the reaction.The left-hand side of Eq. ( 1) follows from the transport equation for a component moving with a fluid (Holzbecher 1998).Dispersion and its effects on the flow and mass transport are fully omitted.However, under high-density differences and high dispersivity, Jamshidzadeh et al. (2013) demonstrated that the dispersive fluid flux considered in the non-reactive saline Elder problem could noticeably influence the shape and number of developed plumes.The strongly nonlinear reaction term on the right side is incorporated to capture collision between the solute and mineral molecule, which leads to the chemical reaction related to consuming the solute and minerals.Since we assume that the surface area for the mineral dissolution is proportional to the volume of the dissoluble minerals (Hinch and Bhatt 1990), the nonlinear term is proportional to the product of the concentrations of the solute and the minerals according to Guldberg and Waage's law (Davis and Davis 2003).
From the mass conservation law for the dissoluble minerals in the solid phase, we have The conservation equation for the reaction products in the liquid phase is given by where D * n [m 2 s −1 ] is the molecular diffusion coefficient for the products in the liquid phase.Note that we do not take into account any effects of different types of intermediates that commonly appear in reaction networks (Davis and Davis 2003).
Darcy's law for porous media flow can be expressed as follows: where * [m 2 ] is the intrinsic permeability, * [kg m −1 s −1 ] is the dynamic viscosity, g * = (0, −g * ) T is the gravitational field and g * [m s −2 ] is the gravitational acceleration.The fluid density * [kg m −3 ] is taken to vary with the concentrations of solute and products, * = * (c * , n * ) , which provides a two-way coupling between the mass transport Eqs. ( 1), (3) and Darcy's Eq. ( 4).The linear approximation of the state equation for the density * is where * c and * n [m 3 mol −1 ] are the concentration expansion coefficients for the solute and the products, respectively, and * 0 ≡ * (c * 0 , n * 0 ) is the reference liquid-phase density at the reference concentrations of the solute and products.
Conservation of fluid mass requires In deriving (1)-( 5), we have applied the Oberbeck-Boussinesq approximation, which amounts to taking the liquid density as constant everywhere in the equations except for the buoyancy term.The validity of the Oberbeck-Boussinesq approximation in the nonreactive case was discussed by Johannsen (2003). ( (3)

Non-dimensionalisation and Streamfunction Formulation
In order to transform Eqs. ( 1) -( 5) to the dimensionless form, we rescale the dimensional spatial and temporal variables as the dimensional concentration fields as and the dimensional velocity and pressure as where Δc * = c * 1 − c * 0 , m * 0 is the uniform concentration of the minerals in the solid and where n * 1 is the reference concentration of the products such that n * 0 < n * 1 .We then obtain a system of dimensionless equations in the form and n = Δc * ∕Δn * are concentration ratios of the solute in the liquid to the minerals in the solid and to the products in the liquid, respectively; the Damköhler number, measuring the solute reaction rate relative to the solute diffusion rate, is defined as n is a molecular diffusion ratio of the solute to the products; p R = p − (Ra c c * 0 ∕Δc * + Ra n n * 0 ∕Δn * ) z is the reduced pressure; the concentration Rayleigh numbers for the solute and the products, respectively, are and ẑ is the upward unit vector.
To describe the flow field, for which the continuity Eq. ( 7e) holds, a streamfunction = (x, t) with u = [− ∕ z, ∕ x] T can be used in place of the velocity variable u (Holzbecher 1998).Introducing , we can rewrite the governing equations as with boundary conditions 1b).In this work, we set c 0 = 0 , c 1 = 1 and n 0 = 0 , which corresponds to zero reference concentrations of the solute and products.
The initial conditions for c, n and are taken to be zero throughout the entire domain, unless stated otherwise.The initial condition for m is identically equal to unity.Note that the zero initial concentrations of the solute and products, ensuring an ideal initial situation for mass transport in the system, can be commonly fulfilled in a laboratory experiment but seldom in a field site.

Numerical Method
Because of the symmetry of the governing Eqs. ( 8) and the boundary conditions (9), we consider only the right half of the model region, i.e. 0 ≤ x ≤ 2 and 0 ≤ z ≤ 1 , with corre- sponding boundary conditions

Pseudospectral Discretization
For spatial discretization, we use a pseudospectral method based on sine and cosine series expansions in the horizontal direction and a Chebyshev expansion in the vertical direction.Expansions for c, m, n and in the x direction are where g k = k ∕2 is a wave number of a horizontal mode k.The sine and cosine basis func- tions for c, n and are chosen such that the expansions for c, n and satisfy the horizontal boundary conditions (10c, d), (10g, h) and (10k, l), respectively.For m, we choose the cosine basis functions, so that each of Eqs.(8a-c) has a definite parity and a problem of parity mixing is avoided (Vasil et al. 2008).
Using a set of new nonlinear variables, the governing Eqs.(8a-c) can be rearranged as follows: where represent the local fluxes of the solute in the horizontal and vertical directions, respectively; represent the local fluxes of the products in the horizontal and vertical directions, respectively, and is the second-order reaction rate.It follows from ( 13), ( 14) and ( 15) that uc and un are odd in the x direction, and wc, wn and r are even.As a result, we expand these variables in the form The sine and cosine series can be constructed using discrete sine and cosine transforms.For example, the coefficients ĉk and ψk can be numerically calculated using the type-II discrete cosine transform (DCT-II) and the type-II discrete sine transform (DST-II), respectively: where are equally spaced collocation points on the interval [0, 2].
Further, we expand these coefficients in the z direction using Chebyshev polynomials (Boyd 2001).For example, Chebyshev expansions for ĉk are defined as for each k, where T l is the l-th Chebyshev polynomial on [−1, 1] and (z) = 2z − 1 .The other coefficients mk , nk , ψk , ûc k , ŵc k , ûn k , ŵn k and rk can be expanded similarly.The Che- byshev coefficients ck,l can then be calculated using DCT-II: for each k, where are Chebyshev points transformed onto the interval [0, 1].The application of the Chebyshev points, which cluster near the ends of the interval, can result in a better approximation of the concentration field c, which is affected by relatively large gradients in the vertical direction near the boundaries of the domain.
We also need to expand partial derivatives of the coefficients ĉk , nk , ψk , ŵc k and ŵn k in the z direction.For example, Chebyshev expansions for  2 ĉk ∕z 2 are given by for each k.Then the Chebyshev coefficients ck,l can be expressed as a linear combination of the coefficients ck,l : where ck = (c k,0 , … , ck,L ) T and ck = (c k,0 , … , ck,L ) T are (L + 1) × 1 vectors of coefficients and D zz is a (L + 1) × (L + 1) second-order differentiation matrix.The matrix D zz can be constructed using relations D zz = D 2 z and D z = 2D , where is the (L + 1) × (L + 1) upper triangular differentiation matrix (Johnson 1996).
Using orthogonal properties of the basis functions, we obtain a system of ordinary differential equations for the time-dependent coefficients in the form ) T , %This expression has been merged to the expression from the previous equation environment.
ũn k = (ũn k,0 , … , ũn k,L ) T , wn k = (wn k,0 , … , wn k,L ) T and rk = (r k,0 , … , rk,L ) T are (L + 1) × 1 vectors of coefficients, and I is the (L + 1) × (L + 1) identity matrix.It should be noted that the unknown coef- ficients ck,l , mk,l , ñk,l , ψk,l , ũc k,l , wc k,l , ũn k,l , wn k,l and rk,l are coupled with each other in the vertical direction whereas in the horizontal direction they can be calculated for each mode k independently.

Time-Stepping
The dynamical system (23), consisting of nonlinear equations for ck , mk and ñk , and of lin- ear equations for ̃ k , can be solved numerically by applying a one-step implicit-explicit Runge-Kutta method.To achieve appropriate convergence behaviour, explicit discretization is used for advection terms and implicit discretization for diffusion terms (Ascher et al. 1997).So, we have with initial values c(0) k , m(0) k and ñ(0) k , where h is the time step and t i = ih is the time discre- tization for i = 0, 1, … , I.
The boundary conditions in the z direction are imposed by incorporating them into the last two rows of differential operators, since the last two rows of D zz contain only zeros.Note that Eqs.(24a, c) give solutions for c(i+1) k and ñ(i+1) k at time t i+1 > 0 , which are subse- quently used for computing a solution for ̃ (i+1) k in (24d).At each time step, the coefficients ũc k,l , wc k,l , ũn k,l , wc k,l and rk,l are recalculated using the discrete sine and cosine transforms and their inverse analogues.
We perform all numerical computations implementing the above approach in Dedalus (Burns et al. 2020).To ensure pointwise convergence of the pseudospectral method, the discontinuity in the step-wise boundary condition (10b) is avoided by approximating c with a smooth function in the form where d is a transition length (van Reeuwijk et al. 2009).As d → 0 , f(x) tends to the origi- nal boundary condition.In this study, we set a value for d to 0.005, which has been shown to be sufficiently small to capture the main characteristics of the classical Elder problem.Although this value of d is too low to achieve pointwise convergence very close to the upper boundary at x = ±1 , the solution converges globally (see van Reeuwijk et al. (2009) for further details).Verification of the numerical code has been done for a non-reactive, purely diffusive case ( Da = 0 , Ra c = 0 and Ra n = 0 ) that can be solved analytically (van Reeuwijk et al. 2009).The agreement between the analytical and pseudospectral steadystate solution on the grid with 256 × 64 collocation points for h = 10 −4 has been achieved at the level of O(10 −5 ).

Results and Discussion
4.1 Weakly Reactive Case without Buoyancy-Driven Convection ( Ra c = 0 and Ra n = 0) When the Rayleigh numbers Ra c and Ra n are set to zero, Eq. (8d) reduces to Laplace's equation with a trivial solution for the streamfunction (no flow).The remaining quantities c, m and n satisfy subject to (10a-h).Note the one-way coupling of c and m to n in Eq. (25c).We will seek solution to this nonlinear reactive system in the limit of small Da and expand the concentration quantities as Note that the leading-order solution m 0 does not change with time, which can be clearly seen by setting Da = 0 in Eq. (25b).Moreover, expressing the series (26b) at t = 0 and using the fact that the initial value of m is independent of Da , we have that m 0 (x, z) ≡ 1 and m 1 (x, z, 0) ≡ 0.
Substituting the series (26) into Eqs.( 25) and comparing coefficients of equal powers of Da at O(Da 0 ) and O(Da 1 ) , we obtain a system of linear differential equations for c 0 , m 0 , n 0 and c 1 , m 1 , n 1 .We then apply the Fourier transform in the x direction and the Laplace trans- form in the time domain and solve the resulting algebraic equations to yield Taking into account the first N + 1 terms in the series (27), we numerically invert the Laplace transform using the Stehfest method (Johansson 2023).The concentration fields c 0 , c 1 and m 1 at time t = 0.5 are shown in Fig. 2. The solution c 0 corresponds to the non-reactive, purely diffusive case.The structure of c 1 indicates that the largest decrease in the solute concentration as compared to the non-reactive case occurs near the point (0, 0.6), arising as a result of solute consumption caused by chemical interactions with the solid minerals.The solution m 0 is equal to unity throughout the domain, corresponding to the non-reactive case with the prescribed initial concentration field m.
The correction m 1 is negative throughout and represents a decrease in the concentration of the minerals located mainly in the region where the concentration field c is highest.The solution n 0 is equal to zero throughout the domain, which corresponds to the non- reactive case with no products generated in the porous medium.Although we do not present the general solution of n 1 due to its complexity, we note from the form of Eqs.(25a,c) that n 1 (x, z, t) = −c 1 (x, z, t) for n = 1 and = 1 .It implies that the decrease in the solute concentration of the weakly reactive case leads to an increase in the concentration field n of a comparable magnitude.

Reactive
Case with Solute-Driven Convection ( Ra c = 400 and Ra n = 0)

Solute Plume Development
For the non-reactive saline Elder problem, the onset of convection occurs at any positive Ra c (van Reeuwijk et al. 2009).Based on the numerical bifurcation analysis with the appropriate selection of initial conditions (Johannsen 2003;van Reeuwijk et al. 2009), three distinct steady-state solutions have been identified in different ranges of Ra c .In particular, the solution featuring one plume exists for 0 < Ra c ≤ 400 , two plumes for 76 ≤ Ra c ≤ 400 and three plumes for 172 ≤ Ra c ≤ 400. Figure 3a, c and e show three qualitatively different solute concentration fields at t = 0.5 for Da = 1 , Ra c = 400 and Ra n = 0 .Values of the other parameters are m = 1 , n = 1 and = 1 .Numerical simulations were carried out using the pseudospectral method with 256 × 64 nodes for h = 10 −4 .The results differ by the number of vertical solute plumes that have developed as a result of convective solute transport.To obtain them, we used three different initial conditions for the solute concentration, as described in the bifurcation analyses of Johannsen (2003) and van Reeuwijk et al. (2009).As an initial condition for the solution featuring one plume (shown in Fig. 3a and denoted as R 1 ), we took the solution for Da = 0 , Ra c = 60 and Ra n = 0 featuring a single plume.The solution featuring two plumes (shown in Fig. 3c, denoted as R 2 ) was calculated using the zero initial condition, and the solution with three plumes (shown in Fig. 3e, denoted as R 3 ) using c(x, z, 0) = [cos( z) + 1]∕2 over the entire domain.
Figure 3b, d and f show a concentration difference c − c NR , where c NR corresponds to the non-reactive case, for all solutions R 1 , R 2 and R 3 , respectively.In all cases, most of the domains are coloured by blue shades, indicating less solute in these regions at t = 0.5 relative to the non-reactive case.The average decrease for R 1 is 1.9% , for R 2 is 1.6% and for R 3 is 1.3% .For solutions R 2 and R 3 , there is a slight symmetrical displacement of the individual plumes compared to the non-reactive case: for R 2 , the roots of the plumes are closer to each other with the maximum absolute difference in the concentration of 7.4% (Fig. 3d); for R 3 , we can observe similar behaviour for the roots of the two outer plumes with the maximum absolute difference of 7.9% (Fig. 3f).In addition, the lower parts of the outer plumes near x = ±0.5 are extended slightly towards the inner plume so that the inner plume is correspondingly narrower.As a result, there are regions with higher concentration (red shades) compared to the non-reactive case.The occurrence of these regions is due to advective mass transport present in the convective system.Note that no regions of increased c appear in Fig. 2b, because the weakly reactive case is purely diffusive and features no fluid flow.
One of the characteristics describing the system behaviour is the vertical solute flux at the top of the domain.Formally, this quantity is defined by and is referred to as the local Sherwood number.The average value of Sh c can be expressed as and is referred to as the average Sherwood number.Since the upper boundary is considered impermeable to the liquid, Sh c measures the average diffusive solute flux across the upper boundary.In Prasad and Simmons (2003), in contrast, the Nusselt rather than Sherwood number is used to quantify the total (advective plus diffusive) solute transport rate relative to the diffusive solute transport rate since the imposed pressure boundary conditions at the top allow for solute advection across the boundary.
Figure 4 compares temporal evolution of Sh c in the time range from 0 to 0.2 for different Da .The average vertical diffusive solute flux at the top for the non-reactive simulation coincides with the results by Prasad and Simmons (2003).They performed calculations up to dimensional time t * = 20 years, which approximately corresponds to dimensionless time t = 0.1 for typical material transport properties.From Fig. 4, we can see that values of Sh c begin to differ from the non-reactive case with increasing Da .The curves for Da = 0 and 0.1 virtually overlap, while those for Da = 1 and 10 vary significantly from the others.The average flux of the solute for Da = 1 decreases more slowly from point D compared to the non-reactive case.The values of Sh c for Da = 10 are larger in the phase A-B and C-E relative to those for the other values of Da , and this difference becomes more pro- nounced behind point D. The explanation for this is the decrease in the solute concentration and increase in the solute transport rate from the source boundary into the domain interior as the chemical reaction rate increases.It occurs only in the diffusive regime D-E because advection is more effective in slowing down the reaction process than diffusion.
In the long term, we might expect the concentration of the minerals to be zero throughout the entire domain since the minerals are gradually being dissolved.At the same time, the system tends towards the non-reactive state (the consumption of the reactive solute is being reduced), thereby effectively reaching a steady state.Therefore, all values of Sh c shown in Fig. 4 are expected to tend to the value of 4 in the long term.This value approximately corresponds to the average diffusive solute flux at the top for the non-reactive steady state

Redistribution of Reaction Products
The choice of the initial condition for c affects not only the development of the solute concentration field but also the concentration field of the products.Since these products result from the chemical reaction between the solute and the minerals, we will assume that at t = 0 they are absent in the liquid phase.
Besides the solutions R 1 , R 2 and R 3 discussed above, we also obtain the corresponding concentration fields for the products, which are shown together with the streamlines of the flow in Fig. 5.All three solutions for the products are affected by the number of the vertical solute plumes in the system.The concentrations of the products are rather low because the mass is absorbed at the upper and lower boundaries.For R 1 , the maximum concentration of the products is 3.9% , for R 2 is 3.3% and for R 3 is just 2.9% .The remaining portion of the products in the system is largely concentrated in the two outer convection eddies.
For a better understanding of the mass discharge in the system, in Fig. 6, we show the local vertical diffusive fluxes of the solute, Sh c , and the products, Sh n , at the top and bottom of the domain at t = 0.5 .The quantity Sh n is defined analogously to Sh c .At the top, the products are absorbed most in the regions where Sh n approaches the smallest (negative) values (Fig. 6b, d and f).In all cases, this occurs at the edges of the solute source near x = ±1 , and for R 2 and R 3 also in the regions between the roots of two adjacent sol- ute plumes.In these regions, there is a higher supply of the solute from the source to the system, especially near x = ±1 (Fig. 6a, c and e), where it is invoked by the discontinu- ous transition of c in the boundary condition.However, both local fluxes at the top, where the plume roots are located, are relatively small in magnitude.This local phenomenon can be explained by the presence of the convective circulation accompanied by vertical solute transport from the source boundary into the domain interior, as shown by the streamlines in Fig. 5. Thus, at the top along the source, the local diffusive flux of the products Sh n is inverse to the local diffusive flux of the solute Sh c .
At the bottom, the solute and products are discharged in the regions below the individual solute plumes, where the downward flow from the top is strongest.Along the lower boundary towards the vertical edges of the domain and in the regions between two adjacent solute plumes, values of Sh c and Sh n tend to zero.The reason for this is the upward flow of fluid with a lower solute concentration from the bottom.To sum up, at the bottom, the local diffusive flux of the products Sh n is proportional to the local diffusive flux of the solute Sh c .The same behaviour also holds for Sh n at the upper boundary outside the solute source.
Following Prasad and Simmons (2003) and Post and Prommer (2007), we define the additional indicator variables The quantity SP c measures the average amount of the reactive solute present in the liquid at a particular time, SP m represents the dissoluble minerals in the solid and SP n represents the reaction products in the liquid.Figure 7a illustrates how SP c varies with time for different Da .In all cases, there is a sharp break in the slope of the individual curves at t ≈ 0.05 indicating a change in the growth rate of the solute saturation.For Da = 0 , when SP c reaches about 25% , the sol- ute saturation slows down until it stabilizes at 36.3% .For Da = 10 , we can see the largest decrease in the saturation rate.This decrease occurs since the diffusive mass transport rate begins to increase (Fig. 4).Values of SP c for Da = 1 and 10 reach about 35% at t = 1 and tend to the value of the non-reactive case with increasing time.For the reactive simulations, SP m decreases with time until all minerals are consumed, and the rate of the decrease is influenced by Da (see Fig. 7b).By t = 1 , this decrease in the average mineral content for Da = 0.1 is relatively small-just about 3% , for Da = 1 is about 27% , and for Da = 10 reaches about 81% .Although the products are not initially present in the liquid, they quickly generate as time proceeds.The average amount of the products at a particular time strongly depends on Da (see Fig. 7c).For Da = 1 , SP n increases with time until it reaches the maximum value of 1.7% at t = 0.33 , then it begins to decline slightly.The quantity SP n for Da = 10 behaves in a similar way, but it grows faster to reach its maximum of 8.7% at t = 0.18 .Subsequently, the products are being removed from the system, and at t = 1 the value of SP n for Da = 10 is nearly identical to that for Da = 1.

Impact of Reaction
Products on Convection ( Ra c = 400 and Ra n ≥ 0) So far, we have considered that the concentration of the products does not affect fluid density in the system ( Ra n = 0 ).However, in general, the fluid density can be affected by the reaction products.An example of this is the chemical reaction between solid calcium carbonate ( CaCO 3 ) and aqueous solution of sodium chloride (NaCl) generating sodium carbonate ( Na 2 CO 3 ) and calcium chloride ( CaCl 2 ) in the liquid phase.This is a typical reaction network that occurs under some natural conditions (Harvey 2000).Such reaction network is nowadays used for industrial production of the sodium carbonate by Solvay process.
In this section, we will analyse buoyancy-driven flow in the porous medium using the streamfunction calculated for reactive simulations ( Da = 1 and Ra c = 400 ) in dependence on Ra n .The temporal development of a maximum value of can be used to characterise the flow regime (Holzbecher 1998).The streamfunction reaches extreme values in the centre of the convection eddies and zero values at the edges of the domain.Therefore, an increase or decrease in the maximum values of at a particular time can be associated with the emergence or disappearance of an eddy in the solution.
Figure 8 illustrates how the maximum values of change with time for different Ra n in the range from 0 to 4000.The curves for Ra n = 0 , 40 and 400 are almost identicaljust the curve for Ra n = 400 is slightly below the others from t ≈ 0.1 .For Ra n = 4000 , the maximum values of are initially greater than those for the other values of Ra n until t ≈ 0.1 , when they become moderately smaller.At t ≈ 0.2 , there is a sharp increase leading to damped periodic oscillations about a nonzero mean.The oscillation period is approximately 0.014 and the frequency is 71 (in dimensionless units).This behaviour is essentially a consequence of a secondary instability in the flow regime due to the gradual decrease in the concentration of the products, as shown in Fig. 7c.The oscillatory convection makes its first appearance at a value of Ra n = 2370.
Figure 9a, c and e show the solute concentration fields for Da = 1 , Ra c = 400 and Ra n = 4000 at t = 0.5 , where the three stable steady states, discussed above, were used as initial conditions for c, respectively.Compared to Fig. 3a, c and e, we observe tongue-like structures in the lower parts of the domains.The associated change in distribution of the products is noticeable in Fig. 9b, d and f.Furthermore, in Fig. 9b and d, we can clearly see how the localized streamlines are deformed by the periodically generated solute tongues.The formation of similar fluid tongues was previously reported in Hele-Shaw cell experiments by Horne and O'Sullivan (1974), who investigated oscillatory convection in a nonreactive porous medium heated from below.

Conclusion
In this paper, we have formulated the two-dimensional porous-medium model including the local chemical interactions between the solute in the pore liquid and the solid space of the porous medium as an extension of the classical Elder problem.For the model in the non-dimensional and streamfunction formulation, we have derived the asymptotic solution to the weakly reactive, purely diffusive case by applying the Fourier and Laplace transforms.Solutions to the reactive case with fully developed convection have been calculated numerically by applying the pseudospectral method for spatial discretization and the implicit-explicit Runge-Kutta method for time-stepping.
Our attention has been further focused on analysing mass distribution in the reactive porous medium in dependence on the Damköhler number Da .For the strongly reactive case with Da = 1 , we have made direct comparisons to the corresponding non-reactive problem (van Reeuwijk et al. 2009) in terms of the solute concentration fields with the different number of vertical plumes.In the regime dominated by the diffusive mass transport, we have found a decrease in the net solute concentration as compared to the non-reactive case.This decrease is more pronounced at higher values of Da .Interestingly, the associated reaction products are largely concentrated in the convection eddies at the side boundaries with the maximum concentration of about 4%.
We have also investigated the impact of the chemical reaction on the flow regimes quantified by the concentration Rayleigh numbers Ra c for the solute and Ra n for the prod- ucts.The flow structure in the reactive convective case with Da = 1 and Ra c = 400 is affected by the reaction products when Ra n is sufficiently high.Numerical simulations with Ra n = 4000 showed the formation of diluted fluid tongues that exhibited damped periodic oscillations about a nonzero mean.material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Fig. 1
Fig. 1 Schematic diagram of a the dimensional and b non-dimensional geometry of the model for buoyancy-driven flow and mass transport in the considered porous medium in the 2D vertical cross-section and where g k = k ∕2 and 2 k = g 2 k + s .The function c 0 represents the Laplace transform of the function c 0 , i.e. c 0 (x, z, s) = ∫ ∞ 0 c 0 (x, z, t)e −st dt .The other functions c 1 and m 1 are defined similarly.

Fig. 2
Fig. 2 Concentration fields c 0 , c 1 and m 1 at time t = 0.5 .Parameter values are m = 1 and N = 500 .Concentration contours are depicted by black solid lines

Fig. 3
Fig. 3 a, c, e Concentration fields c at time t = 0.5 for parameters m = 1 , n = 1 , = 1 , Da = 1 , Ra c = 400 and Ra n = 0 , calculated with three different initial solute concentrations.b, d, f Comparison between the reactive and non-reactive solute concentration fields, respectively.The quantity c NR represents the corresponding solute concentration for the non-reactive simulation with Da = 0 .Concentration contours of c are depicted by black solid lines in all figures

Fig. 4
Fig. 4 Temporal evolution of the average vertical diffusive solute flux Sh c at the top of the domain for the non-reactive simulation with Da = 0 (blue) and reactive simulations with Da = 0.1 (black), 1 (red) and 10 (green).Values of the other parameters are m = 1 , n = 1 , = 1 , Ra c = 400 and Ra n = 0 .Phases of the flow development and solute transport in the considered medium are denoted by letters from A to E. Between points A and B, a boundary layer is formed through diffusion near the solute source at the top.From point B, the onset of instability is exhibited by the formation of vertical solute plumes, which touch the bottom first at point C.During the phase between C and D, the denser part of the generated plumes reaches the bottom.Finally, between D and E, mass transport in the system is dominated by diffusion

Fig. 5
Fig. 5 a, b, c Concentration fields n at time t = 0.5 for the solutions R 1 , R 2 and R 3 , respectively.Parameter values are m = 1 , n = 1 , = 1 , Da = 1 , Ra c = 400 and Ra n = 0 .Streamlines of the velocity field u are depicted by blue solid lines

Fig. 6
Fig. 6 Vertical diffusive fluxes of the solute, Sh c (a, c, e), and the products, Sh n (b, d, f), at the top (black) and bottom (blue) of the domain at time t = 0.5 for the solutions R 1 , R 2 and R 3 , respectively.Parameter values are m = 1 , n = 1 , = 1 , Da = 1 , Ra c = 400 and Ra n = 0

Fig. 7
Fig. 7 Temporal evolution of the average amount of the solute, SP c (a), the minerals, SP m (b), and the products, SP n (c), for the nonreactive simulation with Da = 0 (blue solid) and reactive simulations with Da = 0.1 (orange dashed), 1 (green dash-dotted) and 10 (red dotted).Values of the other parameters are m = 1 , n = 1 , = 1 , Ra c = 400 and Ra n = 0

Fig. 9
Fig. 9 a, c, e Solute concentration fields c at time t = 0.5 for parameters Da = 1 , Ra c = 400 and Ra n = 4000 , calculated with three different initial solute concentrations.Values of the other parameters are m = 1 , n = 1 and = 1 .Concentration contours of c are depicted by black solid lines.Product concentration fields n (b, d, f) correspond to the solute concentration fields (a, c, e), respectively.Streamlines of the velocity field u are depicted by blue solid lines