Anhydrite‐Assisted Hydrothermal Metal Transport to the Ocean Floor—Insights From Thermo‐Hydro‐Chemical Modeling

High‐temperature hydrothermal venting has been discovered on all modern mid‐ocean ridges at all spreading rates. Although significant strides have been made in understanding the underlying processes that shape such systems, several first‐order discrepancies between model predictions and observations remain. One key paradox is that numerical experiments consistently show entrainment of cold ambient seawater in shallow high permeability ocean crust causing a temperature drop that is difficult to reconcile with high vent temperatures. We investigate this conundrum using a thermo‐hydro‐chemical model that couples hydrothermal fluid flow with anhydrite‐ and pyrite‐forming reactions in the shallow subseafloor. The models show that precipitation of anhydrite in warming seawater and in cooling hydrothermal fluids during mixing results in the formation of a chimney‐like subseafloor structure around the upwelling, high‐temperature plume. The establishment of such anhydrite‐sealed zones reduces mixing between the hydrothermal fluid and seawater and results in an increase in vent temperature. Pyrite subsequently precipitates close to the seafloor within the anhydrite chimney. Although anhydrite thus formed may be dissolved when colder seawater circulates through the crust away from the spreading axis, the inside pyrite walls would be preserved as veins in present‐day metal deposits, thereby preserving the history of hydrothermal circulation through shallow oceanic crust.


Introduction
Black smokers on the modern seafloor are the spectacular manifestations of deep hydrothermal processes that vent hot metal-rich fluids into the ocean. The chimney structures at the seafloor are formed as a result of focused fluid flow and mineral precipitation processes. They are primarily composed of sulfide minerals (e.g., pyrite, pyrrhotite, chalcopyrite, and sphalerite), sulfates (anhydrite and barite) and amorphous silica (Haymon, 1983;Tivey, 2007). Many vent systems are associated with polymetallic massive sulfide deposits at and below the seafloor (e.g., Hannington et al., 2005). Around these chimneys, unique chemosynthetic ecosystems exist in the otherwise bleak environment of the deep sea.
Much has been learned about hydrothermal circulation systems from surveys of mid-ocean ridge segments and direct observations of venting at the seafloor (e.g., deMartin et al., 2007;Fornari et al., 2012;Hannington et al., 2005;Tao et al., 2012). However, the processes in the deep reservoir that control the chemistry and physics of the hydrothermal fluids remain largely inaccessible to direct sampling and observations. Here, theoretical work can help to link seafloor observations to the physicochemical regime at depth. Recent advances in hydrothermal flow modeling have revealed the key thermodynamic and fluid-dynamic controls on hydrothermal convection and vent temperatures at oceanic spreading centers. The observed upper limit to black smoker vent temperatures of approx. 400°C has been explained by the thermodynamic properties of water, and high-resolution 3-D flow models have unraveled complex hydrothermal convection patterns at ocean spreading centers (Coumou et al., 2008;Fontaine et al., 2014;Hasenclever et al., 2014;Jupp & Schultz, 2000;Weis et al., 2014). While these insights provide a robust theoretical basis for hydrothermal flow observations at mid-ocean ridges, a number of first order discrepancies between model predictions and observations remain.
One key parameter controlling hydrothermal flow pattern and, thereby, hydrothermal mass and energy fluxes is permeability. Permeability versus depth profiles can be inferred from laboratory and in situ measurements in bore holes and on drill cores (Becker & Fisher, 2000;Fisher, 1998), seismic data (Carlson, 2014;Ingebritsen & Manning, 2010), hydrothermal heat fluxes (Lowell & Germanovich, 1994;Wilcock & McNabb, 1996), and the poro-elastic response to tidal phase shifts (Barreyre et al., 2014(Barreyre et al., , 2018Barreyre & Sohn, 2016;Crone et al., 2011;Xu et al., 2017). All of these data sets point to high permeabilities between 10 −14 and 10 −12 m 2 (and sometimes up to 10 −10 m 2 ) in the oceanic layer 2A decreasing exponentially with depth to values of 10 −16 to 10 −14 m 2 in layer 2B and below-with some estimates based on hydrothermal heat fluxes and tidal phase shifts giving higher values of up to approx. 10 −11 to 10 −10 m 2 (Crone et al., 2011;Wilcock & McNabb, 1996). While the absolute values, the exponential dependence on depth, and the large permeability contrast between layer 2A and 2B appear intuitive and plausible, the inferred values pose a major problem to hydrothermal flow simulations. Numerical experiments consistently show that high permeabilities lead to the increased entrainment of cold ambient seawater which causes a temperature drop that is difficult to reconcile with high vent temperatures (Andersen et al., 2015;Driesner, 2010;Lowell et al., 2003Lowell et al., , 2007. Mineral precipitation reactions are one plausible mechanism that can affect hydrothermal flow and modulate mixing processes. Cann and Strens (1989) formulated the clogged shell model to explain episodic megaplume emissions. Their model aimed at reconciling the required high permeabilities within the hydrothermal circulation zone necessary to sustain megaplume emissions with the much lower baseline flow rates at black smoker systems. Megaplume emissions may occur when tectonic activity and/or hydrofracturing reopen pathways in the upflow zone that were previously clogged by sulfide and quartz precipitation during the waning phases of the previous megaplume event. Fontaine et al. (2007) extended the clogged outer shell model to model the effect of a low permeability shell on vent salinities and exit temperatures. Probably the most spectacular manifestations of reaction-modulated transport in black smoker system are hydrothermal chimneys at the ocean floor. In most cases, they are cemented by and/or enveloped by anhydrite (Hannington et al., 1998;Tivey, 2007). The seafloor chimneys focus the discharge of hot fluids over tens of meters into the cold water column. It is likely that similar mineral precipitations focus the upflow of hydrothermal fluids below the seafloor, for example in "pipe-like" vein networks that tap the geothermal reservoir (Lowell et al., 2003;Tivey, 2007). Large volumes of anhydrite may precipitate near the surface (e.g., in the chimney, mound, and stockwork structures) and are also found in the deeper hydrothermal feeder zones (e.g., Humphris et al., 1995). This is evident in the ODP 158 data from the trans-Atlantic Geotraverse (TAG) hydrothermal mound, which was found to mainly consist of anhydrite, pyrite, silica, and chalcopyrite pointing to the importance of hydrothermal precipitates in focusing hydrothermal upflow towards individual vent sites (Humphris et al., 1995). Anhydrite, which has a retrograde stability, precipitates from heated seawater upon the reaction Ca 2+ + SO 4 2− = CaSO 4 in regions where temperatures reach between~145°C and 300°C (formation window, Figure 1). It is likely to occur in areas where cold seawater is entrained into the hot discharge flow (Bischoff et al., 1996;Fontaine et al., 2001;Kawada & Yoshida, 2010;Lowell & Yao, 2002). However, at temperatures exceeding~300°C, anhydrite is unstable and the evolving hydrothermal fluid becomes more acidic. The remaining SO 4 2− is progressively reduced to aqueous sulfide, and H 2 S becomes the dominant sulfur species in hydrothermal fluids (e.g., Chiba et al., 2001;Foustoukos & Seyfried, 2005;Seyfried & Janecky, 1985). At this stage, the fluid further reacts with the basaltic wall-rock and becomes enriched in, for example, calcium, iron, copper, zinc, lead, and other trace metals as well as H 2 S (e.g., German & Von Damm, 2003;Von Damm, 1995). The hot rising hydrothermal fluid is therefore enriched in metals (incl. Iron, copper, zinc, lead, and other trace metals), H 2 S, and Ca. As the fluid reaches maximum temperatures (~400°C), it becomes increasingly buoyant and starts to ascend from the reaction zone towards the seafloor. Around this upflow zone, additional amounts of anhydrite can precipitate, when the fluid cools to temperatures within the thermal formation window of anhydrite upon mixing with entrained SO 4 2− -rich seawater. These processes are likely to build a chimney-like anhydrite structure around the upflow zone. Finally, sulfide minerals such as pyrite and chalcopyrite precipitate directly from the high-temperature fluids upon cooling below~310°C in response to conductive as well as adiabatic cooling and/or the mixing with seawater ( Figure 1). These processes may be particularly important for enabling hot hydrothermal fluids to reach the seafloor to form metal deposits. The question is how widespread is this process? Under which conditions and over which time scales does anhydrite begin to precipitate and form a seal around the hydrothermal upflow zones? And, what are the implications for large-scale fluid flow and heat removal at mid-ocean ridges? About 90% of the axial discharge at mid-ocean ridge appears to be by diffuse, low-temperature flow, implying widespread shallow mixing and cooling between hot (~350°C) hydrothermal fluids and ambient seawater (Baker et al., 1993;Ginster et al., 1994;Ramondenc et al., 2006;Rona & Trivett, 1992;Schultz et al., 1992;Schultz & Elderfield, 1997). If only 10% of the heat transport at the ridges occurs at black smoker temperatures, then the amount of metal that can be delivered to the seafloor within the neovolcanic zone is much lower than the amount initially mobilized from the reaction zone (Hannington, 2013). Thus, anhydrite precipitation may be a fundamentally important process for metal fluxes to the seafloor.
Here, building on previous modeling work by, for example, Fontaine et al. (2001), Lowell et al. (2003), and Kawada and Yoshida (2010), we investigate the interplay between hydrothermal flow, permeability changes, chemical reactions, and metal transport within a 2-D thermo-hydro-chemical model that simulates hydrothermal flow along with anhydrite and pyrite formation in the subseafloor.

Numerical Model of Hydrothermal Flow
The hydrothermal flow calculations are based on Darcy flow of pure water, in which fluid velocities are proportional to the excess pressure gradient via the fluid's viscosity and the rock's permeability: where v ! is the Darcy velocity, k is permeability, μ f is the fluid's dynamic viscosity, ρ f is the fluid's density, p is pressure, and g ! is gravitational acceleration vector. The fluid and rock properties are referred to with subscripts f and r, respectively. All parameters and their values are listed in Table 1.
Mass conservation is expressed as with ϕ being the rock porosity. Note that compaction phenomena are neglected, and porosity is taken out of the time derivative. Closure of pore space by precipitation reactions is tracked via saturation changes. Substituting Equation 1 into 2 and noting that the fluid's density is a function of temperature T and pressure p yields the pressure equation, ϕρ where α f and β f are the fluid's thermal expansivity and compressibility, respectively. Energy conservation of a single-phase fluid can be expressed using a temperature formulation, with c p being heat capacity, λ r is the bulk thermal conductivity of porous rock, and S f fluid saturation (the remaining volume fraction of the pore space not occupied by precipitated mineral phases). Fluid and rock are assumed to be in local thermal equilibrium (i.e., T = T r = T f ) so that the mixture appears on the left-hand side of Equation 4. Changes in temperature depend on conductive heat transport, advective heat transport by fluid flow, heat generation by internal friction of the fluid, and pressure-volume work. All fluid properties are functions of both pressure and temperature and are evaluated from precalculated lookup tables based on the IAPS-84 formulation of water and steam properties. The tables have been computed using the program PROST 4.1 (Bauer, 1998). The hydrothermal flow model is implemented using a finite element formulation, and the details can be found in Hasenclever et al. (2014) and its supplement.

Chemical Transport
The dissolved aqueous species Ca 2+ and SO 4 2− involved in anhydrite formation are transported by fluid flow. Conservation of each chemical species can be expressed as, where C i is the concentration of each species, which is expressed as mass fraction in unit of kg-ion/kg-solution, S f is fluid saturation, Q i is the source term which represents concentration change when precipitation or dissolution occurs, and D c is the chemical diffusion coefficient. In addition, we track a quantity called pyrite formation potential (see below), which also has units of kg/kg and represents the amount of pyrite that can still precipitate from the hydrothermal fluid. Note that we do not account for a density change in the fluid as a function of chemical composition.

Chemical Reaction Model
When a hydrothermal fluid cools conductively and/or by mixing with ambient seawater, it undergoes a sequence of precipitation reactions. Here, we focus on two hydrothermal minerals, anhydrite, and pyrite, which are likely to have the biggest impact on flow pattern. The temperature-dependent saturation limit of anhydrite in seawater is expressed via the solubility product (e.g., Krauskopf & Bird, 1995). (T) is the apparent solubility product in units of mol 2 /kg 2 (e.g., Kawada & Yoshida, 2010) of anhydrite. To be consistent with units of concentration in Equation 5, we transform the units in Equation 6 to mass fraction where M Ca 2þ ¼ 40; M SO 2− 4 ¼ 96 in units of g/mol are the molecular weights of Ca 2+ and SO 2− 4 , respectively. Here, K a sp T ð Þ is the non-dimensional effective solubility product of anhydrite and we here use the parameterization of Kawada and Yoshida (2010), which is plotted in Figure 2a: 4.1 × 10 −4 kg/kg and 2.7 × 10 −3 kg/kg, respectively. Anhydrite precipitates when seawater is heated beyond~145°C (the intersection of dotted line and solid curve), and the anhydrite mass fraction as a function of temperature is shown by the green line and green axis. In (b), the solid black line denotes the "pyrite formation potential," which is defined as the amount of pyrite that can still precipitate from the fluid upon cooling. The initial value of 67.2 mg/kg was determined by reaction path modeling using the initial fluid composition listed in Table 2. The red line shows the amount of precipitated pyrite.
K a sp T ð Þ ¼ 6:124 × 10 −12 e 5:0594×10 3 Tþ273 ; Precipitation of anhydrite occurs when C Ca 2þ · C SO 2− 4 > K a sp and dissolution happens if the concentrations are lower than the saturation state and anhydrite is present. The related changes in concentration can be expressed as with dC ca 2+ and dC SO4 2− being the concentration changes of Ca 2+ and SO 4 2− so that the fluid composition is back in equilibrium with the effective solubility product. dC Ca 2þ can be calculated from Equation 9 by a quadratic equation with one or two roots. For precipitation, dC Ca 2þ is the negative root, while in the dissolution case it is the positive root. In addition, to keep mass conserved, we compare the theoretical mass change of Ca 2+ in the fluid with mass of Ca 2+ in anhydrite and adopt the minimum of them.
where min is the minimum function.
Implementing pyrite-forming reactions is less straightforward because dissolved iron and sulfur in a hydrothermal fluid are mostly present in form of aqueous complexes. When a hydrothermal fluid approaches metal saturation during the process of cooling and/or mixing with seawater, pyrite precipitates as part of a sulfide and nonsulfide mineral assemblage. Even though pyrite represents the main sink for Fe, it can also be accommodated in other mineral such as chalcopyrite, isocubanite and magnetite. The abundance of these minerals, their mass fractions, and the timing of mineralization are controlled by the source-rock geochemistry, the initial fluid chemistry, and the depositional mechanism, for example, changes in T, pH, redox state (Fuchs et al., 2019). To address the complex conditions of pyrite formation, preceding geochemical model simulations including speciation calculations are required. For these models, we use the thermodynamic constants for aqueous, mineral, and gaseous species from the most recent version of the SUPCRT92 database (Johnson et al., 1992;Shock et al., 1997;Sverjensky et al., 1997) augmented by some additional thermodynamic data ( Table S1 in the supporting information). The B-dot equation, an extension of the Debye-Hückel equation (Helgeson, 1969), was used to calculate activity coefficients for charged aqueous species. An activity coefficient of one was assumed for neutral aqueous species, except nonpolar species for which activity coefficients were used (Drummond & Ohmoto, 1985). In order to calculate the speciation of aqueous and gaseous components in hydrothermal fluid and in seawater during conductive and adiabatic cooling (from vent to ambient temperature), we use the software package Geochemist's Workbench (GWB) 12 (Bethke, 2008), which uses a law of mass action approach to compute of equilibrium state of the system. For this purpose, we have integrated the thermodynamic data of the aqueous and mineral species in the most recent version of the SUPCRT92 database (Johnson et al., 1992;Shock et al., 1997;Sverjensky et al., 1997) into GWB using the software tool DBCreate (Kong et al., 2013). DBCreate generates the required file of equilibrium constants for GWB from the native format in SUPCRT92. This combination of SUPCRT92 and GWB has been successfully used in previous studies (Fuchs et al., 2019;Simmons et al., 2016;Webber et al., 2017) and is a convenient way of performing model calculations under seafloor hydrothermal conditions of up to 350°C and within single-phase fluid states.
The presented simulations have been done in the context of mid-ocean ridge-associated hydrothermal systems, and we have used an average composition of an end-member high temperature fluid from the East Pacific Rise (EPR) as input (Hannington et al., 2016)-see Table 2. Figure 1 presents a reaction-path model showing the mass of mineral precipitation during the conductive cooling of an average EPR fluid of given composition from 350°C to an ambient temperature of 2°C at constant pressure of 500 bar and closed 10.1029/2019JB019035 Journal of Geophysical Research: Solid Earth system conditions. The sulfide minerals chalcopyrite and pyrite are the first to precipitate and therefore likely the ones with the largest impact on high-temperature (>300°C) fluid flow. Only when the fluid reaches lower temperatures (<180°C), the minerals sphalerite, galena, gold, silver, and quartz start to precipitate. The models exhibit that among the sulfide minerals, pyrite is abundant in highest mass fractions; therefore, we here only consider its impact on hydrothermal flow.
The most consistent approach to resolve pyrite precipitation and dissolution in a reaction-transport hydrothermal model would be to track the entire chemical composition of the fluid and to recompute the stable mineral assemblage and fluid composition at every computational point. While this will be hopefully feasible in future reaction-transport models, we feel it is appropriate, for now to explore a simple parameterization designed to single out the effect that pyrite precipitation may have on hydrothermal upflow. Hence, we make the simplifying assumption of a closed system that conductively cools from hydrothermal to seawater temperature. Under these assumptions, pyrite formation is limited by the availability of aqueous iron species (Figure 1). From the amount of pyrite precipitating from the fluid (Figure 1a), we derive a quantity that we call pyrite formation potential (Figure 2b). This quantity is a function of temperature only (for a given starting composition) and for practical purposes, we fit a polynomial to it: K p ¼ −2:4608 × 10 −10 T 5 þ 2:3537 × 10 −7 T 4 − 7:7876 × 10 −5 T 3 þ 1:0235 × 10 −2 T 2 − 0:2092T þ 1:6374: This quantity can be used in a very similar way as the anhydrite solubility product discussed above. If the fluid's pyrite formation potential is higher than its local temperature-dependent value, pyrite precipitation

10.1029/2019JB019035
Journal of Geophysical Research: Solid Earth occurs. If the fluid's value is lower than its local temperature-dependent value and pyrite is already present, dissolution of pyrite occurs.
Changes in the fluid's pyrite formation potential can be calculated by Equation 12, Pyrite will precipitate when dC p is negative, and if dC p is positive, the pyrite will dissolve. To keep mass conservation, dC p should be determined by Equation 13 when in the pyrite dissolution situation.
where dM p is the mass change of pyrite and M p is the mass of pyrite.
After all concentration changes are calculated, the relative saturations can be updated:

Permeability Feedback From Reactions
Fluid-rock interaction has a major influence on both rock and fluid properties (Cann et al., 2015;Lowell et al., 2003), and here we focus on reaction-induced permeability changes and how they affect fluid flow and vent temperatures. The interplay of permeability and fluid flow has been investigated in vent fluid studies (e.g., Shinohara, 2008;Von Damm, 1990), in numerical simulations (e.g., Scott & Driesner, 2018;Symonds et al., 2001), and in laboratory experiments (e.g., Bischoff et al., 1996;Foustoukos & Seyfried, 2007). Numerically, the permeability feedback of reactions is usually based on the Kozeny-Carman correlation (Bear, 1972), for example in CSMP++GEM (Yapparova et al., 2017) and TOUGHREACT (Xu et al., 2006). In our model, the permeability feedback of mineral precipitation is simplified as where k 0 is the initial permeability and S f is the fluid saturation, see Equation 14a. Thus, a closed system of equations is assembled, and the solution strategy is outlined in Figure 3.

Scope and Model Limitations
Fluid-rock interaction and dissolved species transport in submarine hydrothermal systems are highly complex processes. The presented parameterizations are only approximations of natural systems and should be used with these limitations in mind: 1. The model is single-phase, and the chemical reaction parameterizations are not directly transferrable to multiphase situations, where chemical species partition between vapor and liquid phases. 2. The chemical parameterizations are only valid for T < 350°C, and no chemical changes occur until the hydrothermal fluid has cooled from its initial temperature of 400°C down to 350°C. This restriction has been made to stay within the validity range of the underlying thermodynamic calculations. 3. The complexity of sulfide mineralization is not fully captured by the presented parameterization, and it should not be used to "understand" chemical systems. It has been solely derived to study feedbacks between flow and pyrite precipitation reactions. 4. The "pyrite formation potential" is only valid for the assumed initial average EPR fluid composition and closed system conditions. See, for example, Fuchs et al. (2019) for information on how initial fluid composition affect precipitation sequences. 5. We use prescribed initial compositions of seawater and hydrothermal fluids. Element mobilization by fluid-rock interactions is not resolved.
These assumptions allow us to single out the effect that permeability changes induced by mineral precipitation reactions may have on the overall flow patterns. By demonstrating the importance of this feedback, we hope to inspire future work that will make a more rigorous coupling between multiphase transport and a comprehensive thermodynamic model.

Geometry
The model setup is illustrated in Figure 4 and has been designed to be simple yet able to capture the key aspects of hydrothermal transport through a zone of enhanced permeability. The pipe-like inlet in the lower region mimics a focused upflow zone in the deeper crust (e.g., layer 2B) where permeability is sufficiently low (k c = 10 −15 m 2 ) to allow for high-temperature fluid flow. The upper region represents a zone of higher permeability (e.g., extrusive layer 2A), has a thickness h ext = 300 m, and a higher permeability k ext . To avoid boundary effects, the lateral extent of the domain has been set to a large value of 1,200 m. Both regions are discretized with a triangular mesh with the highest resolution of 0.3 m in the central upflow region.

Boundary and Initial Conditions
The different boundary conditions are shown as different colors in Figure 4: sides (black), outlet (magenta), and inlet (red). In correspondence to the heat and chemical transport equations in sections 2.1 and 2.2, boundary conditions for T, P, and C i are applied ( Table 2). The sides (black line in Figure 4) are impermeable  and have zero concentration and temperature gradients. Constant pressure (30 MPa) is assumed at the top boundary (magenta) and at inflow nodes a constant temperature (5°C) and composition (seawater) is assumed. At outflow nodes, the thermal and concentration gradients are zero to allow for free venting.
The bottom of the pipe-like inlet has a constant temperature (400°C) and a prescribed fluid mass influx of 11.5 g/m/s. This mass flux value is consistent with high-temperature (>350°C) venting for permeabilities <10 −14 m 2 (Driesner, 2010, Figure 5a). The injected fluids have a prescribed chemical composition taken from a representative EPR fluid (Hannington et al., 2016), and the values are summarized in Table 2. Note that fluid-rock interactions are not implemented; instead, the influx boundary condition is used to inject fluids with a prescribed chemistry. The modeling domain is initialized to be at a constant temperature of 5°C and seawater composition.

Baseline Simulations
An initial set of four reference simulations with different permeabilities in the top layer is conducted to illustrate the overall fluid dynamic system behavior. Mineral precipitation and dissolution do not occur in these reference simulations. Figure 5 summarizes the results. With increasing permeability, the vertical mass flux at the seafloor increases to values higher than the prescribed influx of 11.5 g/m/s at the basal inlet due to entrainment of colder ambient fluids causing a decrease in vent temperature. Velocities increase and the upflow zone becomes more focused and narrower. High temperature venting (>350°C) only occurs at permeabilities <1 × 10 −14 m 2 , which again demonstrates the fundamental problem of reconciling high-temperature venting with high permeabilities in numerical models.

Impact of Anhydrite and Pyrite Precipitation
The next set of simulations not only is the same as the reference set but also allows for precipitation and dissolution of anhydrite and pyrite. In these simulations, we use the upper limit 4,000 mg/kg Ca 2+ in the hydrothermal fluid (Table 2), which is also the value that Kawada and Yoshida (2010) used in their study; the effects of varying the Ca content of the hydrothermal fluids will be discussed later. Figure 6 shows the temporal evolution of a reactive-transport simulation with k ext = 4 × 10 −14 m 2 (see also Movie S1). Chimney-like structures of anhydrite and pyrite progressively form around the central upflow zone, which become effectively separated from the recharge flow. These sidewalls form in three distinct ways. First, anhydrite starts precipitating from the warming seawater at temperatures higher than~145°C that recharges the central upflow of entrained fluids. This process forms the outermost part of the sidewalls starting at the interface between the two layers and grows towards the surface forming the outer anhydrite wall. The second process forms an inner anhydrite wall due to the cooling and mixing of hot hydrothermal fluids with entrained seawater in the central upflow zone. The upwelling hydrothermal fluid is devoid of sulfate but enriched in calcium due to fluid-rock interactions in the reaction zone (which are not resolved in this study). When these fluids mix with seawater, further anhydrite precipitation occurs around the central upflow zone. Lastly, pyrite forms when the hydrothermal fluids cool below~300°C forming a cap-like pyrite structure.

Journal of Geophysical Research: Solid Earth
After 250 years, mixing and cooling is still intense and pyrite precipitates well below the seafloor (Figures 6a-6c). Over time, more anhydrite forms and mixing is reduced. Consequently, the ascending hydrothermal fluids remain hot enough to carry metals, and pyrite is precipitated further up the upwelling zone ( Figure 6). After 1,000 years, a pyrite deposit begins to form at the seafloor, and after 5,000 years, a stable upflow zone, separated by anhydrite walls from ambient seawater, is established, and pyrite is more effectively transported to the seafloor as the shoaling of the pyrite precipitation zone shows. Figure 7 shows that these precipitations reactions scale positively with the permeability of the upper layer. The more seawater is entrained, the higher is the rate of anhydrite formation. Likewise, the time to create a chimney structure scales directly with permeability. These precipitation reactions result in a strong drop in permeability around the central upflow zone (Figure 7, lowermost panel) causing a significant reduction in mixing between seawater and hydrothermal fluid.
This finding is confirmed in Figure 8, which shows vent temperature evolution and the final discharge mass flux for the different runs. While vent temperatures remain stable in simulations without reactions, they progressively increase in those with reactions. The reason for this is the reduction of seawater entrainment as the upflow zone becomes progressively isolated by precipitation reactions. Total final discharge mass fluxes (15.6-26.5 g/m/s) are significantly higher than the inflow boundary condition (11.5 g/m/s), which shows how significant the entrainment of ambient seawater is over the investigated range of crustal permeabilities. Note that these 2-D values are per meter of ridge axis; to compare to observed/inferred hydrothermal discharge rates, the values need to be multiplied by a representative along-axis distance. Entrainment results in dilution and cooling of the upwelling hydrothermal fluid. Such mixing and entrainment processes become more important with higher permeability, which is reflected in the model with the highest permeability (k = 1 × 10 −13 m 2 ), where the total discharge is more than twice the inflow boundary condition (26.5 vs 11.5 g/m/s).
After 5,000 years, the predicted increases in vent temperature are still moderate (<25°C) yet significant due to the strong temperature dependence of pyrite solubility, which we here for simplicity also treat as a proxy for metal transport (Equation 11). In order to constrain the long-term evolution and maximum plausible increase in vent temperature, we have performed a number of additional model runs. Figure 9 shows the

Journal of Geophysical Research: Solid Earth
long-term evolution of models with permeabilities of 4 × 10 −14 m 2 and 1 × 10 −13 m 2 plus one end-member run that illustrates the maximum plausible increase. Those runs show that vent temperatures continue to increase as the chimney structures become increasingly less permeable. The end-member case has been calculated by taking the precipitation pattern predicted by the 4 × 10 −14 m 2 model run and setting the relative permeability to a very small value, effectively making the chimney structure impermeable. In this limiting case, the maximum increase in vent temperature is 90°C, showing how significant the effect of precipitation reactions can be.

Sensitivity
The presented simulations have shown that the precipitation of pyrite and especially of anhydrite can have a profound effect on predicted flow patterns, flow rates, and thereby vent temperatures. In a final set of simulations, we have explored the sensitivity of the results to variations in fluid chemistry. In the context of this study, the main uncertainty is the Ca concentration in the injected fluid. In natural systems, Ca is mobilized from the host rock at high temperatures within the root (reaction) zone of a black smoker system. The amount of calcium mobilized in this way is a complex function of fluid-rock ratios, fluid chemistry, rock type, and temperature. Here, we do not resolve this complexity but treat the Ca content of the injected hydrothermal fluid as a prescribed model parameter. The observed range of Ca concentrations in EPR vent fluids is (8-4,248 mg/kg), with the average being 625 mg/kg (Hannington et al., 2016;Shanks, 2001). All simulations presented above used a value of 4,000 mg/kg, which is also the value that Kawada and Yoshida (2010) used. Figure

Journal of Geophysical Research: Solid Earth
injected fluid. However, the scaling between the two is small: a 10-fold increase in Ca content of the injected fluid results only in a 30% increase in total anhydrite formed after 3,000 years simulation time. The main reason for this sublinear scaling is that the amount of anhydrite forming is mainly limited by availability of sulfate during mixing with seawater, and the rate of seawater entrainment is primarily a function of permeability ( Figure 10b). Interestingly, the rates of anhydrite formation decrease with increasing simulation time and the amount of anhydrite formed becomes progressively insensitive to permeability. The reason is that upflow becomes increasingly separated from entrained seawater due to the clogging of the pore space.

Discussion
This paper links reaction-path modeling in seafloor hydrothermal systems to heat and fluid flow simulations in the context of subseafloor mineral precipitation. The model runs were constructed to explore the conditions under which anhydrite precipitates in the subseafloor and how that affects flow patterns and heat transport.
The models show that mineral precipitation reactions are significant contributors to isolated and focused hydrothermal upflow. Anhydrite precipitation in the subsurface forms chimney-like structures that are in many ways analogous to hydrothermal chimneys observed at and above the seafloor. These chimney-like structures, sometimes referred to as pipe-like vein networks in the geological literature (Hannington et al., 1998), form in two ways: (1) by direct precipitation from warming seawater and (2) through subsequent mixing with calcium-rich, sulfate-free hydrothermal fluids. These two reactions are primarily limited by the availability of sulfate, which comes exclusively from seawater. As a result, anhydrite formation scales directly with the amount of seawater entrained by the system and therefore the permeability. In the presented simplified simulations, pyrite formation, in turn, is primarily controlled by conductive and adiabatic cooling, which is regarded as the main driver of sulfide mineral precipitation (e.g., Tivey, 1995;Tivey & McDuff, 1990). In natural systems, such processes are more complex, and pyrite mineralization is further influenced by the magmatic contribution of metals (Fe) and H 2 S, the physicochemical regime (e.g., pH, redox state), geochemical gradients in the cooling fluids, seawater mixing, and/or phase separation phenomena. Resolving this complexity is beyond the scope of this study but will hopefully be addressed in future thermo-hydro-chemical models that will couple transport and the thermodynamics of hydrothermal systems.
One important feedback in the presented model system appears to be that the increase in the total mass flux due to entrainment of cold fluids, which initially causes a decrease in vent temperature, also results in an increase in the rate of anhydrite deposition at depth that eventually focuses the upflow, thus increasing  discharge temperatures. Higher temperature venting occurs at the seafloor despite the high bulk permeability of the shallow crust. Time-series show that increased fluid entrainment into the system causes increased subseafloor precipitation and the development of a focused upflow zone (Figure 7). The formation of a pyrite cap is also an important part of this process. Unlike anhydrite, which is ultimately replaced (or dissolved because of its retrograde solubility), a pyrite-sealed zone may be particularly important for focusing fluids near the seafloor over time. In this case, continued hydrothermal upflow (possibly during multiple episodes of venting), could result in overpressuring and fracturing of the relatively impermeable pyrite cap, providing long-lived pathways for high-temperature fluids to reach the seafloor. In this way, anhydrite, and to lesser extents metal sulfides such as pyrite, precipitation ensures that a large proportion of the metals initially mobilized from the deep reaction zones actually make it to the seafloor forming present-day larger deposits.
The above findings are related to the fundamental question of how high-temperature upflow zones become established in otherwise highly permeable oceanic crust. At slow-spreading ridges, where large and deeply penetrating faults are preferential pathways for hydrothermal flow, the proportion of focused flow is high, which can be explained by structural permeability contrasts (e.g., Andersen et al., 2015). Fast spreading ridges are characterized by numerous shallow faults and very high crustal permeabilities down to at least 300 m, corresponding to the top of the dike section resulting in a high proportion of diffuse flow (Elderfield & Schultz, 1996;Fornari et al., 2012). Here, reactions like anhydrite precipitation may play a key role in focusing flow (Lowell et al., 2003Sleep, 1991). Building on those previous studies, the models presented here show that focused flow of high-temperature fluids can be established in highly permeable substrates even in the absence of a controlling structure via the formation of sub-seafloor anhydritechimneys.
While these findings point to interesting interrelations between fluid focusing and mineral precipitation reactions, it is likely that also other processes are important for rapid fluid focusing. In the presented calculations, anhydrite-chimneys form on 100-1,000 year timescales, while individual magmatic events have a recurrence rate of~10-20 years at fast-spreading ridges. One possible way to reconcile these differing timescales is to assume even higher permeabilities in the shallow crust. Lowell et al. (2003) showed that anhydrite precipitation can help establishing a focused flow zone even within years following a magmatic event if crustal permeabilities are sufficiently high (10 −11 to 10 −12 m 2 ). We here explored permeabilities of 10 −13 -10 −14 m 2 , yet our results point in the same direction with the rate of anhydrite precipitation scaling directly with permeability. Furthermore, we have driven the models with a prescribed mass influx of 0.015 kg/s/m, which in combination with the 400°C boundary condition corresponds to a hydrothermal heat flux of 23 MW per kilometer ridge axis. This is approximately the long-term thermal plus magmatic heat flux of a slow-spreading ridge opening at a full rate of 30 mm/year (Sinha & Evans, 2004). Inferred vent field hydrothermal heat fluxes from point measurements and water columns observations can be up to two orders of magnitude higher ranging between 28 MW at Broken Spur to 1,700 MW at TAG along the slow-spreading Mid-Atlantic Ridge and between 50 MW at Salty Dawg and 551 MW at Northern Cleft along the intermediate Juan de Fuca Ridge (Lowell et al., 2013). This shows that hydrothermal heat flux can transiently be much higher than the assumed long-term average heat flux used in the models. Under such transient conditions, chimney structures are likely to form more rapidly thereby reconciling the different time scales. It is also possible that mechanical feedbacks, such as hydrofracturing (Weis et al., 2012), and or fracture flow instead of pervasive Darcy flow (Mezon et al., 2018) help to focus fluid flow on short time scales.
It should also be noted that our reaction-transport model is highly simplified. We have here singled out the effect of anhydrite and pyrite on flow within the upflow zone. Natural systems, precipitation of other mineral phases is also likely to impact on flow. For example, quartz precipitation close to the reaction zone has been linked to heat transfer from the driving heat source into the hydrothermal flow zone (Steele-MacInnis et al., 2012). We have also neglected fluid-rock interaction with the consequence that element mobilization and redeposition cannot be self-consistently resolved. It is our hope that our here presented findings will help to formulate future self-consistent reaction-transport hydrothermal models.

Conclusions
We have shown how mineral precipitation reactions form chimney-like structures in the sub-seafloor through which hydrothermal fluids may eventually discharge without mixing significantly with ambient 10.1029/2019JB019035

Journal of Geophysical Research: Solid Earth
colder ambient fluids. This chemical feedback may help to explain how focused high-temperature flow can be established after a magmatic event, for example, within highly permeable shallow oceanic crust. An important consequence appears to be the increase in the total mass flux with increasing permeability due to entrainment of cold fluids that initially cause not only a decrease in vent temperature but also an increase in anhydrite deposition at depth. Eventually, the upflow zone becomes sealed and discharge velocities and temperatures increase. The result is higher-temperature venting at the seafloor despite the increased permeability below the seafloor. Interestingly, this process is likely to leave little imprint in the geological record. Anhydrite is not preserved as the system cools, leaving behind only sulfide-bearing vein networks and alteration "pipes" with no trace of the original anhydrite (Sleep, 1991). The only evidence that anhydrite might have formed would be the persistence of alteration minerals formed at temperatures that would have caused anhydrite deposition. An interesting question for future investigation is the process of focusing fluid flow in ancient hydrothermal systems where SO 4 2− concentrations were too low (de Ronde et al., 1997) to allow precipitation of anhydrite (e.g., in early Archean oceans or in younger reduced basins). In this case, focusing of hydrothermal discharge could not have been due to sealing by anhydrite, but other minerals such as carbonate and/or quartz must have been involved. Dissolution of anhydrite formed at the ridge as the crust moves away from heat sources will leave behind only the metal sulfide structures which make up present-day subseafloor stockworks.

Data Availability Statement
All data files and processing script are accessible through figshare 10.6084/m9.figshare.10347839