Mineral paragenesis on Mars: The roles of reactive surface area and diffusion.

Geochemical models of secondary mineral precipitation on Mars generally assume semiopen systems (open to the atmosphere but closed at the water-sediment interface) and equilibrium conditions. However, in natural multicomponent systems, the reactive surface area of primary minerals controls the dissolution rate and affects the precipitation sequences of secondary phases, and simultaneously, the transport of dissolved species may occur through the atmosphere-water and water-sediment interfaces. Here we present a suite of geochemical models designed to analyze the formation of secondary minerals in basaltic sediments on Mars, evaluating the role of (i) reactive surface areas and (ii) the transport of ions through a basalt sediment column. We consider fully open conditions, both to the atmosphere and to the sediment, and a kinetic approach for mineral dissolution and precipitation. Our models consider a geochemical scenario constituted by a basin (i.e., a shallow lake) where supersaturation is generated by evaporation/cooling and the starting point is a solution in equilibrium with basaltic sediments. Our results show that cation removal by diffusion, along with the input of atmospheric volatiles and the influence of the reactive surface area of primary minerals, plays a central role in the evolution of the secondary mineral sequences formed. We conclude that precipitation of evaporites finds more restrictions in basaltic sediments of small grain size than in basaltic sediments of greater grain size.

Introduction

Despite several geochemical models of aqueous mineral precipitation on Mars have considered the interactions between the gas phases (H2O, CO2, Cl2, and SO2) and the solution through the water‐atmosphere interface, they generally do not account simultaneously for the transport of ions through the water‐sediment (also referred to as water‐basalt) interface. In closed or semiopen models, the fraction of ions derived from the weathering of primary minerals, either cations (mostly Fe2+, Fe3+, Mg2+, Ca2+, Al3+, and Na+) or anions (particularly SiO3 2− and SiO4 4−), is forced to remain constant in the form of dissolved species or to precipitate as secondary minerals (i.e., clays and evaporites). While closed system models can be useful to describe the mass balance between primary and secondary mineral assemblages, they cannot explain the net balance between ions released after the dissolution of basalt and the total ions incorporated into clays, which results in a residual amount of cations in solution that cannot bind to clays due to the stoichiometric constraints imposed by the limited availability of the silicate anion [Milliken et al., 2009]. Provided the existence of an alternative source of anions, mainly derived from volcanic volatiles [Halevy and Head, 2014], it is expected that these ions would have bonded in evaporites (i.e., sulfates). However, in most localities on Mars, evaporites are commonly absent in clay‐dominated sediments [Bandfield et al., 2003; Bishop et al., 2008; Ehlmann et al., 2008; Osterloo et al., 2008]. As a particular example, evaporites in Gale crater are the result of postdepositional fluid migration, in different late‐stage episodes of fluid flow; therefore, clays and evaporites at Gale have very different depositional histories, occurring spatially closely but with deposition times separated by hundreds of millions of years [Nachon et al., 2014; Rapin et al., 2016]. If clays and evaporites had been formed coevally in the same area on Mars, it should be expected that both secondary minerals would be exposed to similar chemical and physical processes, including dissolution, transport, and reprecipitation in younger terrains [Milliken et al., 2009], and to the effects of dust covering (dust creates an observational bias when analyzing Mars surface composition [Berger et al., 2016] and could even become a relevant input for the mineralogical composition of a given region [Kulowski et al., 2017], affecting equally clays, evaporites, or any other mineral on the surface of Mars).

The reason for this apparent disagreement between geochemical models and direct observational data is twofold. First is a consequence of considering closed or semiopen system conditions in the geochemical modeling, which do not account for simultaneous precipitation and transport processes, particularly through the water‐sediment interface; flow‐through models may be better suited to this end. And second, because equilibrium models do not take into account the role of reactive surface area in determining the differences on crystallization pathways from small‐to‐large sediment particle sizes (kinetic models are needed to accomplish this objective). To help solve this problem, we have applied a reactive‐transport approach to model the formation of mineral sequences known to exist on Mars, as described from lander and orbiter data [e.g., Hamilton and Christensen, 1997; Bandfield et al., 2000; Mustard et al., 2008; McSween et al., 2009; Ehlmann et al., 2011; Vaniman et al., 2014; Rampe et al., 2017; Hurowitz et al., 2017], considering open system conditions both at the atmosphere‐water and water‐rock interfaces and implementing a kinetic approach for the dissolution and precipitation of solid phases. We present a suite of models representing aqueous environments in systems under an early Mars atmosphere, to analyze the role of (i) the reactive surface area of primary minerals; (ii) the reactive transport processes through the basalt interface; and (iii) the cation adsorption by the surface of clays.

Geochemical Setting

Initial Conditions

Geochemical calculations were done with the Phreeqc software [Parkhurst and Appelo, 1999, 2013]. This is a specialized code that applies the mass balance and electroneutrality equations to the simulation of geochemical processes, implementing algorithms for kinetic modeling of dissolution‐precipitation processes, charge balance by sorption on mineral surfaces, and reactive transport along 1‐D sediment columns. Phreeqc is an open source code particularly suitable to do this work because it allows the building of reactive transport models, in particular incorporating kinetic rates and allowing the incorporation of specific pieces of code through its BASIC interface (for a detailed comparison between reactive transport codes, please see Steefel et al. [2015]).

The composition of Martian basalts was obtained from global remote sensing observations [Hamilton and Christensen, 1997; Bandfield et al., 2000; McSween et al., 2009]. Although shergottites and in situ chemical data (MER and MSL) are available, those likely represent local assemblages. Instead, we prefer using global average compositions for the purposes of this work, as we are not attempting to reproduce any particular mineral sequence on a given location on Mars where boundary conditions would be restricted to the already known geochemical information of that area (e.g., the composition of the host‐rock, the pH of the starting fluid, or the influence or not of atmospheric gases). For that, there is already abundant published modeling efforts focused on explaining how local and particular mineral sequences could have formed. As mean values, we used a mineralogy consistent of 65% volume plagioclase (andesine), 27% pyroxenes (orthopyroxenes and clinopyroxenes), 6% olivine (forsterite and fayalite), and 2% of other components. As the starting solution we used water (pH = 7, pe = 4, and T = 273 K) equilibrated with the basalt mineral assemblage, leading to the composition listed in Table 1. This results in a number of ions in solution that could induce the precipitation of a suite of phyllosilicates similar to that reported on Mars [Mustard et al., 2008; Ehlmann et al., 2011], including kaolinite, chlorites, and smectites. The kinetics of dissolution and precipitation of mineral phases was modeled using equation (1).

Table 1

Starting Solution Used in the Modelsa

Element	Molality (mol/kg Water)
Al	1.738 × 10−4
Ca	2.285 × 10−4
Na	1.285 × 10−4
K	1.285 × 10−4
Fe2+	3.671 × 10−5
Mg	4.268 × 10−5
P	4.284 × 10−5
S2−	5.447 × 10−7
Si	1.501 × 10−4
Physicochemical conditions
Temperature (K)	273
pH	11
pe	4

Dissolved ions derived from pure water equilibrated with basalt. Concentrations are expressed as moles of mineral per kg of interstitial water.

In our models the system is open to the atmosphere, and therefore, the loss of water through evaporation is regarded as the main factor inducing supersaturation (as multiple indicators for evaporation have been reported on Mars) [see e.g., Squyres et al., 2004; Sears and Moore, 2005; Tosca et al., 2005; Osterloo et al., 2008; Fairén et al., 2009; Arvidson et al., 2014; Toner et al., 2015; Schwenzer et al., 2016; Hurowitz et al., 2017]; and because the system is likewise open at the water‐sediment interface, reactive transport was also considered (Figure 1). In addition, we considered that dissolved volcanic gases (H2O, CO2, Cl2, and SO2) [e.g., Halevy and Head, 2014] are the main anion source to the aqueous solutions, as well as drivers of redox processes and alkalinity, modifying the dissolution kinetics of primary minerals, which, in turn, are largely dependent of pH, temperature, and the reactive surface. The total reactive surface area of basalt minerals is assumed to be variable depending on the rate of physical weathering of the substrate (see Figure S1 in the supporting information). Kinetic modeling also allowed us to analyze the specific role of the reactive surface area of primary minerals on pH and the subsequent evolution of the system leading to different sequences of secondary mineral precipitation, as well as the evolution of the physicochemical conditions of the system.

Figure 1

Conceptual framework employed in our models to characterize the processes of mineral dissolution/precipitation at the water‐basalt interface of a water body in early Mars. The ionic strength of the upper water layer, in contact with the atmosphere, increases by evaporation and brine rejection, subsequent to ice formation. Denser water accumulates at the water‐basalt interface by the formation of density gradients along the water column. We consider two scenarios (see text) at the water‐sediment interface: (1) if the basaltic basement is not fractured, the absence of diffusion—and diagenesis—restricts the mineral dissolution/precipitation to authigenic processes. In this scenario, we evaluate two conditions: high and low reactive surface of minerals at the interface. And (2) when the solid substrate is fractured at depth, it results in a certain degree of porosity (ϕ), inducing molecular diffusion and the triggering of diagenetic processes.

We analyzed a scenario corresponding to an ocean or lake, assuming that surface water is lost by evaporation/sublimation at different rates. In this situation, a feedback is expected between evaporation and density flux, driving the more saline water to the lower layers [Fairén et al., 2011], in a similar way as is observed in evaporating platforms on Earth. This process produces a “brine reflux” [Jones et al., 2002], in which top waters concentrated by evaporation flow downward, generating density currents and displacing less dense underlying groundwaters (Figure 1). Once high saline waters reach the water‐sediment interface, brines continue to flow into the sediment column, in a process known as “latent reflux” [Jones et al., 2002; Al‐Helal et al., 2012], leading to the precipitation of salts and inducing the persistence of brines in the sediment [Marion et al., 2003; Toner et al., 2015]. At this point, the formation of secondary phases is controlled by diagenetic processes, following a dynamic feedback response between diffusion and reaction. This, in turn, can generate a very complex distribution of diagenetic patterns, in the form of precipitation fronts along the sediment column [Ortoleva, 1994; Berner, 1980], depending on several variables such as sediment porosity, reactive surface, counter diffusion of gases, and adsorption processes.

Our diffusion models are based on a 1‐D sediment column formed by a set of 20 identical cells, 1 cm each. The two upper cells initially contain a solution in equilibrium with basalt, but without solid phases. The next 18 cells down represent a system defined by the interstitial water + mineral phases of the sediment column (which can dissolve or precipitate during the simulation time) (Figure 1). The mean depth of the basal boundary layer in the terrestrial oceanic floor is also of 20 cm [Boudreau and Jørgensen, 2001], and this is the region where early diagenesis is more intense and where the main mineral reactions occur. For modeling purposes, at t = 0, the conditions are the same in all cells, representing the water‐basalt equilibrium mixture before any process of self‐organization (Figure 1). In our models, diffusion starts when a process of evaporation in the upper part of the system increases the ionic concentration of the two upper cells. This situation results in a net vertical flux of ions to the lower cells, triggered by the concentration gradient according to the Fick's laws (equation (10)) [Fick, 1855; Brogioli and Vailati, 2001]. The evolution of the solution through time in the upper layers will be dependent upon the velocity of evaporation, the diffusion coefficient of each ion, and the temporal relationships between evaporation and diffusion processes.

Our models analyze the effect of the reactive surface area on the geochemical evolution of a lake/ocean under different rates of evaporation (assuming a free convection approach) and cooling [Fairén et al., 2011]. The influence of the reactive surface area was tested by using the same interaction with volatile phases in the atmosphere and the same mass amount of basaltic phases in the mineral assemblage. Thus, the only variable parameters in the models are the size of basalt fragments and the generation rate of supersaturation, which is controlled by the rates of evaporation and/or cooling (see Figures 2 and 3 for details about the parameter setting).

Figure 2

(a) Evolution toward alkaline and reducing endpoint conditions when aqueous reaction of small fragments of basalt (cubes, L = 12 μm; S reactive = 480 m2/kgw; ϕ = 0.9) is coupled with a low generation rate of supersaturation induced by simultaneous evaporation (at constant pressure p atm = 0.1 atm) and fast cooling (until 220 K (dT/dt = 27 K/yr)). The remaining mass of water in the figure was calculated from equation (5), with S/V = 0.03. (b) Mineral precipitation sequence with the conditions of Figure 2a. The rapid basalt dissolution offsets the acidification tendency imposed by the atmospheric gases. In the end, the solution remains alkaline and reductive, and only clays and iron oxides precipitate.

Figure 3

(a) The same analysis as in Figure 2a but considering larger basalt fragments (cubes L = 11 cm, S reactive = 0.051 m2/Kgw). In addition, a faster generation rate of supersaturation is induced by fast evaporation—at a constant pressure p atm = 0.07 atm (and slow cooling) until 255 K (dT/dt = 2 K/yr). (b) Mineral precipitation sequence with the conditions of Figure 3a. In such scenario, the slow dissolution of basalt cannot counteract the evolution of the solution toward acidic and oxidizing conditions. Precipitation of clays and iron hydroxides occurs at the initial stages of evaporation, when pH remains in circumneutral values and the temperature is still relatively high. However, the progress of evaporation and the evolution of the system toward acidic conditions induce the precipitation of salts.

Our results are expressed in centimeters after considering time periods ranging from some days to 102 years, so we can show that, even in such short time intervals, the proposed mechanism is able to induce undersaturation in a solution initially in equilibrium with basalt. Indeed, if the initial concentration at the surface and/or the selected time interval is greater, the length of the diffusion gradient will be also greater, but the main rationale and conclusions will remain the same. That is the actual case for the surface of Mars, where the final diffusion processes seem to have been acting on very concentrated ponds generated after vigorous processes of evaporation and for extended time periods [e.g., Tosca et al., 2005; Tosca and McLennan, 2006; Fairén et al., 2009; Fairén, 2010]. But even a short duration of the diffusion gradient would have not precluded the existence of sediment advection or erosion. Our models analyze processes in a lower range of temperature values, relevant for weathering; however, we note that hydrothermal processes would not be incompatible with the specific algorithms employed in our models, as our model conditions are valid up to 573 K.

Kinetics of Mineral Dissolution‐Precipitation

Basalts are volcanic rocks consisting of a number of different minerals and amorphous phases, each one with a specific reactivity and surface area, explaining the extremely complex dissolution behavior of basalts in natural systems [Nesbitt and Wilson, 1992; Sak et al., 2004; Stopar et al., 2006; Hausrath and Brantley, 2010; Gudbrandsson et al., 2011]. Kinetic constraints derived from the three main factors affecting the rate of dissolution of silicates (namely, pH, temperature, and reactive surface area, see Figure S2) can greatly influence the basalt weathering pathways as well as the resulting suite of secondary minerals (Table 2). Geochemical modeling of basalt weathering is usually attained by mean of rate expressions that take into account these dependencies, as for example, the general rate equation of Lasaga [1995], based in the transition state theory, where each term in the sum represents a single reaction mechanism: where dm/dt is the mineral dissolution rate (mol kgwater −1 s−1); SA is the reactive surface area (m2 kgwater −1); K , K , and K are the specific rate constants for dissolution in acidic, neutral, or basic media (mol m−2 s−1); a is the proton's activity; and the exponents n 1 and n 2 are connected with the reaction order, indicating that the reaction is either catalyzed by protons (positive exponent) or inhibited (negative exponent); R (JK−1 mol−1) is the gas constant; T the absolute temperature (K); Ω is the mineral relative supersaturation index; the exponents p and q are empirical constants related to the dependence of dissolution rates on supersaturation; and E , E , and E are the activation energies in acidic, neutral, or basic media, respectively (KJ/mol). The kinetic data used in this investigation (Table 1) are from the compilation by Palandri and Kharaka [2004]. For the precipitation rate constants we applied the principle of reversibility [Lasaga, 1995].

Table 2

Kinetic Data Used in the Modelsa

Primary Minerals: Basalt Phases
	Acid Mechanism			Neutral Mechanism
Mineral Phases	logk	E	n	logk	E
Anorthite CaAl2Si2O8	−3.50	16.6	1.411	−9.12	17.8
Diospside CaMgSi2O6	−6.36	96.1	0.710	−11.11	40.6
Enstatite Mg2Si2O6	−9.02	80.0	0.600	−12.72	80.0
Fayalite Fe2SiO4	−4.80	94.4	–	−12.80	94.4
Forsterite Mg2SiO4	−6.85	67.2	0.470	−10.64	79.0

Rate constant k computed from A and E, 25°C, pH = 0, mole m−2 s−1. E = Arrhenius activation energy, kj mole−1. n = reaction order n respect to H+ (except for calcite, respect to Fe3+). Calcite values (basic mechanism) are referred to the carbonate mechanism (pCO2)

Due to its first‐order dependency in equation (1), the reactive surface area is considered as a major source of variations in timing and extent of dissolution and precipitation rates. In our models, the reactive surface area depends on both the average size of basalt fragments and on the amount of primary minerals in contact with water at the rock surfaces. The relative amount of rock per volume unit, assuming that all pore space is completely filled by water, is usually characterized by the porosity (, defined by where V is the volume of solution involved in voids and V is the total volume (voids + minerals). ranges between 0 and 1 and typically falls within 0.3 to 0.7 in soil sediments [Nimmo, 2004; McSween et al., 2003]. As defined above, porosity () is simply a single‐valued estimation of the amount of space available to fluid within a specific volume of soil. Changes in porosity may be related in a relatively straightforward manner to the volumetric amount of minerals in the sample. However, the amount of mineral surface in contact with the solution also depends on the size and shape of the minerals in the system. Because different size distributions can lead to the same values of , the reactive surface area of minerals is not univocally determined by the porosity. The rate of mineral dissolution may vary considerably even for the same value of and, consequently, changes in the reactive surface area are more difficult to evaluate. Total surface areas are often determined using the adsorption isotherms method (Brunauer‐Emmett‐Teller), or using geometrical approaches to estimate the surface‐to‐volume ratios [White and Peterson, 1990; Brantley et al., 1999; Kump et al., 2000; Cubillas et al., 2005; Hodson, 2006; Zhu, 2009]. In our models, the reactive surface area of the individual minerals depends on both the size of basalt fragments as well as the textural and compositional distributions of minerals at the surface of the rock. An estimation of both values was made using the procedure reported by Sonnenthal et al. [2005] and Van Pham et al. [2012], assigning an average value of reactive surface area to the rock basalt fragments, following:

Here the specific surface of basalt (S sp (m2/g)), which is related to the total volume of rock in the sample through the porosity (), depends on the surface‐to‐volume ratio (A /V ). represents the average density of basalt estimated from its modal composition (3.1 × 106 g/m3). A /V varies with the size and shape of rock fragments in contact with water. The reactive surface area of each individual mineral in the rock was further estimated from the average S sp of basalt by taking into account the molar fraction of each mineral in the surface: where S stands for the reactive surface area of each mineral in the rock; M and n are the molar mass and moles of the mineral, respectively; S sp corresponds to the average specific surface of the basalt fragments, for the size considered; and X is an empirical factor to take into account the effective fraction of the surface that is reactive [Van Pham et al., 2012]. This latter parameter intends to minimize the discrepancies observed between field and laboratory‐based dissolution rates of silicates. These discrepancies arise due to the fact that surface reactive is a dynamic variable, which depends on multiple factors difficult to determine, including roughening, reactive site distribution (e.g., density of the kink sites), and formation of passivating layers among others [see Fischer et al., 2014, and references therein].Values of X between 10−1 and 10−3 have been used in the literature [White and Peterson, 1990; White and Brantley, 2003]. X is highly uncertain and is suggested to vary by various orders of magnitude. Here we used a reduction of one order of magnitude according with the work of Sonnenthal et al. [2005].

Rates and Mechanisms of Supersaturation: Evaporation and Cooling

The rate of generation of supersaturation (when a system attains far from equilibrium conditions) is specific for each phase and is an important factor to take into consideration for the modeling of the precipitation sequences in multicomponent systems, as is the case of solutions in contact with basaltic sediments. Different rates of supersaturation for each possible secondary phase will sequentially modify the composition of the solution and drive the sequence of precipitation to different end‐members. In our models, supersaturation is generated by evaporation and/or simultaneous cooling. In the case of evaporation, supersaturation is induced by removing the solvent at specific rates. Cooling selectively affects the supersaturation of secondary phases as a function of their solubility dependence on temperature and, upon ice formation, by removing reactive water from the system [Fairén et al., 2009]. From the point of view of a liquid‐vapor phase transition, the rate of evaporation in mass terms can be expressed as where M (mol/kgwater) stands for the moles of H2O evaporated; S/V is the surface‐to‐volume ratio (m2/dm3), which depends on the particular geometry and size of the basin; E is the specific rate constant (mol/m2 s); and the saturation term (1 − p/p 0) represents the distance from the equilibrium, where p is the atmospheric pressure and p 0 the saturation pressure of H2O(g) at that temperature [Gibbs, 1878; Becker and Doring, 1935; Appelo and Postma, 1993]. Assuming far from equilibrium conditions, this term tends to 1 and evaporation proceeds continuously. In this case, the mass amount of water evaporated depends on the S/V of the reservoir and on the H2O(g) flow across the air‐water interface, given by E .

This flow of H2O(g) has been characterized for the case of Mars by assuming that evaporation takes place under free convection conditions [Ingersoll, 1970] in a CO2 atmosphere: where (kg/m2 s) represents the Fickian flux of water vapor diffusing across the interface. The parameter 0.17 is a constant for water convection, and Gr is the Grashof number (dimensionless), which characterizes the interplay between buoyancy and viscosity during natural convection [Chevrier et al., 2007a; Sears and Chittenden, 2005; Sears and Moore, 2005]: where L is a characteristic scaling parameter of the interface (usually 1m3); g is the Martian gravity acceleration (3.75 m/s2); Δρ are the density differences between the evaporative surface in the atmosphere, normalized to the average density (dimensionless); and ν is the kinematic viscosity (m2/s).

Assuming a constant value for atmospheric density, the diffusion part represented by J in equation (6) can be linearized to where is the atmospheric density at the water interface (kg/m3); D is the diffusion coefficient of H2O in CO2 (m2/s); M CO2 and M are the molar masses of CO2 and water, respectively (g/mol); p 0( is the saturation pressure of H2O(g) (temperature dependent), which includes the effect of water activity through the dependence of water vapor pressure on the ionic strength; P atm is the atmospheric pressure; and φ (%) is the ambient relative humidity.

The temperature dependence of the specific rate constant for evaporation, under the consideration of a free convection mechanism, is shown in Figure 4a, where E is expressed in mm/h. The remaining fraction of H2O, which is the relevant parameter to evaluate the time dependence of the supersaturation, depends on temperature and pressure through the variations of E with these parameters. It also depends on an additional factor, the surface to volume ratios (S/V) of the aqueous solution at the upper boundary of the sediment column. This variation is illustrated in Figure 4b, which represents the time evolution of the mass of H2O for various S/V ratios (0.3, 0.13 and 0.03), respectively.

Figure 4

(a) Temperature and pressure dependences of the specific rate constant for evaporation in Mars. Calculations are based on free convection models. The red and blue curves show the temperature dependence of E at constant pressure conditions (0.007 atm and 1 atm, respectively). The dotted black lines show the evolution of E for processes involving a continuous decrease on temperature and with the atmospheric pressure varying at different rates (from 1 atm/yr to 0.1 atm/yr). (b) Time evolution of the remaining mass of H2O for various surface to volume ratios of the basin (S/V = 0.3, 0.13, and 0.03, respectively).

Reactive Transport (I): Diffusion Models

Our 1‐D transport models in a sediment column were done with the transport algorithm implemented in the Phreeqc software [Parkhurst and Appelo, 1999]. Vertical diffusion through n cells (represented by n solutions of pore water in the sediment column) was based on part of the general diagenetic equation [Berner, 1980], which includes advection, diffusion, and reaction components to model ionic and water transport into sediments: where C is the concentration of each species in solution (mol/kgwater); F is the flux through the sediment column (mol m−2 s−1); and R is a reaction term that includes the rate at which a component is produced or consumed by different reaction mechanisms (i.e., dissolution, precipitation, and adsorption). The flux term has two parts: first, the advective transport term F adv = vC, which results from the flow of water (v), and second, the diffusive flux, which is governed by Fick's first law: where D (m2 s−1) is the diffusion coefficient of each species in dissolution. Thus, we get the general advection‐diffusion‐reaction equation employed in the calculations:

In the case of basaltic sediments of small size, water can flow through the rock fragments, and the amount of solid phase in the system is greater than the space filled by water. If porosity is such that avoids the flow of water, only the transport of ions occurs (when a concentration gradient exist), in accordance with Fick's laws of diffusion (equation (10)).

A simple criterion to evaluate whether advection is important, compared with diffusion, is the dimensionless Peclet number [Lerman, 1975]: where D is the molecular diffusion coefficient; L is the average distance; and U is the flow velocity of pore water. Diffusion is the dominant process when D > > LU. In our diffusion models, we used porosity values of 0.3 and 0.4, corresponding to weathered basalt or basaltic sands [McWhorter and Sunada, 1977; Freeze and Cherry, 1979], except in Figure 7, where we evaluate the influence of porosity in the transport process. For these porosity values, the advective flow velocities (U) are in the order of 0.001 to 1 cm/yr. And for a value of D  = 100 cm2/yr, advection is only important when migration distances (L) are greater than 1 m [Berner, 1980]. Consequently, diffusion can be considered as the major transport process at the sediment‐water interface during early diagenesis in terrestrial marine sediments [McSween et al., 2003].

Our models allow consideration for both single average diffusion coefficient and multicomponent ionic diffusion transport. In the latter case, we used the values of tracer diffusion coefficients (D ) given in the Phreeqc database. An average value of D  = 10−9 m2 s−1 was employed as default when, for a particular species, D was lacking in the database. The effective diffusion coefficient (D ) in porous media, such as weathered basalts, is estimated from the water diffusion coefficients (D ) and porosity (ϕ), using Archie's law: where n, the cementation exponent, varies from 1.3 for diffusive gas to 5.4 for clays [Navarre‐Sitchler et al., 2009].

Reactive Transport (II): Adsorption Models

In addition to porosity, further chemical aspects like ion exchange and adsorption can greatly influence the transport of ions through the sediment column. Charged species may interact with sediment, in particular with clay minerals. Adsorption onto clay surfaces may occur through the formation of a diffuse layer interface, in which ions are weakly bonded by electrical forces or by chemisorbed complexes [Stumm, 1997]. The way in which surface properties of clays affect ionic diffusion can be represented as a chemical reaction by the adsorption of a mobile ion, A, onto a clay particle [McSween et al., 2003].

To evaluate the effect of adsorption on clay surfaces, we analyzed the reactive surface area in smectites (2:1) able to adsorb ions in its diffuse double layer (DDL). The DDL occurs at the interface between the clay surface and the solution, and it is made up of the permanent negative charge of the clay and the exchangeable cations in the solution that balance the negative charge [Mojid and Cho, 2006]. The net value of the adsorption depends on each cation, on the quantity of surface in weight, and on the number of reactive centers.

The Effect of the Reactive Surface Area on the Precipitation of Aqueous Minerals in an Evaporating Basin

Crystallization and dissolution are far‐from‐equilibrium processes, and the measure of the “distance from the equilibrium” for each mineral is characterized by its saturation ratio. In practice, supersaturation plays the role of a “switch,” allowing or preventing the precipitation of secondary phases on both equilibrium and kinetic modeling. However, in an equilibrium approach, the mass balance is instantaneously restored, and therefore, the time evolution of supersaturation and its influence over the thermodynamic conditions of the system cannot be accounted for. This is a key aspect when we consider a kinetic approach to crystallization. The rate of generation of supersaturation depends on the evolution of the particular parameters inducing far‐from‐equilibrium conditions (i.e., evaporation and/or cooling). In multicomponent solutions, different rates of evolution toward disequilibrium situations can drive the system to completely different end‐points when returning to near equilibrium conditions, depending on the transient values of thermodynamic parameters (pH or pe) and the different crystallization sequences.

We analyze this point by defining a geochemical scenario at the sediment‐water interface, which would mimic the benthic boundary layer in terrestrial oceans. On Earth, two environments are usually defined [Berner, 1980; Boudreau, 1997]. The first is the nondiffusive boundary layer (NDBL), where the mass of water relative to rock (W/R) and the porosity are significantly high (ϕ > 0.7), and where convection is the dominant mechanism of transport. And the second environment is the interstitial water. These two environments may be separated either by an impermeable interface (i.e., pristine basalt or a hardground) with very low porosity, in which case no downward transport is expected or, alternatively, by a gradual decrease of porosity in the sediment column, with values between ϕ = 0.3–0.7, as typically occurs during the early diagenesis of marine sediments. In this latter case, diffusion becomes the dominant mechanism of ionic transport, leading to the definition of the diffusive boundary layer (see Figure S3).

In the NDBL environment, the dissolution rates of primary minerals are mainly controlled by the size and shape of basalt fragments through its effect on the reactive surface area (equation (1)) (additional factors would be the physical properties of basaltic fragments/sediments). Its influence on the evolution of the system was tested by using the same interaction with volatile phases in the atmosphere and the same mass amount of basaltic phases in the mineral assemblage. Thus, the only variable parameters in the models are the generation rate of supersaturation, which is controlled by the rates of evaporation and/or cooling, and the size of basalt fragments. In our models, clast‐size fractions ranged from very fine silt to medium‐sized boulders (Phi 8 to −8, using the classification by Blott and Pye [2012]).

The coupled role of reactive surface area and the evolution of supersaturation, driven by evaporation and cooling, on the pH and the mineral sequence is tested in Figures 2 and 3. In this case, we used a porosity value of ϕ = 0.9 and we did not consider transport. Different rates of evaporation were induced by assuming free convection, according with equation (6), varying the time schedules for the evolution of temperature and/or pressure. Figure 2a shows the system evolution toward alkaline and reducing endpoint conditions when small fragments of basalt (cubes: L = 12 μm; reactive surface area = 480 m2/kgwater) are allowed to react with water, under a low generation rate of supersaturation (compared with Figure 3). Quick basalt dissolution buffers the system with the alkalinity derived from silicates (according to the reaction scheme: proton + silicate mineral ➔ cations + silicic acids), which offsets the acidification tendency imposed by the atmospheric gases [Ghiara et al., 1993; Fairén et al., 2004; Zolotov and Mironenko, 2007; Kato and Yamasato, 2013]. Despite clays incorporate OH− when they precipitate, this process is not fast enough to balance the release of Si by basalt dissolution. In the end, the solution remains alkaline and reduced, and only clays and iron oxides precipitate (Figure 2b).

The opposite situation arises when larger basalt fragments (cubes L = 11 cm) are considered (therefore, the total reactive surface area is much smaller than in Figure 2a: reactive surface area = 0.051 m2/kgwater) and when this condition is coupled with a faster generation rate of supersaturation driven by evaporation, cooling, and reaction with volatiles (Figure 3a). In this case, precipitation of clays and iron hydroxides takes place at the initial stages of evaporation, when pH remains at circumneutral values and the temperature is still relatively high. With progressing evaporation, the system evolves toward acidic conditions, finally leading to the precipitation of evaporites (Figure 3b).

Consequently, the extent of the reactive surface area and the rate of water loss may induce the evolution of the system toward alkaline (high reactive surface area) or acidic values (low reactive surface area), depending on the interplay between the relative rates of acidification (which, in turn, are promoted by atmospheric volatiles) and the OH− delivery (induced by the dissolution of silicate minerals through the formation of silicic acids and often known as “silica alkalinity”).

As a corollary to this analysis, it is important to note that, when the reactive surface area of basalt is large and there is a considerable amount of sediment, the porosity existing along the sediment column triggers an additional process of ionic mass transport. This takes place by diffusion and/or advection, also affecting the kinetic evolution the system. In the following section, we present preliminary models including the reactive transport across basaltic minerals.

Water‐Sediment Diffusion Models

When the basaltic substrate is fractured, the detrital basaltic material shows a degree of porosity that promotes ionic vertical transport through the sediment column. In this scenario, diffusion is a straightforward mechanism of vertical transport into the sediments for the cations remaining in solution after the precipitation of clays. We present here a wide range of possible diffusion processes, with the goal of understanding the interplay between three factors involved in the precipitation path of secondary phases in aqueous environments: (i) the dissolution rate of primary minerals (driven by their reactive surface areas), (ii) the ionic transport by diffusion, and (iii) the reaction with gaseous species.

To track the particular effect of the evaporation rate on the effectiveness of removal of ions contained in surface solutions, we developed a simplified model allowing the downward migration through a gel column, using Ca2+ as tracer. Surface solution is constrained to evaporate at different rates during 500 years (Figure 5). Our results indicate that, under a low evaporation rate (Figure 5a), the concentration gradient disappeared with time because the diffusion velocity is able to decrease the concentration of Ca2+ from the surface layer. Considering an intermediate evaporation rate (Figure 5b), Ca2+ is initially transported to the subsurface, but after 100 years of simulation, a stationary concentration gradient is created. Finally, with a high evaporation rate (Figure 5c), the diffusion of Ca2+ is not enough to remove the supersaturation generated by the rapidly loss of water from the surface layer.

Figure 5

Subsurface migration of Ca2+ under three different evaporation rates (in this case, we use a sequence of linear rates to generate different supersaturation schedules). (a) Slow rate: 0.02 × 10−3 moles/yr. (b) Intermediate rate: 1.8 × 10−3 moles/yr. (c) Fast rate: 40 × 10−3 moles/yr. In the simulations we evaluated the downward migration of Ca2+ through a gel (D  = 1 × 10−9 and ϕ = 0.3) by using an aqueous solution of Ca2+ and SO4 2−.Therefore, the evaporation rate was the only variable in the models that changes the effectiveness of the Ca2+ diffusion. As initial conditions, we started with a concentration gradient of Ca2+ (defined by 1 mM in the upper cell and 0.5 mM in the rest of cells) and a sediment column of 2 m.

The time evolution of the concentration of Ca2+ in the surface layer is shown in Figure 6. The highest rate of cation removal occurs at the slowest evaporation rate. Despite loss of water is the driving force to induce the ionic‐diffusion transport, higher evaporation rates may also induce a level of supersaturation that ionic migration (limited by the effective diffusion coefficient, De) cannot compensate, leading to a lower rate of cation removal. Thus, the rate of cation removal from the surface solution is strongly dependent on the competition between evaporation and diffusion transport: although evaporation increases the ion concentration, diffusion plays an opposite role, decreasing the ion concentration in the surface solution. Consequently, the precipitation of evaporites could be inhibited in a continuous but slow evaporative scenario.

Figure 6

Evolution of the surface concentration of Ca2+ with time (upper cell of the previous models).

As the effective diffusion coefficient is directly proportional to the porosity (equation (13)), it is expected that the main effect of increasing porosity would be a consequential increase of the effectiveness of the downward ionic transport through the sediment column (Figure 7). As is shown in the figure, higher porosities favor faster cation transport into the sediment column, further preventing the precipitation of evaporites near the surface.

Figure 7

Downward diffusion of Ca2+ at different values of porosity, ϕ: (a) 0.1, (b) 0.3, and (c) 0.7. We used the same model approaches described in Figure 5 (evaporation rate = 0.2 × 10−3 mol/yr).

Diffusion Under Instantaneous Precipitation of Secondary Minerals

We use this model to evaluate the net effect of diffusion as a transport mechanism and its ability to remove ions from the water‐sediment interface even in the more adverse scenario (i.e., when precipitation of solid phases occurs instantly). In these simulations, we assume a high reactive surface area of basalt, a slow evaporation rate, and a continuous flux of volatiles from the atmosphere.

Results of this model show that diffusion is an effective mechanism to remove cations from the surface even in short time periods (1 year). The concentration gradient, induced by loss of water in the surface layer, drives the excess of cations (after clay formation) into the sediment column (Figures 8a–8c). Sulfates do not precipitate in this scenario due to the alkalinity of the system (resulting from the high reactive surface area of the basalt). Additionally, the equilibrium model predicts that iron oxyhydroxides would not precipitate together with clays because the content of Fe3+ in solution is very low and falls well between the concentration ranges where hydrolyzed species of Fe+3 are stable in solution. These low Fe+3 concentrations result from the interplay between (i) diffusion that reduces the evapoconcentration of the Fe3+ species in surface waters and (ii) the incorporation of iron into clay mineral lattices, which are very insoluble at neutral pH values (Figure 8d).

Figure 8

Diffusion of cations through a basalt sediment column 20 cm length, during 1 year, and assuming equilibrium conditions for secondary mineral precipitation (D 1 × 10−9; ϕ = 0.3): (a) Ca2+, (b) Mg2+, and (c) Fe3+. (d) Time profile of nontronite‐Mg is shown as an example of the secondary mineral phase formed in the simulation. Model approaches: linear evaporation rate ~ 0.025 mol/d; continue flux of volatiles from the atmosphere (CO2(g) ~ 0.055 Mm/d, Cl2(g) and SO2(g) ~ 0.0237 Mm/d). Dissolution of primary minerals was modeled kinetically, considering an average reactive surface area of 1096 m2/kgw (cubes, L = 5.5 μm). Secondary minerals considered in the model: beidellite, kaolinite, nontronite‐Mg, smectite, clinochlore, gibbsite, goethite, alunite, gypsum, and jarosite. Kaolinite and gibbsite are also formed in the simulation (see Figure 9).

Figure 9

Time evolution of the mineral sequence obtained with equilibrium models: (a) with no diffusion transport and (b) allowing diffusion transport (upper cell of the previous model, Figure 8). The boundary conditions of the model are those described in Figure 8.

To test our results, we analyze our equilibrium model in absence of diffusion, by using the same initial approaches but avoiding the ion transport to subsurface layers. Figure 9 compares the surface mineral assemblages obtained with both models. When vertical transport is not considered, anions derived from volatiles progressively accumulate in solution, and the pH of the system shifts toward more acidic values. This, together with an increase in cation concentration due to evaporation, favors the supersaturation of evaporites (gypsum) and iron oxyhydroxides (goethite). Taken together, our equilibrium models show that diffusion is a qualifying mechanism to inhibit the coetaneous precipitation of clays, evaporites, and iron oxides near the surface.

Results of the Kinetic Model

After testing the effectiveness of surface cations removal by diffusion processes, we have modeled the space‐time variations of the conditions inside a basalt sediment column, where the surface water is in contact with an acidic‐oxidative atmosphere at the upper interface (CO2, SO2, and Cl2) (as on early Mars, see Fairén et al. [2004]) and subject to a continuous loss of water by evaporation. We assume a continuous flux of gases derived from volcanic outgassing through the upper water layer, as these volatiles are the most obvious source of anions that can gradually change Eh and pH conditions. The model represents a 1‐D kinetic system of evaporation + diffusion + reaction (precipitation/dissolution of mineral phases, see Table 2).

Oxidation is mainly driven by chlorine, through transient oxygen reactive species formed at the low pH values induced by the CO2 and SO2 gases [e.g., Fairén et al., 2004; Zelenski and Taran, 2012; Hanley et al., 2012; Levanov et al., 2015]. When the upper acidic and oxidant solution starts to diffuse and mix with the interstitial alkaline water, the pH rapidly decreases with depth. The extent of pH lowering toward acidic values and increase in oxidation potential, on both time and depth, is a function of the rate of diffusion and the buffering capacity of the interstitial pore water, which, in turn, depends again on the dissolution rate of the basalt forming minerals. Therefore, this model considers once again the two different possibilities regarding the reactive surface area of the basaltic minerals, as follows.

First, when the solution is in contact with basaltic sediments of coarser grain size, the total amount of reactive surface area is low and, therefore, the dissolution of the basalt proceeds relatively slowly. In this case, it is expected that the input of volatiles from the atmosphere and the loss of water from evaporation are the main factors inducing variations in pe‐pH conditions and in ionic strength. In this situation the system evolves toward an evaporitic scenario, characterized by extremely acidic and oxidizing conditions. The consequence is the formation of sulfates (i.e., alunite and gypsum), more abundantly in the upper part of the sediment column, while the distribution of clays tends to follow an opposite trend in the column (Figure 10). Note that when our models take the system toward very acidic and oxidizing conditions, our results would predict the formation of thermodynamically related jarosite phases; and, even in the absence of K necessary to form K‐jarosite, other forms of jarosite such as ammonium‐jarosite, natrium‐jarosiste, or even hydrogen‐jarosite would still form under similar conditions.

Figure 10

Depth profiles of the secondary mineral phases formed during aqueous reaction of a basalt column considering gases input from the atmosphere and low reactive surface area of primary minerals (cubes, L = 54 mm; S reactive = 0.11 m2/kgw). Beidellite‐Ca and gibbsite also precipitated. Both phases showed the maximum accumulation at the subsurface (see Figure 12 below). Model conditions: flux of gas input: 0.1 mol/yr of CO2(g) and 0.01 mol/yr of SO2(g) and Cl2(g); length of the sediment column = 20 cm; D  = 1 × 10−9, ϕ = 0.4; simulation time = 100 years; linear evaporation rate = 0.1 mol/yr. Secondary minerals considered in the model: clays (smectite, kaolinite, beidellite, and chlorite), evaporites (jarosite, alunite, gypsum, and calcite), and hydroxides (goethite and gibbsite).

Figure 12

Distribution of the mineral assemblage obtained from the previous models (Figures 10 and 11) after 100 years of simulation: (a) with high reactive surface area and (b) with low reactive surface area. As pH decreases, Fe3+ precipitation leads to the formation of Fe‐clays (nontronite type). When Fe is no longer available and pH continues decreasing, the formation of Al‐clays (beidellite type) is favored.

And second, when the reactive surface area of the basalt increases, the dissolution of basalt occurs faster, resulting in the rapid evolution of the solution toward very alkaline values. In this scenario, sulfates do not form, and clays are the major precipitation product in the surface layer (Figures 11a–11c). The alkaline pH also favors the precipitation of goethite. Goethite's precipitation front is gradually shifted toward the middle point of the column, due to the competition for ferric iron with clay minerals (i.e., nontronite) (Figure 11d).

Figure 11

Depth profiles of the secondary phases formed during aqueous reaction of a basalt column under the same boundary conditions described in Figure 10 but assuming high reactive surface area of primary minerals (cubes, L = 5.5 μm; S reactive = 1096 m2/kgw). Calcite and gibbsite are also formed, showing the maximum accumulation at surface and subsurface layers, respectively.

The comparison of the mineral assemblages obtained in both scenarios after 100 years of simulation is shown in Figure 12. Note that our results are also consistent with the stratigraphic sequence of Noachian phyllosilicates observed in certain locations on Mars: Fe/Mg clays (e.g., nontronite) in the lower layers and Al‐rich clays (e.g., beidellite) in the upper layers [McKeown et al., 2009]. Taken together, these results highlight the importance of considering open kinetic models to reproduce the geochemical conditions on the surface of Mars. As we have shown, cation removal by diffusion, together with the input of volatiles and the influence of the reactive surface area of the basalt, could play a central role in the time evolution of the mineral sequences formed on aqueous environments on early Mars.

Adsorption of Cations to Clays

To make our model more complete, we have also included in our simulations the complementary reservoir of cations that must exist attached to the phyllosilicates, derived from the selective adsorption of cations to the negatively charged surface of clays (smectites (2:1)). Adsorption of cations onto clay surfaces is important when attempting to balance the cation budget in dissolutions and, in particular, is an efficient way to eliminate contaminants [Cama et al., 2005]. During reactive transport, adsorption can act as cooperative mechanism contributing to scavenging ions from solution. Unlike the cation exchange process, which usually does not involve a net decrease of the ions in solution and only affects the rate of ion displacement along the sediment column, the adsorption process via surface complexation [Dzombak and Morel, 1990] is an efficient mechanism to fix ions to specific sites of mineral phases and consequently is an important factor in lowering the ionic concentration.

The effectiveness of the adsorption process to remove cations from the solution, together with the diffusion process, is shown in Figure 13. Adsorption of cations proves to be a very efficient mechanism to immobilize ions, not involving the incorporation of the cations into the mineral lattices. Consequently, these adsorbed cations would be undetectable with the usual orbital spectroscopic techniques, which rely on the recognition and identification of specific functional groups using diagnostic molecular vibrational features: the crystal structure varies from mineral to mineral according to its functional groups and the type of ligands formed, producing different absorptions and making specific mineral identification possible from spectroscopy [Clark, 1999]. The removal of cations from the surface by adsorption processes only was evaluated in Figure 14, where a significant reduction of the surface concentration of dissolved Na+ is observed as the amount of clay increases.

Figure 13

Deep profiles of dissolved cations during 20 days considering adsorption + diffusion processes: (a) Na+, (b) Ca2+, (c) Fe2+, and (d) Fe3+. At the starting point (t = 0) there is a concentration gradient created by a rapid loss of water at the surface (−30.5 mol of H2O). We assume 0.1 surface active sites/g clay and multidiffusion.

Figure 14

Time evolution of Na + concentration at the surface layer considering adsorption onto clay surfaces during 200 days. Four different scenarios for the reactive clay surface are shown, from no ionic absorption (red line) to progressively more surface of negatively charged clays. The lowering of ionic concentration is faster as the quantity of clays increases.

Finally, we analyzed the evolution of cations and anions at the upper layer of a system subject to evaporation, diffusion, and adsorption processes (Figure 15). Cation concentration substantially decreases with time, being more significant for monovalent cations, which are more easily adsorbed onto clays (Figure 15a). In contrast, the concentration of anions increased with time (Figure 15b). In the case of anions, molality increases because the adsorption onto clays is not an efficient process with negative ligands. Therefore, because the loss of water is faster than the diffusion rate of anions, their molality increases. In the end, the net content of anions in the surface solution is progressively lower, but this decrease in concentration is slower than in the case of cations, as the anions are eliminated only by diffusion, with no adsorption to clay minerals.

Figure 15

Evolution of (a) cations and (b) anions molality on the surface during 200 days in a system subjected to evaporation/sublimation. The reactive clay surface has a constant value (1 gr clay each 0.1 mol of surface‐active sites).

Discussion

Authigenic Versus Detrital Sediments in the Modeling of Early Diagenetic Processes

The main goal of our paper was to characterize the formation of different secondary minerals on Mars using geochemical modeling. These minerals would form authigenic sediments (sediments of chemical origin formed “in situ”), not detrital sediments (sediments formed as result of a process of transport and deposition). The only transport mechanisms considered in our models are “molecular transport mechanisms,” which trigger the transport of ions by modifying the supersaturation state of each mineral phase at depth (first Fick's law) and through time (second Fick's law). The word “transport,” as employed here, must not be confused with the macroscopic transport of detrital sediments carried out by the agents considered in basic sedimentology (i.e., rivers, wind, tides, and waves); these, by cessation of mechanical energy plus deposition, led to the formation of detrital sediments with mineral phases that could very well be out of geochemical equilibrium. This would be a different process entirely. Geochemical modeling is an adequate approach to understand the formation of authigenic phases formed by chemical reactions [Berner, 1971, 1980; Boudreau, 1997; Boudreau and Jørgensen, 2001], which represent a fraction (the authigenic fraction) of the whole sedimentary record on Mars. Moreover, these authigenic sediments result from the combination of molecular transport mechanisms (like diffusion and advection), together with chemical reactions. Combined, these processes determine the sequential formation of new mineral phases with depth in the sedimentary column.

We used a sedimentary column because fragments of basalt, despite their igneous origin, generate sediments via transport and deposition and, upon burial, sedimentary columns. The process of chemical alteration involves basalt fragments forming clasts that interact with water (e.g., at the bottom of a lake) and form new authigenic minerals that are sediments [e.g., Schwenzer et al., 2012]. We modeled how pore spaces in the sedimentary column become occupied by authigenic minerals. Consequently, the whole set of sedimentary processes resulting from our models are better understood under the global term of “early diagenetic processes.” Although the term “early diagenesis” and the chemical sedimentary processes associated with it have only recently been applied to sedimentary deposits Mars, some fresh papers do describe the sedimentary context at Gale crater in terms of a lake where diagenetic process took place at the water‐sediment interface, inducing stratification (meaning chemical stratification) of redox gradients and authigenic phases by sequential precipitation [Hurowitz et al., 2017; Rampe et al., 2017; Frydenvang et al., 2017].

Our models showed how, starting with identical primary mineralogy, it is possible to obtain different “suites” of secondary minerals (authigenic sediments) on Mars, depending on the kinetic evolution of the system, driven by two main factors: (1) the rate of generation of supersaturation, controlled by evaporation, cooling, and the mixing with volatiles and (2) the rates of mineral dissolution/precipitation, controlled by the reactive surface of the primary rocks. The primary rocks are either pristine basalts (lavas) or fragmented basalts forming detrital sediments.

Our simulations were strongly dependent on time and on the initial concentration gradient. By diverse combinations of subsurface diffusion versus surface evaporation, and varying the total time evolution of the system, it was possible to modify the spatial relationships between phyllosilicates and evaporites in the uppermost centimeters of the column. Therefore, when substrate porosity enabling diffusion in an open system was considered, the distribution of the mineral phases in each level did not have to mimic the composition of the initial basalt. Diffusion, when taking place in a diagenetic environment, is usually considered as occurring in a dynamic framework with continuous burial and/or denudation of sediments [Berner, 1980]. In this context, our models showed the diffusion behavior in the upper layer, but without considering the deposition of new sediments. This behavior may become stationary over time, in what concerns the migration of cations, but the particular sediments involved could change if processes of sedimentation and erosion would occur (see e.g., in situ determinations of sedimentary rocks exposure by wind‐driven scarp retreat at Gale crater in Farley et al. [2014]).

Setting Our Results in the Context of Previous Models

The mineral assemblages resulting from our models are consistent with results from previous geochemical modeling efforts, studies of laboratory analogs and meteorites, and in situ analyses by rovers. For example, the formation of greenalite (together with other kaolinite‐serpentine, i.e., antigorite, berthierine, and odinite) has been proposed in previous geochemical models of Mars [Chevrier et al., 2007b; Zolotov and Mironenko, 2007]. And although several authors have pointed out how the presence of smectites is incompatible with sulfates from a geochemical point of view [i.e., Chevrier et al., 2007b], our model approach allowed us to explain this apparent contradiction, as well as the formation of mixed facies.

The novelty in our models presented here was the simultaneous use of fully open conditions and a kinetic approach. Previous models, for example, Melwani Daswani et al. [2016], also implemented a 1‐D transport model to understand carbonate formation, or the lack thereof, concluding that evaporation plays a role; however, this work did not consider open conditions at the water‐sediment interface and used only equilibrium (i.e., nonkinetic) conditions. Another example of models that did not consider open conditions at the water‐sediment interface would be Schwenzer et al. [2016], which dealt with the formation of the sulfate veins crosscutting Gale sediments; in addition, our models are not restricted to Gale, as our results are applicable to different Mars localities.

Our work, which considers relatively diluted solutions, is also a departure from previous studies that focused on evaporite mineralogy and brine evolution [e.g., Tosca et al., 2005; Fairén et al., 2009]. In the same vein, some previous investigations of Mars evaporitic precipitation using laboratory analogues had different starting and boundary conditions; for example, Moore et al. [2010] did not include Si in the list of elements forming part of the experimental solutions because their work did not address the weathering of silicates, as their experimental study was exclusively designed to understand the sequence of brine precipitation (i.e., halides, sulfates, and carbonates). Our approach required consideration of the formation/dissolution of silicates to understand the dissolution of basalt. In addition, the time intervals employed in our calculations were not comparable to those used in the laboratory by Moore et al. [2010]. In our analysis, the process of generation of supersaturation is analyzed from the very beginning, starting from just pure water and basalt.

Applicability of Our Models to Mars

Our numerical models allowed us to recreate two different geochemical scenarios relevant to understanding the geological history of Mars. On one side, in a scenario of basaltic sediments of small grain size with a large reactive surface area, silicates were dissolved faster. The geochemical evolution of the system would have been determined by the downward diffusion of cations, which, in this scenario, was a faster process overcoming the progression of evaporation, therefore promoting the maintenance of circum‐neutral pH values and favoring the synthesis of phyllosilicates. Large regions of basaltic sediments of smaller size are expected to have been widespread early in the history of Mars (until the Late Noachian to Early Hesperian), when a more geologically active surface, due to volcanism, impacts, valley formation, weathering, and erosion [e.g., Carr and Head, 2010], would have enhanced mechanical weathering and physical erosion of basalt, producing the clasts with different sizes required to form a “detrital sediment.” This detrital sediment was the ultimate source of ions required to form authigenic sediments, through dissolution precipitation/processes. On the other side, in a scenario of basaltic sediments of greater size with less reactive surface area, the geochemical evolution of the system would have been controlled by the rates of water loss via evaporation, faster than the cation diffusion in this case. This would have promoted the evolution of the system toward acidic conditions, and the slow dissolution of basalt would have resulted in the precipitation of evaporites. This scenario would have been more typical of post‐Noachian periods (i.e., Hesperian to Early Amazonian), when a high average rate of still active volcanism resulted in the formation of new large draping deposits of basaltic materials, while the rates of mechanical weathering and physical/chemical erosion of basalt were already significantly lower [e.g., Carr and Head, 2010].

Our predictive model results are well in line with the reported difference in fracturing of the Martian protolith depending if it was older and altered or young and pristine, as measured by the orbital surface‐penetrating radar Shallow Radar (SHARAD) on board Mars Reconnaissance Orbiter [Seu et al., 2004]. Of course, fractures of the scale resolvable by SHARAD radargrams would not have any effect on reaction rates, but they do provide vital information about the pervasiveness of mechanical weathering on different volcanic sediments and the generation of regolith by physical erosion, which would create a loose surface. In this sense, SHARAD data indicate that young and pristine volcanics must have homogeneous interiors, while the older volcanics become progressively fractured [Stillman and Grimm, 2011]. As such, SHARAD data confirmed that the reactive surface area of Noachian basalts was larger than that of post‐Noachian basalts, therefore offering observation‐based evidence for our geochemical modeling of the different aqueous modification of the volcanic basalts at different time epochs on Mars, due to their different degree of fracturing.

Our results offer an additional, complementary, and not mutually exclusive explanation to that by Bibring et al. [2006] for the transition from a Noachian aqueous mineralogy involving substantial quantities of circum‐neutral liquid water and dominated by the synthesis of clays, to a post‐Noachian geochemistry characterized by the precipitation of evaporites in limited amounts of saltier and acidic liquid water. Moreover, our model also allows for any kind of intermediate state depending on the degree of fracturing of the Martian protolith. Intermediate states would aptly explain the geochemical transition in time observed in Martian mineral sediments. They would also explain particular mixtures of clays and evaporites observed in specific localities: phyllosilicates and evaporites appear sometimes within meters or even in spatial juxtaposition of each other. This has been shown by the MER Spirit in the Columbia Hills at Gusev crater [Wang et al., 2006], by the MER Opportunity on the rim of Endeavour crater at Meridiani Planum [Arvidson et al., 2014], by the MSL Curiosity in Yellowknife Bay at Gale crater [Vaniman et al., 2014], and by MEX‐OMEGA and MRO‐CRISM in several places on Mars, including regional coexposures of phyllosilicates with either carbonates [Ehlmann et al., 2008, 2011; Carter and Poulet, 2012], sulfates [Poulet et al., 2005; Wiseman et al., 2008; Hamilton et al., 2008; Murchie et al., 2009; Wray et al., 2009; Milliken et al., 2014], or chlorides [Murchie et al., 2009] in different locations, and even the presence of sulfate‐bearing materials underlying phyllosilicate‐bearing strata [Wray et al., 2010].

The kilometer‐thick sedimentary sequence observed in the central mound of Gale crater, which beautifully displays the transition from ancient clay sediments to more recent sulfate and oxide assemblages [Milliken et al., 2010, 2014], may serve as a test site for our hypothesis presented here: the clay sediments would have been deposited in ancient lakes covering the crater floor [Grotzinger et al., 2015]; time after, new interactions between water and the central mound of the crater [Fairén et al., 2014] would have resulted in the precipitation of sulfates. This is consistent with the observed transition in sulfate occurrences at Gale: from postdepositional Ca‐sulfate veins ubiquitous in almost all of the clay‐rich bedrock in the crater floor up to the Murray Buttes [Nachon et al., 2014], to abundant Ca‐sulfate cemented sandstones at the Murray Buttes [Newsom et al., 2017], suggesting a possible change in the original mode of deposition that may be interpreted within the framework of the processes modeled in this paper. Indeed, subsequent interactions between water and the central mound could have remobilized/reactivated older deposits. Upcoming data from the Curiosity rover will be useful to test our model results, as the rover exploration has already provided direct evidence for different lakes coming and going over Gale crater during up to tens of millions of years [Grotzinger et al., 2015; Hurowitz et al., 2017].

Conclusions

We have highlighted the contrasting weathering behavior of basalt sediment fragments of different sizes and its influence on the secondary mineral products that may have resulted from aqueous interaction on early Mars. Although the variation in the total reactive surface area of basaltic clasts does not change the pH dependencies for the dissolution of the individual minerals, the release ratios of different cations may be significantly modified. In basaltic sediments of small size, mineral dissolution is counter‐balanced by the precipitation of clays. The simultaneous downward diffusion of cations through the sediment column, together with the selective adsorption of cations to the negatively charged surface of clays, limits the increase in concentration at the water‐sediment interface. This results in preventing the precipitation of salts, even though evaporation is constant in the surface waters in contact with the atmosphere. In contrast, in basaltic sediments of greater size, although mineral dissolution proceeds more slowly, evaporation and long‐term interaction with volatiles (provided by volcanic outgassing) increases the concentration of cations in solution by removing water; this turns the solution toward more acidic and oxidizing conditions. In this case, the absence of ionic diffusion toward the subsurface allows the increase of concentration in the water‐sediment interface, leading to the precipitation of evaporites in the later stages of the process. Our models describing the subsurface diffusion of cations may help opening a new venue in the analysis of the geochemical processes generating supersaturation and precipitation in aqueous solutions on early Mars, beyond surface evaporation. Our results, implying no widespread formation of evaporites together with clays during the Noachian, contribute to explain why no global exposed subsurface reservoir of evaporite minerals in cross sections in craters, canyon walls, or fractures has been detected on Mars by orbital means. Taken together, these results highlight the importance of considering open kinetic models to reproduce the geochemical conditions on the surface of early Mars. As we have shown, cation removal by diffusion, along with the input of atmospheric volatiles and the influence of the reactive surface area of the basalt, could play a central role in the time evolution of the mineral sequences formed in aqueous environments on early Mars.

Supporting Information S1

Click here for additional data file.