Understanding controls on hydrothermal dolomitisation: insights from 3D Reactive Transport Modelling of geothermal convection. Solid Earth,

. The dominant paradigm for petrogenesis of high-temperature fault-controlled dolomite, widely known as “hy-drothermal dolomite” (HTD), invokes upwelling of hot ﬂuid along faulted and fractured conduits from a deep overpressured aquifer. However, this model has several inherent ambiguities with respect to ﬂuid sources and their dolomitisation potential, as well as mechanisms for delivering enough of these reactive ﬂuids to form substantial volumes of dolomite. Here, we use generic 2D and 3D reactive transport simulations of a single transmissive fault system to evaluate an alternative conceptual model whereby dolomitisation is driven by seawater being drawn down into the subsurface and heated. We examine the evolution of ﬂuid chemistry and the distribution of diagenetic alteration, including predictions of the rate, distribution, and temperature of HTD formation, and consider the possible contribution of this process to the Mg budget of the world’s oceans. The simulations suggest that it is possible for convection of seawater along the fault damage zone to form massive dolomite bodies that extend hundreds of

Abstract. The dominant paradigm for petrogenesis of hightemperature fault-controlled dolomite, widely known as "hydrothermal dolomite" (HTD), invokes upwelling of hot fluid along faulted and fractured conduits from a deep overpressured aquifer. However, this model has several inherent ambiguities with respect to fluid sources and their dolomitisation potential, as well as mechanisms for delivering enough of these reactive fluids to form substantial volumes of dolomite. Here, we use generic 2D and 3D reactive transport simulations of a single transmissive fault system to evaluate an alternative conceptual model whereby dolomitisation is driven by seawater being drawn down into the subsurface and heated. We examine the evolution of fluid chemistry and the distribution of diagenetic alteration, including predictions of the rate, distribution, and temperature of HTD formation, and consider the possible contribution of this process to the Mg budget of the world's oceans.
The simulations suggest that it is possible for convection of seawater along the fault damage zone to form massive dolomite bodies that extend hundreds of metres vertically and along the fault within a timescale of a few tens of thousands of years, with no significant alteration of the country rock. Dolomitisation occurs as a gradient reaction by replacement of host limestones and minor dolomite cementation, and it results in the discharge of Mg 2+ -poor, Ca 2+ -rich fluids to the sea floor. Fluids sourced from the basement contribute to the transport of heat that is key for overcoming kinetic limitations to dolomitisation, but the entrained seawater provides the Mg 2+ to drive the reaction. Dolomite fronts are sharper on the "up-flow" margin where Mg 2+ -rich fluids first reach the threshold temperature for dolomitisation, and the "down-flow" dolomite front tends to be broader as the fluid is depleted in Mg 2+ by prior dolomitisation. The model demonstrates spatial contrasts in the temperature of dolomitisation and the relative contribution of seawater and basement-derived fluids which are also commonly observed in natural fault-controlled dolomites. In the past, such variations have been interpreted in terms of major shifts in the system driving dolomitisation. Our simulations demonstrate that such changes may also be a product of emergent behaviour within a relatively stable system, with areas that are dolomitised more slowly recording the effect of changes in fluid flow, heat, and solute transport that occur in response to diagenetic permeability modification.
Overall, our models robustly demonstrate that hightemperature fault-controlled dolomite bodies can form from mixed convection and act as a sink for Mg in the circulating seawaters. In addition, comparison of our 3D simulations with simplifications to 2D indicate that 2D models misrepresent critical aspects of the system. This has important implications for modelling of systems ranging from geothermal resources and mineralisation to carbonate diagenesis, including hydrothermal karstification and ore genesis as well as dolomitisation.
of limestone, and/or overprinting of dolomitic country rocks (Gregg, 2004;Hollis et al., 2017;López-Horgue et al., 2010;Machel, 2004). The economic importance of HTD often lies in the associated hydrothermal dissolution and/or mineralisation, in particular Mississippi Valley-type Pb-Zn ore deposits (Davies and Smith, 2006;Qing and Mountjoy, 1994). More controversially, fluids thought to have formed faultrelated HTDs have also been interpreted as being responsible for the formation or overprinting of large volumes of matrix-replacement dolomite that form important hydrocarbon reservoirs and infer basin-wide hydrothermal activity (Al-Aasm, 2003;Davies and Smith, 2006;Sharp et al., 2010;Weissenberger et al., 2006). Machel and Lonnee (2002) were particularly critical of this suggestion and argued strongly against the use of the term hydrothermal (as distinct from geothermal) without evidence that the dolomites formed at temperatures of at least 5-10 • C higher than the country rock. This would require very rapid rates of both fluid flow and precipitation, and the physics of heat transport mean that it is practically impossible to form HTD away from highpermeability conduits. The requisite temperature data to define HTD sensu stricto is rarely presented; thus, this distinction, whilst technically correct, is often ignored. Rather, the term HTD has become widely applied to describe hightemperature fault-related dolomites.
Despite a substantial body of previous work describing HTD geobody distributions and dimensions in outcrops, along with petrographical and geochemical characterisation, important questions remain regarding the process of formation. These include the origin and chemistry of the dolomitising fluids, the volumes of fluids that can be released from deep sources, and the mechanisms controlling the distribution and rate of fluid flow and, thus, the temperature and rate of dolomitisation. Over-pressured aquifers exist at depth as a result of rapid burial and/or orogenic compression, and tectonic activity can create permeable pathways for the release of such fluids (Ge and Garven, 1994;Oliver, 1986;Whitaker et al., 2004). However, these processes can release, at the most, only a single pore-volume; thus, dolomitisation requires significant focusing of flow from large reservoirs (Consonni et al., 2010;Frazer et al., 2014;Kaufman, 1994). Potential dolomitising fluids are hypothesised as brines or saline connate waters, and this is supported by high-salinity fluid inclusions in the dolomite cements in some cases (Al-Aasm, 2003;Boni et al., 2000;Davies and Smith, 2006;Luczaj et al., 2006). However, a review of the thermodynamic potential of formation waters to drive dolomitisation suggests that fluids hosted in deep aquifers are close to dolomite-calcite equilibrium and, thus, have limited dolomitisation potential.
Given the ambiguities in existing conceptual models for HTD formation, we evaluate an alternative model for the formation of high-temperature dolomites by "top-down" circulation of Mg 2+ -rich fluids from an effectively infinite reservoir of seawater (or modified seawater) overlying a perme-able fault damage zone. In high-enthalpy areas, heating of groundwater at depth increases buoyancy and, where the critical Rayleigh number is exceeded (Nield, 1968), drives free convection that can entrain water that is cooler and therefore more dense. Controls on this flow system have been analysed numerically by means of 2D models of heat and fluid flow (Person and Garven, 1994;Simms and Garven, 2004), and they include sediment and fault permeabilities, spacing between faults, lithologic heterogeneity, and basal heat flow. These studies show convection cells that are characteristically laterally elongated, with an aspect ratio of 1.8, and inflow over broad areas but discharge focused over narrower zones. The power of this approach was demonstrated by simulations coupling large-scale free convection of groundwater flow with gravity-driven flow in the Dead Sea rift, Israel (Gvirtzman et al., 1997a, b), that replicate both thermal anomalies and the composition of springs that are mixtures of shallow groundwater with hot brine from a deep aquifer.
An increasing number of studies have suggested faultrelated dolomitisation by seawater or seawater-derived fluids (Carmichael and Ferry, 2008;Corbella et al., 2014;Haeri-Ardakani et al., 2013;Hollis et al., 2017;Wilson et al., 2007). For example, in the Hammam Faraun Fault (HFF) area, Egypt, strontium isotope and rare earth element (REE) signatures suggest that massive dolostone bodies limited to the damage zone of the fault formed from Oligo-Miocene seawater during development of the Suez Rift (Hirani et al., 2018;Hollis et al., 2017). We have previously suggested a conceptual model for the formation of these dolostones whereby Mg 2+ -rich seawater could have been drawn down into a convection cell along the fault, heated, and then ascended to form massive dolostone (Hollis et al., 2017). Testing this hypothesis by reference to the behaviour of modern analogues is challenging using direct observation of hydrological and chemical processes, and process-based modelling provides an alternative means for rigorous evaluation in such scenarios.
Following the early work of Wilson et al. (2001) and Jones and Xiao (2005), over the past decade process-based numerical modelling of heat and fluid flow coupled with the simulation of water-rock reactions using reactive transport models (RTMs) have contributed to the understanding of shallow, relatively low-temperature dolomitisation (Al-Helal et al., 2012;Garcia-Fresca, 2009;Whitaker and Xiao, 2010;Xiao et al., 2013). RTM simulations of HTD have focused on understanding controls on the nature of dolomitisation fronts, which are characteristically sharp in fault-related dolomites, and on the development of fault-associated stratabound dolostone bodies Xiao et al., 2013;Yapparova et al., 2017). These studies demonstrate the importance of the geometry of the fault, termination of the fault against a sealing layer, and distribution of permeability and precursor mineralogy in the country rock, but they are all predicated on a supply of hot, Mg 2+ -rich fluid at the base of the fault. Consonni et al. (2018) highlight some of the challenges in moving beyond these outcrop-scale models of HTD that arise from uncertainties in fluid composition and volumes, and they contrast these with the RTM simulation of reflux and thermal convection systems for which fluid flow and chemistry can be better constrained from an understanding of the geometry of the carbonate platform, relative sea-level, and climate. Where fluids are derived from compaction, the geometry of the system can help constrain fluid volumes, but specifying fluid chemistries is more challenging -even when the mineralogy of the compacting formation is known (Consonni et al., 2010;Frazer et al., 2014). This points to the need for a better understanding of deep aquifer and reservoir fluid chemistry to identify controls on dolomitisation potential among different lithological and geological settings.
Simulations of free convection within a magmatic hydrothermal system volcanic caldera system exemplified how faults can either act as conduits for focusing the discharge of high-temperature fluids or as conduits for the influx of cool shallow fluids from a shallow reservoir (Jasim et al., 2015). Recently, the study of natural fault-controlled dolomite bodies have provided evidence for the mixing of seawater and mantle-derived fluids, during dolomitisation, implying that hot, deeply sourced brines at least provide a thermal drive for dolomitisation (Breislin et al., 2020;Koeshidayatullah et al., 2020a). The mixing that occurs between such fluids expelled from a fault and host-rock pore waters of marine salinity was invoked by Smith et al. (2008) in order to explain HTD in the Trenton and Black River reservoirs of eastern North America, following 1D RTM studies by Salas et al. (2007). However, in a limestone aquifer separated from overlying seawater by a confining layer, a single pore-volume of seawater would clearly be insufficient to provide the Mg 2+ required for the scale of hydrothermal dolomitisation described by Smith et al. (2008).
The objective of this study is to use RTM to evaluate whether, for a highly simplified generic single transmissive faulted system, physical and chemical constraints might allow the formation of significant volumes of dolomite by entrainment of seawater into a geothermal convection system. It is not our objective to reproduce the development of a real dolostone body; rather, we consider a generic scenario with a simple geometry under conditions that are clearly defined.
The key question addressed is whether the proposed flow system could account for the volumes and geometries of dolostone geobodies and temperature characteristics observed in fault-related HTDs, given heat and solute transport, chemical thermodynamics, and kinetic constraints. We examine the evolution of fluid chemistry and the distribution of diagenetic alteration, including predictions of the rate, distribution, and temperature of dolomite formation, and the relative contribution of seawater and basement-derived fluids along with sensitivity analysis of certain model parameters. Additionally, we compare the results of our 3D simulations with those from 2D sections perpendicular to the plane of the fault (echoing this key element of many previous 2D models of HTD) as well as 2D models oriented along the plane of the fault.

Methods
Fluid flow and reactive transport simulations were performed using TOUGHREACT (Sonnenthal et al., 2017;Xu et al., 2011), which simulates multiphase fluid flow, heat, and solute transport with physical and chemical heterogeneity and incorporates feedbacks between diagenesis and evolving porosity and permeability (Xu et al., 2004). A number of generic simulations were run to evaluate the potential of geothermal convection of seawater to form fault-related dolomites within a single transmissive fault system open to the sea surface. Fluid flow and heat transport were simulated in 2D sections oriented both parallel and perpendicular to the plane of the fault and compared with that in a 3D domain to constrain an initial condition for the reactive transport simulations. The 3D domain has a dimension of 2.7 km along the plane of the fault, 2.5 km perpendicular to the fault, and is 5 km deep (Fig. 1). A symmetry boundary along the centre of the fault damage zone (FDZ) has been assumed, and the half model is discretised into 4158 active grid blocks, 150 m wide and 250 m deep, which vary in width from 10 m within the FDZ to 500 m at the boundaries. The simulated sequence comprises a 3000 m thick carbonate unit with a permeability of 2 × 10 −14 m 2 (≈ 20 mD) that is sandwiched between two 1000 m thick lower-permeability (10 −15 m 2 ≈ 1 mD) layers, herein termed "caprock" and "basement" respectively as generic terms indicating rocks of lower permeability relative to our reactive aquifer. An 80 m wide permeable FDZ extends vertically through the carbonate unit and the overlying layer to the sea floor. Matrix cells are isotropic, but within the FDZ vertical permeability is 10 times higher than horizontal permeability and the porosity is 3 % higher than that of the aquifer. The properties of all units are listed in Table 1. A full 3D flow model (2.7 × 2.5 × 5 km) has validated the assumed symmetry, and, as such, the presence or absence of a low-permeability fault gouge in the core of the fault would not significantly affect the result.
The two 2D models are based on slices through the 3D model, both perpendicular to the plane of the fault and along the FDZ (Fig. 2). The fault-perpendicular 2D model is 2.5 km wide and includes an 80 m wide FDZ from the surface to a depth of 4 km within a model that extends to a total depth of 5 km. The fault-parallel 2D model represents an 80 m wide FDZ that is one cell thick and extends 2.7 km laterally and to a depth of 4 km within a model that extends to a total depth of 5 km. Both models are discretised into grid blocks of identical dimensions and with identical rock properties to those in the 3D model; a single porosity-single permeability approach was used in all cases.  In all three models, the side boundaries are no-flow, and the bottom and top boundaries are constant temperature and pressure. Infinite volume cells at atmospheric pressure and 25 • C are implemented at the top boundary to represent a shallow overlying body of seawater. The basal cells are at 250 • C, giving an initial conductive temperature gradient of 45 • C km −1 , and a hydrostatic pressure gradient is imposed. The flux of fluids into the base of the FDZ is determined by the permeability of the 1 km thick basement. Initial parameters and boundary conditions for all model runs are summarised in Fig. 1. The resulting fluid flux is expressed as a Darcy velocity (in m yr −1 ), also known as Darcy flux or specific discharge, and is calculated by multiplying the pore velocity by porosity where all porosity is assumed to be effective. We compare the total mass flow, differentiating seawater and basement fluid drawn into the model, and the discharge of mixed fluid to the sea floor via the 80 m wide FDZ and from the matrix. We model solute transport by diffusion as well as advection, with the same diffusion coefficient (10 −9 m 2 s −1 ) for all species.
The baseline RTM was constructed using modern openocean seawater (Nordstrom et al., 1979) as the initial fluid throughout the reservoir and caprock. Within the basement the fluid is high-temperature (250 • C) modern seawater at equilibrium with calcite and dolomite and, thus, has a higher Ca 2+ concentration, a lower Mg 2+ concentration and alkalinity, and is slightly acidic ( Table 2). Geochemical calculations include 10 primary aqueous species and two minerals, calcite and dolomite. We use the concentration of Br − to trace the contribution of fluids from the basement, hereinafter expressed as the basal fluid fraction. The upper 1 km (matrix and FDZ) and basal 1 km (basement) are assumed to be unreactive. Initial mineralogy throughout the interval, 1-4 km depth, is 99 % calcite, with 1 % dolomite as a seed. As the rate of dolomite precipitation is much slower than that of calcite dissolution, calcite was modelled as an equilibrium thermodynamic mineral (fluid and solid equilibrating within a single time step). Thus, the rate of dolomitisation is determined by the kinetics of dolomite precipitation. The kinetic rate of dolomite precipitation was constrained by the laboratory experiments of Arvidson and Mackenzie (1999) in which dolomite was precipitated at 115-196 • C. Thus, we include the effects of temperature on both thermodynamics and kinetics of dolomitisation, as well as the evolving fluid composition and mineralogy. The baseline simulation assumes an effective reactive surface area (RSA) of 1000 cm 2 g −1 , representative of 25 µm diameter idealised grains ), but does not consider changes in the RSA due to diagenetic changes in texture.
Porosity changes in the matrix and fractures are directly tied to the volume changes as a result of mineral precipitation and dissolution, and the model estimates changes in permeability from these. TOUGHREACT tracks diagenetic changes in both porosity and permeability that result from changes in mineral volume fractions due to water-rock interaction, and it subsequently returns this feedback to modify fluid flow. The porosity-permeability relationship for the matrix is approximated by Eq. (1), the Kozeny-Carman equation (Bear, 1972), which is one of the simplest and most widely used models for porosity-permeability relationship (Costa, 2006;Ehrenberg et al., 2006): where k i and ϕ i are the initial permeability and porosity respectively. Within the fracture zone, we represent changes in permeability using the cubic law relation (Steefel and Lasaga, 1994), Eq. (2): In two additional 3D simulations, we explore the influence of heat flux on temperature, fluid flux, and dolomitisation by reducing the permeability of the basement from 10 −15 m 2 (≈ 1 mD) to either 5 × 10 −16 m 2 (≈ 0.5 mD) or 1 × 10 −16 m 2 (≈ 0.1 mD). This limits the effective flux of hot fluids from the basement whilst maintaining the same temperature at the base of the model. We further investigate the effect of basal temperature by specifying a lower temperature (200 • C) at the base of the model, giving an initial geothermal gradient of 35 • C km −1 (Fig. 1) which is comparable to the average thermal gradient of 31.3 and 32.0 • C km −1 in the onshore Netherlands region (Bonté et al., 2012) and geothermal areas in Germany (Agemar et al., 2012) respectively. Finally, we evaluate the sensitivity of the system to the permeability of the caprock, reducing this to 0.5 and 0.1 mD from 1 mD in the baseline model. Post-processing data analysis includes Mg 2+ -flux calculation and differentiation between dolomite cementation and replacement dolomite. The Mg flux is calculated as the product of the aqueous concentration of Mg 2+ (mol kg −1 ) and mass fluid flux (kg s −1 m −2 ), is reported in millimoles per second per square metre (mmol s −1 m −2 ), and provides both a better understanding of source-sink behaviour and a mechanism to track Mg 2+ circulation and exchange. The changes in the total volume of dolomite formed by the replacement of calcite (as commonly described by the stoichiometric equation, Reaction R1) is differentiated from dolomite cementa-tion (primary precipitation of dolomite, Reaction R2) using the reduction in the volume of calcite and adjusting for the molar volume difference between calcite (36.93 cm 3 mol −1 ) and dolomite (64.37 cm 3 mol −1 ).
This ignores the effect of any changes in the volume of calcite that may occur independent of dolomitisation. Because of the retrograde solubility of calcite, this calculation yields a maximum estimate of the fraction of dolomite that forms by replacement in the areas of increasing temperature in the direction of flow and a minimum estimate in areas of decreasing temperature.

Baseline 2D and 3D simulations of geothermal convection: fluid flow and heat transfer
The 2D and 3D simulations show the circulation of fluids by geothermal convection within the fault damage zone (FDZ) and the combined effects of conduction and advection of heat for a very simple planar fault (Figs. 2 and 3). Fluid circulation is driven by differences in fluid density generated by temperature contrasts between the surface and basement, mediated by the spatial distribution of permeability. The fault-perpendicular 2D model shows the FDZ acting as a conduit for advective discharge of hot fluid sourced largely from the basement, with a vertical Darcy velocity within the FDZ of up to 159 m yr −1 , some 1-2 orders of magnitude higher than rates in the matrix (Fig. 2a, b). Convection and heat transfer between the matrix and the FDZ mainly occur in the aquifer zone (1-4 km depth). A combination of conduction and the ascending flow of buoyant hot fluid from the basement (average Darcy velocity of 1.1 m yr −1 ) elevates the temperature throughout the FDZ (185-246 • C) and in the lower part of the aquifer (235-245 • C at > 2.75 km). The total discharge across the 80 m wide FDZ at 30 kyr is 3.44 × 10 3 m 3 yr −1 per metre length of fault, which is equivalent to 9.30 × 10 6 m 3 yr −1 for the entire 2.7 km length of the model. Water from the basement represents 97 % of this discharge, as the influx of seawater from the upper boundary through the matrix is limited by the low permeability of the caprock (with a Darcy velocity for inflow of only 0.01-0.04 m yr −1 ) and contributes < 25 % of fluid within the matrix (Fig. 2b, c). The fault-parallel 2D model shows convection within the FDZ, with flow rates approximately 1 order of magnitude slower those within the FDZ in the fault-perpendicular 2D model (Fig. 2e, f). Several narrow convection cells develop which draw cold seawater downward into the FDZ across a broad zone. This circulates to the base of the aquifer, with Darcy velocities for the descending limb of up to 6 m yr −1 , resulting in widespread cooling within the FDZ, and in the ascending limb of up to 16.5 m yr −1 at temperatures of up to 72 • C (Fig. 2e, f). The total flow of seawater drawn into the FDZ from the sea along the 2.7 km length of the model is 3.14 × 10 5 m 3 yr −1 , with only 15 % of this being fluid from the basement. This 2D simulation ignores fluid exchange with the matrix but allows lateral circulation along the plane of the fault and gives much lower fluid fluxes and temperatures in the system than the simulation oriented perpendicular to the fault.
The 3D simulation clearly demonstrates the importance of considering both flow along the plane of the FDZ and exchange with the surrounding matrix (Fig. 3). Fluid flux is strongly focused within the FDZ which accounts for 88 % of the total fluid circulation within the system within only 3 % of the total modelled volume. Rates of fluid flow within the FDZ in the 3D simulation are almost 1 order of magnitude higher than those of the fault-parallel 2D case. Darcy velocities reach a maximum of 280 m yr −1 within the FDZ at a depth of 1.2-1.5 km and decrease with depth below this. However, temperature within the FDZ is much higher than in the fault-parallel 2D simulation despite the 3D model predicting an influx of cold seawater within the FDZ 5 times higher (1.61 × 10 6 m 3 yr −1 ) than that in the fault-parallel 2D model. This apparent contradiction arises from the inability of the 2D fault-parallel model to account for advective transfer of heat between the FDZ and the matrix. In the 3D model, seawater at 25 • C that is drawn into the FDZ results in cooling that extends to the top of the basement and laterally into the surrounding matrix at a depth of 2.5-4.0 km. High temperature at the base of the model leads to upward flow of fluid through the basement both within the FDZ and, at lower rates, throughout the aquifer (with a Darcy velocity of up to 33 and 18 m yr −1 respectively). The lower permeability of the caprock forces ascending fluid within the aquifer matrix back into the FDZ at depths shallower than 2.5 km, with a total influx from the upper aquifer to the FDZ of 1.11 × 10 7 m 3 yr −1 at temperatures that range from 115 to 222 • C. This leads to a major contribution of basement water in mixed fluid circulated in most parts of the system, with the exception of the downwelling core (Fig. 3).

Baseline 2D and 3D reactive transport simulations of dolomitisation
The low convective flux of Mg 2+ -rich fluid (seawater) means that, despite the high temperature (> 185 • C), dolomitisation in the 2D, fault-perpendicular model is very limited, with replacement of < 0.05 volume fraction of the calcite within 30 kyr. Minor dolomitisation occurs where surface-derived fluids enter the reactive aquifer ("D1" in Fig. 2d) and also adjacent to the FDZ at 2.3-3 km depth ("D2" in Fig. 2d). For the 2D fault-parallel model, an even smaller amount (< 0.005 volume fraction) of dolomite forms at the base of the ascending limbs of the free-convection cells close to the basement. Although there is a high convective flux of Mg 2+ -rich fluid, the associated cooling (< 72 • C) means that dolomitisation is limited by sluggish kinetics. This illustrates the importance of having both a sufficient flux of Mg 2+ and temperatures hot enough to overcome kinetic barriers to dolomitisation. The 3D simulations show how a free-convection system developed within the FDZ of a simple planar fault could result in the flux of Mg 2+ -rich seawater across strong thermal gradients and drive dolomitisation. Reactions are largely confined to within the plane of the fault (Fig. 3), although convection also drives some exchange of fluids between the FDZ and the matrix. Flow rates within the matrix are 1-2 orders of magnitude less than those within the FDZ, yet flow from the upper part of the aquifer into the fault (at 1-2.1 km depth) is critical to maintaining temperatures within the FDZ that are kinetically favourable to dolomitisation. This "top-down" circulation of seawater drives the replacement of calcite by dolomite, reflected in a steep decline in Mg 2+ along the flow path in areas of active dolomitisation within the FDZ, mirrored by an increase in Ca 2+ (Fig. 4). The activity of these ions is a complex function of the effect of water-rock interaction and that of mixing between seawater and basement fluid.
Replacement of calcite by dolomite leads to an increase in porosity of up to 7 %, controlled by the difference in molar volume between dolomite and calcite (Fig. 5). Dolomitisation initiates at a depth of 1-2 km in the FDZ, either side of the descending limb of the convection cell, where the temperature exceeds 75 • C. In these areas the precursor calcite is completely replaced by dolomite within 9 kyr, with a commensurate increase in porosity. Replacement dolomitisation occurs as a gradient reaction, driven by Mg 2+ -rich fluids flowing from cooler to warmer areas. There is minimal alteration of the surrounding matrix, even at very high temperatures (> 200 • C in the rising limb of the convection cell), because there is insufficient flux of Mg 2+ -rich fluid. The dolomite front extends laterally over time, reaching a limit of 1 km from the centre of the downwelling limb. At depths > 2 km within the descending limb of the convection cell, dolomitisation initially occurs at a slower rate during the period when dolomitisation "up-flow" at shallower depths consumes Mg 2+ . After 15 kyr, the fluids circulating to depth are less depleted in Mg 2+ , and complete dolomitisation gradually extends to the base of the aquifer. In addition, there is late-stage replacement of calcite in the core of the downwelling limb where low temperatures initially limited reaction rate. Figure 6 illustrates the importance of temperature and fluid flux in controlling patterns and rates of dolomitisation at contrasting locations. There are sharp lateral differences in dolomite abundance within the shallow FDZ (1-2 km depth), with a change from 8 % to 70 % dolomite between adjacent grid cells (150 m model resolution). Where dolomite first develops, temperatures are > 100 • C (Location B), and dolomitisation occurs at a rate of up to 10 % kyr −1 at simulation times of 3-10 kyr following an initial 1 kyr induction period of slower dolomitisation. The temperature is slightly higher 0.6 km further along the flow path at Location C, but here the dolomitisation rate is no more than 1.4 % kyr −1 within the entire simulation period, reflecting the relatively low Mg 2+flux. Dolomitisation occurs most rapidly at Location B, reflecting the ample supply of Mg 2+ from seawater, paired with a sufficiently high temperature (Fig. 5). Dolomitisation within the downwelling zone at the same depth (Location A) is very limited for the first 13 kyr, and only initiates after temperatures start to increase in response to the influx of hot waters from the matrix (which leads to a doubling of heat flux from 13 to 30 kyr). Dolomite abundance decreases more gradually with depth (Fig. 6), and further insights into the behaviour of the system are provided by comparing locations at different depths in the downwelling core. Here the higher temperature at depth results in more dolomite forming in the early stages of the simulation (compare locations A, D, and E). However, dolomitisation occurs more slowly at depth, limited by the lower flux of Mg 2+ .

Role of porosity-permeability feedbacks in controlling patterns of dolomitisation
Most dolomite forms by replacement of calcite, with the associated increase in porosity moderated only slightly by precipitation of a minor amount of dolomite cement. Once all calcite has been replaced, dolomite formation continues but at a very slow rate (0.015 % kyr −1 ). Dolomite cement represents a minor percentage (< 8 %) of the total dolomite at shallow depth (locations A and B) and decreases with depth (to 5 % at Location E) due to the lower Mg 2+ -flux. The baseline 3D simulation incorporates feedbacks between diagenetically controlled changes in porosity, estimated permeability, heat and solute transport, and reactions. An increase in porosity of 5.4 % which results from the complete replacement of calcite by dolomite increases the permeability by a factor of 2. One impact of this in the upper aquifer is that the laterstage dolomitisation occurs within the core of the descending limb of the convection cells due to the increased transfer of heat from the matrix. In contrast, in a 3D simulation with no porosity-permeability feedback (Fig. 7) the distribution of temperature and fluid flux do not evolve with time. The total discharge of hot fluid from the FDZ to the ocean floor amounts to 7.63×10 6 m 3 yr −1 (reduced by 9 % relative to the baseline simulation with the porosity-permeability feedback incorporated). In absolute terms, the volume of recirculated seawater is higher, but this is outweighed by the reduction in fluids from the basement which represent 70 % of the total discharge (compared with 81 % in the baseline simulation). The result is that the FDZ in the shallow aquifer remains un-dolomitised within the descending limb of the convection cell in the simulation in which there is no increase in permeability driven by replacement dolomitisation (or reduction due to dolomite cements).

Sensitivity of dolomitisation by geothermal convection of seawater to basement and caprock properties
Reducing permeability of the basement from 10 −15 m 2 (≈ 1 mD) to 5 × 10 −16 m 2 (≈ 0.5 mD) and 1 × 10 −16 m 2 (≈ 0.1 mD) significantly reduces the flux of both fluid and heat from the basement, with a drop in the maximum Darcy velocity within the FDZ from 280 m yr −1 to 155 and 64 m yr −1 respectively. This results in a rather broader zone of descending flow and a more focused discharge zone (Fig. 8). The total discharge of hot fluid from the FDZ to the ocean floor is 3.64×10 6 and 8.28×10 5 m 3 yr −1 respectively, or 44 % and 10 % that of the baseline simulation. Halving the permeability of basement, the influx of seawater into the FDZ is 16 % higher, but there is a substantial (72 %) decrease in total volumetric fluid flux from the basement (Fig. 8). For a basement that is 10 times less permeable, there is a 58 % reduction in the seawater descending into the FDZ and a 97 % reduction in the basal fluid flux. Temperatures in the fault plane are lower due to reduced heat transfer from the basement to the aquifer and, subsequently, from the matrix to the fault zone. The temperature within the FDZ is consequently lower, with a broad (600 m) zone of cool (< 80 • C) waters extending to the base of the aquifer and laterally in the FDZ within the lower part of the aquifer. There is also a reduced heat transfer from the basement into the aquifer matrix and, subsequently, to the FDZ in the upper part of the aquifer. The result of these differences is a considerably smaller volume of dolomite in the FDZ compared with the 3D baseline simulation, with only partial dolomitisation below 1750 m depth and a less sharp dolomite front in the zone of ascending flow than seen in the baseline model. There is also a marked difference in dolomite distribution, with dolomite forming in the ascending limb of the convection cell where the tempera-ture exceeds 100 • C; however, as in all 3D simulations, reactions are constrained to within the FDZ. Dolomitisation occurs at temperatures lower than in the baseline simulation, ranging from 155 • C down to a minimum of 73 • C where dolomitisation rates are low. With a 90 % lower permeability basement, the very low basal heat flux and temperature significantly limit dolomitisation within FDZ (Fig. 8c).
Reducing the temperature at the base of the model from 250 • C (baseline model) to 200 • C decreases the flux of both fluid and heat from the basement (Fig. 8d). As in the lower-permeability basement models, descending flow occurs within a broader zone and discharge is more focused than in the baseline simulation. The total discharge of hot fluid from the FDZ to the ocean floor is 55 % that of the baseline simulation (4.56×10 6 m 3 yr −3 ), with a reduction in total volumetric fluid flux from the basement of 54 % (Fig. 9). This reduction in heat transfer from the basement to the aquifer and, subsequently, from the matrix to the FDZ results in a cooler FDZ. Dolomitisation within the FDZ is very limited below 75 • C, generating an un-dolomitised core in the zone of descending seawater, and rarely exceeds 50 % dolomite below 2000 km depth. The rate of dolomitisation is lower and dolomite fronts are more diffuse both laterally and vertically along the fault plane. Whilst the basement properties influenced dolomitisation via control on the transfer of both heat and fluid, the sensitivity of the system to caprock permeability seems negligible. There is no significant effect on fluid circulation, temperature, or dolomite distribution by reducing the permeability of caprock from 10 −15 m 2 in the baseline model to 5 × 10 −16 m 2 or 1 × 10 −16 m 2 (Figs. 9, 10).

Discussion
It has been more than 35 years since Land (1985) recognised that seawater is the only widely available fluid with the capacity to form extensive bodies of replacement dolomite, although reaction kinetics can limit the rate of dolomitisation at near-surface temperatures (Arvidson and Mackenzie, 1999). Recent geochemical studies (Ryb and Eiler, 2018;Shalev et al., 2019) have presented isotopic data consistent with largescale dolomitisation by convection of seawater through carbonate sediments. In contrast, a major focus of sedimentological research over the last decade has been to understand the genesis of dolostones developed along extensional or transtensional fault systems. Many of these dolomites form at high temperatures and, thus, are generally assumed to be formed from fluids that originate at depth. Release and upward escape of fluids from deep aquifers has become the dominant paradigm (Davies and Smith, 2006), despite considerable debate over the physical and chemical basis for this model and the extent to which it can account for dolomitisation of surrounding limestone (Machel, 2004;Machel and Lonnee, 2002;Robertson et al., 2015).

Comparing RTM predictions with fault-related dolomites
Whilst recognising the simplifications inherent in this narrow subset of possible 3D models, our aim is simply to evaluate the potential of the new conceptual model for fault-related high-temperature dolomitisation and capture the essence of such a system. Thus, it is important to compare the simulation results with available data describing these dolostones at outcrops. Data extracted from RTM simulations include the dimensions and position of dolomite geobodies within the permeability field, the extent to which alteration extends beyond the FDZ into the matrix, and, less commonly, the nature of reaction fronts. Outcrops provide a snapshot of dolomitisation at the time that reactions ceased, but models can provide time series maps of the development of simulated dia- genetic geobodies. Rarely explored is the potential of RTM simulations to provide metrics that relate to the properties of the dolomitising fluids (temperature, salinity, and elemental concentrations and ratios) and can be directly compared with spatial variations in the geochemical and isotopic characteristics of the dolomites. In this example, we have illustrated the further insight offered from tracking changes in temperature through the progress of the reaction at specific locations, and the difficulties in interpreting evidence of such changes from the rock record alone.

Dolomite geobody dimensions
For a FDZ with vertical and along-fault permeability 2 orders of magnitude greater than that of the country rock, our simulations suggest dolomitisation would be constrained almost entirely within the plane of the fault. In the matrix, reactions are limited not by temperature, but by the low flux of Mg 2+ . Note that, sensu stricto, dolomite in the FDZ could be classified as hydrofrigid rather than hydrothermal, as the temperatures at which it forms are lower -rather than higher -than that of the country rock (Machel and Lonnee, 2002). As the free-convection system draws cool seawater down from the surface, the transfer of heat from the matrix to the FDZ contributes to dolomitisation. Here, we do not investigate systems with sufficient matrix permeability to give the potential for significant fluid flow between the matrix and the FDZ, and the rates of diffusion are insufficient to allow for dolomitisation. However, even in such scenarios, with a transmissive FDZ providing direct connectivity between overlying seawater and the geothermal heat source at depth, the drive for dolomitisation of rocks more than a few tens of metres from the FDZ is likely limited. An exception could be where a permeable matrix may allow hydraulic continuity between two or more faults, such as across a relay ramp, allowing for the development of a single-pass convection system. The baseline simulation suggests that, for a high basal heat flux, vertically extensive dolomite bodies could form and, over time, extend laterally, merging to form a single body that might extend for > 1 km along the plane of the fault. In comparison, with a lower heat flux, dolomite bodies form at a shallower depth, and these are more limited vertically and extend only a few hundred metres along the fault. This pattern of dolomitisation focused in zones within the uppermost part of the limestones sequence beneath a barrier or baffle to flow is consistent with observations in the subsurface by Smith et al. (2008) in the Trenton and Black River formations. Comparing our generic 3D simulations with the faultcontrolled dolostone bodies of the Hammam Faraun Fault described by Hollis et al. (2017) and Hirani et al. (2018) shows marked geometric similarities. The dolostone bodies are constrained to the FDZ of the Hammam Faraun Fault, postdating earlier strata-bound dolostones within Eocene limestone country rock. Whilst the outcrop provides limited information about the vertical extent of the dolostone bodies (which exceed the 75-80 m height of the outcrop), they extend for up to 500 m along the plane of the fault. Similarly, faultrelated dolostones in the Ranero area, northwest Spain, that extend > 1 km and are defined by sharp contact within the host Aptian-Albian limestone Shah et al., 2012).

Dolomitisation timescales
Our simulations suggest that, for the combination of heat flux and permeability geometry that has been simulated, complete replacement of limestones can occur within the FDZ within a timescale of thousands to a few tens of thousands of years (Fig. 4). The unlimited supply of Mg 2+ from seawater, coupled with a high geothermal heat flux in tectonically active or rifting settings, enables much more rapid reactions than suggested by previous RTM simulations of faultcontrolled diagenesis that propose timescales ranging from 60 kyr to a few million years, even assuming more favourable kinetics (Abarca et al., 2019;Consonni et al., 2018;Corbella et al., 2014). Our model simulates a single event lasting a few tens of thousands of years within a fault that retains its transmissivity, but many outcrop studies suggest multiple episodes of fluid flow, each with the potential to recrystallise the previous dolomite. For example, in Hammam Faraun, Sr-isotope dating suggests that the massive dolomite formed   over a 10 Myr period but likely within multiple episodes during which fault zone permeability permitted more active circulation (Hirani et al., 2018;Hollis et al., 2017). Fault transmissivity is dynamic, and the movement of large quantities of fluids within faults during earthquakes via "seismic pumping" has been linked to submarine hydrothermal mineralisation that is episodic on timescales ranging from a minutes to tens of thousands of years and is controlled by volcanic and tectonic processes (Brantut, 2020;Sibson, 2001;Sibson et al., 1975). Data on permeability changes over longer timescales are sparse (Ingebritsen and Gleeson, 2017), but the reconstruction of episodic hydrothermal flow along the Malpais Fault in Nevada suggests either intermittent fluid pulses that lasted ≈ 1 kyr or continuous flow that shifted in location every 1 kyr (Louis et al., 2019). Such timescales are not incompatible with those of dolomitisation suggested by our generic simulations.

Characterisation of dolomite fronts
Our rather simplistic simulations are not designed to evaluate the characteristics of realistic dolomite fronts within the FDZ, given the large size of the model cells both vertically and along the fault as well as the associated requirement to use effective permeabilities for cells which are known to be both highly heterogeneous and temporally variable. It would be rewarding to further explore suggestions from these simulations that at early times, and with low heat flux, large areas of the fault-bound geobody remain only partially dolomitised, particularly at depth. The model also raises important questions about contrasts between dolomite fronts in different areas of the flow system depending on whether reactions are flux or reaction rate limited. We predict that the character of the dolomite front varies between the margins that receive unmodified seawater (here termed the "up-flow" margins) and those where fluids flow out of the dolomite body ("down-flow margins"). Sharper dolomite fronts (complete alteration across a single model cell) occur in the "up-flow" dolomite front where dolomitisation is kinetically limited (Mg 2+ -rich fluids first reach the threshold temperature for dolomitisation). In comparison, the "down-flow" dolomite front tends to be broader, as rates here are flux-limited as a result of the depletion of Mg 2+ by upstream dolomitisation. Reaction fronts perpendicular to the plane of the fault, where the model resolution is also higher, are within a single 10 m wide cell. This is in the same order of magnitude as those predicted by 2D simulations of HTDs formed by flow along through-going faults (Yapparova et al., 2017) as well as many field examples with sharp dolomitisation fronts (Grandia et al., 2003;Nurkhanuly and Dix, 2014;. Detailed analysis of dolomitisation fronts within fault-controlled dolomite bodies in nature are remarkably understudied, however, and are likely complex; for example, Koeshidayatullah et al. (2020b) show that reaction fronts can back-step in time, as progressive phases of recrystallisation result in porosity reduction (over-dolomitisation), restricting the extent to which subsequent fluids can flux away from the fault. An obvious benefit of flow simulations such as these is that they can provide information on the controls on dolomite body termination during the earliest phases of fluid flux, before the porosity-permeability feedback becomes an important control on the extent of fluid migration and reaction.

Linking RTM output to dolomite geochemistry and isotopic characteristics
In addition to evaluating the concept of fault-related dolomitisation by hydrothermal seawater circulation, we hope to illustrate the insights offered by a more forensic analysis of simulation data. Specifically, changes in the temperature and chemistry of fluids during the dolomitisation process can be extracted from RTM output for a more meaningful comparison between simulations and paragenetic sequences recorded in the rock record. The contrast between two simple simulations performed at higher heat flux (baseline) and lower heat flux (reduced basement permeability) demonstrates the sensitivity of geobody size and location within the convection cell, as well as the robustness of other characteristics. The temperature of dolomitisation across much of the dolomite body in the baseline simulation is relatively high and tightly constrained (120-150 • C) at depth (1-2 km). This is slightly higher but comparable to the range for the spatially more restricted dolomite bodies that form in the lower heat flux simulation (100-135 • C), reflecting the dominant role of temperature in controlling reaction rate. However, there are also contrasts between the core and the more marginal areas of the dolomite body, which are clear from the baseline simulation. Within the "up-flow" part of the dolomite body, dolomitisation initially occurs at considerably lower temperatures (80-100 • C) (e.g. locations D and E, Fig. 11), with a step change to higher temperature as the percentage of dolomite increases. Similarly, there is a sudden increase in the dolomitisation temperature in the "down-flow" margins, although here it is from ≈ 160 to ≈ 190 • C (e.g. Location C, Fig. 11). These changes occur at the same time, reflecting shifts in patterns of heat and solute transport that occur in response to diagenetic changes in porosity and permeability.
The incorporation of an unreactive tracer (Br − ) to our RTM allows us to evaluate the relative contribution of basal fluids and seawater (Figs. 2, 3). Chemical differences in the two fluids (Table 2)  areas. The association between basally derived fluids, identified from distinctive geochemical and/or isotopic signatures, and high-temperature fault-related dolomite is one commonly observed in the rock record and is generally interpreted to indicate that there is a chemical drive for dolomitisation that is associated with the input of fluid from a deep aquifer (Al-Aasm, 2003;Boni et al., 2000;Davies and Smith, 2006;Luczaj et al., 2006). The RTM clearly calls this assumption into question, as the source of Mg 2+ for dolomitisation is from the fraction of the fluid derived from the surface seawater. Instead, the basal fluid controls heat transport within the aquifer and the FDZ by advection, and thus temperature distribution. Therefore, the basal fluid plays a critical role in providing the heat required for relatively rapid formation of dolomite but is not critical to the supply of Mg 2+ for dolomitisation. The geochemical and isotopic signatures of the resulting dolomites would reflect the component of basement-derived fluids. We might expect, for instance, to see mixing trends between a 87 Sr / 86 Sr signature of downwelling fluids reflective of the seawater values, ∼ 0.7068-0.70907 (Burke et al., 1982), and more radiogenic basinderived fluids. Fluid inclusion and isotope analysis from several studies suggest that the dolomitising fluid might be formed by an evolved seawater or its interaction with rocks and/or other fluids Jacquemyn et al., 2014;López-Horgue et al., 2010).
When the porosity-permeability feedback is not implemented in the model, the temperature and basal fluid fraction remain constant through the simulation. With this feedback implemented, the evolution of thermal and chemical conditions for dolomitisation during reaction progress is a clear example of emergent behaviour, arising from the dynamic nature of the coupled system. However, because of spatial differences in the dolomitisation rate, this behaviour would be recorded at different stages within the paragenetic sequence of the dolomites at different locations in a natural dolomite. Thus, simulations help to explain complexities in defining a unique paragenetic sequence within some large, natural dolomite bodies.
We suggest that these thermal and geochemical trends are likely to be robust, irrespective of uncertainties in dolomite kinetics, although assumptions about reaction rates will affect the absolute temperatures (and, thus, basal fluid fractions) recorded in the simulated dolomites, as will differences in the composition of fluids modelled. Even relatively minor evaporative concentration of the overlying seawater will significantly lower the temperature of dolomitisation (Al-Helal et al., 2012;. Thus, clumped isotope analyses suggesting temperatures < 100 • C for the massive dolostone developed in the Hammam Faraun Fault (Hirani et al., 2018) could suggest formation from mesohaline seawater; this is not unreasonable as dolomitisation has been dated to the Miocene, when the Gulf of Suez was undergoing desiccation and evaporite precipitation (Hollis et al., 2017). Temperatures predicted from our simple models are consistent with fluid inclusion homogenisation temperatures of most dolomite formed in the Trenton and Black River groups (85-160 • C; Smith et al., 2008) and are more similar to those reported for fault-related dolomites in the Ranero area in northern Spain (120-200 • C; Shah et al., 2010) which are interpreted to form during multiphase dolomitisation, potentially at variable salinities of up to 22 wt % NaCl (López-Horgue et al., 2010). Numerous other studies combining petrographic, isotopic, and fluid inclusion analysis support the suggestion that HTDs with multiple populations of temperature and fluid composition are likely to be the norm. For example, Luczaj et al. (2006) report fluid inclusion temperatures of 120-150 • C in the HTD of the Devonian Dundee Formation of the Michigan Basin, with an average 21 • C difference between the interior and outer zone of individual dolomite crystals reflecting the episodic faultcontrolled transport of fluids and heat which is interpreted to originate from the deeper basin. Koeshidayatullah et al. (2020a) demonstrated complex geochemical proxies within the Mount Whyte Formation in the West Canada Sedimentary Basin, suggestive of mixing between seawater and deep crustal fluids, resulting in dolomite precipitated at high temperature (> 200 • C) from very high-salinity, metal-enriched fluids at shallow depths. This highlights the need for further work to evaluate the impact of changing composition of the end-member fluids, as well as the effect of solute concentration on density-driven circulation, which is driven solely by differences in fluid temperature in this simulation.

Model limitations and new insights from simulations of dolomitisation by hydrothermal convection of seawater
Systems that are governed by interactions between heat and solute transport and water-rock interaction are challenging to predict or quantify a priori (Ingebritsen et al., 2010). Nu-merical simulations that represent the coupled nature of these processes can reveal emergent behaviour that is governed by relationships between these separate components. Thus, they offer a route to evaluate conceptual models of complex systems. However, model results are inevitably affected by the chosen geometry and permeability of the system which controls solute flux and heat transport and, thus, the rate and distribution of reactions.
Simulations of structurally complex reservoirs are particularly challenging as they comprise discrete structures with properties that contrast strongly with those of the country rock and range over vastly different length scales (Matthäi et al., 2007). Variations in permeability and fault zone thickness are fundamental controls on fluid circulation and subsequent mineralisation distribution (Guillou-Frottier et al., 2020;Harcouët-Menou et al., 2009), whereas the effect on hydrothermal convection of other parameters, such as porosity, is considerably less significant (Gow et al., 2002;Kühn et al., 2006;Zhao et al., 2003). However, our models do not attempt to capture the complexities of flow within FDZs, for which the porous media formulation employed by TOUGHREACT is inappropriate. We also do not consider how the permeability of the FDZ may change due to geomechanical responses to variations in the stress field that are localised in time and space nor the physical and chemical effects of changes in state of the fluid, for example, due to boiling. Rather, we simply demonstrate that, where major discontinuities provide high-permeability pathways for flow, it is possible for interaction between heat supplied from the basement and Mg 2+ -rich seawaters from above to form hightemperature fault-bound dolomite bodies in high-enthalpy settings. The magnitude and anisotropy of permeability, as well as the effective reactivity, of the country rock will also influence the extent to which this fault-associated convection system could form dolomites in the adjacent matrix or, alternatively, might drive recrystallisation of a dolomite precursor formed during earlier diagenesis.
Although this model deals solely with seawater and an equivalent salinity basal fluid, there is a clear opportunity now to evaluate the reactions that occur within convective systems involving basal fluids of different geochemical characteristics and densities as well as to investigate alternative surface water compositions. More complex mineral assemblages should also be investigated, both in terms of hostrock mineralogy and diagenetic minerals that may form. Our simple simulations incorporating the effect of replacement dolomitisation on permeability demonstrate the importance of better understanding the impact of the full range of both physical and chemical controls on the evolving permeability. They also offer three substantive lessons with implications that extend beyond the contentious realm of genesis of HTD.
Firstly, although our generic models are very simple, they robustly demonstrate that, under the conditions simulated, high-temperature fault-controlled dolostone bodies can form by seawater convection. Mg 2+ -rich seawater is drawn down-ward into the permeable damage zone around a fault that extends to shallow depth, mixing with hot buoyant fluids escaping from the basement. Dolomitisation occurs within the FDZ as a gradient reaction, with the replacement of host limestones and minor dolomite cementation as circulating fluids flow across isotherms. Fluids from the basement play an important role in transporting heat into the FDZ, but the Mg 2+ to drive dolomitisation comes from the entrained seawater even though this is < 50 % and locally < 20 % of the fluid within the FDZ. Thus, our reactive transport simulations reconcile the apparently contradictory observations that on the one hand identify geochemical and isotopic signatures for basement fluids in many high-temperature dolomites, and on the other hand indicate that almost all basement fluids are depleted in Mg 2+ and, therefore, do not have the potential for dolomitisation.
Secondly, comparison of our 3D baseline model with simplifications in 2D, both perpendicular and parallel to the plane of the FDZ, illustrate how 2D models fail to represent critical aspects of system behaviour ( Fig. 2a-h). The alignment of convective flow with fault planes has long been recognised Smith, 1995, 1996;Murphy, 1979) and is supported by recent 3D thermal models of the Upper Rhine Graben (Armandine Les Landes et al., 2019) that show good agreement between simulated and predicted thermal anomalies. Our RTM simulations show that 2D models oriented perpendicular to the fault cannot represent the convection that occurs along the plane of the fault within the permeable damage zone, which is key to generating massive fault-constrained dolostone bodies. It is likely that prior simulations using such 2D models oriented perpendicular to the plane of the fault Frixa et al., 2014;Jones and Xiao, 2013;Xiao et al., 2013) misrepresent these important aspects of reactions driven by the circulation of fluid along the fault plane. Therefore, they might not explain the diagenetic evolution of fault-guided diagenetic reactions well, specifically overestimating the influence of fault-controlled fluid convection control on the adjacent matrix. Circulation within the FDZ is seen in the 2D model oriented along the plane of the fault (Fig. 2e-h). However, because this model fails to account for the effect of heat transfer between the matrix and the FDZ, the system appears to have been too cool for significant dolomitisation in this case. Our clear conclusion is that numerical analysis of systems where faults impact permeability must be conducted in 3D, and this has important implications for modelling a range of systems, from geothermal resources and mineralisation fluid to carbonate diagenesis, including dolomitisation and hydrothermal karstification.
Finally, dolomitisation associated with the hydrothermal and geothermal circulation of seawater has been suggested as an important but poorly constrained sink for Mg in seawater (Shalev et al., 2019). Thus, these processes have fundamental implications for understanding the global budgets of Mg, Ca, and C, which are linked by the common con-trols of weathering, volcanism, and carbonate precipitation (Arvidson et al., 2011;Elderfield, 2010;Holland, 2005). Combining recent magnesium isotope ( 26 Mg) data for lowtemperature hydrothermal fluids with evidence of a constant 26 Mg over the past 20 Myr, Shalev et al. (2019) suggest that there is a significant unexplained flux of Mg from the modern ocean (estimated at 1.5-2.9 Tmol yr −1 ). Mg mass balance calculations for our simple generic models suggest that the dolomitisation that results from this mixed convection system within permeable fault zones may provide a sink for 2.0 × 10 7 -2.6 × 10 7 mol Mg yr −1 for each kilometre of fault. At such rates, in order to account for all the cryptic dolomite identified by Mg-isotope studies, such a system would need to be established along total fault length of 60 × 10 3 -150 × 10 3 km. For comparison, the East African Rift system extends some 6.4 × 10 3 km, with extension accommodated by numerous intra-rift faults in addition to those bordering the rift (Shillington et al., 2020, and references therein). This suggests that Mg drawdown by fault-related dolomitisation in areas of very high localised heat flux may be significant, but it is likely to be subsidiary to dolomitisation driven by geothermal ("Kohout") convection within platform carbonates. The absence of permeable conduits or high basal heat flux mean much lower rates of dolomitisation in these systems (Whitaker and Xiao, 2010), but this is amply compensated for by the vast area of major carbonate platform systems globally (currently some 800 000 km 2 of lowlatitude oceans; Milliman, 1993). Further resolution of these fluxes is key to critically assessing the role of dolomitisation relative to weathering and seafloor spreading in determining secular changes in the Mg/Ca of global ocean (Hardie, 1996;Holland and Zimmermann, 2000;Shalev et al., 2019).

Conclusions
Hydrothermal dolomitisation has been synonymous with dolomitisation driven by release of over-pressured fluids from depth, but our simulations offer another potential interpretation. An alternative concept considers that hightemperature dolomite could replace a precursor limestone where seawater is drawn down into a convection system that is developed within a steep high-permeability conduit during periods of high heat flux. Such a scenario may be most likely developed in extensional basins, when faults breach the seafloor and geothermal gradients are elevated. Evaluating this concept using 2D and 3D reactive transport simulations shows that when a deep-seated permeable fault extends to the sea floor, the influx of high-Mg 2+ / Ca 2+ seawaters can drive replacement of calcite by dolomite with only minor dolomite cementation. The rate of advective transfer of heat by basally derived fluids relative to that of cold seawater from above will determine the distribution of dolomite within the flow field and the size of the dolomite geobodies that form within the convective cell. This new model provides one so-lution to the Mg 2+ mass balance problems inherent in HTD formation by release of hot over-pressured fluids from a deep reservoir. It may also explain at least some of the missing Mg 2+ from the world's oceans. The new insights and key conclusions from our study include the follows: -RTM simulations suggest that convection of seawater that leads to a mixture of hot basal fluids and seawater along the FDZ can form sizable (10 2 -10 3 m) faultbound dolomite bodies within timescales of thousands of years in high-enthalpy settings.
-Reactions are largely confined to within the FDZ, and dolomitisation occurs at temperatures of 80-160 • C and from fluids that comprise a significant (> 0.5, and locally > 0.8) fraction of basally derived fluids.
-Modification of porosity and permeability during diagenesis can lead to systematic changes in the temperature and fraction of basal fluid in the hydrothermal mixture as progress of the reaction that may be recorded in the paragenetic sequence.
-A permeable fault damage zone provides a corridor along which fluids will preferentially flow. As this is not represented in 2D simulations perpendicular to the plane of the fault such model construction is fundamentally inappropriate and could result in misleading suggestions (e.g. about the extent of the associated alteration of the country rock).
-Heat transfer between the fault zone and country rock plays an important role in contributing to dolomitisation within the FDZ which cannot be observed using 2D simulations along the plane of the fault.
3D RTM simulations can enhance understanding of the mechanisms and controls on dolomitisation within faulted and fractured systems, and the geometry and spatial distribution of the resulting high-temperature dolomite geobodies. RTMs offer insights into temporal evolution of reactions, and an opportunity to examine not only the geometry of diagenetic geobodies and the nature of reaction fronts but also a broader range of metrics that can be directly related to isotopic and geochemical characteristics of diagenetic products. Critically, numerical analysis of systems where faults impact permeability must be conducted in 3D, and this has important implications for the modelling of systems ranging from geothermal resources to ore formation and carbonate diagenesis, including hydrothermal karstification as well as dolomitisation.
Data availability. In our experiments we used TOUGHREACT, a numerical simulation program developed by the Lawrence Berkeley National Laboratory (LBNL) which is available at https://tough.lbl. gov/software/toughreact (Sonnenthal et al., 2017, last access: 5 December 2020. The output files from TOUGHREACT are visualised using trexplot, a Python script developed by Hamish Robertson which can be found at https://github.com/hammytheham/trexplot (Robertson andBenjakul, 2019, last access: 5 December 2020). The output data from all simulation models in this article are available at https://doi.org/10.5281/zenodo.4312236 , last access: 8 December 2020.
Author contributions. FFW and RB wrote the article and designed and developed the conceptual model and experiments. Model construction, simulations, and analyses were conducted by RB under the supervision of FFW. ELS provided technical support with the TOUGHREACT simulator. HAR offered advice on the initial model construction and developed a Python script, trexplot, for TOUGHREACT data visualisation. CH, ELS, and HAR gave advice, feedback, and comments on the paper.
Competing interests. The authors declare that they have no conflict of interest.