Quantifying Water Diffusivity and Metamorphic Reaction Rates Within Mountain Belts, and Their Implications for the Rheology of Cratons

The distribution of rheologically strong cratons, and their weakening by metamorphic hydration reactions, is of fundamental importance for understanding first‐order strength contrasts within the crust and the resulting controls on the tectonic evolution of the continents. In this study, the Douglas Harbor structural window within the Paleoproterozoic Trans‐Hudson orogen of Canada is used to study the hydration of the footwall Archean Superior craton basement by water released from the overlying Paleoproterozoic Cape Smith thrust‐fold belt. Phase equilibria modeling is applied to quantify the Archean and Paleoproterozoic metamorphic conditions, and to determine the effect of hydration on basement mineralogy. The amount of structurally bound water within the basement is calculated and shown to decrease as a function of distance below the basal décollement of the thrust‐fold belt. Applying a reactive fluid transport model to these results, the rate coefficient for fluid‐rock reaction is constrained to be 10−19 mol−1m3s−1 , and the diffusivity of water through the grain boundary network to be 10−9 m2s−1 at the ambient metamorphic conditions of 570 ° C and 7.5 kbar. This newly documented rate of water diffusion is three orders of magnitude slower than thermal diffusion, implying that hydration by diffusion may be the rate‐limiting factor in the weakening of cratons, and therefore plays an important role in their geological persistence. This conclusion is consistent with field observations that Paleoproterozoic strain in the Douglas Harbor structural window is restricted to hydrated portions of the Archean Superior craton basement.

2 of 24 hydrous minerals can be present provided they are not interconnected (e.g., Menegon et al., 2011). Such conditions, which we hereafter refer to as "anhydrous," can be most easily achieved for large volumes of crust by partial melting during high-grade metamorphic events (e.g., due to crustal thickening), which concentrates water and other volatiles into the melt that is emplaced at shallow crustal levels, leaving behind an anhydrous residue (e.g., Jackson et al., 2008;White & Powell, 2002). Following erosion, what remains is the lower ∼ 40 km of the crust, representing the anhydrous residue following melt extraction (e.g., McKenzie & Priestley, 2016).
Although the formation of strong anhydrous crust in cratons is becoming increasingly well understood, what is less clear is whether this crust is ever weakened or destroyed, and if so by what mechanisms. Some attention has focused on earthquakes as a mechanism for rheological change. For example, influx of fluids associated with earthquake slip can result in metamorphic transitions from anhydrous metastable granulite to hydrous eclogite, which field relations show is significantly weaker than the surrounding granulite (Austrheim, 1986;Austrheim & Boundy, 1994;Jamtveit et al., 2018). Alternatively, the extremely small grain size of pseudotachylytes formed by earthquake slip can lead to effective diffusion creep in even anhydrous lithologies (e.g., Menegon et al., 2017). Such mechanisms are feasible agents of rheological change, but are only relevant to brittle portions of the crust breaking in earthquakes, and their importance depends upon the volume and connectivity of the regions undergoing seismic strain. In some locations large volumes of pseudotachylytes can be generated in zones that span multiple structural levels (Menegon et al., 2021), whereas in other locations pseudotachylytes are sparse, or the deformation is dominated by ductile processes (e.g., Sibson & Toy, 2006).
In this study, we investigate an alternative and under-appreciated mechanism for rheological weakening of cratons, which could potentially be important at a regional scale and under different geological conditions: the diffusion of water through anhydrous rocks, causing hydration reactions and the development of a weaker mineral assemblage. Such a mechanism is likely to be dominant in regions that are undergoing ductile deformation, in contrast to the earthquake-related processes described above. For this process to be important, large volumes of hydrous rocks would need to be juxtaposed against anhydrous crust, which is a geologically common situation whenever volatile-rich thrust-fold belts and accreted island arcs are emplaced onto the margins of cratons. Present day examples include the Andes overriding the Brazilian shield and the Himalayas being thrust onto the Indian shield (e.g., Barazangi & Ni, 1982;Molnar & Lyon-Caen, 1988;Nabelek et al., 2009;Yin & Harrison, 2000). Compilations of earthquake depth estimates and thermal models suggest that in such settings ductile regions may represent most of the thickness of the orogenic crust, once the underthrusting plate has passed the down-dip seismogenic limit on the range-front thrust fault (e.g., Craig et al., 2012;Cattin & Avouac, 2000). The effectiveness of water transport along grain boundaries in this setting depends upon the rate of fluid influx into the anhydrous crust, the rate of motion through the grain boundary network, and the rate of reaction with the surrounding minerals. Addressing this topic allows us to quantify the poorly known rates of fluid diffusion and metamorphic reactions, which are of general geological significance (e.g., Ague, 2003;Baxter, 2003;Walther & Wood, 1984). The wider implications of understanding this topic arise from the first-order role that strength contrasts play in controlling the geological and tectonic evolution of the continents.
Our study focusses on the Douglas Harbor structural window, which exposes Archean Superior craton basement beneath Paleoproterozoic cover units, primarily composed of the Cape Smith thrust-fold belt (CSB) (Figure 1a). The cover units were thrust over the margin of the Superior craton along a basal décollement during the Himalaya-Tibet-scale Trans-Hudson orogen (THO), which formed following terminal collision of the Superior and Churchill plates at ca. 1.83 Ga (Hoffman, 1988;St-Onge et al., 2006). Late-orogenic, long-wavelength folding and subsequent erosion created the window, exposing a continuum of structural levels including the basement-cover interface (Figure 1b; Lucas, 1989). Here, basement is observed to transition structurally upwards both in terms of metamorphic grade and fabric toward the décollement (over 10-20 km of structural depth), from anhydrous Archean granulite-facies rocks with an Archean gneissosity, to hydrous and recrystallized Paleoproterozoic amphibolite-facies rocks that exhibit a penetrative mineral foliation broadly sub-parallel to the décollement (St-Onge & Ijewliw, 1996). Previous workers have ascribed this relationship to water infiltration and reaction during the propagation of the overlying thrust-fold belt . Dating of accessory minerals (primarily titanite) that are petrographically tied to WHYTE ET AL.
10.1029/2021GC009988 4 of 24 the Paleoproterozoic metamorphism shows that this process occurred over an interval of ∼ 30 Myr at ca. 1,815-1,785 Ma (Lucas, 1989;Scott & St-Onge, 1995;St-Onge et al., 1999. Structurally beneath the zone of hydration and deformation, the Superior craton basement preserves Archean fabrics, and lacks evidence of Proterozoic faulting or ductile deformation . Deformation along the overlying décollement is predominantly ductile (Lucas, 1990), implying that the hydration and deformation occurred down-dip of the seismogenic portion of the rangefront thrusts. This setting is analogous to underthrusting Indian basement material in the mid-crust ∼ 100-250 km north of the range-front in the present-day Himalayan-Tibetan orogen (e.g., Cattin & Avouac, 2000;Craig et al., 2020). As such, the Douglas Harbor structural window presents an ideal natural laboratory to study the effects of weakening cratons by water infiltration in a ductile mid-crustal setting.
By analyzing and modeling the amount of structurally bound water as a function of structural depth beneath the thrust-fold belt, the rate of reaction with the surrounding rocks is constrained. Additionally, because we conclude that diffusion is the dominant mechanism of fluid transport in this setting (as described in Section 6), we are also able to estimate the diffusivity of water through the grain boundary network. These results allow an assessment of the relative contributions of changes in temperature and hydration to the rheology of cratons as they underthrust mountain ranges. Our constraints on the rate of reaction provide a new observation-based estimate of this poorly constrained parameter (e.g., Baxter, 2003).
We begin by describing the geological setting of the Douglas Harbor structural window. Petrographic observations are then integrated with phase equilibria modeling to place constraints on the metamorphic pressure-temperature ( -) conditions and hydration state of the analyzed samples. Using this information, we then describe numerical models of reactive fluid transport that allow us to constrain the rate of water movement and the metamorphic reaction rate. Finally, the implications of our results for the controls on the large-scale deformation of the continental crust in convergent settings are discussed.

Geological Setting
The Ungava Peninsula, northern Quebec, exposes the southern extent of the THO, which is a Paleoproterozoic orogenic belt that extends for over 4,600 km across North America and was formed by early microcontinent and island arc terrane accretions at ca. 1.88-1.865 Ga followed by collision between the lower Superior and upper Churchill plates at ca. 1.83 Ga (Hoffman, 1988;St-Onge et al., 2006;Weller & St-Onge, 2017). In the study region (Figure 1a), the lower plate comprises Archean basement gneiss (Superior craton) and Paleoproterozoic cover (Povungnituk and Chukotat groups), and the upper plate comprises an ophiolite (Watts Group), fore-arc clastic apron (Spartan Group), and a magmatic arc terrane (Parent Group and Narsajuac arc) (St-Onge et al., 2009).
The Superior craton basement is dominantly composed of tonalitic to dioritic orthogneiss, which was subject to granulite-facies metamorphism ( ≥ 800 • C and ≤ 5 kbar) coeval with plutonism and deformation at ca. 2.7 Ga . Sedimentary and volcanic strata (Povungnituk and Chukotat groups) accumulated along the Superior craton margin during rifting from 2.04 to 1.87 Ga (Parrish, 1989;St-Onge et al., 1992), while the ophiolitic Watts Group at ca. 2.00 Ga (Parrish, 1989;Scott et al., 1992) and the Narsajuac arc from 1.86 to 1.82 Ga  formed outboard of the Superior craton margin. Terminal collision of the Superior and Churchill plates was accompanied by obduction, accretion of the allochthonous upper plate units, and imbrication of the parautochthonous cover units, which were thrust

of 24
southwards over the Superior craton margin, forming the CSB that marks the southern front of the THO in northern Quebec (Figure 1b).
During the THO collision, the CSB experienced regional metamorphism lasting ∼30 Myr (1. 82-1.79 Ga;St-Onge et al., 1999). Metamorphic conditions in the CSB vary from greenschist-to amphibolite-facies conditions and increase northwards from 400 • C to 575 • C and 6.3-9.0 kbar, consistent with the interpretation of the CSB as a southwards-tapering thrust wedge (Lucas, 1989). The Superior craton margin, CSB and accreted Narsajuac arc terrane were subject to post-orogenic folding (ca. 1.76-1.74 Ga), producing map-scale, orogen-parallel, basement-cored anticlines. Subsequent erosion of these folds has exposed two primary structural windows into the Superior craton basement in northern Quebec (Figure 1a): the western Kovik structural window and the eastern Douglas Harbor structural window, described below.

Douglas Harbor Structural Window
The Douglas Harbor structural window mostly exposes Superior craton basement, and includes several elongate klippen that correspond to tight, upright fold hinges cored by Paleoproterozoic cover units (Figure 1c). These folds have a lobate-cuspate morphology, characteristic of folding in the presence of significant strength contrasts across the basement-cover interface (Harris et al., 1987;J. Ramsay, 1967). A basal décollement is present at the basement-cover interface, and comprises a ductile shear zone (with a southwards direction of hangingwall transport) that increases in width from 3 m along the southern margin of the structural window to at least 400 m along the northern margin (Figure 1b; Lucas, 1990;. Within the window, the pre-existing Archean gneissosity is progressively transposed into a new Proterozoic foliation sub-parallel to the basal décollement at structural depths of up to 10 km below the shear zone . This deformation represents the creation of a pervasive fabric parallel to the overlying décollement, and is not expressed as distinct shear zones. The Archean mineralogy of the crystalline basement is observed to be modified up to 20 km below the Proterozoic basal shear zone, such that St-Onge and  divided the Superior craton basement into four mineral zones, here listed from deep to shallow structural level: orthopyroxene-clinopyroxene, coronitic garnet, hornblende-epidote, and muscovite-biotite. The orthopyroxene-clinopyroxene mineral zone was interpreted to represent unaltered basement, with Archean fabrics and assemblages consistent with metamorphic conditions of ≥ 800 • C and ≤ 5 kbar, whereas the remaining three mineral zones were interpreted as products of Paleoproterozoic metamorphism. The order of the mineral zones is consistent across the window, but their structural depths relative to the basal décollement are laterally variable (primarily due to variation in protolith lithology, as described below). The coronitic garnet zone is present at structural depths between a minimum of 1 km and a maximum of 20 km and generally retains Archean fabrics, but the precursor Archean granulite-facies assemblage has been modified by the development of garnet coronae around plagioclase grains and the conversion of some pyroxene to amphibole. The corona texture formed over -conditions spanning from 7.7 kbar and 640 • C at the margins of the window, to 9.8 kbar and 700 • C in the core of the antiform, consistent with deeper structural levels being exposed in the center of the fold. The -conditions are markedly higher pressure and lower temperature than the Archean metamorphism, and are interpreted to represent a reaction front for Paleoproterozoic overprinting under fluid-limited conditions (St-Onge & Ijewliw, 1996). The hornblende-epidote mineral zone, characterized by a well-developed Paleoproterozoic foliation defined by the alignment of amphibole, biotite, and epidote, is typically present within 5 km structural depth of the décollement (but locally present at structural depths up to 10 km), and was interpreted to represent increasing levels of hydration and overprinting proximal to the overlying thrust-fold belt. Finally, the muscovite-biotite mineral zone is present at similar structural levels to the hornblende-epidote zone, but is restricted to felsic lithologies, and is characterized by a very well developed, penetrative Paleoproterozoic schistosity defined by the growth and alignment of muscovite and biotite. The distribution of hydration can therefore be seen to have controlled the distribution of the Proterozoic deformation, with pervasive re-orientation of the structural fabric only occurring where the dominant lithologies have become sufficiently hydrated to contain a through-running framework of amphibole, mica, or both.
St-Onge and  interpreted the "hydrous-side up" distribution of the mineral zones to have formed as a result of downward fluid infiltration driven by thermal and/or chemical potential gradients, with generally pervasive fluid transport in the basement contrasting with focused fluid flow along the basal shear zone.

of 24
The source of the fluid was considered to have been the dehydrating metamorphosed clastic rocks in the overlying thrust belt. Based on analysis of the major element compositions of the mineral zones, St-Onge and  suggested that the mineral zones were hydrated without significant changes to the major element chemistry, with the exception of overprinted basement lithologies within the muscovite-biotite mineral zone, which were metasomatized by the addition of sodium and potassium during hydration.

Sample Selection
This study uses thin sections of samples collected from the Douglas Harbor structural window (Figure 1c) by St-Onge and  and St-Onge and Ijewliw (1996). Of the samples analyzed in those studies, we focus only on samples from the northern part of the structural window. We limited our sample selection to this geographical extent based on two observations. First, because the peak Paleoproterozoic metamorphic conditions in the CSB increase from south to north (St-Onge & Lucas, 1991), this spatially restricted subset of samples likely experienced a narrow range of peak Paleoproterozoic -conditions. Second, the geometrically simple and well-constrained nature of the basement-cover interface in the northern part of the structural window allows determination of the structural distances of the samples beneath that interface, which are used as an input for the modeling presented below. We determined the bulk composition of each sample (see below) and divided the samples into three categories: mafic (45-52 wt% SiO2 ), intermediate (52-63 wt% SiO2 ), and felsic ( 63 wt% SiO2 ) (Le Bas et al., 1986). Two samples were excluded due to their anomalous composition, and the resultant data set comprises 22 samples (12 mafic, 6 intermediate, and 4 felsic), which span a range of structural depths (ca. 16 km) below the décollement, and include all four mineral zones described above (Figure 1c).
The ultimate aim of this study is to understand the hydration of the Superior craton basement by water influx from the overlying CSB, for which we need to calculate the water content of the samples as a function of structural depth beneath the CSB. To achieve this aim, first we describe the petrography of the study samples, which we use to determine the H 2 O (the molar content of H2O ) of each sample. Then we undertake phase equilibria modeling to determine and conditions, and the relationships between assemblage, bulk composition and H 2 O . Subsequently, we describe how we calculated the structural depth of each sample. We then combine these sources of information in a reactive fluid transport model for the region.

Petrography
Mineral proportions (Table S1 in Supporting Information S1), mineral compositions (e.g., Table S2 in Supporting Information S1), and whole-rock data ( Table S3 in Supporting Information S1) were acquired for each study sample ( Figure 1c). Representative photomicrographs and phase maps are shown in Figure 2.

Methods
The petrographic methods we used depended on whether the sample thin section had a cover slip, and whether previous whole-rock and/or mineral compositions were available. Mineral proportions were acquired by two methods. For thin sections without a cover slip, the samples were analyzed by Quantitative Evaluation of Materials by Scanning Electron Microscopy (QEMSCAN) using the Quanta650F Scanning Electron Microscope (SEM) at the University of Cambridge, using 10-15 m pixel sizes. Cropped areas of a sub-set of these QEMSCAN images are shown in Figure 2. Thin sections with a cover slip were analyzed by point-counting scanned images of the thin sections using the software JMicroVision. Mineral compositional data were also acquired from two sources. Where available, the electron microprobe data from St-Onge and  and St-Onge and Ijewliw (1996) were used. This data set was supplemented using a Cameca SX-100 electron microprobe at the University of Cambridge, with analyses carried out using a 20 keV acceleration voltage, 20 nA probe current and a 1 m spot size. Representative mineral compositional data for the samples analyzed in the subsequent pseudosection modeling are given in Table S2 in Supporting Information S1. Mineral abbreviations follow Whitney and Evans (2010), and mineral cation proportions were calculated using AX (Holland, 2009; accessed January 2020), which provides an estimate of Fe 3+ where stoichiometric criteria can be applied. Anhydrous mineral compositions were calculated to standard numbers of oxygen per formula unit (pfu), whereas biotite analyses were calculated to a total of 11, amphibole to 23 8 of 24 and epidote to 25 oxygen pfu. Finally, whole-rock data were acquired from two sources, either using the reported X-ray fluorescence data from St-Onge and , or from combining the mineral proportions with the mineral compositions to generate a thin-section-specific bulk composition. Where both sources were available for a given sample, no significant (all within 1 wt.%) or systematic differences were observed, and the thin-section-specific bulk composition is listed in Table S3 in Supporting Information S1 and used in subsequent pseudosection modeling.

Orthopyroxene-Clinopyroxene Zone
Sample D192, described in detail below, is representative of the three samples in the orthopyroxene-clinopyroxene zone (Figure 1c). This sample comprises plagioclase (41 vol.%), clinopyroxene (26 vol.%), amphibole (18 vol.%), orthopyroxene (11 vol.%) and ilmenite (2 vol.%), with accessory biotite, epidote, quartz and pyrite (Table S1 in Supporting Information S1). The major phases form a medium-grained ( ∼ 0.5-1 mm), granoblastic, polygonal aggregate (Figure 2a), typical of equilibration under granulite-facies conditions. Within this aggregate, amphibole occurs in two forms: (a) as platy, khaki green crystals in the matrix (6 vol.%); and (b) as green-blue needles decorating grain boundaries and fractures (12 vol.%; Figure 2b), with the latter commonly intergrown with quartz. Both forms have a similar chemistry, clustering on the edenite-pargasite boundary (Table S2 in Supporting Information S1). The platy amphiboles are interpreted to be part of the peak assemblage, owing to their 120 • grain boundaries and polygonal form, whereas the rim amphiboles are interpreted as the products of retrogression. Plagioclase is present as unzoned grains with a relatively calcic composition (bytownite; Table S2 in Supporting Information S1). Minor epidote is present as anhedral grains, commonly along fractures within plagioclase, and minor biotite is present as ragged laths along some grain boundaries; both minor phases are interpreted as retrograde. Overall, sample D192 is interpreted to have had a peak Archean assemblage defined by orthopyroxene, clinopyroxene (platy), amphibole, plagioclase feldspar, and ilmenite. The minor retrogression observed in this sample is typical of the Archean Superior craton basement and is interpreted to be Archean in age , as opposed to the distinctive textural and mineralogical overprints described below that are Paleoproterozoic in age.

Coronitic Garnet Zone
Garnet-poor sample D158 and garnet-rich sample S132 are representative of the textural variations within the coronitic garnet zone (Figures 2c-2f). Both samples dominantly comprise a matrix of plagioclase, amphibole (pargasite), clinopyroxene, and garnet. Plagioclase grains are optically and chemically zoned, with labradorite interiors surrounded by thin andesine rims, and variably surrounded by garnet coronae (Figures 2c-2f;St-Onge & Ijewliw, 1996). Orthopyroxene forms ragged grains, interpreted as relict products of Archean granulite-facies metamorphism ( Figure 2d; St-Onge & Ijewliw, 1996). Both ilmenite and rutile are observed in the matrix; in places, rutile is observed to mantle ilmenite, suggesting that it is the later formed stable Ti-bearing phase. The conversion of background Archean basement (orthopyroxene-clinopyroxene zone) to this coronitic garnet mineral assemblage can be described by the reaction Pl + Opx = Grt + Cpx + Qz, which is a typical expression of increasing pressure in mafic rocks in a simple chemical system (Eskola, 1915). Consistent with this interpretation, previous studies have shown that the Paleoproterozoic coronitic garnet assemblage developed at higher pressures and lower temperatures than Figure 2. Petrography of the orthopyroxene-clinopyroxene (a and b), coronitic garnet (c-f), and hornblende-epidote (g and h) mineral zones. Images are QEMSCAN phase maps (a, c, e, g) or photomicrographs in plane-polarized light (b, f, h) and cross-polarized light (d). Mineral abbreviations follow Whitney and Evans (2010). (a) Sample D192 matrix, comprising granular orthopyroxene, clinopyroxene and plagioclase, with minor platy amphibole and ilmenite. (b) Close-up of (a), showing two generations of amphibole: platy, khaki green amphibole, which forms 120 • grain boundaries with adjacent pyroxene and plagioclase and is interpreted as a peak Archean phase; and blue-green, fine-grained amphibole developed along grain boundaries and fractures within grains, which is interpreted as Archean retrograde in origin; see text for details. (c) Sample D158 matrix, showing anhedral garnet (orange), which in places forms corona structures mantling plagioclase. Plagioclase forms zoned grains, with bytownite cores (medium blue) rimmed by thin margins of andesine (light blue). (d) Close-up of (c), highlighting an anhedral orthopyroxene grain that is interpreted as relict. (e) Sample S132 matrix, showing a higher proportion of garnet, and thicker sodic rims on plagioclase, than sample D158. (f) Close-up of (e), showing widespread development of garnet coronae around plagioclase. (g) Sample L056A matrix. The preferred alignment of amphibole, biotite and epidote forms a penetrative foliation. (h) Close-up of (g), highlighting titanite surrounding rutile. 9 of 24 the Archean metamorphism, with the corona texture thought to be indicative of fluid-limited conditions (St-Onge & Ijewliw, 1996;.

Muscovite-Biotite Zone
Two samples from this zone (L154 & S108; Figure 1c) are texturally similar to the hornblende-epidote zone, but have a fabric defined by the preferred alignment of muscovite and biotite rather than amphibole (Table  S1 in Supporting Information S1).

Phase Equilibria Modeling
Phase equilibria (pseudosection) modeling was applied to determine the -conditions of the Archean granulite-facies metamorphism and Paleoproterozoic overprinting (Figures 3a and 3b), using Archean orthopyroxene-clinopyroxene sample D192 (Figures 2a and 2b) and Paleoproterozoic hornblende-epidote sample L056A (Figures 2g and 2h), respectively. In addition, to analyze the assemblage variation observed across the Douglas Harbor structural window as a function of water content, -H 2 O pseudosections were constructed for all samples (Figures 3c, 3d, 4, and S1 in Supporting Information S1).

Methods
Pseudosections were constructed using Theriak-Domino (De Capitani & Petrakakis, 2010) with the thermodynamic data set ds-62 of Holland and Powell (2011) Green et al. (2016). For modeling granulite-facies Archean metamorphic assemblages, the "augite" clinopyroxene model was used, whereas for modeling lower-temperature, amphibolite-facies Paleoproterozoic assemblages, the "omphacite" clinopyroxene model was applied (Green et al., 2016). All model bulk compositions are provided in Table S4 in Supporting Information S1, and were calculated in two ways, depending on whether the whole-rock compositions (Table S3 in Supporting Information S1) were acquired by XRF or thin section-specific techniques. For XRF-derived compositions, the Ca content of the samples was modified by assuming that all P2O5 was hosted by apatite (e.g., Weller et al., 2013), and the Fe 3+ was set at 0.1, as this value provided the best fit to the observed assemblages (particularly oxides present) based on -Fe 3+ pseudosection analysis. For thin section-derived compositions, which were calculated by combining mineral proportions with compositions, these compositions were recalculated to discount minerals external to the model system (e.g., pyrite, apatite). In addition, for sample D192, to better recreate the effective bulk composition at peak conditions, we substituted the Archean retrograde amphibole rims (12.1 vol.%) for orthopyroxene and clinopyroxene in a 2:1 (Opx:Cpx) ratio based on petrographic observations. Figure 3a shows a -pseudosection for sample D192. The interpreted peak assemblage is shown in red, indicating Archean metamorphic conditions of ∼ 800 • C and ∼ 4 kbar. The peak assemblage field is bounded at higher temperatures by the appearance of silicate melt, at lower temperatures by the appearance of biotite, Where present, zero isomode field boundaries are colored for orthopyroxene (dark blue), garnet (red), clinopyroxene (dark green), H2 O (light blue), epidote (light green), and titanite (orange). (a) -pseudosection for orthopyroxene-clinopyroxene mineral zone sample D192. The inferred equilibrium assemblage is highlighted in red text, suggesting Archean peak -conditions of ∼ 4 kbar and 800 • C. (b) -pseudosection for hornblende-epidote mineral zone sample L056A. The inferred equilibrium assemblage field is highlighted in red text, and contoured for epidote modal percent. The upper part of the field provides the best fit to observation, suggesting Paleoproterozoic re-equilibration -conditions of ∼ 7.5 kbar and 570 • C. (c) -H 2 O pseudosection at 7.5 kbar for sample D192, showing that most mineral assemblage changes observed in the Douglas Harbor structural window can be attained by increasing water contents at this pressure; see text for details. (d) Equivalent -H 2 O pseudosection for sample L056A, showing that this sample likely equilibrated (red text) under watersaturated conditions; see text for details. 11 of 24 and at higher pressures by the appearance of quartz. Amphibole is stable to temperatures of over 900 • C, consistent with our inference that the platy amphibole was part of the peak assemblage. Garnet is only stable at pressures greater than 5.5 kbar for temperatures over 600 • C (red line, Figure 3a), consistent with the suggestion that samples in the coronitic garnet mineral zone (with a similar bulk composition) experienced higher-pressure conditions during Paleoproterozoic reworking . Likewise, orthopyroxene is only stable at pressures below 7 kbar for temperatures over 600 • C (dark blue line, Figure 3a), consistent with the interpretation that samples that lack (or contain relict) orthopyroxene were overprinted by higher pressure/lower temperature metamorphic conditions. Our result of ∼ 800 • C and ∼ 4 kbar is consistent with the previous estimate of 3.5 ± 2 kbar and ≥ 800 • C by St-Onge and , and consolidates the -conditions of granulite-facies metamorphism of the northern Superior craton (which were coeval with plutonism and deformation at 2.7 Ga; . Uncertainties related to the thermodynamic data used for our -modeling can be calculated, but these underestimate the total uncertainty, which is likely to be at least ± 0.5 kbar and 50 • C (2 ; Palin et al., 2016;Powell & Holland, 2008).

P-T Conditions of Archean and Paleoproterozoic Metamorphism
To constrain the -conditions of Paleoproterozoic hydration during the emplacement of the CSB, we constructed a pseudosection for the hornblende-epidote mineral zone sample L056A described above (Figure 3b). The interpreted equilibrium assemblage (shown in red) occupies a narrow temperature range of 520-570 • C. At lower temperatures, albite is predicted to join the assemblage, and at higher temperatures clinopyroxene is stabilized. Pressure is less well constrained, ranging from 4.5 to 7.5 kbar across the identified equilibrium field. Consideration of the variation of epidote mode through this field (contoured on Figure 3b), which varies from 0 to 4.5 vol.%, compared with the observed epidote mode of 5.1 vol.% (Table  S1 in Supporting Information S1), suggests that the sample most likely experienced pressures in the upper part of the possible range. Mineral mode agreement to within 1 vol.% is considered a reasonable fit, given the accuracy of such estimations in pseudosections (e.g., Lanari & Duesterhoeft, 2018;Waters, 2019). As such, we suggest Paleoproterozoic metamorphic conditions of 570 ± 50 • C and 7.5 ± 0.5 kbar for the 12 of 24 hornblende-epidote mineral zone. This result lies at the low-grade end of the previously reported range of -conditions across the window (7.7-9.8 kbar, 590-720 • C), which were calculated based on samples from the structurally deeper coronitic garnet mineral zone (St-Onge & Ijewliw, 1996). Accordingly, our data are consistent with, and further corroborate, the inference that -conditions in the basement decrease toward the margins of the window (St-Onge & Ijewliw, 1996). Furthermore, our -result is most applicable to the ambient metamorphic conditions of basement hydration.
Highlighted on Figure 3b is the H2O -in line (light blue), which is situated at slightly lower temperatures than the L056A assemblage field (red text), and separates the pseudosection into water-saturated (at higher temperatures, including the observed assemblage) and water undersaturated conditions. These calculations support the interpretation that sample L056A was fully saturated with structurally bound water and formed in the presence of a free fluid phase.

T-MH 2 O Conditions of the Paleoproterozoic Overprint
The effect of water content was further explored by constructing -H2O pseudosections for all samples (Figures 3c, 3d, 4, and S1 in Supporting Information S1). The pseudosections are calculated at a pressure of 7.5 kbar and a temperature range of 500-600 • C, reflecting the conditions of the hornblende-epidote mineral zone defined above, with H2O varying from 0.1 to 7.6 mol.% to encompass near-anhydrous to fully saturated conditions. As discussed below, the -H2O pseudosections show that all four mineral zones in the Douglas Harbor structural window can be accounted for by variations in water content (as a function of proximity to the décollement), the difference in pressure between the Archean and Paleoproterozoic metamorphism, and the range of mafic-felsic bulk compositions. Figure 3c shows the effects of varying water content on stable mineral assemblages in mafic sample D192. For the considered -conditions, the H2O -in line (light blue) indicates that water-saturation is reached at ∼ 5-6 mol.%. Had sample D192 reached equilibrium at these -conditions with its calculated water content of 1.4 mol.% (Table S5 in Supporting Information S1), the equilibrium assemblage (Cpx-Grt-Hbl-Pl-Qz-Rt ± Ttn) would be similar to that observed in the coronitic garnet mineral zone. In the hypothetical scenario that this sample experienced an influx of water during Paleoproterozoic metamorphism, the following assemblage changes would take place in order of increasing H2O : titanite (orange line on Figure 3c) replaces rutile, epidote-in (light green), garnet-out (red), free water-in (light blue) and, below ∼ 510 • C, clinopyroxene-out (dark green).
Notably, the hypothetical assemblage evolution described above closely matches all assemblage changes observed in the Douglas Harbor structural window bar the muscovite-biotite zone: growth of garnet and rutile and breakdown of orthopyroxene due to an increase in pressure at water-undersaturated conditions, and then sequential transformation of garnet-clinopyroxene-rutile-bearing amphibolite-facies assemblages to titanite-epidote amphibolite-facies assemblages as a function of increasing water content. Furthermore, these results support the suggestion that the growth of titanite was directly related to Paleoproterozoic metamorphism between 1,814 and 1,789 Ma (Lucas, 1989;Scott & St-Onge, 1995), and match the petrographic observation of rutile being replaced by titanite in sample L056A (Figure 2h). Figure 3d shows the equivalent analysis for the intermediate-composition sample L056A. Although the topology of the pseudosection is different in detail, similar key assemblage changes are observed from low to high water contents: titanite (orange line on Figure 3d) replaces rutile as the Ti-bearing phase, garnet (red) and clinopyroxene (dark green) become unstable, and epidote (light green) and free water (light blue) join the assemblage, with water-saturation reached at c. 4-5 mol.%, which overlaps the calculated water content of 4.46 mol. % for this sample (Table S5 in Supporting Information S1).
A similar analysis for felsic sample L158, which was collected from the orthopyroxene-clinopyroxene zone, shows that if this sample had experienced water infiltration at the same Paleoproterozoic metamorphic conditions, it would have also recorded the same key assemblage changes as noted above ( Figure S1 in Supporting Information S1). In addition, muscovite is also stabilized at high water contents, yielding a "muscovite-biotite" zone assemblage. This result indicates that hydration of felsic lithologies by pure H2O is sufficient to produce the muscovite-biotite assemblage, and that metasomatic effects are not necessarily required to explain the presence of this lithology in the Douglas Harbor structural window (St-Onge & 13 of 24 . This result also accounts for the co-location of hornblende-epidote and muscovite-biotite zone samples, which is likely controlled in part by variations in basement bulk composition. In the reactive fluid transport models described below, the water content required to saturate the samples will be utilized. To further explore the effect of the bulk compositional variation present in the Douglas Harbor structural window on the position of the water saturation line during the overprinting Paleoproterozoic metamorphic -conditions, the position of the H2O -in lines for all samples are shown in Figure 4, with the dotted lines representing felsic samples, dashed lines showing intermediate samples, and the solid lines marking mafic samples. Most of the samples cluster in the range ∼ 4-6 mol.%, and are relatively insensitive to temperature variations (at these pressure conditions). Some overlap is observed between the lithologies, although the more felsic samples on average require less H2O to be saturated. How these results are incorporated into the modeling is discussed below.

Hydration State as a Function of Structural Depth
For the reactive fluid transport modeling described below, the water content of the samples needs to be quantified as a function of structural distance (i.e., paleodepth) beneath the base of the CSB. The amount of structurally bound water in the samples (in mol.%) was calculated by combining the estimates of the mineral proportions in each sample (Table S1 in Supporting Information S1) with the compositions of the minerals (e.g., Table S2 in Supporting Information S1).
To calculate the structural depth of the samples below the basement-cover interface, the fold geometry was reconstructed (as shown in Figure S2 in Supporting Information S1), making use of the following observations. First, the large-scale structure of the Douglas Harbor structural window is an anticline. At the margins, Paleoproterozoic foliations within the cover units, and those in the upper structural levels of the basement, are broadly parallel to the interface and thus give local estimates of the décollement orientation, indicating typical values of 30-35 • (Figure 1c). Second, within the window, klippen of the overlying Paleoproterozoic cover sequences are observed, defining near-isoclinal synclines. Where numerous foliation measurements are available within a klippe (e.g., inset, Figure 1c), the deformation fabrics indicate that the limbs dip at ∼ 60-80 • toward the hinge of the klippe. These steep dips are typical for the cuspate component of folding generated by buckling of the interface between strong (i.e., basement) and weak (i.e., cover) lithologies, which leads to an overall lobate-cuspate fold geometry where the cusps point into the stronger material (Harris et al., 1987;J. Ramsay, 1967). An average value of 70 • was used for all klippen syncline limbs in the study area. Third, locally cylindrical fold geometries were reconstructed using these constraints and structural distances of the samples below the base of the CSB were calculated (as shown in Figure S2 in Supporting Information S1). Finally, possible errors in these estimates were assessed by calculating the maximum and minimum plausible structural depth for each sample using an angular geometry, and cylindrical geometries with artificially lowered limb angles (by 10 • ), respectively.
The results of this analysis are shown in Table S5 in Supporting Information S1 and Figure 5. The water content of the samples varies from 0.19 to 5.25 mol.%. The range in structural depth of our sample suite is 0.4-16.3 km. Estimated errors for the structural depth are in the range 0.1 km for the shallowest samples to 4.8 km for the deepest ones (Table S5 in Supporting Information S1).
The results are broadly consistent with the documented mineral zone distribution of St-Onge and , with the orthopyroxene-clinopyroxene samples generally more distal from the décollement and containing low water contents, and displaying no spatial overlap with the hornblende-epidote and muscovite-biotite zone samples, which contain higher water contents. However, the coronitic garnet samples span the full range of structural levels, and the samples with both the lowest (T129) and highest (S107) water contents belong to this zone. This analysis indicates that the four-fold classification of St-Onge and  cannot be uniquely correlated with water content, as bulk composition variation within each zone exerts a strong mineralogical control. For example, sample S107 has the second most mafic composition and as such stabilizes abundant amphibole (Table S1 in Supporting Information S1). More generally, the water contents correlate strongly with the amount of amphibole in each sample (which varies from 0 to 60.2 vol.%; Table S1 in Supporting Information S1), with the exception of samples rich in biotite and epidote.
Overall, the data show a decrease in structurally bound water content with increasing distance from the décollement ( Figure 5). In agreement with the suggestion of St-Onge and , the most hydrated samples occur closest to the overlying thrust belt, with those more than 6 km into the basement having a water content indistinguishable from Archean basement samples, assumed to be unmodified during the THO (Figures 2a and 2b).

Reactive Fluid Transport Modeling
Having established the variation of structurally bound water with structural depth (Figure 5), the amount of water required to fully saturate the samples (Figure 4), and the pre-existing water content in the Archean basement (Table S5 in Supporting Information S1), reactive fluid transport modeling can be used to investigate the rate of propagation of the water through the basement underlying the thrust-fold belt, and the rate of fluid-rock reaction. For this modeling, the mol.% data were converted to mol m −3 using a calculated density for each sample based on its mineralogy (typically around 3,000 kg m −3 ).
We first need to establish what effects we need to include in the model, from the perspective of fluid transport mechanisms. Field observations (St-Onge & Ijewliw, 1996; and our petrographic observations above indicate that fluid infiltration occurred through the rock matrix, rather than being concentrated into vein systems (except for isolated locations immediately adjacent to the basement-cover interface; e.g., Lucas, 1990). This observation shows that we should model fluid movement through the rock matrix, rather than in discrete vein systems. This approach is supported by the pattern of water content versus structural depth (Figure 5), which implies that the fluid was transported through the rock matrix from the overlying thrust belt, as such a consistent pattern over a scale of kilometers would not be expected if the fluid budget were controlled by local proximity to veins or other fluid sources. One approach is to model diffusion of the fluid through the grain boundary network. Such a mechanism explains why the fluid moved downwards during hydration of the basement: diffusion driven by chemical potential gradients could result in the downwards motion of the water into the basement from the overlying CSB.
An alternative possibility is that the downwards motion of the fluid could be by advective flow, driven by gradients in the fluid pressure, linked to the thermal and mineralogical evolution of the thrust belt and underlying basement. Even in low-porosity lithologies, such a process might be possible in the presence of reaction-induced nanoscale porosity (Plümper et al., 2017). However, Hanson (1992) and Lyubetskaya  Table S5 in Supporting Information S1. The felsic samples are excluded from the reactive fluid transport modeling; see text for details.

of 24
and Ague (2009) used models of metamorphic reaction and fluid advection to show that the dominant advective flow in collisional belts is upwards, and that downward advection could only occur under special circumstances, often as a transient and short-lived feature. Specifically, advection in the up-temperature direction in the models of Lyubetskaya and Ague (2009) occurred in the initial 1-3 Myr of mountain-building, which in their models features an inverted geotherm. Such model predictions are inconsistent with the ∼ 30 Myr duration of basement hydration and re-equilibration in the Douglas Harbor Window, mostly constrained by the dating of titanite (Scott & St-Onge, 1995;St-Onge et al., 2006). Additionally, comparison of the metamorphic temperatures in the CSB and the underlying basement (St-Onge & Ijewliw, 1996;Lucas, 1991, 1995 andSection 4) shows that the geotherm was not inverted during the hydration of the basement, with the basement samples recording temperatures that are similar to, or larger than, those in the CSB. Downwards advection driven by metamorphic reactions occurred on the margins of the orogens modeled by Lyubetskaya and Ague (2009), but only at lower temperatures and shallower structural levels than those experienced in the Douglas Harbor Window, or during the late stages of exhumation.
The further possibility of hydrodynamic and electrokinetic fluid transport through nanopores was quantified by Plümper et al. (2017), and found to be ≤10 −15 m 3 s −1 for feldspathic rocks. Such a rate would correspond to ∼ 1 m 3 of water, per square meter of decollement surface, passing from the thrust belt into the underlying basement during the 30 Myr duration of hydration. This value is insufficient to account for the ∼ 70 m 3 of water required, per vertical column, to reproduce the pattern seen in Figure 5. We therefore proceed with a model in which the fluid transport is by diffusion through the grain boundary network.

Model Setup
In this section, we describe a numerical model we have used to estimate the fluid mobility and the hydration reaction rates that give rise to the variations of fluid concentration with structural depth below the basement-cover interface. Due to the limited along-or across-strike variability in the properties of the lithologies observed (as discussed above), we model the evolution of water only as a function of depth beneath the thrust interface, which is the direction of greatest variability in water content. We use a model for diffusive transport in the presence of metamorphic hydration reactions that consume water from the grain boundary network, and result in it being structurally bound within hydrous minerals (e.g., amphibole). Our model is based upon that developed by Rudge et al. (2010) to study the serpentinization and carbonation of peridotite, which is a geodynamically equivalent scenario of transport of a reacting fluid through a rock mass. In this approach, the spatial and temporal evolution of the grain-boundary water is given by where is the concentration of free water, is time, is distance below the basal décollement (i.e., the source of free water), is the fluid diffusivity through the grain boundary network, and is the rate of consumption by metamorphic hydration reactions. We solved the diffusive term using a Crank-Nicholson (joint implicit-explicit) finite difference approach (Press et al., 1992).
The reaction rate ( , with units of mol m −3 s −1 ) varies in space and time, and is parametrized using the approach of Rudge et al. (2010), in which the rate of reaction depends upon the local concentrations of the grain-boundary water ( , in mol m −3 ), the amount of additional water that could be structurally bound into the local mineral assemblage through hydration reactions (denoted by , also measured in mol m −3 ), and a rate constant ( , measured in mol −1 m 3 s −1 ), such that ( ) = As the reaction progresses, at each point in the model we track the concentration of free and structurally bound water, and how they evolve through time. The structurally bound water content ( ) is equivalent to our petrological estimates on Figure 5.
As the initial setup for the model, we impose the appearance of the overlying thrust-fold belt with a water content of on the upper boundary of the model domain, and calculate how the water diffuses into the 16 of 24 underlying basement, reacting as it does so. The free parameters in the model are , the diffusivity ( ), and the reaction rate coefficient ( ).
There are also some fixed parameters in the models. We have produced results using our best estimates of these parameters, and also run models using maximum and minimum plausible values, to investigate the effects on the results. For the water content at full saturation, based on the dominantly mixed intermediate and mafic lithologies in the field area, we have used the mean of the highest water content intermediate and mafic samples (2,541 mol/ m 3 ; D039 and S107; Table S5 in Supporting Information S1). We have also run models using just the highest intermediate (2,321 mol/ m 3 ) and mafic (2,760 mol/ m 3 ) values. We impose an initial concentration of structurally bound water in the basement of 500 mol/ m 3 , based on the mean of the values observed in the unaffected granulite-facies samples distant from the basal décollement (at distances greater than 10 km in Figure 5; Table S5 in Supporting Information S1). We have also explored using values of the means of only the intermediate (222 mol/ m 3 ) and only the mafic (639 mol/ m 3 ) lithologies in this distance range. We run the model for 30 Myr, based on the duration of metamorphism estimated by the compilation of St-Onge et al. (2006) (their " M2b " phase of re-equilibration of the basement and metamorphism of the CSB; 1,815 ± 1 Ma-1,785 ± 4 Ma). This duration is mainly constrained by dating of titanite, a mineral found only in re-equilibrated basement units and which the petrological modeling above suggests was the direct result of the water influx modeled here (orange line on Figure 3). The implication of the 30 Myr duration is that following this time, temperatures reduced to a sufficient level that reactions ceased, despite the presence of fluid. We also explore the effect of using durations of 20 and 40 Myr, which covers a range of twice the uncertainty in the dating reported by St-Onge et al. (2006). When comparing our model results to the petrological observations, we make use of our preferred structural depths of the samples below the base of the CSB, and also consider the minimum and maximum plausible depths, as discussed above.
We assume that the thrust-fold belt has a high enough water content that the value on the top surface of our model remains at through time, which is equivalent to assuming that the water consumed by hydration of the basement is small compared with the water content in the thrust-fold belt (present as both pore fluid and that generated by prograde metamorphism of the cover units). This assumption is justified based on: (a) the large thickness of the cover sequences ( 20 km, based on the metamorphic pressure estimates), including clastic sedimentary rocks (Figure 1) that would have been undergoing prograde dehydration reactions; and (b) the small depth-extent of the basement hydration compared with the thickness of the overlying thrust-fold belt, highlighting the small total amount of water that diffused into the basement rocks compared to that available in the overlying thrust-fold belt. For example, the ∼ 70 m 3 of water input to the basement, per vertical column, required to reproduce the pattern shown on Figure 5 is small compared to that contained in a 20 km thickness of thrust belt with a similar hydration state as the most hydrous samples within the basement (900 m 3 ), and the pore fluid released by metasedimentary rocks (e.g., 500 m 3 for 5 km of sediments with an initial porosity of 10%). The presence of syn-tectonic quartz veins along the décollement (Lucas, 1990; supports the inference of large volumes of water being present in the overlying fold belt, in addition to hydrous lithologies such as the abundant pelite, semipelite and micaceous psammite lithostratigraphic units of the lower Povungnituk Group (e.g., St-Onge & Lucas, 1990).
It is not known whether the water input from the thrust belt would have been continual through time, or limited to short-duration events due to metamorphic reactions within the belt. However, metamorphic dehydration reactions are known to release water over a wide range of -conditions due to principally being controlled by continuous reactions (e.g., Nicoli & Dyck, 2018;Ramsay, 1973). Coupled to the dynamics of thrust-fold belt evolution, in which rocks pass through a spatially and temporally varying temperature field (e.g., Cattin & Avouac, 2000;Craig et al., 2020), such a situation would be expected to result in continual fluid availability above the décollement. The spread of metamorphic ages throughout the ∼ 30 Myr duration of hydration (Scott & St-Onge, 1995) supports this view, and implies that the fluid input was relatively continual through time. Fluctuations on timescales that are short compared with the diffusive timescale would have little effect on the overall results, in which case the values of that we estimate would represent the temporal average. We show below that the estimates of our main parameters of interest ( and ) are relatively insensitive to the value of . Our assumption of constant through time is therefore justified, and has little effect upon our results.
By modeling a range of diffusivities ( ), reaction rate coefficients ( ), and thrust belt water contents ( ), we are able to use the variation of structurally bound water as a function of depth below the basement-cover interface to place bounds on these values. We will first describe the behavior of the model, and then the results of inversions to obtain the best-fitting combination of parameters, and explore the trade-offs between them. Figure 6a. The free water is initially zero everywhere except on the boundary with the overlying Paleoproterozoic thrust-fold belt (i.e., the left side of the graphs), representing the lack of free water in the Archean granulite-facies basement. The level of structurally bound water is initially non-zero, set by the level of existing hydrous phases in the Archean granulite-facies assemblages (i.e., the highplaty amphiboles). As free water diffuses into the basement rocks (dotted lines), it reacts with the anhydrous assemblage, and the level of structurally bound water increases (solid lines). Once the structurally bound water has reached the saturation level for a given composition (i.e., "full humidity," in the sense of not being able to contain any further structurally bound water, and any additional fluid being present as a free phase; 2,500 mol m −3 in the example in Figure 6a), it stops increasing at those locations. In this situation, the concentration of free water trends toward a linear diffusive gradient in the fully-saturated regions (i.e., at distances less than 1 km at 30 Myr in the figure). Where the rocks are undersaturated with respect to structurally bound water, that concentration, and that of the free fluid, are non-linear, representing the dual influences of diffusion and reaction.

An example model result is shown in
The main influence of diffusivity on the model results is to control how far into the crystalline basement the free water penetrates over the timespan of the model, and therefore control the depth-extent of the resulting structurally bound water. Figure 6b shows the final profiles of structurally bound water after 30 Myr for models where the diffusivity has been varied. This figure demonstrates the non-linear relationship between 18 of 24 the diffusivity and the distance over which the hydration reactions occur, due to the characteristic timescales and lengthscales of diffusion being related by = 2 ∕ 2 , where is timescale and is lengthscale, meaning that is proportional to √ .
Changing the reaction rate coefficient has the effect of changing the gradient of the structurally bound water content. If the reaction rate is high, all free water is consumed as soon as it encounters undersaturated rock, and a sharp hydration front is formed (Figure 6c). For low reaction rates, little free water is consumed, so for a given diffusivity it propagates further into the basement, producing a wide zone of low-concentration hydration (Figure 6c).
The water content in the overlying thrust-fold belt controls the model behavior by changing the free water concentrations (i.e., the y-intercept of the dotted lines on Figure 6a). The effect on the level of structurally bound water is due to this value controlling the concentration of free fluid at a given location within the basement, if the other parameters are unchanged. The form of Equation 2 means that higher concentrations of free fluid result in higher reaction rates and therefore higher values of structurally bound water (up to the saturation value).
The models in Figure 6 show that where sufficient fluid input and reaction has occurred to generate a domain of fully saturated rocks adjacent to the décollement, the pattern of structurally bound water has a sigmoidal shape. The form of the reaction rate parameterization means that the sigmoid is not symmetric about the mean value. Where there has been insufficient fluid input to generate a region of fully saturated rocks, the concentration of structurally bound water instead has a gradient adjacent to the décollement.
In order to investigate the ranges of diffusivity, reaction rate coefficient, and that are able to reproduce the structurally bound water distribution shown on Figure 5, we performed a series of grid-searches through the parameter space, in which we have varied the above parameters, and examined the fit with the petrological and structural results described above. We initially describe the results obtained using our best estimates of the fixed parameters described above (duration, structural depth of samples, fluid-saturated water content, and pre-existing hydration level), and then the effects of varying these parameters. In these inversions we excluded the four felsic samples with anomalously low concentrations of water required to attain saturation at the -conditions documented for the Paleoproterozoic metamorphism (L154, D016, M045, T129; Figure 4). We instead use the 18 intermediate and mafic samples with relatively consistent calculated water concentrations required to attain saturation, and which are thought to be lithologically representative of the Douglas Harbor structural window as a whole. Figure 7 shows the best-fitting model, along with graphs that illustrate the misfit between the model predictions and the petrological and structural results as a function of the combination of free parameters. Our main parameters of interest, the diffusivity and reaction rate coefficient, are well-constrained. The diffusivity has a value of 10 −9 m 2 s −1 , and the reaction rate coefficient a value of 10 −19 mol −1 m 3 s −1 . We discuss below the wider importance of these values. The models are not very sensitive to the water content in the overlying thrust belt, which only shows a minor trade-off with the other model parameters (Figure 7).
The supplemental information (Figures S3-S10 in Supporting Information S1) contains equivalent results to Figure 7, produced using alternative values for the fixed parameters described above. The results are summarized on Figure 8, which shows the best-fitting diffusivity and reaction rate for each inversion, and the region in which any model fits the observations to within an RMS misfit of 400 mol m −3 . One inversion, based on an intermediate-composition background hydration state, failed to produce a good fit to our data. The poor fit is accompanied by a higher diffusivity, as the model compensates for the unrealistically low background hydration state by including water propagation 10 km from the fold belt, which in turn fails to reproduce the sharp decline in concentration values at distances 4 km. The remaining inversions show best-fitting models that cluster in a tight parameter space. The models do not lie symmetrically to either side of our preferred model, due to the influence of tradeoffs with the third free parameter, . However, based on the distribution of parameter values that result in a good fit to the data (the gray shading on Figure 8), we can conclude that the diffusivity lies in the range 1 × 10 −9 -1 × 10 −8 m 2 s −1 , and the reaction rate coefficient lies in the range 2 × 10 −20 -3 × 10 −18 mol −1 m 3 s −1 , regardless of the values chosen for the fixed parameters.

Discussion
By combining petrography and phase equilibria modeling, we have quantified the water content within the Archean Superior craton crystalline basement rocks underlying the Paleoproterozoic Cape Smith thrustfold belt, and by modeling this distribution have obtained estimates of water diffusivity and hydration reaction rates during Paleoproterozoic metamorphic overprinting of the basement.
Despite much attention given to the importance of fluid-rock interaction in metamorphism and for rheology, to the best of our knowledge there are no existing estimates from a natural setting for the rate of diffusion of water along grain boundaries in an initially anhydrous assemblage. Our estimated rate of fluid diffusion ( 10 −9 m 2 s −1 ) is slightly lower than the lab experiments of water infiltration into quartzite of Nakamura and Figure 7. Best-fitting model using our preferred parameter values (see text for details). The upper panel shows the data used in the inversions as red circles, and the best-fitting model (solid lines are the structurally bound water at times of 0.1, 1, 5, 10, 20, and 30 Myr, and dotted lines the free water at the same times). The lower panels show the RMS misfit between the model and the data as a function of the varied parameters. The best-fit has a misfit of 375 mol m −3 , and values over 700 are colored red. The best-fitting diffusivity is 4 × 10 −9 m 2 s −1 , reaction rate coefficient is 2 × 10 −19 mol −1 m 3 s −1 , and is 2,500 mol m −3 .

of 24
Watson (2001) (2 × 10 −8 m 2 s −1 ), which may be a result of the experiments being conducted at higher temperatures than were experienced in our study region (600-900 • C vs. 570 • C). Our estimated diffusivity is also less than the rate of cation diffusion through an intergranular fluid film (e.g., 10 −8 m 2 s −1 ; Rubie, 1986), implying that in water-limited systems, the rate of chemical transfer may be limited by the rate of fluid influx, rather than the rate of cation motion through that fluid once it is present.
The thermal diffusivity of silicate rocks is ∼10 −6 m 2 s −1 (e.g., Miao et al., 2014;Whittington et al., 2009), which is considerably greater than our estimated rate of fluid diffusion. A tectonic implication of this appreciable difference in diffusion rates is that where grain-boundary diffusion is the dominant control on water transport, as in the Douglas Harbor structural window, anhydrous crust underthrusting the margin of a mountain range will heat up faster than it becomes hydrated by fluids derived from the overlying thrust-fold belt (Figure 9). The significance of this comparison is that both temperature and hydration state have an important influence on crustal rheology, with the difference in viscosity between "dry" and "wet" lithologies being orders of magnitude (e.g., Mackwell et al., 1998;Rybacki & Dresen, 2004). This strength contrast would require a temperature increase of hundreds of degrees to achieve by thermal means alone. Therefore, in anhydrous crust underthrusting a mountain belt, two effects with different timescales will affect the rheology. First, thermal diffusion will occur, with a timescale that depends on the thermal diffusivity. Second, the underthrusting plate will become hydrated, at a rate that depends on the (lower) chemical diffusivity that governs the water transport into the rocks (Figure 9). Exactly when during this process the underthrusting plate becomes weak enough to deform will depend on the temperature, mineral assemblage, and the forces being exerted upon it. In the case of the Douglas Harbor structural window, the deformation as defined by the presence of overprinting Paleoproterozoic fabrics is strongly correlated with hydrous mineral Figure 8. The effects of changing the values of parameters fixed in the inversions. The red star shows the best-fitting diffusivity and reaction rate coefficient for our preferred model, as shown in Figure 7. The upwards-pointing triangles show the effects of increasing each parameter (color-coded according to the key), and the downwards-pointing triangles the effects of decreasing that parameter. See text for description of the ranges the parameters are varied within. The gray-shaded region shows the parameter space over which any model fits the data with an RMS misfit of 400 mol m −3 .

of 24
assemblages (St-Onge & Ijewliw, 1996;; see, for example, sample L056A above). Therefore, we conclude that in at least some orogens, and at some pressure-temperature conditions, the weakening of the underthrusting plate is controlled by the slow rate of hydration, rather than the more rapid rate of thermal diffusion, implying that temperature increases alone are insufficient to weaken the rocks sufficiently to deform. This point has considerable implications for our understanding of the rheological evolution of the lower crust in orogenic belts.
To compare our estimated reaction rate coefficient of ∼ 2 × 10 −19 mol −1 m 3 s −1 with other estimates, we convert to the units of g cm −2 yr −1 used by Baxter (2003). We assume that the water-consuming reactions occur on the surface of grains of 0.5 mm diameter (based on the petrography described above), and convert our reaction rate of water into a reaction rate of solid minerals using a weight percentage of water in amphibole of 2%. To make our values directly comparable to those of Baxter (2003), we make the same assumption that 20% of the rock is composed of the reacting minerals, which is equivalent to the ∼ 20% increase in amphibole proportion between sample D192 and L056A (Figure 2). Our resulting reaction rate is 3.2 × 10 −8 g cm −2 yr −1 . This value is significantly lower than laboratory-based estimates, but in agreement with other geologically constrained values as compiled by Baxter (2003). Along with Hetnyi et al. (2021), our study provides an important additional constraint on this value, because of the natural setting we have used being characterized by significantly longer lengthscales than the previous estimates.
A striking feature of the above analysis is that the penetration of hydration and deformation into anhydrous crust on timescales of tens of millions of years is on the order of kilometers. This result is independent of whether our choice of a diffusive model was correct-the dating of the duration of hydration (St-Onge et al., 2006), and the distribution of hydrous lithologies, show this to be the lengthscale of water-rock reaction over 30 Myr regardless of the underlying mechanism. The implication of this result is that over the lifetimes of most orogenic belts, it is unlikely that the entire 30-40 km thick crust of a craton underthrusting the margins of a convergent thrust-fold belt will be hydrated in the ductile regime, and become deformable. Although faulting-based mechanisms of weakening have been proposed (e.g., Austrheim & Boundy, 1994;Figure 9. Schematic cartoon (vertically exaggerated) illustrating a range of factors affecting the rheology of an anhydrous craton underthrusting hydrous crust (such as imbricated cover units and accreted terranes). Blue arrows indicate areas of H2O diffusion from the dehydrating overriding wedge into the anhydrous craton below. Thermal diffusion (white arrows) operates along thermal gradients from warmer (lightly shaded) to cooler (dark shaded) portions of the lithosphere. The brittle-ductile transition is shown as a green dashed line, taken schematically from Craig et al. (2012). The Douglas Harbor structural window (purple dashed box) falls within the ductile part of the mountain belt, where the deformation is controlled primarily by the degree of hydration. Jackson et al., 2004;Jamtveit et al., 2018;Menegon et al., 2017), the spatial and temporal effectiveness of such weakening is variable. In cratonic plate interiors, the entire thickness of the crust can be seismogenic (e.g., Jackson et al., 2008), and large volumes of pseudotachylyte can be generated that will have important rheological effects (e.g., Menegon et al., 2021). However, within active mountain belts, the efficiency of this mechanism is more limited, because the heating of the crust means that a smaller volume of material is cool enough to break in earthquakes, which are limited to sparse events where the underthrusting plate bends near the range-front, and close to the Moho beneath the mountain belt (e.g., Craig et al., 2012) as indicated on Figure 9. Mountain-building therefore does not represent an effective mechanism for weakening and reworking anhydrous cratonic crust. This trait may in part explain the geological persistence of cratons since the Archean, as well as the apparent spatial continuity of present-day regions of thick sub-cratonic lithospheric mantle in the Permian (McKenzie et al., 2015), which implies significant reworking and removal of the intervening mobile belts and oceanic lithosphere, but minimal cratonic strain.

Conclusions
By combining field observations, petrology, structural constraints and numerical modeling, we have investigated the hydration and re-equilibration of the anhydrous Archean Superior craton basement that was overthrust by CSB during the Paleoproterozoic Trans-Hudson orogeny. We show that the variation in mineral assemblages in the Douglas Harbor window can be principally accounted for by the difference in pressure between Archean and Paleoproterozoic metamorphism, and by the variation in water content as a function of proximity to the basal décollement. We have been able to constrain both the rate of water diffusion through a grain boundary network ( 10 −9 m 2 s −1 ), and also the rate coefficient for metamorphic hydration reactions ( 10 −19 mol −1 m 3 s −1 ) converting anhydrous granulite-facies mineral assemblages to hydrous amphibolite-facies assemblages. Our results demonstrate that hydration occurs over a fundamentally different timescale than crustal heating during mountain-building, slower by three orders of magnitude, and that this longer timeframe may explain the geological longevity and apparent indestructibility of cratonic continental crust.

Data Availability Statement
This work is based upon samples originally described by St-Onge and Lucas (1995) and St-Onge and Ijewliw (1996). All of our new analyses are presented in the tables in Supporting Information S1, and are also hosted on the Apollo repository (https://www.repository.cam.ac.uk; https://doi.org/10.17863/CAM.76855).