Nonlinear Formulation of Multicomponent Reactive Transport With Species‐Specific Dispersion Properties

The modeling of reactive transport through porous media is a challenging numerical problem. Methods of solution have leveraged the stoichiometry of chemical reactions to address the transport of multiple aqueous species by expressing them in terms of an equivalent, linearly independent variable (component). This approach effectively decouples advection‐dispersion transport from the source terms associated with equilibrium reactions. A common assumption found in the literature is that all species disperse with the same transport coefficients. Recent experimental studies have discussed that this is not necessarily the case, particularly for transverse mixing, which is limited by the species‐specific molecular diffusion. This article presents a formulation of multicomponent reactive transport that takes into account the differences in dispersion coefficients. These differences lead to a nonlinear transport equation for the components, from where an expression for evaluating reaction rates is derived. It is demonstrated that this expression simplifies to the well‐known equations assuming the same dispersion for all species. Numerical simulations of a binary chemical system under diffusion‐ and advection‐dominated transport conditions are used to evaluate the influence that differential transport coefficients have upon the output of chemical reactions. Results indicate that differences in transport coefficients are particularly relevant when the chemical signature of the input solutions is not strongly dominated by one of the species in the component. Unexpectedly, this opens the possibility to mineral dissolution coexisting with precipitation during the mixing of two waters in equilibrium. This phenomenon can be explained by nonlinear mixing processes proportional to the differences in transport coefficients.


Introduction
Modeling reactive transport in heterogeneous porous media, with accurate predictions of mixing and chemical reactions, is a complex numerical task (Dentz et al., 2011;Kräutle & Knabner, 2005;Steefel & MacQuarrie, 1996;Valocchi et al., 2019).Reactive transport models involve coupling the chemical interaction between several species and compounds, which may eventually distribute in multiple phases, with transport equations for each of them.As a result, computational costs escalate rapidly with the complexity of chemical systems and the scale of resolution of the medium, as well as chemical heterogeneities.A common approach used to overcome this computational load is to decrease the total number of equations that need to be solved by creating a reduced-order basis of conservative components, typically defined as a linear combination of the chemical species based on the stoichiometry of reactions (e.g., De Simoni et al., 2007;Lichtner, 1985;Molins et al., 2004;Saaltink et al., 1998;Sanchez-Vila et al., 2010;Steefel & MacQuarrie, 1996).In practical applications, common simplifying assumptions have considered that reactions occur sufficiently fast with respect to the time scales of transport (local equilibrium), and that all chemical substances disperse with the same transport coefficients.Such conditions give several advantages for reactive transport simulations and their interpretation.Namely, (a) the linear combination naturally leads to a smaller number of transport equations expressed in terms of an equivalent variable (also known as component); (b) it removes the explicit computation of reaction-related source terms for aqueous species, decoupling the transport equations from chemistry; and (c) it allows obtaining a spatio-temporal description of reaction rates from the distribution of the equivalent variable and chemical speciation.Methods following these assumptions have been employed, for example, in the interpretation of transport processes involving mineral precipitation/dissolution (e.g., Guadagnini et al., 2009;Katz et al., 2011;Rezaei et al., 2005), and in understanding the influence of aquifer heterogeneities on mixing-controlled reactions (e.g., Fernàndez-Garcia et al., 2008;Luo et al., 2008).Similar formulations have decoupled reactive transport in terms of the mixing fraction of waters with a known chemical signature (e.g., Cirpka & Valocchi, 2007;De Simoni et al., 2007;Ginn et al., 2017).Furthermore, studies have also discussed techniques for the simulation of kinetic reactions (e.g., Molins et al., 2004;Saaltink & Rodríguez-Escales, 2022;Sanchez-Vila et al., 2007, 2010) and the application to reactive transport in multicontinuum media (e.g., Donado et al., 2009;Willmann et al., 2010).
In most applications of reduced-order reactive transport models based on conservative components, the assumption that species disperse with the same transport coefficients has consistently remained.Most likely, this interpretation stems from the common practice of considering mechanical dispersion as a linear function of the velocity, with dispersivity (a length scale characteristic of the porous medium) as the proportionality constant (e.g., Bear & Cheng, 2010;Freeze & Cherry, 1979).Furthermore, early models of mechanical dispersion have considered a species-specific dependence that vanishes for sufficiently large Péclet numbers (e.g., Bear & Bachmat, 1967, 1990).In this context, classical models of mechanical dispersion become independent of speciesspecific properties (e.g., diffusion).However, this interpretation is in conflict with recent experimental evidence, and numerical studies (both at the pore and continuum scales) contending that, even in the limit of advectiondominated transport, classical dispersion models cannot accurately describe the transverse mixing of species with different molecular diffusion coefficients (e.g., Chiogna et al., 2010;Hochstetler et al., 2013;Olsson & Grathwohl, 2007;Rolle et al., 2012).Diffusion becomes a limiting factor of the extent of mixing, most notably, in the transverse direction.An accurate characterization of transverse dispersion is crucial for reactive transport.Several studies have discussed that transverse mixing is the process that ultimately controls the extent of reactions and dilution in porous media (e.g., Chiogna et al., 2011;Cirpka et al., 1999;Fernàndez-Garcia et al., 2008;Klenk & Grathwohl, 2002;Rolle, Chiogna, et al., 2013).In addition, it is remarked that diffusion-limited transport is likely to occur in many aquifers due to the large spatial variability of their hydraulic properties.Most aquifers include lenses of clay or other low-permeability materials where solutes are primarily transported through diffusion (e.g., Carrera et al., 1998;Haggerty & Gorelick, 1995), which is inherently a property specific to each chemical species.This is why an accurate representation of low-permeability inclusions in aquifers is considered essential to predict the tailing of concentrations and to assess the time needed for remediation (e.g., de Barros et al., 2013;Haggerty & Gorelick, 1995).Diffusion-limited transport is also clearly the case in fractured systems represented by double porosity media (e.g., Martinez-Landa et al., 2012) and in the unsaturated zone (e.g., Arora et al., 2019;Haberer et al., 2015;Molins et al., 2010;Neale et al., 2000) with transport coefficients influenced by the heterogeneous distribution of residual water contents.This underscores the necessity for a formulation of reactive transport based on components that is compatible with scenarios in which the transport coefficients of various species cannot be assumed to be identical.Ideally, one could achieve a more comprehensive formulation that reconciles discrepancies in interpreting dispersion properties.This formulation should also solve the problem using a reduced set of transport equations, while providing the means to assess reaction-related source terms by considering different transport coefficients.In this context, some studies have already discussed this issue to a certain extend.For example, Liu et al. (2011) evaluated multiple alternatives for the transport of components with different diffusion coefficients, also considering effects related to the charge interaction between species, which becomes necessary for charge balance under certain conditions (Muniruzzaman et al., 2014;Muniruzzaman & Rolle, 2017).More recently, Sole-Mari et al. (2021) discussed a formulation of reactive transport based on components and Random Walk Particle Tracking, with species-specific dispersion coefficients.Nevertheless, none of the aforementioned studies provided expressions to evaluate the impact of neglecting differences in transport coefficients on reaction-related source terms, similar to the work of De Simoni et al. (2005).These authors presented a framework to evaluate source terms from the transport of aqueous components, allowing for example, to determine the amount of precipitation (dissolution) of a mineral.However, the analysis was performed under the assumption that species disperse with the same transport coefficients.This motivates a fundamental analysis of the reactive transport formulation based on components that aims to discuss its applicability to scenarios in which chemical species disperse with different transport coefficients.The analysis should provide an explicit expression for evaluating reaction-related source terms while accounting for these differences.
This article provides a new nonlinear formulation of the reaction rates based on a reduced-order basis of conservative components for species-specific diffusion and dispersion coefficients.At this stage, the formulation does not contemplate the electrochemical coupling of species.So, it shall be considered to be representative of those scenarios were chemical species disperse with properties closer to their liberated state (self-dispersion, refer to Liu et al., 2011;Muniruzzaman et al., 2014;Muniruzzaman & Rolle, 2017;Rolle, Muniruzzaman, et al., 2013).Similarly, the formulation is discussed considering only local equilibrium reactions.This implies that deviations from the equilibrium state induced by transport are instantaneously corrected by reactions (Engdahl et al., 2013).The nonlinear formulation preserves the basic idea of an equivalent problem expressed as a linear combination of species, which is solely determined by the stoichiometry of the chemical system, independently of dispersion properties.However, a new interpretation of the formulation becomes necessary for addressing the case of species-specific transport properties.In these conditions, the transport coefficients of the equivalent variable become nonlinear functions of the variable itself (as in Fujita, 1952;Philip, 1960).This transformation arises from the evaluation of terms related to chemical speciation while still excluding the explicit consideration of source terms from the transport equations.The general expression obtained is shown to reduce to the well-known equations for equal transport coefficients for all chemical species.The analysis reveals that the influence of nonlinear processes is directly proportional to the differences in transport coefficients, enabling the evaluation of the impact that differential dispersion coefficients have upon the estimation of reactions.
The structure of this article is as follows.Section 2 provides the nonlinear formulation of multicomponent reactive transport with species-specific dispersion and a generalized expression for evaluating source terms.In Section 3, the nonlinear formulation is applied to a binary chemical system.Numerical simulations performed in both diffusion-and advection-dominated conditions are shown to evaluate the influence of differential transport coefficients on reactions.Conclusions of the analysis are summarized in Section 4.

Transport Equations
In multicomponent reactive transport, it is common practice to express the set of transport equations for all chemical species in compact form, employing a matrix-vector notation.The system of transport equations for the concentration of all N c species in saturated porous media is expressed as (De Simoni et al., 2007;Saaltink et al., 1998;Sanchez-Vila et al., 2010;Sole-Mari et al., 2021) where m = (m 1 ,…,m N c ) T is the vector grouping the mass per unit volume of porous medium for all species, T is a vector of scalar linear transport operators for each species, and T is a vector comprising source/sink terms for each species, in this case, due to equilibrium reactions.The vector of concentrations can be organized in subvectors that group the concentrations of constant activity species (e.g., minerals, pure phases), and the concentrations of species dissolved in the aqueous phase, such that c = (c κ ,c α ) T , where the subscript κ stand for constant activity species, and α for the species dissolved in the aqueous phase.For the latter, the mass per unit volume is defined as m α = ϕc α , being ϕ the medium porosity.For constant activity species, the mass per unit volume is equivalent to the concentration, that is, m κ = c κ .In the context of advection-dispersion, the linear transport operator for a mobile species is defined as (2) In the previous expression, q is the specific discharge, which is taken to be the same for all mobile species, and D i is a hydrodynamic dispersion tensor specific for the species i.In the context of this work, constant activity species indicated with subscript κ are considered to be only immobile minerals, hence, the associated entries in the linear transport operator are L κ = 0.

Classical Definition of Components in Equilibrium
The system of reactive transport equations, expressed in terms of the concentration of chemical species, can be transformed into an equivalent transport problem by taking advantage of the additional equations provided by the set of N r equilibrium reactions.These are expressed by the mass action law, summarized in the expression where S is the matrix of stoichiometric coefficients with dimensions N r × N c .This matrix contains positive coefficients for the products of reactions and negative for the reactants (Lichtner, 1985;Saaltink et al., 1998); loga = (loga 1 ,…,loga N c ) T is a vector containing the logarithm of the activity of species and logK = (logK 1 ,…,logK N r ) T contains the logarithm of equilibrium constants for all reactions.The activity of a species c i is defined as a i = γ i (c)c i , where γ i is an activity coefficient.Constant activity species (pure phases) are by convention considered to be of unit activity.Following De Simoni et al. (2007), the vector of concentrations is further subdivided as c = (c κ ,c p ,c s ) T where the subscripts p, s stand for primary and secondary aqueous species, respectively.Thus, the stoichiometry matrix can be written as where S κ are the stoichiometric coefficients for constant activity species, S p for the aqueous primary species, and I N r is the identity matrix with dimensions of N r .The method based on components postulates that instead of solving N c transport equations, including the evaluation of source terms due to chemical reactions, one can express the reactive transport problem in terms of N u = N c − N r − N κ linearly independent aqueous components u (Saaltink et al., 1998;Steefel & MacQuarrie, 1996), with N κ the number of constant activity species.Components are defined as linear combinations of the vector of concentrations u = U α c, where U α is known as the component matrix for the aqueous species.This matrix is of dimensions N u × N c , and is part of the component matrix U employed to linearly combine the transport equations of all chemical species.The matrix is not unique and several methodologies for its construction have been proposed in the literature in order to facilitate a straightforward interpretation of chemical processes (e.g., De Simoni et al., 2005, 2007;Molins et al., 2004;Saaltink et al., 1998).
The main idea behind the method is that the product between a component matrix U and the vector of sources/ sinks F, retains the source terms only for the equations of constant activity species c κ and removes them for the linear combination of the transport equations of aqueous species c α (transport equations for the aqueous components u).This feature allows to solve the advection-dispersion of components without explicitly calculating reaction-related source/sink terms, which is particularly convenient for problems dealing with immobile minerals.For equilibrium reactions, source/sink terms are expressed as the product between the transposed stoichiometry matrix and the vector of reaction rates F = S T r.This allows to define a convenient form of the component matrix (De Simoni et al., 2007) where U κ is the matrix for constant activity species.Notice that there are as many components N u as aqueous primary species c p , and as many aqueous secondary species c s as equilibrium reactions N r .It can be easily verified that the product between U (Equation 5) and S T (Equation 4) removes source/sink terms from the transport equations of components, whereas these are explicitly retained in the equations for constant activity species.It is noteworthy to remark that common practice considers a component matrix for a given chemical system to be independent of space and time.The latter is satisfied as long as the set of chemical reactions remains consistently active within the domain of interest.However, this may not always be the case and studies have discussed alternatives for dealing with situations where the species assemblage changes, for example, due to the formation or depletion of minerals (e.g., Ginn et al., 2017;Saaltink et al., 1998).The latter is not considered at this stage, preserving in subsequent analyses the assumption that the component matrix remains constant.Once this matrix is defined, the classical formulation of the problem continues with the linear combination of transport equations, typically under the assumption that dispersion is the same for all aqueous species.This means that all the transport coefficients of the linear transport operators become the same L(c).Therefore, the multiplication of the transport equations (Equation 1) by the matrix U (Equation 5) leads to for j = 1, …, N u .In this last expression, u j is the jth-component defined as u j = ∑ i U α ji c i , with U α ji being an element from the aqueous component matrix U α .Once the transport of components is solved, there is a final Water Resources Research 10.1029/2023WR036358 PÉREZ-ILLANES ET AL. necessary step in order to obtain the magnitude of concentrations, considering also the mass action law (Equation 3).This last step can be done, for example, with the Newton-Raphson method and is of relatively low cost in comparison to solving the complete set of transport equations (Molins et al., 2004).

Formulation for Species-Specific Dispersion
As seen previously, the definition of components as a linear combination of chemical species, is independent of dispersion properties, and is defined only by the stoichiometry of the chemical system.This implies that the linear combination of equations can also be applied to the case with species-specific dispersion, leading to the following transport equation for the components The difference between Equations 7 and 8 is the dispersion term: applying the component matrix leads to a linear combination of dispersive fluxes (similar to Lichtner, 1985) rather than concentrations.This is a direct consequence of considering different dispersion coefficients, which cannot be defined as a common factor that ultimately enables the linear combination of concentrations.However, in a strict sense, expressing the set of transport equations in terms of the equivalent variable u is a change of variable, meaning that the concentrations of species become a function of the new variable u, that is, c = c(u) (De Simoni et al., 2005;Sanchez-Vila et al., 2010).Therefore, 2005), the concentrations are also considered to be affected by the variability of equilibrium parameters K, whose magnitude may depend on temperature, for example.For the purposes of this work, K is considered to be homogeneous and constant such that the expression for the concentration gradients is well defined.Replacing into Equation 8 leads to The term pre-multiplying the gradient of components can be understood as an equivalent dispersion coefficient, defined as which determines the dispersive transport of the component j due to the spatial gradients of the component k.This term is a nonlinear function of u because of its explicit dependence on the derivatives from the change of variable (Jacobian).This dependence is a direct consequence of hydrodynamic dispersion coefficients being different for each species.To prove this statement, consider the definition of components ji ∇c i , assuming that the component matrix is spatially uniform.Application of the chain rule gives and is zero otherwise.Hence, in the case that dispersion is the same for all species, that is, D i = D, the equivalent dispersion (Equation 10) adopts this same value, and the cross-component dependence is lost, meaning that the dispersive transport of one component is not any more affected by the gradients of other components.Defining J u (c) as the Jacobian matrix (containing the derivatives of concentrations with respect to components), the previous result is more succinctly expressed as U α J u (c) = I N u , which is an identity of the change of variables.

An Expression for Source/Sink Terms
In this section, an expression for source/sink terms is derived from the linear combination of equations and the chain rule.For this purpose, the transport equations with species-specific dispersion are simplified by considering a fully saturated, undeformable porous medium without external flow sources.Furthermore, the influence of temporal changes in porosity due to reactions of precipitation/dissolution is considered to be of low relative magnitude in comparison to the other transport processes.This allows to rewrite the transport equation for a species as

Water Resources Research
Application of the chain rule for the derivatives of concentrations leads to The simplifying assumptions for the properties of the porous medium are also made for Equation 9, which is then replaced into the squared brackets of the previous expression (Equation 12), leading to an expression for the source term Appendix A shows that this expression reduces to the form presented in De Simoni et al. ( 2005) if the same dispersion coefficients are considered for all species.Notice in Equation 13that dispersion coefficients of all species play a role in the calculation of the source/sink terms.Therefore, this expression can be used to study the effect of simplifying assumptions for dispersion on reaction rates.From Equation 13, by only considering the transport equations referring to secondary species and bearing in mind the considered definition of the stoichiometry matrix (Equation 4), it is straightforward to obtain an expression for the reaction rate r i , related to a secondary species i where D s i is the dispersion tensor of the species c s i .By noticing that the expression for reaction rates (Equation 14) can also be written as The previous equation implicitly considers that the equivalent dispersion is a spatially symmetric tensor, which is true because this term is the weighted sum of the symmetric dispersion tensors of the chemical species.This expression shows more clearly two major contributions to the rates of reaction.The first term encompasses the nonlinear mixing of all components driven by the equivalent component dispersion.The second term illustrates local nonlinear effects proportional to the difference between the flux driven by the dispersion of the secondary species and the fluxes explained by the equivalent component dispersion.Notice that in case that all species share the same dispersion properties, Equation 16readily reduces to the well-known formula for reaction rates (De Simoni et al., 2005).Considering that the spatial gradients of the components are zero at the boundaries of a domain Ω, the domain-integrated form of Equation 16shows that changes in the global product of reaction are explained by the integral nonlinear mixing of the components, driven by the equivalent dispersion.That is, The previous case can be considered to be representative, for example, of scenarios where mixing-controlled reactions are occurring sufficiently far away from the domain boundaries, or also for a fully impermeable domain.

Results and Discussion
This section presents results from numerical simulations considering different scenarios of reactive transport with species-specific dispersion.First, the nonlinear formulation is applied to a binary chemical system, which is consistently employed for subsequent simulations.Nonlinear dispersion is simulated by means of a finitedifferences code.Details about the solver and validation test cases are shown in Appendix B.

System Definition
Simulations discussed throughout this section consider a simple binary chemical system in local equilibrium, which undergoes a reversible mineral precipitation/dissolution reaction, as follows where p and s are the aqueous primary and secondary species, and κ is a mineral.The reaction rate is expressed as r.Aqueous species are considered to posses different transport coefficients, indicated accordingly in subsequent sections.Further simplifications consider that the product of reaction is ubiquitous and activity coefficients for the aqueous species are considered unitary.In this case, the mass action law (Equation 3) simplifies to c p c s = K, where K is the solubility equilibrium constant.For this setup, the stoichiometry matrix is given by S = (1, 1, 1), which leads to an aqueous components matrix written as U α = (0,1, 1).This means that the chemical system can be solved in terms of a single conservative component u = c p c s .The main advantage of the binary system is that one can derive explicit expressions for the Jacobian of concentrations with respect to components.These are obtained by calculating the derivative of the solution to the system of equations c p c s = K and u = c p c s for constant K. Concentrations obtained for this system are expressed in terms of a non- From these expressions, it is possible to analyze the ratio cp / cs (Figure 1), which provides some additional meaning to the different values of the component ũ, in this case, related to the purity of a solution.For example, for values of ũ ≈ 5, the ratio is cp / cs ∼ 100.Hence, although the solution satisfies the equilibrium reaction (Equation 18), it is almost exclusively composed of the primary species cp .Similarly, for values of ũ ≈ 0.5, the ratio is cp / cs ∼ 3. So, although the primary species is dominant, close to one-quarter of the solution is composed of the secondary species cs .From the definition of the component u, derivatives of primary and secondary species are related by ∂c p /∂u = 1 + ∂c s /∂u.It is straightforward to verify that this relation meets the identity product U α J u (c) = 1, as discussed previously.Derivatives of the expression for the secondary species (Equation 20) with respect to the component are given by Water Resources Research Notice that Equation 21 satisfies 1 < ∂ cs /∂ ũ < 0, which is verified by taking the limit of ũ toward negative and positive infinity.The second derivative (Equation 22) is always positive in this case.These remarks should be kept in mind for the interpretation of subsequent expressions for the reaction rates.

Equivalent Dispersion for the Component
Recall the definition of the equivalent dispersion for components (Equation 10).The binary chemical system leads to a single component, hence, there is only one equivalent dispersion given by where D p , D s are the dispersion tensors for the primary and secondary species, respectively.By introducing the relationship between derivatives of concentrations with respect to the components, it becomes possible to rewrite the equivalent dispersion as a function of the derivative of the secondary species and after replacing the expression for derivatives (Equation 21), it is obtained Notice that when ũ = 0, that is cp = cs , dispersion for the component is the average between the dispersion of the primary and secondary species (refer to Figure 2).In the vicinity of values ũ ∈ [ 5,5], the equivalent dispersion transitions nonlinearly between the characteristic dispersion of each species.Taking limits shows that as ũ grows toward positive infinity (c p ≫ cs ) , the dispersion for the component tends to D u → D p , and similarly, when ũ goes toward negative infinity (c p ≪ cs ) , the dispersion for the component tends to D u → D s .This behavior for the equivalent dispersion is in line with remarks from the work of Liu et al. (2011), where the approach of dispersing the component with the transport coefficients of the dominant species was demonstrated to be an effective alternative for simulating transient problems.One important aspect to highlight about Equation 25 is that dispersion coefficients are inherently positive by definition.As a result, the equivalent dispersion also remains positive, regardless of the component's value.The fact that the component is defined as the difference between concentrations of the primary and secondary species reveals that the relative magnitude of concentrations also plays a role in the transport of the component.This role is amplified by the differences in dispersion coefficients.For this particular case, both effects are grouped in the same nonlinear term (first term in the right hand side of Equation 25).One consequence of the nonlinearity of dispersion is that the temporal evolution of a given initial condition for the component is not trivial, and requires keeping track of the component magnitudes at every time step.

Reaction Rate
Starting from the generalized expression for reaction rates (Equation 14), and neglecting strong spatial changes in porosity for simplicity, reaction rates for the binary chemical system are given by Expanding the first term of the right hand side leads to From here, it is straightforward to see that in cases where the dispersion coefficients are equal (D u = D s ), the result from De Simoni et al. ( 2005) is recovered, stating that the rate of reaction is solely determined by the product between a term characteristic of the chemistry (∂ 2 c s /∂u 2 ), and the mixing of the component by hydrodynamic dispersion (∇ T uD s ∇u).Introducing the definition of equivalent dispersion D u (Equation 24) leads to This result demonstrates more clearly that beyond the mixing of components driven by the dispersion of the secondary species, there exist additional nonlinear mixing processes that impact reaction rates.These processes are directly proportional to the differences between dispersion coefficients.While evaluating this expression is not straightforward for the general case where dispersion coefficients vary spatially, some simplified scenarios can be considered.These scenarios involve spatially uniform parameters for each species.

One-Dimensional Diffusive Mixing
A natural advantage of the nonlinear formulation for components, is that the method can be applied to a more extensive range of transport conditions.Particularly, mixing driven solely by diffusion is probably a well-known limit scenario where the transport coefficients are expected to be different for each chemical species.This becomes especially important in aquifers with lenses of clay or other low-permeability materials, which is a common scenario.This section discusses the diffusion-based mixing of two solutions in a one-dimensional homogeneous porous medium column (length L = 1 [m], centered at x = 0).The initial condition is such that a sharp gradient of the component is imposed at the center of the column, simulating a situation where the left-half is dominated by the primary species (c l p > c l s , ũ l > 0), and the right-half of the column is dominated by the secondary species (c r s > cr p , ũr < 0).Effective molecular diffusion coefficients (D p , D s ) are considered to be spatially uniform.Individual values are systematically modified such that the average between them remains constant D = 1 × 10 9 [m 2 / s]) and the ratio D p / D s < 1 adopts different magnitudes.Numerical simulations are performed until reaching a final time of T = 30 [day].For all combinations of diffusion parameters, the same time step is used, determined as δ t = C T δ 2 x / D m , where C T = 0.1, the cell size is δ x = L/500, and D m = 1 × 10 8 [m 2 / s] is a large diffusion scale.The influence of different diffusion coefficients is evaluated by analyzing the concentration profiles at the end of the simulation and the corresponding reaction rate.This is obtained from Equation 28 for a one-dimensional domain, and written in terms of non-dimensional concentrations While interpreting numerical results, all spatial derivatives in the previous expression are evaluated using finitedifferences, once a distribution of the component is obtained from the nonlinear dispersion solver.For validation purposes, the method to compute numerical derivatives is applied to the case in which diffusion coefficients are equal.This case possesses known analytical expressions, as shown in Appendix C. The conditions of equal diffusion (D p / D s = 1) are taken as a baseline scenario for comparing the effect of nonlinearities on the expected output of the transport model.The analytical solution of this case is employed to specify the initial condition for the component, defined as anti-symmetric with respect to the center of the column ũr = ũ l ) and with a slight smoothing (Figure 3a), such that the numerical estimation of the spatial derivatives of the component at the beginning of the simulation is in agreement with the analytical expressions considering the average diffusion coefficient D) .
Concentration profiles from simulations driven by a strong initial gradient of the component ũ l = 5, ũr = 5) illustrate the influence of considering different transport coefficients (Figure 3).In particular, the distribution of the component displays clear deviations with respect to the scenario of equal diffusion, for ratios of diffusion below D p / D s < 0.5.Increasing this difference results in the left-half of the column tending to maintain high values of the component, which is a consequence of the lower diffusion coefficient.This translates in steep gradients which ultimately impacts the reaction rate (Equation 29).In this regard, the transition zone from high to low component values exhibits the major discrepancies with respect to the reference solution.
Notice that the point of zero component concentration ( ũ = 0) shifts toward the region with lower diffusion coefficients as the differences in transport coefficients increase.In the reference case, this point consistently remains at the center of the column.Evaluating the reaction rates provides a further confirmation of this effect (Figure 4a).In this regard, two remarkable aspects are observed: (a) the magnitude of the peak reaction rate is directly proportional to the difference in diffusion coefficients, and (b) the location of the peak reaction rate is consistently shifted from the center of the column toward the region of low molecular diffusion.This result is in agreement with similar studies found in the literature contending that for a binary reaction with initially separated reactants and symmetric concentrations c l p = cr s ) , the reaction front would remain at the center of the column as long as the diffusive fluxes coming from each side of the column are balanced, which occurs when diffusion coefficients are equal (e.g., Jiang & Ebner, 1990;Koza & Taitelbaum, 1996).In the case of different diffusion coefficients and the symmetric initial condition, the reaction front will displace toward the direction of the solute with smaller diffusion (D p in this case).This effect reveals some interesting consequences for the cumulative (time-integrated) reaction product in each cell δ mκ (Figure 4b).In the base simulation with equal diffusion (D p / D s = 1) , the location of the peak reaction rate remains consistently at the center of the column, which ultimately leads to a strong accumulation of the mineral in this place.However, the differences in diffusion produce less concentrated product, with a higher spatial extension toward the region of low diffusion.From the initial state of the system, the peak reaction rate consistently shifts toward the region of lower diffusion.As a result of this effect, the timeintegrated reaction product does not remain concentrated in a fixed location, leading to reduced mineral accumulation at the center of the column.
From the preceding results, it is evident that all combinations of diffusion coefficients result in a predominant mineral precipitation (δ mκ > 0) .This observation aligns with findings in the literature, which suggest that mineral precipitation is the anticipated outcome of a reversible binary reaction with uniform equilibrium constant K and unitary activity coefficients (e.g., De Simoni et al., 2005;Rubin, 1983).However, these remarks consider that the input solutions are either mixed in a batch reactor (fully mixed) or that they are mixing at the same rate (same dispersion).In the case of different transport coefficients, the latter is not the case.Notice in Equation 29that not all terms contribute positively to the total reaction rate.In particular, the second spatial derivative of the component ∂ 2 ũ/∂x 2 ) and the Jacobian derivative (∂ cs /∂ ũ) , are potential sources of negative contributions to the reaction rate for a given difference in diffusion coefficients.If local negative contributions are higher in magnitude than the mixing term related to the secondary species (first term in right hand side of Equation 29), it would be possible to also observe mineral dissolution.To explore this idea, an additional simulation is conducted by reducing the magnitude of the initial condition for the component to ũ l = 0.5 and ũr = 0.5.Within this range of values, the terms related to the chemistry of the problem display a stronger variability with respect to changes in the magnitude of the component.This controls the equivalent component dispersion and Jacobian derivatives, which ultimately dictates how terms related to the differences in transport coefficients contribute to the reaction rates.Both the reaction rates (Figure 5a) and the time-integrated reaction product (Figure 5b) indicate clearly that in cases where D p / D s < 0.25 dissolution is occurring in areas of higher diffusion.In this regard, dissolution is more noticeable in this problem setup, primarily because the amount of precipitation is considerably smaller compared to the previous case.It was also verified for the problem with a strong initial gradient ( ũ = ±5) that some dissolution can be observed in the region of higher diffusion.Nevertheless, this dissolution was negligible in comparison to the amount of precipitation.
As a final addition to these findings, the domain-integrated reaction product Δ mκ is discussed.This quantity is estimated from the time-cumulative spatial distribution of reaction product at the end of the simulation (T = 30 [day]).In this context, results indicate that, in the simulation with stronger component gradients ( ũ = ±5) , the total amount of precipitated product decreases slightly with increasing differences in diffusion (refer to Table 1).When D p / D s = 0.1, the domain-integrated precipitation Δ + mκ ) is about 98% of that achieved in the scenario with equal diffusion coefficients.More clear differences are observed with the initial condition of smaller magnitude ( ũ = ±0.5) .In this case, the total precipitation increases with increasing differences in diffusion with a maximum discrepancy of 74% with respect to the case employing the average diffusion.Regardless of these differences, from a net perspective, the cumulative outcome of the problem is mineral precipitation.As remarked previously, in a problem without gradients of the component at the boundaries, the domain-integrated reaction (Equation 17) will be ultimately controlled by the mixing of the components driven by the equivalent component dispersion (Equation 25 in this problem).Because this dispersion is positive, it is expected that from a global perspective the net product of reaction is precipitation, regardless of the presence of local dissolution patterns.

A Conceptual Interpretation of the Phenomenon
Expanding upon the previous results, this section provides a comprehensive interpretation of the intricate process of mixing solutions characterized by distinct transport coefficients.This interpretation provides an explanation for the mineral dissolution observed in the numerical test cases.As remarked previously, studies have considered that the expected output from a binary chemical reaction in equilibrium is mineral precipitation (assuming unitary activity coefficients).This result has been discussed with an emphasis on its occurrence during mixing in batch reactors (Rubin, 1983), or in cases where the input solutions mix with the same dispersion coefficients (De Simoni et al., 2005).However, simulations showed that the mixing of solutions with species-specific dispersion opens the possibility to mineral dissolution.
To further explore this effect, consider the previous one-dimensional diffusive test case involving the mixing of two distinct waters initially placed into two segments of a homogeneous column.Each water contains two aqueous species (c p , c s ) in equilibrium.The primary species c p is predominantly present in the left-half of the column (x < 0), while the secondary species c s dominates the right-half (x > 0).Based on this setup, the curve of chemical equilibrium, which dictates that c p c s = K in the space of concentrations, is compared with the concentrations of the mixture (as in De Simoni et al., 2005;Sanchez-Vila et al., 2007, 2010).The analytical solution for the species concentrations at a given time t, considering only the mixing of the two waters without any reactions, can be expressed as follows where ( This is illustrated in Figure 6, which shows the parametric curve estimated for different ratios of the diffusion coefficients, considering both the strong ( ũ = ±5, Figure 6a) and the mild ( ũ = ±0.5, Figure 6b) component gradient.Notice that the figure that in cases where both species mix with the same rate (same dispersion), the mixing curve is a straight line, connecting the two chemical signatures, from where the only possible output to reach equilibrium is precipitation.The differences in diffusion give rise to a nonlinear mixing curve with deviations from the straight line proportional to the difference in diffusion coefficients.These curves illustrate the interaction between the multiple factors determining whether the mixture will precipitate or dissolve at a given point x of the column.For example, in the case of a strong initial gradient ( ũ = ±5, as shown in Figure 6a), differences in diffusion coefficients have a limited potential to induce dissolution.In all cases, the mixing curves predominantly remain above the equilibrium line.In cases with strong differences in diffusion coefficients (D p / D s = 0.1) , it is observed that the mixing curve falls below the equilibrium line around the point where the secondary species dominates (c s ≈ 5) .This observation aligns with the small dissolution observed in the simulations of the one-dimensional test case (as well as in the two-dimensional test case discussed in the following).A more interesting perspective arises in the case of a mild initial gradient of the component ( ũ = ±0.5, as depicted in Figure 6b).In this scenario, curves corresponding to diffusion ratios D p / D s < 0.25 distinctly lie below the equilibrium curve, particularly in the region characterized by a chemical signature dominated by the secondary species.These mixing regions are where mineral dissolution originates.Moreover, in concurrence with numerical simulations, dissolution is anticipated as a possible outcome in scenarios where the concentration of the species with the higher diffusion coefficient dominates (as in the case of cs here).Finally, notice that the nonlinear nature of the mixing curves illustrates the escalated reaction rates in areas where the primary species holds dominance, especially when differences in diffusion coefficients are important.

Two-Dimensional Flow-Through System
A second numerical test case considers the application of the nonlinear dispersion model to study different scenarios of reactive transport in a rectangular, homogeneous flow-through system with uniform velocity oriented along the x-coordinate, from left to right.The objective of these simulations is to evaluate the influence that different transport coefficients can have upon the output of chemical reactions, for different scenarios of advectiondominated transport.Dimensions of the aquifer model are L × W = 35 × 8 [cm].Similar setups have been a popular choice in the literature when discussing transverse dispersion through granular media (e.g., Chiogna et al., 2010;Klenk & Grathwohl, 2002;Olsson & Grathwohl, 2007).Based on the results of experiments and numerical simulations, studies have suggested a parameterization of transverse dispersion that takes into account the influence of species-specific diffusion coefficients, properties of the porous medium, and incomplete mixing.For a species i, a proposed model of transverse dispersion follows the expression (Chiogna et al., 2010;Hochstetler et al., 2013;Rolle et al., 2012) where D aq i is the aqueous molecular diffusion, and the product ϕD aq i is an approximated form of pore-diffusion.Pe = vd/ D aq i is the grain Péclet number with v being the average flow velocity, and d the characteristic grain size.Notice that the Péclet number is also species-specific because of the specific molecular diffusion.δ is a parameter of the porous medium representing the ratio between the length of pore channels and their hydraulic radius, and β is a nonlinear exponent representing the effects of incomplete mixing.Although the dispersion model is empirical, it can be seen as an adapted form of a theoretically derived expression (in Bear & Bachmat, 1967).At this stage, the model (Equation 35) will not be questioned and is taken as a reference formulation that has been experimentally validated with parameters δ ≈ 6 and β ≈ 0.5 for spherical particles of diameter d < 2 [mm].The aquifer porosity and grain diameter are both considered spatially uniform with values of ϕ = 0.35 and d = 1.5 [mm], respectively.Similar to the previous test case, the average value of diffusion is fixed at D = 1 × 10 9 [m 2 / s] .In this case, the ratio of coefficients is modified, adopting the values D p / D s ∈ {0.1,0.5,1}.Simulations consider different magnitudes of the uniform flow velocity v ∈ {0.5, 1, 10, 20} [m/day], preserving in all the configurations of diffusion a grain Péclet number Pe > 1 (Table 2).The main feature of the transverse dispersion model (Equation 35), which makes it interesting for the purposes of this work, is that it preserves differences in the magnitude of dispersion coefficients for different species in the whole range of flow velocities (refer to Figure 7).This is a direct consequence of the species-specific diffusion, based on the argument that for the same groundwater flow conditions, the extent of transverse mixing is limited by diffusion (e.g., Rolle et al., 2012).The ratio of transverse dispersion coefficients for species with different diffusion coefficients is a function of the flow velocity, which equals the ratio of diffusion coefficients only in the case that velocity is close to zero (Figure 7b).For intermediate velocities, a transition zone exists where the coefficients approach in magnitude (reference values in Table 2).For very large velocities, an asymptotic difference in dispersion is preserved, proportional to the ratio of molecular diffusion coefficients.It is important to remark that these differences are also influenced by the choice of dispersion parameters D,δ,d,β) , such that the values in Table 2 shall be considered to be specific to this test case.The longitudinal dispersion coefficient is simply taken as D L i = ϕD aq i + vd/2.The adopted model conditions lead to spatially uniform dispersion parameters.Hence, the following expression for reaction rates is obtained from Equation 28 where D L p,s and D T p,s stand for the longitudinal and transverse dispersion coefficients for the primary and secondary species, respectively.From this expression, a more simplified form can be obtained by considering a stationary continuous injection, in which both the longitudinal gradients and the variations in longitudinal dispersion have a negligible influence on reaction rates.This leads to Notice that this expression shares the same form as the one obtained for the problem of diffusive mixing (Equation 29), but now takes into account the differences in transverse dispersion coefficients, as well as the gradients of the component in the transverse direction.The results are discussed in terms of reaction rates normalized by the system pore-volume, that is, r = , where the pore-volume is defined as P V = L/v.
A first exercise with the numerical model (Setup 1) contemplates that the upper half of the aquifer inlet is injecting a solution where the concentration of the secondary species is dominant (c s > cp ) , while the lower half injects with the opposite conditions (c p > cs ) .As in the previous test, two different magnitudes of the components are considered, representing strong differences ( ũ = ±5) and mild differences ( ũ = ±0.5) in the composition of the input solutions.The aquifer is initially saturated with a solution where the primary species is dominant ( ũ > 0) and the simulation progresses until steady-state conditions are reached.A mixing-zone develops at the center of the transverse direction, leading to expected reactions in this area.The stationary distribution of the component is then used to calculate reaction rates using Equation 37. Transverse reaction rate profiles under the simulated conditions are depicted in Figure 8, and the domain-integrated reaction product in one pore-volume is summarized in Table 3.These results provide straightforward insights about the influence of different diffusion coefficients in this problem.
The discussion is first focused on the case with a strong difference in chemical signature ( ũ = ±5; Figures 8a, 8c, 8e, and 8g).In particular, both the spatial distribution of reaction rates and the integrated-reaction product show a great similarity between simulations performed with the averaged diffusion and D p / D s = 0.5.For all the considered velocities, differences in total precipitation remained always below 1% compared to the equal diffusion case.More pronounced differences arise when considering D p / D s = 0.1.Similar to the  Water Resources Research 10.1029/2023WR036358 one-dimensional test case, the reaction peak is slightly shifted from the center toward the region where the primary species dominates (lower diffusion).This shift, however, is small given the conditions of this problem.The shift is more pronounced in the case with lowest flow velocity (v 1 = 0.5 [m/day]; Figure 8a), reaching an extent of approximately two grain diameters (∼2d) near the system outlet.In terms of the total precipitation in one pore-volume Δ + mκ ) , simulations with the highest diffusion difference exhibit larger discrepancies compared to the case of equal diffusion, with deviations of up to 6% in the case with the highest flow velocity (v 4 = 20 [m/ day]).The results from simulations with a mild difference in the composition of inflow solutions ( ũ = ±0.5)were more sensitive to differences in diffusion.In these cases, the total precipitation obtained from simulations with D p / D s = 0.5 is on average around ≈2.5% higher than in the case with equal diffusion.However, when considering D p / D s = 0.1, this difference consistently exceeds 10% compared to the reference model.Furthermore, a clear presence of mineral dissolution throughout the system is observed, concentrated in the region of higher diffusion (upper half of the aquifer, where cs > cp ).
In the second transport problem setup (Setup 2), the continuous injection of a solution where the secondary species dominates (c s > cp ) is considered.Injection occurs through the central part of the domain with a width of w = 1 [cm].The rest of the inlet introduces a solution with the opposite conditions (c p > cs ) .As before, the system is initially dominated by the primary species ( ũ > 0) and evolves until steady-state conditions are reached.In this case, strong transverse gradients develop in the surroundings of the injection.Similar to the previous setup, when considering a strong difference in the chemical signature of inflow solutions (Figures 9a,9c,9e,and 9g; Table 4), simulations with a difference in diffusion of D p / D s = 0.5 generate almost no impact in the total precipitated product with respect to the case of equal diffusion.The case with stronger diffusion differences leads to a symmetric shifting of the peak reaction rates, moving away from the center toward the ambient solution of lower diffusion.This effect is more clear in the case of highest velocity (Figure 9g), where the reaction profile tends to maintain a high relative magnitude across the longitudinal extent of the system.Probably, the more evident consequence of differences in diffusion is the occurrence of a dissolution pattern with peak dissolution rates occurring at the center of the injection.When considering mild differences in the component ( ũ = ±0.5;Figures 9b, 9d, 9f, and 9h), dissolution is of high relative magnitude with respect to precipitation, with advection playing a role in maintaining the presence of dissolution throughout the system.In this sense, the magnitude of the total dissolved mass in one pore-volume (Δ mκ ) is always above 20% of the total precipitated mass Δ + mκ ) for cases with the highest difference in diffusion (D p / D s = 0.1) .The maximum influence occurs when the highest flow velocity is considered (Figure 9h), which shows a relative dissolution magnitude of approximately ≈30% compared to the precipitated mass.This effect is interesting because it ultimately means that when attempting to estimate the total amount of precipitation, one must consider it from a net perspective, that is, the total precipitation product minus the total dissolution.Similar dissolution patterns of smaller magnitude are also observed in simulations with diffusion ratios of D p / D s = 0.5.In this case, the magnitude of the total dissolved mass is on average ≈6% of the total precipitated mass.
In summary, the findings from these simulations suggest that, under certain conditions, the impact of different transport coefficients might not significantly affect the outcomes of chemical reactions.Notably, simulations involving the mixing of solutions where one species strongly dominated the composition ( ũ = ±5) generally demonstrated more resilience to modest differences in diffusion coefficients (D p / D s = 0.5) across all considered flow velocities.The output of reactions presented much more sensitivity to differences in diffusion in cases where the chemical signature of the input solutions was not significantly dominated by one species ( ũ = ±0.5) .For these scenarios, not only the quantitative metrics of the problem were affected, but also in qualitative terms.Most importantly, mineral dissolution arises as an output of relevant magnitude with respect to precipitation, regardless of the considered flow velocity.

Conclusions
This article discussed a nonlinear formulation for multicomponent reactive transport with species-specific diffusion/dispersion coefficients.The analysis was motivated by a recurrent assumption found in formulations of reactive transport based on components, contending that species disperse with the same transport properties through porous media.Scenarios where the diffusion/dispersion coefficients cannot be assumed to be the same for all species are easily found in problems of transport through porous media, for example, for diffusion-dominated regimes occurring in the unsaturated zone and low-permeability materials, or for the modeling of transverse dispersion where the specific molecular diffusion of a given species limits the extent of transverse mixing.Introducing species-specific dispersion coefficients into the formulation of reactive transport based on components still allows to reduce the number of transport equations and eliminates the need for explicitly evaluating reaction-related source terms in the transport of aqueous species.However, the transport equations for the components become nonlinear and coupled to each other.Components are now dispersed by a nonlinear equivalent dispersion that is a function of the chemical speciation.The analysis has enabled the derivation of a new expression for evaluating reaction-related source terms, taking into account the influence of all the different transport coefficients.It was shown that this expression reduces to the well-known equations obtained when Water assuming the same dispersion for all species.The methodology was applied to study different scenarios of reactive transport, considering a binary chemical system that emulates reversible mineral precipitation/dissolution.These illustrated the occurrence of nonlinear mixing processes in direct proportion to the differences in transport coefficients, and influenced by the differences in the chemical signature of the mixing solutions.These nonlinear effects can produce dissolution coexisting with precipitation for a equilibrium binary reaction in dilute systems.Only mineral precipitation would be predicted if assuming that all aqueous species disperse with the same coefficients.A one-dimensional test case of diffusion-dominated transport illustrated the consequences of nonlinear mixing on the spatial distribution of reaction rates and the cumulative product of reaction.In configurations emulating initially separated species with symmetric concentrations, the peak of the reaction rate displaced toward the solution where the species of lower diffusion is dominant.This effect is contrary to the case of equal diffusion coefficients, where the peak reaction rate is always located at the center of the mixing region.A numerical flow-through system was employed to examine the relevance of different dispersion coefficients in the output of reactions under advection-dominated transport.A parameterization of transverse dispersion recently discussed in the literature was employed, which preserves differences in the dispersion coefficients of different chemical species throughout the whole range of velocities.Two scenarios of transport illustrated that the output of reactions presented little sensitivity to mild differences in diffusion in cases that the chemical signature of solutions was strongly dominated by one of the species in the component.As the signature of the mixing solutions became more similar, the reactionrelated source terms became more sensitive to differences in diffusion coefficients, leading to the occurrence of dissolution regions within the system, and influencing the spatial distribution of precipitation.In the long term, these effects have the potential to influence the formation of chemically-driven heterogeneities in the porous medium, leading to different preferential structures in comparison to a case where all the dispersion coefficients are considered the same.The new nonlinear formulation provides a framework to evaluate the impact that differences in diffusion/dispersion coefficients can exert upon the output of reactive transport models.The results discussed in this work allow to expand the applicability of components-based formulations of reactive transport, and could serve as a base to extend existing models considering slow kinetics and mass transfer processes in multicontinuum porous media.

Appendix A: Source Terms With the Same Dispersion Properties
Considering a unique dispersion tensor D in the generalized expression for source terms (Equation 13) leads to where spatial changes in porosity are considered to be of small magnitude for simplicity.Second term in the right hand side of the previous equation can be expanded to obtain ∂ 2 c i ∂u j ∂u k ∇ T u j D∇u k ∑ j ∂c i ∂u j ∇ ⋅ (D∇u j ).

(A2)
Simplifying, leads to an expression for the source term which follows the same form derived in De Simoni et al. (2005).

Appendix B: Nonlinear Dispersion Solver
Nonlinear dispersion is solved with a finite-differences code.Concentration-dependent dispersion is evaluated by a lagged one time step approach (Özişik et al., 2017).The implementation is based on the transport module of MODFLOW (Langevin et al., 2022), employing the MODFLOW-API (Hughes et al., 2022) to update in real-time the memory addresses of the dispersion model.The concentration-dependent solver is benchmarked with the synthetic problem from Fujita (1952), which considers a one-dimensional semi-infinite domain with prescribed unitary concentration in the left boundary and nonlinear diffusion coefficient evaluated as (1 αc) 2 , (B1) where D 0 is a scale of diffusion, c is a normalized concentration, and α ∈ [0, 1] is a parameter controlling the change in diffusion with respect to changes in the concentration.Higher values of α lead to stronger differences in diffusion between regions of high-concentration and low-concentration.Time step for simulations is obtained from the most restrictive condition of diffusion where δ x is the grid size, and the value of C T = 0.1 is adopted.Different values of α expose the solver to different conditions of diffusion variability (Figure B1a).The highest considered value leads to a difference of two orders of magnitude in the diffusion coefficients near the vicinity of the boundary condition.Concentration profiles are shown in Figure B1b in terms of the self-similar variable η = x/2 ̅̅̅̅̅̅̅ ̅ D 0 t √ , displaying a good agreement with the analytical solution for all the considered values of α.
For problems of advection-dispersion, the species-specific transverse dispersion model (Equation 35) is implemented by establishing an analogy with the simple linear dispersion model D T c = ϕD aq c + α T v, which leads to the following expression for the transverse dispersivity Implementation is validated with the analytical solution for a continuous injection into a two-dimensional homogeneous system (Srinivasan et al., 2007), for different magnitudes of velocity and diffusion coefficients (Figure B2), employing the dispersion parameters discussed in the main text.In all cases, calculating the transverse dispersivity with expression (Equation B3) leads to numerical results in agreement with the analytical model.
Numerical approximation of derivatives displays a good agreement with the previous expressions (refer to Figure C1), which gives confidence in the numerical estimates of terms involving spatial derivatives while computing reaction rates.

Figure 1 .
Figure 1.Ratio between the concentrations of primary and secondary species (c p / cs : black, left axis; cs / cp : light-gray, right axis) as a function of the component ũ.Vertical lines indicate ũ = ±0.5 and ũ = ±5.

Figure 2 .
Figure 2. Equivalent component dispersion D u as a function of the nondimensional component ũ, for different order relations between coefficients, considering the same difference.Vertical lines indicate ũ = ±0.5 and ũ = ±5.

Figure 3 .
Figure 3. Simulated profiles for the one-dimensional column at time T = 30 [day], solved with nonlinear component diffusion, for different ratios of D p / D s .Panel (a) presents the distribution of the component, and panels (b) and (c) the corresponding concentrations for the primary cp and secondary cs species.In all panels, the initial distribution is shown as a solid gray line.

Figure 4 .
Figure 4. Product of reaction for the simulation of diffusive mixing in onedimensional column, for different ratios of molecular diffusion.Panel (a) presents the product between the rate obtained from the profiles at the end of the simulation (T = 30 [day]) and the simulation time step δ t .Panel (b) is the time-integrated product of reaction.

)
being c l p,s and c r p,s the initial concentrations of the primary and secondary species on the left-and right-half of the column, respectively.c m p and c m s denote the species concentration resulting only from the mixing of the two waters.This system of equations defines a parametric curve in the c p -c s plane, which determines c m p and c m s of each point of the curve as functions of x, where x ∈ [ L/2, L/2].In general, the mixing curve is a nonlinear function that passes through the points c l p ,K/ c l p ) and K/ c r s ,c r s ) .When D p = D s , the mixing curve is linear and α p = α s = α, representing a straight line passing through these two points, characterized by the following parametric form (c p ,c s ) = c r p ,c r s )

Figure 5 .
Figure 5. Results of diffusive mixing with initial condition ũ = ±0.5.Panel (a) is the reaction product in one time step at the end of the simulation (T = 30 [day]).Panel (b) is the time-integrated product of reaction.

Figure 6 .
Figure 6.Mixing curves for the problem of one-dimensional diffusive transport, considering different ratios of the diffusion coefficients.Panel (a) is the case with a strong difference in the initial condition for the component ( ũ = ±5) , and panel (b) with the mild difference in the initial condition ( ũ = ±0.5).Equilibrium curve is the black-dotted line.

Figure 8 .
Figure 8. Profiles of reaction rates for the problem of transverse mixing in homogeneous flow-through system (Setup 1).First column group results with different velocities for the strong gradient of component ( ũ = ±5; a, c, e, g) , and the second column for the mild gradient ( ũ = ±0.5;b, d, f , h).Spatial coordinates are normalized by the respective domain size.

Figure 9 .
Figure 9. Profiles of reaction per pore-volume, for the problem of centered continuous injection (Setup 2).First column group results for the strong difference in component ( ũ = ±5; a, c, e, g) , and second column for the mild difference in component ( ũ = ±0.5;b, d, f , h).Rows group the results for a given flow velocity, indicated accordingly.Horizontal gray lines mark the width of the injection.Spatial coordinates are normalized by the respective domain size.

Figure B1 .
Figure B1.Analytical nonlinear diffusion from Fujita (1952).(a) Diffusion coefficient as a function of c for different values of α, (b) concentration comparison between numerical (scatter) and analytical (solid lines) in terms of the self-similar variable η.Numerical results are plotted for simulation time T = 30 [day], considering D 0 = 1 × 10 10 [m 2 / s].

Figure C1 .
Figure C1.Comparison between numerical and analytical derivatives for the solution of a one-dimensional diffusive sharp initial condition.Panel (a) first derivative, and panel (b) second derivative.Parameters follow ũ l = 1, ũ l = 1, L = 1 [m], D = 1 × 10 9 [m 2 / s] .Curves are shown for two characteristic times indicated in panel (a).The time t = 2 [hr] corresponds to the initial condition for the problem of one-dimensional diffusional mixing.

Table 1
Integrated Product of Reaction for the Problem of One-Dimensional Diffusive Mixing at Time T = 30 [day] + mκ and Δ mκ are the domain integrated products for precipitation and dissolution, respectively.Average diffusion is D = 1 × 10 9 [m 2 / s].

Table 2
Characteristic Values of the Péclet Number for a Given Flow Velocity, and Average DiffusionD = 1 × 10 9 [m 2 / s]

Table 3
Integrated Product of Reaction in One Pore-Volume for the Problem of Transverse Mixing in Homogeneous Note.Integrated precipitation Δ + mκ ) and dissolution (Δ mκ ) are presented for a given flow velocity, component gradient and ratio of diffusion coefficients.

Table 4
Integrated Product of Reaction in One Pore-Volume for the Problem of Centered Injection in Homogeneous Flow-Through System (Setup 2) Note.Integrated precipitation Δ + mκ ) and dissolution (Δ mκ ) are presented for a given flow velocity, component gradient and ratio of diffusion coefficients.