Multiple Kinetic Parameterization in a Reactive Transport Model Using the Exchange Monte Carlo Method

1 Department of Solid Earth Geochemistry, Japan Agency for Marine-Earth Science and Technology (JAMSTEC), Yokosuka 237-0061, Japan 2 Department of Environmental Studies for Advanced Society, Graduate School of Environmental Studies, Tohoku University, Sendai 980-8579, Japan; okamoto@mail.kankyo.tohoku.ac.jp (A.O.); tsuchiya@mail.kankyo.tohoku.ac.jp (N.T.) * Correspondence: royanagi@jamstec.go.jp; Tel.: +81-46-867-9667


Introduction
The reactive transport model is an essential tool to understand the changes in physical, chemical, and hydrological properties of surface and subsurface environments during water-rock interaction, including weathering, hydrothermal alteration, and metasomatism [1][2][3].The model also plays important roles in predicting the evolution of the system in applications to geothermal exploration, disposal of nuclear waste, and carbon capture and storage [4][5][6][7].The validity and prediction performance of the reactive transport model highly depend on the accuracy of kinetic parameters, such as the dissolution/precipitation rate constants of minerals and the diffusivity of aqueous species, as well as reactive surface areas of the minerals.
Many laboratory experiments were conducted to estimate the dissolution-and precipitation-rate constants of major rock-forming minerals and the diffusivity of aqueous species [8][9][10][11][12][13][14][15][16][17].Nevertheless, the available data of such kinetic parameters are still limited compared to those on the solubility of minerals and equilibrium constants for the dissociation/association of aqueous species [18,19].Water-rock interaction in nature commonly involves several kinds of minerals, and numerous reactions happen simultaneously.Even for monomineralic systems, such as feldspar-water or olivine-water systems in batch experiments, coupled reactions could occur, marked by the dissolution of the primary mineral and the formation of a secondary mineral [9,[20][21][22][23].Accordingly, an effective methodology is required to simultaneously estimate the reaction-rate constants of several minerals and/or diffusivity of elements from complex water-rock systems.
However, multiple parameterizations from a coupled transport and reaction system, pertaining to the dissolution and precipitation of many minerals and the transport of aqueous species (ions and complexes), raise a difficult geochemical inverse problem due to computational difficulties.This is especially true when multiple parameters are to be estimated; estimation sometimes fails to converge to the global minima due to convergence to the local minima through error minimization [24,25].Moreover, the data used for parameterization may be "noisy" due to analytical uncertainties, and therefore, the fitted parameters would also contain errors.Such errors on rate constants are typically ignored to avoid the complex computations required to handle error propagation.Therefore, it would be ideal to develop a versatile method that can extract multiple kinetic parameters and associated errors from noisy laboratory datasets without becoming trapped in local minima.
In this study, we present a novel method to estimate not only multiple reaction-rate constants, but also the diffusivity of aqueous species using only the dataset of spatial mineral distribution, and by combining the reactive transport model and the exchange Monte Carlo (EMC) method [26][27][28].The latter is well-known as an improvement on Markov chain Monte Carlo (MCMC) methods [29].The EMC method is able to estimate a probability distribution on each parameter, which could be used for estimating the uncertainty associated with the parameters more efficiently than the MCMC methods [30], without becoming trapped in local minima [31].In the field of Earth science, the EMC method has been applied to problems pertaining to geophysics [30,32] and rheology [33], and, to the authors' best knowledge, this is the first time that the EMC method has been coupled to equations of a geochemical kinetic model in order to estimate kinetic parameters.The proposed method was tested by examining multiple parameterizations of a synthetized serpentinization system.Using our kinetics-based and data-driven approach, we derived several rate constants for surface reactions, and the diffusivity of aqueous species of primary and secondary minerals was estimated from noisy, observed abundance data.

General Framework
The mass-conservation equation for a chemical species i in solution, and in a reactive transport process within a fluid-saturated porous medium, is generally written as follows [34]: where ϕ, C i , J i , and R chem are the porosity (-), concentration of i (mol/cm 3 solution), total transport flux of aqueous species i, and the rate of gain or loss of chemical species i by surface reactions (mol/cm 3 solution/s), respectively.We consider a situation wherein either spatial variations in minerals or ionic concentrations are observable from experiments or natural settings, and the values of kinetic parameters (such as rate constants of surface reactions) are unknown.Here, a set of N unknown kinetic parameters (Θ) is expressed as Θ = {P 1 , P 2 , . . ., P N }.By running the reactive transport model with Θ, spatial variations in minerals and ionic concentrations can be forward-calculated.Given that the pressure-temperature (P-T) condition and timescales are constrained, Θ is inversely determined from maximum-likelihood estimation by minimizing the difference between the calculated and observable data (E(Θ)) as follows: i denotes each of the M different dependent variables.y i,j , y i,j (Θ), and σ i are the abundance of i at the j-th distance obtained from observable data, the abundance of i at the j-th distance calculated from Θ, and the standard deviation of the observable data for i, respectively.N obs is the number of observations in each dataset.Therefore, the value of E(Θ) denotes the average weighted square error per data point.
Finding a global minimum for E(Θ) and not becoming trapped in a local minimum is the objective of maximum-likelihood estimation.Becoming trapped in a local minimum is a common problem for traditional nonlinear optimization algorithms [24,25], preventing determination of the parameters.In this study, the EMC method [26][27][28]35,36], which is also called parallel tempering [37], is used as an iterative minimization routine to overcome this problem.
As stated previously, the EMC approach is a Markov chain Monte Carlo method [29] and provides more advantages over other approaches as it is able to: (1) obtain samples much more effectively, (2) overcome the trapping problem, and (3) estimate each conditional probability density functions for unknown parameters, which allows the determination of uncertainties that arise from measurement errors in datasets.
The flowchart for estimating Θ using the EMC method is shown in Figure 1.Firstly, candidate values for the new parameters were obtained by adding a random number (n r = [-1,1]) to one of the unknown parameters, and by using log 10 (P i,can ) = log 10 (P i ) + n r , where P i and P i,can are the current unknown parameter and the candidate unknown parameter, respectively.Therefore, the candidate parameter vector is Θ can = {P 1,can , P 2 , . . ., P N } (Figure 1a).E(Θ can ) is then calculated from the candidate values for the parameters, and ∆E is obtained from ∆E = E(Θ) − E(Θ can ), where E(Θ) and E(Θ can ) are the values of E calculated from Θ (={P 1 , P 2 , . . ., P N }) and Θ can (={P 1,can , P 2 , . . ., P N }), respectively.The candidate vector is accepted if the probability is min(1, exp(-∆Eβ n )).β n is the "inverse temperature" on replica n, which controls the acceptance or rejection updates of the candidate values.At low β n , the candidate value is easily updated when ∆E > 0, but this is not the case at high β n .The parameter vector after this candidate update attempt is {P 1,new , P 2 , . . ., P N }, where P 1,new = P 1,can when the candidate value is accepted.However, P 1,new = P 1 is when the candidate value is rejected.After one local update, the next unknown parameter is changed in a similar way.
One Monte Carlo step (MCS) is defined as N trials (i.e., candidate values for N unknown parameters were examined) of this local update (Figure 1b).At the end of each MCS, replicates of the system, which were simultaneously and independently simulated at different β n , were exchanged with probability min(1, exp(-∆E∆β n )) (Figure 1b).This allows the local minima to be overcome, because a system at low β n can provide new local optimizers for a system at high β n , thereby allowing tunneling between the local minima and improving convergence to a global optimum.
When we implement the EMC method, the settings of the value of ∆β n and number of replicas are important since these values have a close relation to the exchange ratio (i.e., the acceptance ratio of the exchange process) [38], which are directly related to the efficiency and computational cost of parameter search.Since a large ∆β n (i.e., a small number of replicas) would result in a low average exchange ratio, ∆β n should not be too large.On the other hand, in order to make the average exchange ratio high, ∆β n has to be small (i.e., a large number of replicas is required), although the associated computational cost is huge.The number of replicas to be used in the EMC method is unclear and completely problem-dependent; they are chosen by trial and error.In a verification of the method described later, we used 20 replicas since >20 replicas are typically used [39,40].The range of β was set from 20 to 350, and ∆β n between each replica was constructed as a geometrical progression, as it would improve sampling efficiency [28].As discussed later, we show that the EMC method can find the true values, given the setting of ∆β n and number of replicas, as described above.Note that the exchange between the replicas occurred frequently during optimization, and the EMC method provided the true answer of multiple unknown parameters.
parameter and the candidate unknown parameter, respectively.Therefore, the candidate parameter vector is Θcan = {P1,can, P2, …, PN} (Figure 1a).E(Θcan) is then calculated from the candidate values for the parameters, and ΔE is obtained from ΔE = E(Θ) − E(Θcan), where E(Θ) and E(Θcan) are the values of E calculated from Θ (={P1, P2, …, PN}) and Θcan (={P1,can, P2, …, PN}), respectively.The candidate vector is accepted if the probability is min(1, exp(-ΔEβn)).βn is the "inverse temperature" on replica n, which controls the acceptance or rejection updates of the candidate values.At low βn, the candidate value is easily updated when ΔE > 0, but this is not the case at high βn.The parameter vector after this candidate update attempt is {P1,new, P2, …, PN}, where P1,new = P1,can when the candidate value is accepted.However, P1,new = P1 is when the candidate value is rejected.After one local update, the next unknown parameter is changed in a similar way.

Model System
As an example for the application of the proposed methodology, we consider here a coupled diffusion and reaction in the olivine-quartz-H 2 O system (Figure 2).This system was chosen for the following reasons: (1) The diffusional metasomatic zoning of talc and serpentine zones between quartz and olivine was classically modeled [41,42].This system is relatively simple, but it involves the dissolution/precipitation processes of several minerals as well as element diffusion; accordingly, it serves as a good example of reactive transport processes in water-rock interactions.(2) The progress of the metasomatic reaction in the olivine-quartz-H 2 O system was experimentally investigated [43], which is a useful when considering a realistic model.(3) The reaction rate is critical to understanding the various processes in subduction zones, including metamorphism, earthquakes, and volcanism [44,45].
The model considers silica metasomatism involving the dissolution-precipitation of talc, serpentine (lizardite), and brucite after olivine in hydrothermal batch experiments.This scenario is modeled as a coupled process involving the diffusion of an aqueous species and surface reactions at 300 • C at a vapor-saturated pressure (P sat ) of 8.58 MPa in a porous medium.The porous medium was initially composed of 63 vol.% olivine (forsterite) grains and 37 vol.%fluid-saturated porosity (Figure 1).Porosity was relatively higher than that of metamorphic rocks, although the high porosity used in this study could be found at surface and subsurface water-rock hydrothermal environments, as well as in hydrothermal experiments with mineral powder.At distance (x) = 0, a source of aqueous silica (i.e., quartz) was present, and the aqueous silica in the solution was transported by diffusion into the porous olivine zone (Figure 2).Therefore, silica concentration changed in time and space due to diffusion and surface reactions.In response to the local silica concentration, several surface reactions could take place simultaneously.

Kinetic Forward Model
The surface reactions in the olivine zone are characterized by the coupled dissolution of olivine and the precipitation of secondary minerals (i.e., talc, lizardite, and brucite).Three aqueous species (i.e., SiO2(aq), Mg 2+ , and H + ) are expected to be involved in surface reactions in the system.Of course, these species could be implemented in the reactive transport model and they would normally be considered.However, several numerical models suggest that the serpentinization system could be modeled with an Mg-fixed frame [46,47].Therefore, in this study, we simplified the reactions and assumed that the surface reactions in the olivine zone could be expressed by the following overall reactions: For the overall reactions involving olivine (i.e., R1, R2, and R4), only the forward reaction was considered since olivine would not be stable under the given P-T conditions.For the overall reactions between secondary minerals (i.e., R3 and R5), both forward and backward reactions were considered.
When the activities of the minerals and water were unity, we considered the first-order equations for the olivine-hosted region to be as follows [34,48]: where kn, Ai, and C n SiO2(aq),eq are the reaction-rate constant (mol SiO2(aq)/cm 2 mineral/s) for the overall

Kinetic Forward Model
The surface reactions in the olivine zone are characterized by the coupled dissolution of olivine and the precipitation of secondary minerals (i.e., talc, lizardite, and brucite).Three aqueous species (i.e., SiO 2(aq) , Mg 2+ , and H + ) are expected to be involved in surface reactions in the system.Of course, these species could be implemented in the reactive transport model and they would normally be considered.However, several numerical models suggest that the serpentinization system could be modeled with an Mg-fixed frame [46,47].Therefore, in this study, we simplified the reactions and assumed that the surface reactions in the olivine zone could be expressed by the following overall reactions: For the overall reactions involving olivine (i.e., R1, R2, and R4), only the forward reaction was considered since olivine would not be stable under the given P-T conditions.For the overall reactions between secondary minerals (i.e., R3 and R5), both forward and backward reactions were considered.
When the activities of the minerals and water were unity, we considered the first-order equations for the olivine-hosted region to be as follows [34,48]: where k n , A i , and C n SiO2(aq),eq are the reaction-rate constant (mol SiO 2(aq) /cm 2 mineral/s) for the overall reaction n, the bulk surface area (cm 2 mineral/cm 3 rock) of mineral i, and the equilibrium silica concentration (mol SiO 2(aq) /cm 3 solution) for reaction n, respectively.For the reactions where both the forward and backward reactions were considered (i.e., R3 and R5), the rate constant for reaction n was expressed as k n+ for the forward reaction and k n-for the backward reaction (Figure 2).The C n SiO2(aq),eq values for R1-R5 were obtained as follows: (1) the equilibrium silica molarities (mol/kg solution) for the overall reactions R1-R5 are calculated from logK obtained from SUPCRTBL [49], assuming that the activities of water and minerals were one; and (2) the molarities are converted to volumetric concentrations (mol SiO 2(aq) /cm 3 solution) by extrapolation of the partial molar volume of SiO 2(aq) (-9.1 cm 3 /mol at 300 • C and P sat ; [50]), as in Watson et al. [16].Log K and C n SiO2(aq),eq values obtained for R1-R5 are summarized in Table 1.Given that the overall reactions involve elementary reactions for dissolution and precipitation, A i represents the surface area of the reactant or product for each overall reaction.In the present study, we introduce effective rate constant k n [34]: Therefore, k n has the unit of s -1 , and k n values for R1-R5 (k R1 , k R2 , k R3+ , k R3-, k R4 , k R5+ , and k R5-) are the fitting parameters in these calculations.
Porosity is updated after calculation at each time step, and it is calculated from the sum of the mineral volume fractions as follows: where m i and V i indicate the amount of mineral (mol/cm 3 rock) and the molar volume of minerals (cm 3 mineral/mol), respectively, which are sourced from Reference [51].N i is the number of minerals in the model.In our calculations, forsterite, talc, lizardite, and brucite were considered; thus, N i = 4.Given that the model considers an Mg-fixed system, the reactive transport model can be written as follows: Minerals 2018, 8, 579 The diffusional transport flux of SiO 2(aq) in a porous medium is given as follows [48,52]: F ∇C SiO 2(aq) (7) where D SiO2(aq) is the diffusion coefficient of SiO 2(aq) in pure solution.Formation factor (F) is defined as in Reference [52], where based on Archie's law [53], where m is the cementation exponent that typically ranges between 1.3 and 2.5 [52].m is one of the fitting parameters in our calculations.Moreover, given that D SiO2(aq) has only been determined for limited P-T conditions [16,[54][55][56], and not for the P-T conditions of the present experiment, D SiO2(aq) is also one of the fitting parameters in our calculations.Consequently, we considered a one-dimensional model for SiO 2(aq) involving diffusion and surface reactions in the olivine region, as follows: This reaction-diffusion equation is solved numerically using the implicit-difference method from x = 0 to 7.0 mm, and the matrix equations are solved using the Thomas algorithm [34].The value of C SiO2(aq) at t = 0 (initial conditions) was uniformly set to 7.12 × 10 -8 mol/cm 3 , similar to the value of C SiO2(aq) obtained after a 24 h run in hydrothermal serpentinization experiments on olivine and NaCl-NaHCO 3 fluid [20].The value of C SiO2(aq) at x = 0 (boundary conditions) was set to 7.41 × 10 -6 mol/cm 3 , which is a quartz-saturated SiO 2(aq) concentration at 300 • C and vapor-saturated pressure.Given that our model does not consider the nucleation of secondary minerals, we added a small amount (1.0 × 10 -7 vol.%) of seed lizardite, talc, and brucite at t = 0.

Synthetic Datasets
The value of Θ estimated by fitting of the reactive transport model contains the following nine parameters (N = 9): To verify our approach, we applied our method to noisy synthetic mineral-abundance datasets.We considered a situation in which data on the number of primary (Forsterite; Figure 3a) and secondary minerals (talc, lizardite, and brucite; Figure 3b-d), and porosity (Figure 3e), as a function of distance, were available (M = 5).Each dataset contained 70 data points (N obs = 70), and the distributions were obtained from the reactive transport model (Equation ( 9)) after 180 h of reaction time when using the following randomly generated parameters: log 10 (D SiO2(aq) ) = -3.0;log 10 (k R1 ) = -4.0;log 10 (k R2 ) = -3.0;log 10 (k R3+ ) = -6.0;log 10 (k R3− ) = -8.0;log 10 (k R4 ) = -1.0;log 10 (k R5+ ) = -2.0;log 10 (k R5-) = -3.0;and log 10 (m) = -0.47.To improve the efficiency of the parameter sampling using the EMC method, the parameter search range was limited from 0 to -10 (in log units) for D SiO2(aq) and k n , and 0 to 0.7 (log units) for m.The observation noise of σ y = 10 -1.5 , which equates to a 6.3% relative deviation of the maximum amount of each mineral, was added to the mineral and porosity data (Figure 3a-e).The noisy datasets (i.e., red data points in Figure 3a-e) were then analyzed using our method to observe its effectiveness in extracting the known reaction-rate constants.

Extraction of Reaction Kinetic Parameters from Noisy Datasets with the Reactive Transport Model
Figure 4 shows probability-distribution histograms for the fitting parameters as the average ± one standard deviation (1σ) of the sampled parameters obtained using the EMC method.To compare each parameter, the bin ranges for DSiO2(aq) and k′n were uniformly set to 0.05 log units, whereas the ranges for m were set to 0.01 log units.
For all estimated parameters, each true value was within the average value (±1σ) of the sampled parameters (Figure 4).However, the histogram shapes were different for some parameters.The histograms of the sampled parameters for DSiO2(aq), k′R1, k′R2, k′R4, k′R5+, and m appeared to be normally distributed, whereas the histograms of the other parameters (k′R3+, k′R3-, and k′R5-) had a uniform distribution with single peaks at around -10 and/or 0 (Figure 4).These single peaks would be an artificial effect due to a boundary condition of the parameter search ranges, which were limited to 0 to -10 in log units.

Extraction of Reaction Kinetic Parameters from Noisy Datasets with the Reactive Transport Model
Figure 4 shows probability-distribution histograms for the fitting parameters as the average ± one standard deviation (1σ) of the sampled parameters obtained using the EMC method.To compare each parameter, the bin ranges for D SiO2(aq) and k n were uniformly set to 0.05 log units, whereas the ranges for m were set to 0.01 log units.
constant for R2 (k′R2), which controls the amount of lizardite (Figure 1), is estimated with low error (log10(k′2) = -2.99 ± 0.05; Figure 4).This is due to the fact that lizardite (0-10 mmol/cm 3 rock; Figure 3c) was the most abundant secondary mineral (cf.0-0.15 mmol/cm 3 rock for talc (Figure 3b); 0.013 mmol/cm 3 rock for brucite (Figure 3d)) and that this reaction also controlled the amount of observed olivine and porosity (Figure 3a,e).This suggests that the main reaction controlling the observed mineral distribution can be qualitatively constrained using its estimated parameter errors.3a-e).Each gray region represents the average ± one standard deviation of the parameters sampled using the EMC method.Dotted red lines show the true values of each parameter described in Section 3. To allow comparison of each parameter, the bin ranges for DSiO2(aq) and k′n were uniformly set to 0.05 log units, whereas the range for m was set to 0.01 log units.
In contrast, for parameters with uniform distributions (k′R3+, k′R3-, and k′R5-; Figure 4), the proposed method was unable to estimate the true values.This is because the parameters have a broad minimum, and therefore, any value could result in a good fit.This reflects the fact that these reactions (i.e., forward and backward reactions for R3, and backward reaction for R5) do not affect the amounts of different minerals present, and these reactions could be excluded from the model for simplicity.In fact, the backward reaction of R5 with a reaction-rate constant of k′R5only occurs when CSiO2(aq) < C R5 SiO2(aq),eq, which was never the case in our calculations (Figure 5).Therefore, the proposed method could also exclude this reaction path (i.e., k′R3+, k′R3-, and k′R5-could be set to zero), as it is not important in reconstructing the observed mineral distributions.3a-e).Each gray region represents the average ± one standard deviation of the parameters sampled using the EMC method.Dotted red lines show the true values of each parameter described in Section 3. To allow comparison of each parameter, the bin ranges for D SiO2(aq) and k n were uniformly set to 0.05 log units, whereas the range for m was set to 0.01 log units.
For all estimated parameters, each true value was within the average value (±1σ) of the sampled parameters (Figure 4).However, the histogram shapes were different for some parameters.The histograms of the sampled parameters for D SiO2(aq) , k R1 , k R2 , k R4 , k R5+ , and m appeared to be normally distributed, whereas the histograms of the other parameters (k R3+ , k R3-, and k R5-) had a uniform distribution with single peaks at around −10 and/or 0 (Figure 4).These single peaks would be an artificial effect due to a boundary condition of the parameter search ranges, which were limited to 0 to −10 in log units.
For normally distributed parameters (D SiO2(aq) , k R1 , k R2 , k R4 , and k R5+ ; Figure 4), each true value was within the average value (±1σ) of the sampled parameters, indicating that the true parameters were successfully estimated despite errors from the noisy datasets.In particular, the reaction rate constant for R2 (k R2 ), which controls the amount of lizardite (Figure 1), is estimated with low error (log 10 (k 2 ) = -2.99 ± 0.05; Figure 4).This is due to the fact that lizardite (0-10 mmol/cm 3 rock; Figure 3c) was the most abundant secondary mineral (cf.0-0.15 mmol/cm 3 rock for talc (Figure 3b); 0.013 mmol/cm 3 rock for brucite (Figure 3d)) and that this reaction also controlled the amount of observed olivine and porosity (Figure 3a,e).This suggests that the main reaction controlling the observed mineral distribution can be qualitatively constrained using its estimated parameter errors.In contrast, for parameters with uniform distributions (k R3+ , k R3-, and k R5-; Figure 4), the proposed method was unable to estimate the true values.This is because the parameters have a broad minimum, and therefore, any value could result in a good fit.This reflects the fact that these reactions (i.e., forward and backward reactions for R3, and backward reaction for R5) do not affect the amounts of different minerals present, and these reactions could be excluded from the model for simplicity.In fact, the backward reaction of R5 with a reaction-rate constant of k R5-only occurs when C SiO2(aq) < C R5 SiO2(aq),eq , which was never the case in our calculations (Figure 5).Therefore, the proposed method could also exclude this reaction path (i.e., k R3+ , k R3-, and k R5-could be set to zero), as it is not important in reconstructing the observed mineral distributions.The reaction paths identified as not being important are typically ignored in other studies, such as those using the forward-modeling data approach.However, such arbitrary exclusion of reactions may be based on subjective and/or erroneous assumptions.In contrast, our proposed method can objectively identify and exclude these unimportant reaction paths by producing a distribution of the most probable values for each parameter.8), with the parameters given in Section 3).aSiO2(aq) (and thus CSiO2(aq)) is always larger than C R5 SiO2(aq),eq.

Dependence of Estimation Accuracy on Observation Noise
To investigate how observation noise affects the estimation performance of our proposed method, we examined the estimation accuracy for fitted data with different noise levels of σy = 10 -2.0 and 10 -1.2 , which corresponded to 2.0% and 12.6% deviations, respectively, relative to the maximum amount of each mineral.The fitted noisy datasets are shown in Supplementary Figure S1.
Figure 6 shows the noise dependency of the estimated values.Parameters were successfully estimated even when the observed data had noise levels of σy = 10 -1.2 .The errors were low when the observable noise was low.For example, log10(DSiO2(aq)) values were estimated to be -2.96± 0.06 at σy = 10 -2.0 , -2.96 ± 0.24 at σy = 10 -1.5 , and -2.80 ± 0.26 at σy = 10 -1.2 (Figure 6).The reaction paths identified as not being important are typically ignored in other studies, such as those using the forward-modeling data approach.However, such arbitrary exclusion of reactions may be based on subjective and/or erroneous assumptions.In contrast, our proposed method can objectively identify and exclude these unimportant reaction paths by producing a distribution of the most probable values for each parameter.

Dependence of Estimation Accuracy on Observation Noise
To investigate how observation noise affects the estimation performance of our proposed method, we examined the estimation accuracy for fitted data with different noise levels of σ y = 10 -2.0 and 10 -1.2 , which corresponded to 2.0% and 12.6% deviations, respectively, relative to the maximum amount of each mineral.The fitted noisy datasets are shown in Supplementary Figure S1.
From these results, it is evident that our proposed method is robust to noise, with the parameter accuracy for several reaction-rate constants and silica diffusivity being excellent over a wide range of noise intensities.Moreover, our method identified and excluded noneffective parameters from the model.From these results, it is evident that our proposed method is robust to noise, with the parameter accuracy for several reaction-rate constants and silica diffusivity being excellent over a wide range of noise intensities.Moreover, our method identified and excluded noneffective parameters from the model.

Summary and Outlook
We developed a probabilistic parameter-estimation method based on the reactive transport model using the EMC method.Accuracy tests on synthetic datasets demonstrated that the method can successfully estimate reaction rate constants and silica diffusivity as a probability distribution, rather than a single value, even when the fitted dataset is noisy.The proposed method has the following advantages: (1) multiple rate constants can be estimated from noisy datasets; (2) the uncertainties of rate constants can be determined; and (3) the main reaction paths controlling the observed mineral abundances can be identified, and unimportant reaction paths can be objectively identified and excluded, rather than being arbitrarily ignored.
In the present study, small seeds of lizardite, talc, and brucite were added at t = 0, since our model did not consider the nucleation of secondary minerals.When therate constants from a real dataset to be estimated, the seed volume would affect the model estimates.Therefore, it is required to estimate the seed volume of each mineral or consider the nucleation theory in the reactive transport model.If some unknown parameters appear in the model (e.g., mean surface energy of secondary mineral), parameters could be also estimated by the EMC method.In such cases, an accuracy test must be done with synthetic datasets to ensure the EMC method could find true answers.
When the number of observable data points (Nobs) decreases, the error on parameters estimated by the proposed methodology would increase, as suggested by previous studies on parameter estimation [57,58].Although the effects of sparseness in the dataset on the estimation error were not

Summary and Outlook
We developed a probabilistic parameter-estimation method based on the reactive transport model using the EMC method.Accuracy tests on synthetic datasets demonstrated that the method can successfully estimate reaction rate constants and silica diffusivity as a probability distribution, rather than a single value, even when the fitted dataset is noisy.The proposed method has the following advantages: (1) multiple rate constants can be estimated from noisy datasets; (2) the uncertainties of rate constants can be determined; and (3) the main reaction paths controlling the observed mineral abundances can be identified, and unimportant reaction paths can be objectively identified and excluded, rather than being arbitrarily ignored.
In the present study, small seeds of lizardite, talc, and brucite were added at t = 0, since our model did not consider the nucleation of secondary minerals.When therate constants from a real dataset to be estimated, the seed volume would affect the model estimates.Therefore, it is required to estimate the seed volume of each mineral or consider the nucleation theory in the reactive transport model.If some unknown parameters appear in the model (e.g., mean surface energy of secondary mineral), parameters could be also estimated by the EMC method.In such cases, an accuracy test must be done with synthetic datasets to ensure the EMC method could find true answers.
When the number of observable data points (N obs ) decreases, the error on parameters estimated by the proposed methodology would increase, as suggested by previous studies on parameter estimation [57,58].Although the effects of sparseness in the dataset on the estimation error were not tested, in this study's example at least, the effect of the amount of data was not large since the spatial patterns of individual minerals are rather smooth.
The result of parameter fitting may be poor despite the application of the EMC method for one of the two following reasons: trapping in local minima or an unrealistic forward model.One of the strongest advantages of using the proposed method is that it can overcome local minima problems.Therefore, in this case, we conclude that the latter caused bad fitting, and other forward models should be similarly tested.Thus, it would be necessary to conduct iteration tests on many forward models and check the fitting results to objectively identify the most realistic one.This method is called model selection in the field of information science [59], and it will likely be a useful concept to objectively evaluate forward models.
To date, many reaction-rate constants of mineral dissolution and precipitation have been estimated from laboratory experiments and natural settings.However, even in simple dissolution reactions involving pure phases, many discrepancies are observed among datasets collected under similar conditions.Recently, Fischer et al. [60] suggested that a probabilistic approach must be incorporated in the mineral dissolution process for quantitative analysis of the reactivity of the materials.Probabilistic components causing spatial-temporal heterogeneity in the dissolution rate are potentially present in the dissolution process (e.g., the complex interaction of the areal distributions and concentrations of reactive sites, defects, and grain boundaries [60][61][62]).When the rate of the overall reaction is to be estimated, potential probabilistic components would increase since the overall reaction involves complex and dynamically interacting processes.Therefore, we suggest that it is more reasonable to estimate each rate constant in the overall reaction as a probability distribution, and the EMC used in this study could be a strong tool for estimating the probability distribution of multiple kinetic parameters We applied the reactive transport and the EMC method to an artificial system of silica metasomatism and serpentinization.Typically, the datasets that could be used for kinetic inversion from the reactive transport system relate to (1) spatial variation in minerals or (2) spatial variation in solution chemistry.We showed that multiple parameters could be estimated from spatial variation in minerals.Maher et al. [10] applied the reactive transport model to U isotopes and pore-fluid chemistry in marine sediments, and successfully estimated the dissolution rate of plagioclase by fitting a model to noisy natural chemical datasets.Given that our method can extract multiple independent variables as a probability distribution from natural datasets, it could facilitate modeling of more complex reactive transport systems.Thus, improved prediction of a variety of geological phenomena and physical and chemical responses during water-rock interaction in various settings, such as weathering processes, geothermal exploration, disposal of nuclear waste, and carbon capture and storage, is possible over a wide range of spatial and temporal scales.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2075-163X/8/12/579/s1. Figure S1: S Synthetic datasets with noise levels of (a) σ y = 10 -2.0 and (b) σ y = 10 -1.2 used in the accuracy tests.Gray lines show the values predicted using the reactive transport model (Equation ( 9), with the parameters described in Section 3) and the red data points are the noisy data.

Figure 1 .
Figure 1.Schematic flow chart of the exchange Monte Carlo method for four replicate systems.(a) Flow chart of local update.(b) The local update was repeated after the exchange process.See Section 2.1 for details.MCS: Monte Carlo step.

Figure 1 .
Figure 1.Schematic flow chart of the exchange Monte Carlo method for four replicate systems.(a) Flow chart of local update.(b) The local update was repeated after the exchange process.See Section 2.1 for details.MCS: Monte Carlo step.

Minerals 2018, 8 , 15 Figure 2 .
Figure 2. Schematic illustration of the reactive transport model considered in this study.At x = 0, SiO2(aq) was transported by diffusion at a rate of DSiO2(aq) through the porous media of powdered forsterite, and secondary-mineral (talc, lizardite, and brucite) was produced.The gray bold line shows the overall reaction between two minerals.See Section 2.1 for details.

Figure 2 .
Figure 2. Schematic illustration of the reactive transport model considered in this study.At x = 0, SiO 2(aq) was transported by diffusion at a rate of D SiO2(aq) through the porous media of powdered forsterite, and secondary-mineral (talc, lizardite, and brucite) was produced.The gray bold line shows the overall reaction between two minerals.See Section 2.1 for details.

Figure 3 .
Figure 3. Synthetic datasets with a noise level of σy = 10 -1.5 used in the accuracy tests.(a-e) Plots with gray lines show true values predicted using the reactive transport model (Equation (9), with the parameters given in Section 3) for (a) forsterite, (b) talc, (c) lizardite, (d) brucite, and (e) pores.Red data points are synthetic noisy data.

Figure 3 .
Figure 3. Synthetic datasets with a noise level of σ y = 10 -1.5 used in the accuracy tests.(a-e) Plots with gray lines show true values predicted using the reactive transport model (Equation (9), with the parameters given in Section 3) for (a) forsterite, (b) talc, (c) lizardite, (d) brucite, and (e) pores.Red data points are synthetic noisy data.

Figure 4 .
Figure 4. Histograms showing estimated candidate values sampled using the exchange Monte Carlo method from the noisy datasets (σy = 10 -1.5 ; Figure3a-e).Each gray region represents the average ± one standard deviation of the parameters sampled using the EMC method.Dotted red lines show the true values of each parameter described in Section 3. To allow comparison of each parameter, the bin ranges for DSiO2(aq) and k′n were uniformly set to 0.05 log units, whereas the range for m was set to 0.01 log units.

Figure 4 .
Figure 4. Histograms showing estimated candidate values sampled using the exchange Monte Carlo method from the noisy datasets (σ y = 10 -1.5 ; Figure3a-e).Each gray region represents the average ± one standard deviation of the parameters sampled using the EMC method.Dotted red lines show the true values of each parameter described in Section 3. To allow comparison of each parameter, the bin ranges for D SiO2(aq) and k n were uniformly set to 0.05 log units, whereas the range for m was set to 0.01 log units.

Figure 6 .
Figure 6.Noise effects on the accuracy of the estimated parameters.The dotted red lines show the true values of each parameter described in Section 3. The parameter errors are low when the observable noise (σy) is low.

Figure 6 .
Figure 6.Noise effects on the accuracy of the estimated parameters.The dotted red lines show the true values of each parameter described in Section 3. The parameter errors are low when the observable noise (σ y ) is low.

Author
Contributions: R.O. conducted the overall study, and N.T. and A.O. supervised it.R.O. wrote the original draft, and all authors commented on the manuscript.Funding: This work was supported by JSPS KAKENHI (Grant Nos.18J01649, 16H06347, and 17H02981), ERI JURP (Grant No. 2018-B-01), and JST PRESTO (Grant No. JPMJPR1676).

Table 1 .
Equilibrium constants of the overall reaction considered in the reactive transport model.
[49]ue at 300 • C and P sat obtained from SUPCRTBL[49].ξ Calculated from log K at 300 • C and P sat (see main text for details).