Evolution of Fault-Zone Hydromechanical Properties in Response to Different Cementation Processes

Progressive cementation and sealing of fault-localized fractures impact crustal mass transport and the recovery of fault strength following slip events. We use discrete fracture network (DFN) models to examine how fracture sealing during end-member cementation mechanisms (i.e., reactionversus transported-limited cementation) influences the partitioning of fluid flow through time. DfnWorks was used to generate randomized fracture networks parameterized with fracture orientation data compiled from field studies. Single-phase flow simulations were performed for each network over a series of timesteps, and network parameters were modified to reflect progressive fracture sealing consistent with either reactionor transport-limited crystal growth. Results show that when fracture cementation is reaction-limited, fluid flow becomes progressively channelized into a smaller number of fractures with larger apertures. When fracture cementation is transport-limited, fluid flow experiences progressive dechannelization, becoming more homogeneously distributed throughout the fracture network. These behaviors are observed regardless of the DFN parameterization, suggesting that the effect is an intrinsic component of all fracture networks subjected to the end-member cementation mechanisms. These results have first-order implications for the spatial distribution of fluid flow in fractured rocks and recovery of permeability and strength during fault/fracture healing in the immediate aftermath of fault slip.


Introduction
Brittle fractures in fault zones exert a first-order control on the hydrological and mechanical properties of the upper crust. Understanding their distribution, properties, and evolution is essential for a variety of geological applications. For example, the formation of brittle fractures results in dramatic increases in the permeability of fault zones relative to protolith [1], which may degrade the integrity of CO 2 or other waste-storage reservoirs [2][3][4][5][6]. Conversely, increased fracture-related permeability along faults may also provide pathways for the migration and accumulation of hydrocarbons [7][8][9][10]. Multikilometer fluid flow along such pathways may also be a major driver of active geothermal fields [11,12], leading in some cases to the localization of significant deposits of economic minerals (e.g., [13,14]). From a mechanical perspective, the localization of fractures within fault zones has been argued to partially mediate radiated seismic energy and long-period ground motions during earthquakes [15][16][17][18][19]. Thus, the occurrence and longevity of fault-related fractures is an active and long-lived area of research with substantial societal implications. Research in this area is complicated; however, by the fact that fault-zone fracture networks are not static structures in the crust. Rather, a variety of rock record, geophysical, and experimental studies indicate that fracture systems are dynamic, repeatedly opening and closing throughout fault history [1,15,16,[20][21][22][23][24].
One common mechanism by which fractures become closed (hereafter "sealed") is cementation, whereby minerals precipitating from fluids progressively occlude fracture porosity (i.e., vein formation). The rate of fracture cementation is controlled by the relative rates of mass transport and chemical reaction on mineral surfaces. Reaction-limited cementation describes a scenario within which the rates of mass transport in the fluid phase outpace the rates of reaction on mineral surfaces, resulting in fracture sealing rates that are invariant with respect to the rates of solute transport ( Figure 1). Thus, reaction rates and the spatial distribution of reaction products should be relatively homogeneous over length scales that do not result in significant variations in temperature/kinetics (i.e., a few 10's or perhaps 100's of meters in the crust). Reaction-limited transport in both porous media and fracture networks in the crust has been inferred by a variety of experimental, numerical, and rockrecord studies [25][26][27][28][29][30][31][32] and may be particularly applicable during cementation of fracture networks when the rates of fluid advection are modest or where the precipitation reaction is kinetically limited. Thus, in a fractured but homogeneous rock mass or in a rock mass where mineralogical variations (and associated surface chemistry heterogeneities) occur over smaller spatial scales than those of fracture length, the rate of aperture decrease and fracture sealing during cementation should be approximately constant through space and time.
Transport-limited cementation describes a scenario within which the rate of reaction on mineral surfaces is proportional to the rate of solute transport in the fluid phase ( Figure 1). Stated alternatively, transport-limited cementation occurs when a mineral reaction is not kinetically inhibited, and the reactant is effectively consumed as quickly as it is supplied. In transport-limited systems, the rates of reaction and thus the spatial distribution of reaction products are often heterogeneous and fundamentally reflective of variations in the rates of solute transport [33]. This effect was demonstrated in a fault-localized fracture networks by Williams et al. [34], who used geochronology of syntectonic calcite veins to demonstrate that fracture cementation rates varied by over an order of magnitude (from 50 to 800 μm/kyr) during fluid advection in the Loma Blanca fault zone, New Mexico, USA. These veins were localized over an along-strike distance of less than~40 m in the fault damage zone. The estimated rates of cementation and sealing were shown to be directly proportional to width of the examined cement-growth increments, which for continuous growth histories provides a minimum estimate of the fracture aperture at the time of cementation (i.e., fractures with larger apertures experience faster rates of calcite growth and sealing). Assuming that (1) the transmissivity of individual fractures scales as a cube of their aperture [35] and (2) and fluid flow preferentially localizes within the most transmissive fractures in an interconnected network (e.g., "flow channeling," [36][37][38][39][40]). Williams et al. [34] interpreted these observations to record transport-limited fracture cementation during cementation of the Loma Blanca fault. They also proposed a general conceptual model, wherein if solute transport is advection-dominated, and the fractures with the largest apertures (i.e., most transmissive fractures) maintain higher flow rates; then, the rates of aperture decrease, and sealing will be variable throughout the network and proportional to the aperture width. In this way, fractures with larger apertures would be expected to seal more quickly than smaller fractures, an inference that contrasts starkly with previous studies examining cementation in reaction-limited systems [26,41]. These inferences have been demonstrated repeatedly for dissolution processes occurring within porous media [30,[42][43][44]. Similar outcomes are expected for precipitation reactions [31].
In summary, reaction-and transport-limited cementation may be considered end-member regimes contributing to sealing of fault-zone fracture networks. The effects that these end-member regimes have on the evolution of network properties, however, are largely unknown.
In this paper, we employ a computationally simple modeling strategy to quantify the effects of reaction-limited (constant precipitation rate) and transportation-limited (simulated here as an aperture-dependent precipitation rate) mineral growth and cementation in discrete fracture network models parameterized with geometric data from real fault zones. We find that these two end-member mechanisms of mineral precipitation have opposite effects on the distribution of fluid flow and cementation in fault-zone fracture networks. Specifically, when mineral growth in fractures is reaction limited, the locus of fluid flow and cementation becomes increasingly channelized within discrete segments of the fracture networks through time (e.g., [45]). When mineral growth is transport limited, the locus of fluid flow and cementation becomes increasingly homogeneous and widespread in the fracture networks through time. These results have important implications for the evolution of hydrological and mechanical properties of fault damage zones and, more generally, of fractured media in the brittle crust.

Methods
In the last two decades, significant progress in coupling fluid flow, transport, stress, and chemical reaction equations has been achieved by thermal-hydrologic-mechanical-chemical modeling techniques (THMC modeling, [46][47][48][49][50]). These models combine of a variety of different physical modeling approaches with finite-element mesh generation to simulate complex processes occurring at the continuum scale in dualporosity media. THMC models have been successfully used for several multi-and single-phase applications in fractured and porous media, such as CO 2 and radioactive waste storage and geothermal energy production [51][52][53][54]. The complexity of THMC models, however, also results in several important limitations. For example, THMC models are computationally intensive, which often necessitates their application to individual discrete fractures rather than networks of fractures (e.g., [55]). The results of models   Lithosphere considering the effects of mineral precipitation and/or dissolution are also highly dependent on the specific geochemistry (i.e., fluid and rock composition, temperature, pressure, and hydraulic gradient) of the modeled system. Recognizing these limitations, we opted for a simplified approach to quantifying the first-order, network-scale effects of reactionand transport-limited cementation in fault-zone fracture networks in the most general sense. We use dfnWorks to simulate the effects of reaction-and transported-limited cementation on the hydrologic properties of fracture network during sealing. dfnWorks is an open source tool developed by Los Alamos National Lab that generates discrete fracture network (DFN) models and simulates fluid flow and mass transport. The workflow is fully described by Hyman et al. [56]. A three-dimensional DFN is generated based on specified geometric parameters (e.g., fracture intensity, orientation, length, and aperture) using the feature rejection methodology described by Hyman et al. [57]. These initial DFNs are then converted to a high-resolution computational mesh using the parallelprocessing code LaGriT [58]. High-resolution computational meshes are then converted to an input mesh for the flow solver PFLOTRAN [59], which computes a steady-state pressure solution and associated Darcy flow velocities. We extend the dfnWorks approach described above to allow fracture apertures (and thereby fracture permeability) to be updated between model timesteps, thus simulating the effects of progressive fracture sealing in an otherwise static network (i.e., no new fractures are formed between timesteps, nor are the original fractures permitted to growth in length or aperture). For each analyzed fracture network, aperture values were updated using either a constant (i.e., simulating the reaction-limited case) or aperture-dependent (i.e., simulating the transport-limited case, discussed more below) rate between timesteps.
We argue that our modeling approach, although highly simplified relative to the more common THMC approach, is still capable to capturing the first-order, network-scale effects of reaction-and transport-limited cementation processes in complex and densely populated fracture networks. Moreover, we highlight several advantages of our methodology when compared to the THMC modeling: (i) It is independent of the exact geochemistry of the system. Quartz, carbonates, oxides, and phyllosilicates are relatively common fracture cementing phases in fault zones [34,[60][61][62][63][64][65][66]. Each of these precipitation reactions, however, has different kinetic parameters and therefore geochemical behaviors. Our approach is nonspecific with respect to the mineralogy of the cementing phase and instead describes network-scale flow associated with cementation in the most general sense (ii) It makes no assumption nor predictions about the exact timescales of the network's modification and sealing, which are subject to large uncertainties. Fully coupled models specific to a particular cementing phase, for example, have to rely on laboratory-derived kinetics, which frequently overpredict the rate of many geochemical reactions that occur in nature (e.g., [67][68][69][70]) (iii) It allows us to use empirical, rock-record estimates as an exemplar form of the relationship between fracture aperture and sealing rate in a natural, transport-limited fault system (e.g., [34]) (iv) It focuses on easily quantified, network-scale fluid partitioning effects that occur within realistic fracture networks in a computationally efficient way 2.1. Network Parameterization. To ensure that the orientation of fractures within our DFNs is realistic, we utilize data from three previous publications that quantified fracture orientations in fault damage zones: the Stillwater normal fault (Dixie valley, Nevada; [71]), the Hammam Faraun normal fault (Suez Rift, Egypt; [72]), and the oblique-slip Alpine fault (New Zealand; [73]). Fracture orientation data from the Stillwater and Hamman Faraun normal faults were collected in field studies, whereas the data from the Alpine fault are the product of borehole geophysical logging. Although many previous papers have described the average orientation of fault-zone fractures, relatively few provide the large, detailed data sets available within the sources above. Thus, the data sets were selected primarily on the basis of their availability. Notably, the data sets we use are purely geometric and do not include sequential fracture opening through time as might be established via crosscutting relationships or vein microstructures, and as such our models cannot simulate mechanical interactions between fracture formation/growth and cementation (we refer the reader to a recent review which provides a detailed description of these issues [27]). The adoption of these orientations to parameterize our DFNs effectively assumes that all of the measured fractures were active simultaneously, perhaps after a singular slip event. We recognize that this assumption is simplistic with respect to geological realities, and that many of the measured fractures likely formed episodically throughout the slip history of each fault. Here, we take the employed geometries as simply representative of the stress fields and kinematics of their formation, which provide a basis for parameterizing only the orientation of fractures within our DFNs. Parameterization of our DFNs began with identifying fracture sets with distinct orientations within stereonet data from our three exemplar fault zones (Figure 2(c)). The orientation and dispersion of each set were estimated and used for input to dfnGen, the initial fracture network generation routine of dfnWorks. This process was repeated iteratively until dfnGen consistently produced randomized DFNs (Figure 2(a)) whose fracture orientations were qualitatively similar to those of the field data (Figures 2(b) and 2(c)). The details of geometric parameters for each fault are included in Table S1 of the Supplementary material document. Discrete fractures in each network were modeled as rectangles with a lognormal distribution of 3 Lithosphere lengths. The fracture intensity of each network, commonly described as p32, is defined as follows: where S f is the surface area occupied by the fractures and V is the total volume of the domain. Generally, a larger fracture intensity results in a larger number of discrete fractures when other parameters are held constant, although in some circumstances, a smaller number of fractures with larger lengths/surface areas may yield similar results. We note that all of the DFNs produced in this way are unique, being the product of random generation given a fixed set of "target" parameters (i.e., fracture orientation, dispersion, length distribution, and intensity).
A preponderance of field studies has confirmed the fundamental scaling relationship between the aperture and length of isolated fractures [41,[74][75][76][77], consistent with theoretical predictions based on linear-elastic fracture mechanics [78][79][80]. Thus, the initial aperture of each fracture is determined from its length by a power-law relationship: where a is the aperture, l is the half length of the fracture, and c and b are constants. Here, a value of 0.006 is assigned 4 Lithosphere to c, and a value of 0.5 is assigned to b to target aperture values ranging from 3 mm to 15 mm. This aperture range was targeted to coincide with that documented at the Loma Blanca fault [34,81], from which our aperture-dependent sealing rate relationship was derived (see next section). Fracture permeability (k) is then estimated as function of the fracture aperture (a), in accordance with the "parallel plate" model commonly used in fluid-flow modeling [35][36][37]82]: The rock matrix is treated as impermeable media in our models.

Simulating Reaction-and Transport-Limited Fracture
Cementation. Single-phase fluid flow is simulated imposing a pressure gradient from West to East across the model domain. No flow conditions characterize the front, back, top, and bottom faces of the domains. We simulated the effects of reaction-and transport-limited fracture cementation by modifying the aperture (and thereby permeability) of individual fractures in our randomized networks between flow simulations over a series of discrete timesteps. The initial state of each network is determined by the methods described in the previous section. In subsequent timesteps, we simulate the effects of fracture cementation by reducing aperture values between fluid-flow simulations. In models simulating reaction-limited cementation, fracture apertures are decreased at a constant rate at each timestep. In this way, the rate of closure of each fracture is independent of its aperture. To maintain a constant network geometry, the amount of aperture "loss" per timestep is imposed such that the fracture with the smallest aperture in the network remains open following the last timestep. In models simulating transport-limited cementation, we employ an aperturedependent closure rate of the same form documented by Williams et al. [34] for the Loma Blanca fault.
where r is the rate in mm/kyr (thousands of years) and a is the aperture in millimeters. Thus, the amount of aperture "loss" in each fracture depends on the width of its aperture in the previous timestep, where fractures with larger apertures seal faster. We use the aperture-dependent rate from Williams et al. [34] here as it provides us with a rate relationship derived from a natural system, and we are unaware of any similar data which relate sealing rate to fluid-flow velocity (which would more accurately simulate the transport-limited case, discussed more in Section 4.1 below). When applied to millimeter-scale apertures typical in our generated DFNs, Equation (4) implies simulation timescales on the order of several thousand years, consistent with that documented for the Loma Blanca fault. We note, however, that the equation is meant only to provide a generalized form of aperture-dependent cementation during transportlimited sealing. We make no argument here regarding the absolute timescales of fracture cementation, which are almost certainly variable between faults and indeed within individual faults as a function of depth, lithology, temperature, local geochemistry, and prevailing hydraulic gradients. Moreover, we recognize the potential for mechanical interaction between progressive cementation and ongoing deformation as outlined in Laubach et al. [27] and elsewhere.

Quantification of Network-Level Effects during
Cementation. We employ a "fracture participation" ratio to quantify the effects of reaction-and transport-limited cementation on network-scale fluid-flow properties during our simulations. Fracture participation (or flow channeling density indicator as defined in Maillot et al. [83]) is defined as the ratio between the "active" p32 and the overall p32 (Equation (1)) of the network during flow [84].
where fp is the dimensionless fracture participation ratio, S is the surface area of each fracture (m 2 ), and Q is the flow rate in each fracture (m 3 /s). Given a constant network geometry, decreasing fracture participation implies that fewer fractures are participating in fluid flow, and therefore, that flow is more channelized. Conversely, increasing fracture participation implies that more fractures are participating in fluid flow, and therefore, that flow is more homogeneously distributed throughout the fracture network.

Results
We quantify changes in fracture participation in model runs simulating reaction-and transport-limited cementation in randomized representations of our exemplar faults over equal model durations. These observations allow us to identify a striking difference between the two end-member mechanisms of fracture closure. When the rate of fracture closure is constant (simulating reaction-limited cementation), fracture participation decreases over time (Figure 3(a)). When the rate of fracture closure is dependent on the width of the aperture in the previous timestep (simulating transport-limited cementation), fracture participation increases over time (Figure 3(b)). Thus, fluid flow in reaction-limited simulations becomes more channelized over time, whereas fluid flow in transport-limited simulations becomes more homogeneously distributed over time. Moreover, we observe similarly striking differences in the fracture aperture distributions prior to and following simulated cementation (Figure 3(c)). In simulations of reactionlimited cementation, the fracture aperture distributions shift uniformly to smaller values, reflecting the constant rate of cementation regardless of initial aperture width. In transport-limited cementations, we observe that fracture apertures effectively converge on a single value with time. These trends are observed in models representing each of our exemplar faults, although the magnitude of the effect is variable between them. We hypothesize that these effects are intrinsic to reaction-and transport-limited cementation 5 Lithosphere in fractured rocks, rather than emergent properties of the specific network geometries we examined. To test this hypothesis and examine the effect of specific parameters on overall flow behavior, we repeat our analysis while systematically varying individual parameters of our randomized fracture networks.

Reproducibility.
To assess the reproducibility of our initial observations, we constructed four additional network permutations for each of the three case studies. These permutations used identical input parameters to the original case shown in Figure 3 but utilized different random seed numbers, resulting in new fracture networks satisfying the original input parameters. Stated alternatively, the distribution of fracture orientations is similar between the permutations, but the position, length, and aperture of individual fractures change. These permutation tests yielded effectively identical behavior to that observed for our original model (Figure 4), demonstrating that the difference in fracture participation trends between reaction-and transport-limited simulations is a reproducible, intrinsic property of the examined network geometries when exposed to the end-member cementation mechanisms.

Fracture Orientation.
Although fracture permeability is known to be a major control on fluid-flow partitioning in fracture networks, the orientation of individual fractures relative to the direction of the maximum hydraulic gradient is also important [85][86][87]. We investigate the potential for fracture orientation effects on fluid flow in our simulation in two ways. First, we conducted several trials where the dispersion in fracture orientation was sufficiently large as to yield randomly oriented fractures ( Figure 5). These trials indicate that the differences in fracture participation trends between models simulating reaction-and transportlimited cementation persist even in highly irregular fracture networks.
Second, we plot the fracture orientation in our fieldparameterized models ( Figure 2) as a function of the relative Darcy flux resulting from simulations shown in Figure 3 at 6 Lithosphere the initial timestep. In Figures 6(a)-6(c), the average of the aperture values (in millimeters) for each fracture set is also reported as superimposed numbers. The results show that the relative Darcy flux is largely independent of fracture orientation in our simulations. Rather, the highest fluxes appear to be correlated with larger aperture fractures ( Figure 6(d)). This correlation is likely induced by the modalities in which the networks are built and how the aperture and permeability values are defined. As mentioned in Section 2.1, permeability for each fracture is directly proportional to the fracture aperture. Fractures with larger aperture are characterized by higher permeabilities and a higher Darcy flux (m/s) in our simulations. Thus, it is unlikely that the observed differences in fracture participation trends between models simulating reaction-and transport-limited cementation can be explained by emergent properties associated with orientation effects in a fracture network.

Fracture Intensity.
We assess the effect of changes in fracture intensity (p32, Equation (1)) on the observed trends in fracture participation by systematically varying the target intensity during network generation and rerunning our flow simulations. Fracture intensity was increased or decreased by an equal factor for each set in each test case between flow simulations while other factors (e.g., size of the model domain as well as the minimum, maximum, and standard deviation of the fracture length distribution) were held constant. The results of this test show that the overall trends in fracture participation between models simulating reactionand transport-limited cementation persist regardless of the fracture intensity of the fracture network ( Figure 7). The results also show that the magnitude of the deviation in fracture participation between simulating reaction-and transport-limited simulations increases with increasing fracture intensity values. We attribute this effect to the fact that the networks generally contain more fractures as the fracture intensity is increased, leading to increased potential for changes in flow channelization relative to the initial state.
3.4. Initial Aperture Distribution. We assess the potential effects of the initial distribution of apertures in the fracture networks by systematically varying the distribution parameters while holding all other factors constant. To achieve different aperture distributions without changing the distribution of fracture length or fracture intensity, we modify the scaling constants (b and c) that describe the relationship between fracture aperture and length (Equation (2)). We accomplish this by first varying the parameter b while c is    Lithosphere and transport-limited cementation simulations increases. These observations are most readily explained by the aperture-dependent nature of fracture permeability. The contribution of variations in fracture aperture to flow channelization is minor when those apertures are relatively uniform between individual fractures. Changes in fracture participation through time are similarly suppressed in models simulating transport-limited cementation, as more uniform aperture distributions lead to uniform sealing rates between fractures (Equation (4)).
We similarly examine the effect of systematic changes in the c parameter of the aperture-length scaling (Equation (2)) in our simulated fractures by holding the b value constant. As the c value is increased, both the mean and variation of the initial aperture distribution increase (Figure 9). The effect of these changes on initial fracture participation and its change through time in our models, however, is negligible. Only a slight difference in the magnitude of decrease in fracture participation during reaction-limited simulations is observed, where larger c values lead to smaller changes in fracture participation. We attribute this effect to the constant rate of aperture decrease in reaction-limited simulations and the increasing mean aperture of the initial distributions with increasing c value; the percent occlusion and permeability change (and therefore changes in flow channelization) experienced by each fracture in reaction-limited simulations is smaller where initial apertures are larger.

Micron-Scale Fracture
Apertures. The millimeter-scale fracture apertures utilized in the models occur locally within the upper most brittle crust (e.g., <10 km [81,[88][89][90][91]). At greater depths, however, fracture apertures are often submillimeter in scale (e.g., [92][93][94]), although exceptionally large fractures may occur [95,96]. To assess the potential portability of our results to fluid flow and cementation in fracture networks at greater depth in the crust, we conducted a subset of simulations employing smaller aperture values (50-200 microns). All other parameters were kept identical to the original case shown in Figure 3. Even at these smaller aperture values, we observe the same general deviation in fracture participation between models simulating reactionand transport-limited cementation (see Figure S1 in the Supplementary material document). Thus, we tentatively propose that our results may also be applicable to fluid flow and fracture sealing in faults deep within the brittle crust, although we stress that this interpretation should be treated with caution, and additional research will be required to fully validate this claim.
3.6. Flow Direction. As a final test, we also ran our simulations employing a vertical hydraulic gradient (from bottom 9 Lithosphere to top of the cubic domain rather than the original West to East). We observe the same general behavior as documented above in models simulating both reaction-and transportlimited cementation ( Figure S2 in the Supplementary material document).

Model Limitations.
In this study, we simulate reactionand transport-limited cementation processes by imposing a constant or aperture-dependent fracture sealing rate, respectively. In real systems, however, fracture sealing rates in a transport-limited regime would be proportional to rates of fluid/mass flux, which for a given hydraulic gradient are dependent on fracture aperture, orientation, and degree of interconnectivity within the network (which may be influenced by length). Ultimately there are no empirical data available that describe the rates of calcite (or any other phase) precipitation as a function of flow rate in natural systems. Hence, it is reasonable to utilize aperture-dependent growth rates derived from geochronology of calcite veins [34] to simulate a transport-limited cementation process in our models. Although this approach neglects the potential role(s) of flow channeling due to fracture orientation in real systems, we argue that the effects on our results are limited. As shown in Figure 6, fracture apertures are the dominant control on flow rates at model initiation. The eventual plateau in fracture participation observed in all models simulating transport-limited cementation at later timesteps is a result of progressive aperture normalization (Figure 3(c)). It is conceivable, however, that fracture participation in real systems will continue to evolve following normalization of fracture apertures when variations in orientation become the dominant control on rates of fluid flow. Alternatively, formation of new fractures (or growth of old fractures) may also exert a significant effect on the evolving properties of the fracture network (e.g., [27]). Thus, we make no specific comment on the rate or magnitude of changes in fracture participation in natural systems based on our models. We contend, however, that simplifications in our models relative to natural systems have no impact on the more general results we demonstrate: fluid flow can be expected to become less channelized through time in systems experiencing transport-limited cementation and sealing, but more channelized in systems experiencing reaction-limited cementation and sealing. As such, understanding the geochemical character of fracture systems is an important component of predicting the nature of future changes in fluid flow through fracture networks.

Implications for the Hydrological and Mechanical
Properties of Fault Zones. This study demonstrates that  10 Lithosphere variations in the mechanism of fracture cementation strongly affect the properties of evolving fracture networks.
In general, quantitative works simulating fluid flow through fracture networks assume that fracture properties (e.g., aperture, length, intersections, and abundance) are constant over time [97][98][99][100]. This assumption is likely only valid over very short timescales or in settings where the geochemical component of fluid-rock interaction can be neglected. Over longer timescales or in reaction environments, however, both dissolution and precipitation of mineral phases occur, and the configuration of fracture networks will necessarily change over time. Abundant veins observed in exhumed fault zones and other fractured rocks provide a clear example of coupled geochemical and transport processes that drastically modify fracture properties [23,93,[101][102][103][104].
The abundance of veins relative to dissolution features in exhumed systems may arguably be considered evidence that precipitation is a more common process affecting fracture   (2)), while c value is held constant (=0.006) for all the simulations. (b) Fracture participation behavior calculated for each model at each timestep. If the aperture distribution is narrow, initial fracture participation at timestep 0 is higher and the magnitude of the deviation between the reaction-and transport-limited cementation simulations is smaller, compared to wider aperture distributions. 11 Lithosphere networks. It is important to recognize, however, that fracture apertures are not the only property of a network that will be affected by cementation. Although modeling studies (included this one) typically treat fracture aperture as being constant across an individual fracture's length, aperture profiles in a real rock mass are necessarily elliptical, decreasing to zero towards the fracture tips [105,106]. Thus, it is likely that progressive cementation of fractures will also lead to the destruction of fracture porosity and complete sealing from the tips inward, resulting in an effective decrease in fracture length [27,107] and potentially elimination of fracture interactions depending on their location during cementation.
This work shows that changes in fracture network configuration (at least in the form of aperture values in an otherwise static network) are strongly influenced by cementation, and that manner of this change is dependent on whether the cementation process is reaction-or transportlimited. When cementation is reaction-limited, the rate of   12 Lithosphere cementation (and aperture reduction) is constant between fractures that are not concurrently experiencing widening. At the onset of fluid flow, the most transmissive (widest aperture) fractures that are favorably aligned within the flow direction will accommodate the majority of fluid flow and mass transport. As reaction-limited cementation proceeds, fractures with narrower apertures will necessarily close first (e.g., Figure 3(c)), resulting in fluid flow that becomes progressively more channelized into the "backbone network" [45]. As a result, the propensity for strongly channelized fluid flow through a backbone network increases through time during reaction-limited cementation, as represented by decreasing fracture participation observed in all of our models regardless of their parameterization (Figures 3-5 and 7-9). When cementation in a fracture network is transportlimited, the rates of cementation are proportional to fluid flow and mass transport and thus are variable both within and between individual fractures. More specifically, rates of cementation (and aperture closure) are highest in the most transmissive (widest aperture) fractures in the network but decrease through time as those fractures seal and become less transmissive. Similar to the reaction-limited case, the most transmissive (widest aperture) fractures that are favorably aligned within the flow direction will accommodate the majority of fluid and mass transport at the onset of fluid flow. Higher rates of fluid flow and mass transport in these backbone networks then lead to preferential cementation and sealing, leading to fluid flow that is gradually redistributed into a greater number of originally narrower aperture fractures. As a result, the propensity for strongly channelized fluid flow through a backbone network decreases through time during reaction-limited cementation, as represented by increasing fracture participation observed in all of our models regardless of their parameterization (Figures 3-5 and 7-9). Thus, the identification or prediction of the locus of fluid flow over time is considerably more challenging if cementation is transport-limited, as it is strongly dependent on not only network geometry but also the time since the onset of cementation. Collectively, our results show that mineral precipitation and its associated kinetics have important implications for fluid flow through fracture networks (fault-localized or otherwise). Applied DFN studies should consider these factors in scenarios where the geological setting and relevant timescales permit precipitation reactions. These may include examination of petroleum migration (long timescales, moderate reactivity), geothermal energy (short timescales, high reactivity), and storage of carbon or nuclear waste (short-to-long timescales, low-to-moderate reactivity).
Mineral precipitation within fracture networks also influences the mechanical properties of a rock mass, particularly in fault zones (e.g., [71,81,93,[108][109][110][111]). The seismic cycle, for example, is typified by dramatic, coseismic reductions in the strength of fault rocks during slip (e.g., [112]), which is then progressively recovered during the interseismic period [22,[113][114][115]. This process (known as fault "healing") has important implications for the magnitude and recurrence of earthquakes and is generally thought to be dependent in part on the cementation and sealing of slip-related fractures [22,108,109,116,117]. Unlike permeability reduction, which is continuous (but not linear) during cementation of fractures, it is likely that strength recovery exhibits a more discrete behavior. For example, a small amount of mineral precipitation along fracture walls results in aperture reduction and an associated decrease in permeability but is unlikely to greatly affect rock strength as long as the fracture remains open. Strength recovery can therefore be described as a more binary process that only occurs appreciably once a fracture is completely sealed with cements. As such, our model results have first-order implications for the recovery of strength and stiffness in fault zones following slip. Results suggest that during a reaction-limited cementation process, for example, the smallest aperture fractures in the network rapidly seal. Some amount of strength recovery therefore likely begins relatively shortly after fault slip and the onset of reaction-limited cementation. During transport-limited cementation, however, the apertures and transmissivity of larger fractures in the backbone network must first be reduced before smaller fractures receive sufficient fluid flux to begin sealing. This effect is such that closure of any individual fracture in the network may be delayed until fracture apertures normalize toward similar values. Thus, it is possible that the postslip recovery of fault strength following slip is substantially delayed when cementation is transport-limited when compared to cementation that is reaction-limited.
Our models are, by definition, simplified relative to the geological reality of fault-zone fracture systems. For example, many fault zones exhibit systematic spatial variations in the aperture, orientation, length, and abundance of fractures. Fracture aperture values may increase with proximity to fault cores [118] or other structural features such as relays [81]. They may also be strongly affected by the fracture toughness and elastic modulus of the protolith [119], which likely varies both along strike and down dip in large faults. Those the same variations in protolith also lead to variations in mineral composition and distribution which may affect cementation patterns (e.g., [111,120,121]). Variations in damage-zone width and asymmetry of internal structures represent additional complications not accounted for by our models. For example, our models incorporate fractures that are generally comparable in length to the model domain. Some damage zones, however, extend for 100's of meters perpendicular from the fault core(s) and thus contain fractures whose lengths are short relative to the overall fractured domain. Moreover, pressure gradients (and consequently flow rates) may be subject to variations that are extrinsic with respect to the network/area of interest. Resultant variations in the rate of mass transport may have direct consequences on the ratio between the relative rates of mass transport and surface reaction and hence on the mechanism of mineral precipitation (i.e., reaction-limited or transport limited). Such variations may even be driven by differential rates of cementation and permeability reduction upstream of an area of interest in advective systems. Collectively, these arguments suggest that any real, complex fault zone may be typified by both reaction-and transport-limited cementation 13 Lithosphere processes that vary both in space and geologic time. Oversimplification of reactive fluid transport in fracture networks (fault-localized or otherwise) and ignoring the temporal component of this problem can lead to erroneous estimation of the evolution of geometry and fundamental hydromechanical of fracture networks.

Conclusions
We quantify the effects of reaction-limited (constant precipitation rate) and transportation-limited (aperture-dependent precipitation rate) mineral growth and cementation in fault zones. Models were constructed using the software dfnWorks and were parameterized with fracture geometry data from real fault damage zones. We find that these two endmember mechanisms of mineral precipitation have opposite effects on the distribution of fluid flow and cementation in fault-zone fracture networks. Specifically, when mineral growth in fractures is reaction limited, the locus of fluid flow and cementation becomes increasingly channelized within discrete segments of the fracture networks through time. When mineral growth is transport limited, the locus of fluid flow and cementation becomes increasingly homogeneous and widespread in the fracture networks through time. We show that this effect is persistent regardless of overall network geometry and flow direction, although the magnitude of the effect is dependent on the fracture intensity (p32, Equation (1)) and the initial fracture-aperture distribution. This observation conflicts with previous efforts in forward modeling of fluid flow through fracture networks, which typically assume static network properties. We argue that future DFN studies should consider these factors in scenarios where the geological setting and relevant timescales permit precipitation of mineral cements during fluid movement. This may include applied settings where the operative timescales of relevance are long (e.g., petroleum migration and nuclear waste storage) or the rates of reaction are relatively rapid (e.g., injection of CO 2 -rich brines and geothermal energy production).
Our results also have important implications for fault "healing" or the progressive recovery of fault-zone strength following slip events. In models simulating reaction-limited cementation, for example, the smallest aperture fractures in the network experience closure relatively early after the onset of cementation. Some amount of strength recovery therefore likely begins relatively shortly after fault slip and subsequent postseismic fluid flow. During transport-limited cementation, however, the apertures and transmissivity of larger fractures in the backbone network must first be reduced before smaller fractures receive sufficient fluid flux to begin sealing. This effect is such that complete closure (and associated strength recovery) of any individual fracture in the network may be delayed until fracture apertures normalize toward similar values.
In summary, our results show that mineral precipitation and its associated kinetics have important implications for the evolution of hydromechanical properties of fault zones in the brittle crust. The significant geologic complexities of these structures may ultimately be such that reaction-and transport-limited cementation represent a continuum of behaviors that are both variably operative within individual fault zones and vary as a function of space and time. We argue that a more accurate understanding of fracture systems in the crust (fault-localized or otherwise) will rely on a more complete accounting for the myriad of geochemical processes that may dictate fracture network properties over decadal-to-millennial timescales.

Data Availability
Data sets for this research including generation files for creating the discrete fracture networks of Figure 2 and python scripts for running the simulations with dfnWorks are available at https://osf.io/t6wvk/.

Additional Points
Key Points.

Conflicts of Interest
The authors declare that they have no conflicts of interest.