Thibaut Barreyre
email: thibaut.barreyre@uib.no

Jean-Arthur Olive

Timothy J Crone

Robert A Sohn

Depth-Dependent Permeability and Heat Output at Basalt-Hosted Hydrothermal Systems Across Mid-Ocean Ridge Spreading Rates

The permeability of the oceanic crust exerts a primary influence on the vigor of hydrothermal circulation at mid-ocean ridges, but it is a difficult to measure parameter that varies with time, space, and geological setting. Here we develop an analytical model for the poroelastic response of hydrothermal exit-fluid velocities and temperatures to ocean tidal loading in a two-layered medium to constrain the discharge zone permeability of each layer. The top layer, corresponding to extrusive lithologies (e.g., seismic layer 2A) overlies a lower permeability layer, corresponding to intrusive lithologies (e.g., layer 2B). We apply the model to three basalt-hosted hydrothermal fields (i.e., Lucky Strike, Main Endeavour and 9846 0 N L-vent) for which the seismic stratigraphy is well-established, and for which robust exit-fluid temperature data are available. We find that the poroelastic response to tidal loading is primarily controlled by layer 2A permeability, which is about 3 orders of magnitude higher for the Lucky Strike site ($10 210 m 2 ) than the 9846 0 N L-vent site ($10 213 m 2 ). By contrast, layer 2B permeability does not exert a strong control on the poroelastic response to tidal loading, yet strongly modulates the heat output of hydrothermal discharge zones. Taking these constraints into account, we estimate a plausible range of layer 2B permeability between $10 215 m 2 and an upper-bound value of $10 214 (9846 0 N L-vent) to $10 212 m 2 (Lucky Strike). These permeability structures reconcile the short-term response and long-term thermal output of hydrothermal sites, and provide new insights into the links between permeability and tectono-magmatic processes along the global mid-ocean ridge.

Introduction

Circulation of hydrothermal fluids through young oceanic crust at mid-ocean ridges (MORs) accounts for up to $10% of Earth's internal heat loss [START_REF] Elderfield | Mid-ocean ridge hydrothermal fluxes and the chemical composition of the ocean[END_REF][START_REF] Sclater | The heat flow through oceanic and continental crust and the heat loss of the earth[END_REF][START_REF] Stein | A model for the global variation in oceanic depth and heatflow with lithospheric age[END_REF][START_REF] Williams | Heat loss from the Earth: New estimate[END_REF], controls the thermo-mechanical state and degree of hydration of newly formed oceanic lithosphere (e.g., [START_REF] Emmanuel | Suppression and stimulation of seafloor hydrothermal convection by exothermic mineral hydration[END_REF][START_REF] Iyer | Feedbacks between mantle hydration and hydrothermal convection at ocean spreading centers[END_REF], profoundly affects the chemistry of the oceanic lithosphere and global ocean (e.g., [START_REF] Edmond | Ridge crest hydrothermal activity and the balances of the major and minor elements in the ocean: The Galapagos data[END_REF][START_REF] Elderfield | Mid-ocean ridge hydrothermal fluxes and the chemical composition of the ocean[END_REF]Humphris et al., 1995;[START_REF] Wolery | Hydrothermal circulation and geochemical flux at mid-ocean ridges[END_REF], and provides a unique habitat for complex chemosynthetic ecosystems (e.g., [START_REF] Kelley | Volcanoes, fluids, and life at Mid-Ocean Ridge spreading centers[END_REF][START_REF] Lutz | Interrelationships between vent fluid chemistry, temperature, seismic activity, and biological community structure at a mussel-dominated, deep-sea hydrothermal vent along the East Pacific Rise[END_REF][START_REF] Marcon | Megafaunal distribution and assessment of total methane and sulfide consumption by mussel beds at Menez Gwen hydrothermal vent, based on geo-referenced photomosaics[END_REF][START_REF] Nees | Hydrothermal vent mussel habitat chemistry, preand post-eruption at 9o50 North on the East Pacific Rise[END_REF][START_REF] Shank | Temporal and spatial patterns of biological community development at nascent deep-sea hydrothermal vents along the East Pacific Rise[END_REF]. All of these hydrothermal systems processes are modulated by the permeability structure of the elastic matrix hosting fluid flow (e.g., [START_REF] Coumou | The structure and dynamics of mid-ocean ridge hydrothermal systems[END_REF][START_REF] Driesner | The interplay of permeability and fluid properties as a first order control of heat transport, venting temperatures and venting salinities at mid-ocean ridge hydrothermal systems[END_REF], which is a heterogeneous tensor field that varies in time and space, locally, and with geological context, globally. Constraining the matrix permeability structure is a fundamental problem for virtually all fields of hydrothermal research, but placing meaningful constraints on the parameters is technically difficult and often extremely costly (e.g., ocean drilling), which helps explain why permeability is probably the most poorly constrained hydrologic parameter for the oceanic crust [START_REF] Fisher | Rates and patterns of fluid circulation[END_REF].

Direct in situ measurements (e.g., borehole packers) and indirect estimates and inferences (e.g., borehole temperatures, seafloor heat flow) for oceanic crust permeability are highly variable and span many orders of magnitude, with values as high as 10 211 m 2 near the seafloor to values lower than 10 218 m 2 at the base of the upper crust (e.g., [START_REF] Becker | New evidence for age variation and scale effects of permeabilities of young oceanic crust from borehole thermal and pressure measurements[END_REF]Becker & Fisher, 2008;[START_REF] Davis | Thermally stimulated ''runaway'' downhole flow in a superhydrostatic ocean crustal borehole: Observations, simulations, and inferences regarding crustal permeability[END_REF][START_REF] Fisher | Permeability within basaltic oceanic crust[END_REF][START_REF] Fisher | Borehole-to-borehole hydrologic response across 2.4 km in the upper oceanic crust: Implications for crustal-scale properties[END_REF][START_REF] Winslow | Characterizing borehole fluid flow and formation permeability in the ocean crust using linked analytic models and Markov chain Monte Carlo analysis[END_REF]. The decreasing trend with depth is thought to reflect the closure of cracks and pores due to increasing lithostatic pressure and changes in crustal lithology. Although ocean crust permeability measurements and estimates from direct and indirect in situ methods vary widely, convection models of hydrothermal circulation at MORs suggest that the ''effective permeability'' of all discharge zones is in the range of 10 215 -10 212 m 2 (e.g., [START_REF] Coumou | The structure and dynamics of mid-ocean ridge hydrothermal systems[END_REF][START_REF] Driesner | The interplay of permeability and fluid properties as a first order control of heat transport, venting temperatures and venting salinities at mid-ocean ridge hydrothermal systems[END_REF][START_REF] Hasenclever | Hybrid shallow on-axis and deep off-axis hydrothermal circulation at fast-spreading ridges[END_REF][START_REF] Lowell | Hydrothermal processes at mid-ocean ridges: Results from scale analysis and single-pass models[END_REF][START_REF] Lowell | Characteristics of magma-driven hydrothermal systems at oceanic spreading centers[END_REF]. The effective permeability represents the permeability that a single layer stretching between the magma chamber and the seafloor would have to have in order to match the thermo-chemical output observed at a hydrothermal field. The effective permeability thus represents an average in some sense over the entire layered crustal section between the magma chamber and the seafloor, immediately below the venting sites. It is noteworthy that the range of effective permeabilities estimated by [START_REF] Lowell | Characteristics of magma-driven hydrothermal systems at oceanic spreading centers[END_REF] is considerably narrower than the variability documented by in situ measurements and indirect estimates. Therefore, the applicability of borehole measurements carried out mostly off-axis and at the sample scale is questionable for crustal-scale convection on-axis, where very few subseafloor in situ measurements have been made [START_REF] Wilcock | Geophysical constraints on the sub-seafloor environment near mid-ocean ridges[END_REF]).

These considerations have motivated a variety of innovative approaches for constraining matrix permeability over meaningful length scales and in young crust (e.g., [START_REF] Crone | Permeability structure of young ocean crust from poroelastically triggered earthquakes[END_REF][START_REF] Davis | Thermally stimulated ''runaway'' downhole flow in a superhydrostatic ocean crustal borehole: Observations, simulations, and inferences regarding crustal permeability[END_REF][START_REF] Lowell | Hydrothermal processes at mid-ocean ridges: Results from scale analysis and single-pass models[END_REF][START_REF] Theissen-Krah | Coupled mechanical and hydrothermal modeling of crustal accretion at intermediate to Fast Spreading Ridges[END_REF][START_REF] Wilcock | Geophysical constraints on the sub-seafloor environment near mid-ocean ridges[END_REF][START_REF] Wilcock | Estimates of crustal permeability on the Endeavour segment of the Juan de Fuca mid-ocean ridge[END_REF][START_REF] Winslow | Characterizing borehole fluid flow and formation permeability in the ocean crust using linked analytic models and Markov chain Monte Carlo analysis[END_REF], including the use of vent exit-fluid temperature and ocean tide (bottom pressure) time series data to model the poroelastic response of the matrix to tidal loading [START_REF] Barreyre | Poroelastic response of mid-ocean ridge hydrothermal systems to ocean tidal loading: Implications for shallow permeability structure[END_REF]Barreyre et al., 2014a). These sorts of passive techniques that do not require sampling rocks from below the seafloor provide a cost-effective way to constrain permeability, but to this point it has been difficult to separate variability due to methods from that due to site-specific geology and hydrology.

In this paper, we investigate the permeability structure of the discharge zone of hydrothermal systems from a variety of geological settings on the global MOR system using two approaches. First, we develop a multilayer model for the poroelastic response to ocean tidal loading. We focus on the case of two layers of decreasing permeability with depth, which is consistent with the lithostratigraphy of volcanic crust formed at MORs. This approach allows us to compare the poroelastic response from sites with known and variable lithostratigraphies, and allows us to constrain the depth intervals over which the poroelastic response is controlled. Second, we attempt to reconcile our estimates of permeability layering with the observed heat output of hydrothermal sites using an analytical framework derived from [START_REF] Driesner | The interplay of permeability and fluid properties as a first order control of heat transport, venting temperatures and venting salinities at mid-ocean ridge hydrothermal systems[END_REF]. Our analysis shows that the poroelastic response to ocean tidal loading is controlled by the permeability structure of extrusive layer 2A, which varies considerably across sites along the global MOR, while the effective permeability of the system and the efficiency of convective heat transfer is primarily controlled by the permeability structure of intrusive layer 2B. Taken together, these results improve our understanding of how the permeability structure for deep-sea hydrothermal systems varies across MOR spreading rates, and how it modulates the vigor of hydrothermal exchanges.

A Model for Tidal Modulation of Hydrothermal Systems With Depth-Dependent Permeability

We begin by presenting a 1-D model aimed at predicting the perturbations of the pressure, temperature, and velocity of fluids moving through a layered oceanic upper crust subjected to periodic fluctuations in seafloor pressure, imposed by tidal cycles (Figure 1). Our model is a generalization of the homogeneous permeability model of Jupp and Schulz (2004), and it extends the preliminary approach of [START_REF] Xu | A preliminary 1-D model investigation of tidal variations of temperature and chlorinity at the Grotto Mound, Endeavour Segment[END_REF], who derived two-layer solutions for the phase lag of vent temperatures relative to tides but did not completely analyze the system behavior and solution space.

Fundamental Concepts and Equations of Poroelasticity

When a porous medium (e.g., the oceanic crust) is placed under an external load (e.g., ocean tides), the resulting stress is borne partly by the solid matrix and partly by the interstitial fluid, where it manifests as a change in fluid pressure, often referred to as a pore pressure perturbation. The partitioning of the total stress between the fluid and the solid matrix is a function of the elastic properties of the system and the porosity. It is thus possible for a locally uniform seafloor load to produce a time-periodic, spatially variable pressure perturbation field when there are gradients in these properties. The associated pressure gradients can drive interstitial fluid flow and perturb the background flow pattern (e.g., a hydrothermal convection cell). By perturbing the advection of mass and heat in regions of background fluid upwelling, the tidal forcing can, in principle, produce a phase lag between the loading function (seafloor pressure) and the velocity and temperature of the exit-fluids, which depends upon the poroelastic parameters of the system (e.g., [START_REF] Crone | Modeling the effects of tidal loading on mid-ocean ridge hydrothermal systems[END_REF][START_REF] Jupp | A poroelastic model for the tidal modulation of seafloor hydrothermal systems[END_REF][START_REF] Wang | Theory for the propagation of tidally induced pore pressure variations in layered subseafloor formations[END_REF].

The equations governing the temporal (t) and spatial (z) evolution of pore pressure perturbation (p, departure from a reference pressure) within a porous, fluid-saturated domain consisting of N layers of uniform physical properties [START_REF] Biot | General theory of three-dimensional consolidation[END_REF][START_REF] Van Der Kamp | Theory of Earth tide and barometric effects in porous formations with compressible grains[END_REF][START_REF] Wang | Theory for the propagation of tidally induced pore pressure variations in layered subseafloor formations[END_REF], and subjected to variations in surface load r B t ð Þ is:

@ 2 p j @z 2 5 1 g j @p j @t 2c j @r B @t ; (1)
where the notation for this and subsequent equations is given in Table 1, and subscript j (j51; 2; . . . ; N 2 1; N) refers to the properties of the jth layer, numbered downward. Parameters g and c are the hydraulic diffusivity and loading efficiency, respectively. The hydraulic diffusivity is defined as:

g5 k lS ; ( 2
)
where l is the dynamic viscosity of the pore fluid, k is the permeability of the porous medium, and S is the uniaxial storage compressibility, defined as:

S5 K 21 2K 21 s À Á 12 4w 3 1/ K 21 f 2K 21 s À Á ; ( 3
)
where K is the matrix drained bulk modulus, K s is the solid grain bulk modulus, K f is the fluid bulk modulus, / is the porosity of the medium, and w is the dimensionless poroelastic stress coefficient [START_REF] Detournay | Fundamentals of poroelasticity[END_REF]. Detailed expressions for drained parameters K, determined by Gassmann's equation [START_REF] Gassmann | € Uber die elastizit€ at por€ oser medien[END_REF], and for dimensionless coefficient w are given in the Appendix A.

The loading efficiency c is defined as:

c5 b 11m ð Þ 3 12m ð Þ22ab 122m ð Þ ; ( 4
)
where m is the drained Poisson's ratio, a is the Biot-Willis parameter, and b is the Skempton's coefficient, all defined in the Appendix A. [START_REF] Wang | Theory for the propagation of tidally induced pore pressure variations in layered subseafloor formations[END_REF] rewrite equation (1) in terms of an instantaneous pore pressure perturbation component p i j and a diffusive pore pressure perturbation component p d j , which emphasizes the role c and g play in characterizing the solution: where

p i j 5c j r B ; (5) @p d j @t 2g j @ 2 p d j @z 2 50; (6)
p j 5p d j 1p i j : (7)
Note that equation ( 6) has the form of a linear diffusion equation with diffusivity a piece-wise constant function of depth.

Analytical Solutions for 1-D Multilayer (and Multifrequency) Model

In the following subsections, we outline the general solution for each parameter of interest in layer j, namely: pore pressure perturbation, vertical velocity, and temperature. These solutions rely on generic coefficients (e.g., a j , b j , A j , B j , and C j ) that are determined by applying appropriate boundary conditions (BCs) to equation (6).

Pore Pressure

Tidal loading at the seafloor can be described as a periodic function of time: where p T and x are the amplitude and angular frequency of the tidal loading (e.g., corresponding to a 12 and 24 h period for the M2 and K1 tides, respectively). With this type of harmonic forcing, the solution for pore pressure perturbation p j z; t ð Þ À Á in a multilayer system can be written: p j z; t ð Þ5 a j e l:p:z=Dj 1b j e 2l:p:z=Dj 1c j h i p T e ixt ;

where l is a complex constant equal to 11i, and D is the skin depth defined as

D j 5 ffiffiffiffiffiffi ffi pg j f r ; (10)
and a j and b j denote 2N complex coefficients. These are determined using the Fourier-transformed BCs, which involves solving a system of 2N complex algebraic equations. In the simplest case with N 5 2 (Figure 1), BCs for the pore pressure perturbation solution (equation ( 9)) are:

1. p d 1 0; t ð Þ5 12c 1 ð Þr B , 2. Continuity of pressure across layers: p d 1 2h 1 ; t ð Þ1c 1 r B 5p d 2 2h 1 ; t ð Þ1 c 2 r B ,
3. Continuity of velocities across layers:

k 1 l 1 @p d 1 @z j z52h1 5 k 2 l 2 @p d 2 @z j z52h1 4. @p d 2 @z j z52h2 50,
where h 1 and h 2 denote the depth to the bottom of layers 1 and 2, respectively. The corresponding system of complex algebraic equations (combination of BCs and equation ( 9)) writes: :

1 1 0 0 e 2lph 1 =D 1 e lph 1 =D 1 2e 2lph 1 =D 2 2e lph 1 =D 2 lp D 1 k 1 l 1 e 2lph 1 =D 1 2 lp D 1 k 1 l 1 e lph 1 =D 1 2 lp D 2 k 2 l 2 e 2lph 1 =D 2 lp D 2 k 2 l 2 e lph 1 =D 2 0 0 e lph 2 =D
The complex coefficients a 1 ; b 1 ; a 2 ; and b 2 are then determined by solving the two complex algebraic equations systems above using MATLAB's Symbolic Math Toolbox.

Velocity

The perturbation in vertical fluid velocity u j z; t ð Þ À Á due to tidal loading can be modeled as a Darcy velocity driven by gradients in pore pressure perturbation (p j ):

u j z; t ð Þ52 k j l @p j @z : (11)
The analytical solution for velocity in a multilayer system is obtained through a spatial derivative of equation (9): u j z; t ð Þ52 k j lp lD j a j e l:p:z=D j 2b j e 2l:p:z=D j p T e ixt :

(12)

Note that the velocity perturbation at the seafloor (z50) is of particular interest since it can be measured by seafloor instruments such as a direct flow meter (e.g., [START_REF] Germanovich | Measuring fluid flow and heat output in seafloor hydrothermal environments[END_REF] or camera systems that use image analysis methods (e.g., [START_REF] Crone | Flow rate perturbations in a black smoker hydrothermal vent in response to a midocean ridge earthquake swarm[END_REF][START_REF] Mittelstaedt | A noninvasive method for measuring the velocity of diffuse hydrothermal flow by tracking moving refractive index anomalies[END_REF][START_REF] Mittelstaedt | Quantifying diffuse and discrete venting at the Tour Eiffel vent site, Lucky Strike hydrothermal field[END_REF][START_REF] Mittelstaedt | Diffuse venting at the ASHES hydrothermal field: Heat flux and tidally modulated flow variability derived from in situ time-series measurements[END_REF]. It is expressed as:

u 1 0; t ð Þ52 k 1 lp lD 1 a 1 2b 1 ð Þ p T e ixt : (13)
The phase lag of the velocity perturbation at the seafloor u 1 0; t ð Þ ð Þrelative to the tidal loading at the seafloor r B t ð Þ ð Þis given by:

Du u 5arg u 1 0; t ð Þ ð Þ 2arg r B t ð Þ ð Þ62np ½ ; (14)
where n can be any integer.

Temperature

So far, we have considered the pore pressure change p j z; t ð Þ À Á and velocity perturbation u j z; t ð Þ À Á due to tidal loading. However, recording time series of these quantities at the seafloor (i.e., in hydrothermal vent exit fluids) is technically challenging. By contrast, it is much simpler to measure temperature fluctuations in fluids discharging at the seafloor DT 1 0; t ð Þ ð Þ(e.g., [START_REF] Barreyre | Poroelastic response of mid-ocean ridge hydrothermal systems to ocean tidal loading: Implications for shallow permeability structure[END_REF]Barreyre et al., 2014a[START_REF] Barreyre | Temporal variability and tidal modulation of hydrothermal exit-fluid temperatures at the Lucky Strike deep-sea vent field, Mid-Atlantic Ridge[END_REF][START_REF] Fornari | Time-series temperature measurements at hightemperature hydrothermal vents, East Pacific Rise 9849 0 -51 0 N: Evidence for monitoring a crustal cracking event[END_REF][START_REF] Larson | Parameters of subsurface brines and hydrothermal processes 12-15 months after the 1999 magmatic event at the Main Endeavor Field as inferred from in situ time series measurements of chloride and temperature[END_REF][START_REF] Lilley | Magmatic events can produce rapid changes in hydrothermal vent chemistry[END_REF][START_REF] Scheirer | Temperature variations at diffuse and focused flow hydrothermal vent sites along the northern East Pacific Rise[END_REF][START_REF] Sohn | Stochastic analysis of exit fluid temperature records from the active TAG hydrothermal mound (Mid-Atlantic Ridge, 268N): 1. Modes of variability and implications for subsurface flow[END_REF]Tivey et al., 2002). We therefore complement our previous solutions with an analytical expression for temperature fluctuations due to tidal loading DT j z; t ð Þ À Á .

Our assumption of 1-D geometry means that we must ignore lateral flow and instead consider the idealized case of uniform vertical flow across the permeable layers, which is an acceptable approximation for a region of hydrothermal upwelling that underlies seafloor vent fields. We define a characteristic vertical velocity U j and temperature T j for the background upwelling flow that is perturbed by tidal loading, and introduce corresponding perturbation quantities (u j and DT j , respectively). We further assume a linear background (steady) temperature T j 5T 0 2C j z, with C the background temperature gradient [START_REF] Jupp | A poroelastic model for the tidal modulation of seafloor hydrothermal systems[END_REF], and approximate the steady upwelling velocity as U j 5

kj l j q sw 2q f ð Þ g
, where q sw is the density of background seawater and q f the hydrothermal fluid density.

Assuming that tidally induced changes in the fluid density and adiabatic cooling are negligible, the advection of temperature under tidal loading is expressed as:

@ T j 1DT j À Á @t 1r U j 1u j À Á@ T j 1DT j À Á @z 52rC j U j ; (15)
where the right side represents a 1-D approximation of the rate at which upwelling fluids are cooled conductively and adiabatically [START_REF] Jupp | A poroelastic model for the tidal modulation of seafloor hydrothermal systems[END_REF]. This equation can be linearized and rewritten as:

@DT j @z 1 1 rU j @DT j @t 5 C j U j u j : (16)
Assuming that the velocity perturbation solution is given by equation ( 12), we solve equation ( 16) and express the solution for temperature perturbation as follows:

DT j z; t ð Þ5C j e 2 ix rU j z e ixt 1A j e l:p:z D j 1B j e 2l:p:z

D j ; (17)
where A j and B j are two complex coefficients written as functions of the complex coefficients a j and b j defined earlier, and of various medium/fluid parameters:

A j 52 C j k j lp U j l j D j 1 lp Dj 1 ix rUj p T e ixt a j
and

B j 51 C j k j lp U j l j D j 1 2lp Dj 1 ix rUj p T e ixt b j :
The third coefficient C j is determined using the following BCs for temperature perturbation (equation ( 17)):

• DT 1 2h 1 ; t ð Þ5DT 2 2h 1 ; t ð Þ • DT 2 2h 2 ; t ð Þ50:
For the two-layer case with an impermeable lower boundary at z 5 -h 2 , the C coefficients are:

C 1 5 2e ix rU 2 h12h2 ð Þ A 2 e 2lph 2 D 2 1B 2 e lph 2 D 2 1A 2 e 2lph 1 D 2 1B 2 e lph 1 D 2 2A 1 e 2lph 1 D 1 2B 1 e lph 1 D 1 h i 3e 2 ix rU 1 h1 e 2ixt
and

C 2 52e 2 ix rU 2 h 2 A 2 e 2lph 2 D 2 1B 2 e lph 2 D 2 e 2ixt :
The two-layer temperature perturbation at the seafloor i:e:;

DT 1 0; t ð Þ5C 1 1A 1 1B 1 ð
Þ is then obtained by inserting expressions of C 1 , A 1 , and B 1 in equation ( 17), where z 5 0. This yields:

DT 1 0; t ð Þ5 C 1 2 C 1 k 1 lp U 1 l 1 D 1 1 lp D1 1 ix rU1 a 1 1 1 2lp D1 1 ix rU1 b 1 ! " # p T e ixt : (18)
Further, the phase lag of the temperature fluctuation at the seafloor

DT 1 0; t ð Þ ð Þrelative to the ocean tide r B t ð Þ ð
Þis given by:

Du T 5arg DT 1 0; t ð Þ ð Þ 2arg r B t ð Þ ð Þ 62np ½ ; (19)
where n can be any integer. This solution for phase lag is of particular interest and will be used in the remainder of this study to infer permeability layering in real systems, using time series of seafloor vent temperatures. Before doing so, however, we briefly describe how this phase lag varies as a function of layer geometry and permeability. In the following, we restrict ourselves to the two-layer case, which is appropriate for typical upper oceanic crust comprising an upper layer (2A, j 5 1 in Figure 1) of extrusives and a lower layer (2B, j 5 2 in Figure 1) of sheeted dikes. Figure 2b provides an example solution of phase lag Du T as a function of k 2A and k 2B , and compares it to the single layer solution (Figure 2a).

Effect of Heterogeneous Permeability

In order to gain insight into the underlying physics of poroelastic modulation due to periodic loading in a two-layer system, we define four distinct system behaviors (regimes in Figure 2c), which depend on both the permeability (k 2A , k 2B ) and geometry of the layers (thicknesses H 2A , H 2B ), as follows:

1. Regime 1 (R1) corresponds to low permeabilities in layer 2B, and a wide range of permeabilities in layer 2A. Low permeability makes layer 2B behave as an impermeable boundary at the base of 2A, which effectively decouples the two layers. This regime is thus analogous to a shallow single-layer system where the poroelastic behavior is entirely controlled by the layer 2A thickness and permeability (Figures 2a and2b). 2. Regime 2 (R2): Wide, intermediate range of permeabilities for both layer 2A and 2B. This transitional regime describes a system where the overall poroelastic behavior reflects strong coupling between the parameters of both layers, and maps in (k 2A , k 2B ) space (Figure 2c) as a ''funnel-shaped'' pattern in phase lag contours. Layer geometry exerts a strong modulation on this pattern (not shown in Figure 2c). The lower the ratio of H 2A /H 2B (i.e., H 1 /H 2 in Figure 2c), the more the funnel-shaped phase lag pattern is pushed toward greater layer 2A permeabilities, and vice versa. Thicker H 2B (i.e., H 2 in Figure 2c) shifts the funnel-shaped pattern toward greater layer 2B permeabilities, and vice versa. 3. Regime 3 (R3): Layer 2A and 2B have relatively similar, high permeability. This regime describes a uniform system where the poroelastic behavior is therefore controlled by the uniform permeability of layer 2A and 2B, as well as the overall thickness of the system, following the single layer solutions of Jupp and Schulz (2004). 4. Regime 4 (R4): A wide range of permeabilities for Layer 2B, and high permeabilities for layer 2A. Layer 2A has such a high permeability that it poses no resistance to the tidal modulation (signals diffuse through it very efficiently). This manifests as a decoupling of the two layers where a ''ghost'' layer 2A has very little influence. This regime describes a system where the poroelastic behavior is therefore controlled by the permeability and thickness of layer 2B.

This basic exploration of parameter space demonstrates that combining knowledge of the crustal structure (e.g., H 2A and H 2B from seismic profiles) with phase lag estimated between time series of in situ exit-fluid temperature and tidal pressure can place strong constraints on the poroelastic regime experienced by the layered system, and therefore on its first-order permeability structure.

While permeability controls the temperature-tide phase lag to first order, Figures 2d and2e further explore the sensitivity of phase lag to the geometry of the system (i.e., layer 2A and 2B thickness-Figure 2d) and the thermodynamic properties of the fluids (primarily modulated by the temperature and pressure of the fluid-Figure 2e, Table 2). This is done over a large range of layer 2A and 2B permeability spanning 6 orders of magnitude. In Figures 2d and2e, we specifically outline the combinations of (k 2A , k 2B ) that can explain two phase lag values: 2108 (red) and 1608 (blue), which represent two commonly observed values of phase lag reported by [START_REF] Barreyre | Poroelastic response of mid-ocean ridge hydrothermal systems to ocean tidal loading: Implications for shallow permeability structure[END_REF]. Different line shading and shapes represent different model assumptions (e.g., different fluid temperature/pressure, or layer 2A/2B thickness), which are detailed in the caption of Figures 2d and2e. Figure 2d shows the relative importance of the geometry of the system in controlling the phase lag under fixed fluid temperature and pressure. In particular, the thickness of layer 2A exerts an important control on phase lag, and increasing H 2A from 100 to 600 m (with every other parameter fixed) can significantly change the values of (k 2A , k 2B ) required to explain a given phase lag, primarily toward greater values of k 2A (by up to 2 orders of magnitude) with small changes in k 2B . Figure 2e illustrates the effect of fluid properties on phase lag for a fixed system geometry. The combinations of (k 2A , k 2B ) that are required to explain a given phase lag appear less sensitive to changes in fluid properties, as they change by less than 1 order of magnitude, when pressure-temperature conditions are varied within a reasonable range. Greater sensitivity is observed for temperatures approaching 4008C, which reflects the high sensitivity of fluid properties to changes in temperature and pressure near the triple point of water.

In summary, layer permeability and thickness, particularly in layer 2A (i.e., H 2A and k 2A ), control the temperature phase lag to first order, while fluid properties exert a second-order control. This exploration of parameter space gives us confidence that if layer geometry is well constrained by seismological techniques, measurements of temperature-tide phase lags can be used to constrain combinations of layers 2A and 2B permeabilities within $1 order of magnitude, and that greater phase lags generally provide a more robust constraint on k 2B .

3. Available Data Sets of Exit-Fluid Temperature and Pressure Lags

Obtaining Estimates of Temperature-Tide Phase Lags

We use the phase lag estimates for exit-fluid high temperature versus tidal pressure published by [START_REF] Barreyre | Poroelastic response of mid-ocean ridge hydrothermal systems to ocean tidal loading: Implications for shallow permeability structure[END_REF] as the basis for our analyses. These phase lag estimates represent the most robust values obtained from a comprehensive analysis of the complete set of publicly available (Marine Geoscience Data System, PANGAEA, European Multidisciplinary Seafloor and water column Observatory, and Ocean Networks Canada) exit-fluid high-temperature time series data for vents at the Lucky Strike Hydrothermal Field (LSHF), the Main Endeavour Field (MEF), and the East Pacific Rise (EPR) 9846 0 N sites (e.g., [START_REF] Barreyre | Temporal variability and tidal modulation of hydrothermal exit-fluid temperatures at the Lucky Strike deep-sea vent field, Mid-Atlantic Ridge[END_REF][START_REF] Fornari | Time-series temperature measurements at hightemperature hydrothermal vents, East Pacific Rise 9849 0 -51 0 N: Evidence for monitoring a crustal cracking event[END_REF][START_REF] Larson | Parameters of subsurface brines and hydrothermal processes 12-15 months after the 1999 magmatic event at the Main Endeavor Field as inferred from in situ time series measurements of chloride and temperature[END_REF][START_REF] Scheirer | Temperature variations at diffuse and focused flow hydrothermal vent sites along the northern East Pacific Rise[END_REF]Tivey et al., 2002). These hydrothermal sites are located on MORs spanning a large range of spreading rates, from slow to fast (MAR: 2.2 cm yr 21 [START_REF] Cannat | Mid-Atlantic Ridge-Azores hotspot interactions: Along-axis migration of a hotspot-derived magmatic pulse 14 to 4 ma ago[END_REF]; JdFR: 6 cm yr 21 [START_REF] Riddihough | Recent movements of the Juan de Fuca Plate System[END_REF]; EPR: 11 cm yr 21 [START_REF] Carbotte | Comparison of sea floor tectonic fabric created at intermediate, fast and superfast spreading ridges: Influence of spreading rate, plate motions, and ridge segmentation on fault patterns[END_REF]).

Coherency and phase lag between tidal pressure at the seafloor (i.e., data generated with the GOT4.7 global ocean tide model [START_REF] Ray | Precise comparisons of bottom-pressure and altimetric ocean tides[END_REF] as described in [START_REF] Barreyre | Poroelastic response of mid-ocean ridge hydrothermal systems to ocean tidal loading: Implications for shallow permeability structure[END_REF], and exit-fluid high temperatures are estimated by applying multitaper [START_REF] Thomson | Spectrum estimation and harmonic analysis[END_REF] cross-spectral methods with adaptive weighting (Percival & Walden, 1993) to time series data (see details in [START_REF] Barreyre | Poroelastic response of mid-ocean ridge hydrothermal systems to ocean tidal loading: Implications for shallow permeability structure[END_REF]. After applying stringent 2004). e [START_REF] Vera | The structure of 0-to 0.2-m.y.-old oceanic crust at 98N on the East Pacific Rise from expanded spread profiles[END_REF]. f [START_REF] Luyendyk | On-bottom gravity profile across the East Pacific Rise crest at 218 north[END_REF]. g [START_REF] Becker | Large-scale electrical resistivity and bulk porosity of the oceanic crust, Deep Sea Drilling Project Hole 504B. Costa Rica Rift[END_REF]. h [START_REF] Carmichael | Handbook of physical properties of rocks[END_REF] and [START_REF] Christensen | Sea floor spreading, progressive alteration of layer 2 basalts, and associated changes in seismic velocities[END_REF]. i [START_REF] Pros | The velocity of elastic waves in diabase and greywacke under pressure up to 4 kilobars, Stud[END_REF]. j Computed from the equations of state for pressure at 330 bar and hydrothermal fluid temperature at 3508C.

Geochemistry, Geophysics, Geosystems 10.1002/2017GC007152 criteria for phase analysis (e.g., coherency c 2 ð Þ ! 0:85) to maximize the signal-to-noise ratio, [START_REF] Barreyre | Poroelastic response of mid-ocean ridge hydrothermal systems to ocean tidal loading: Implications for shallow permeability structure[END_REF] provided a highly coherent and stable data set of phase lag estimates consisting of two records ($2.4 years of recording) for the EPR-9846 0 N, 30 records ($27 years of recording) for the LSHF, and two records ($0.8 years of recording) for the MEF. The corresponding average values of phase lag and the associated uncertainty for each site/time series are listed in Table 3.

3.2. Background Information on the Hydrothermal Sites From Our Compilation 3.2.1. East Pacific Rise 9850 0 N Field (EPR-9850 0 N) The 9850 0 N hydrothermal field is located on a volcanically active segment of the fast-spreading EPR (11 cm yr 21 , [START_REF] Carbotte | Comparison of sea floor tectonic fabric created at intermediate, fast and superfast spreading ridges: Influence of spreading rate, plate motions, and ridge segmentation on fault patterns[END_REF]) that has experienced recent volcanic/magmatic activity in 1991 [START_REF] Haymon | Volcanic eruption of the mid-ocean ridge along the East Pacific Rise crest at 9845-52 0 N: Direct submersible observations of seafloor phenomena associated with an eruption event in April, 1991[END_REF][START_REF] Rubin | 210 Po-210 Pb dating of recent volcanic eruptions on the sea floor[END_REF]), 2005-2006[START_REF] Tolstoy | A sea-floor spreading event captured by seismometers[END_REF], and possibly in 1995 [START_REF] Germanovich | Magmatic origin of hydrothermal response to earthquake swarms: Constraints from heat flow and geochemical data[END_REF]. Because a considerable amount of seismic work has been carried out at the EPR (e.g., [START_REF] Detrick | Multi-channel seismic imaging of a crustal magma chamber along the East Pacific rise[END_REF][START_REF] Sohn | Fine-scale seismic structure of the shallow volcanic crust on the East Pacific Rise at 9850 0 N[END_REF][START_REF] Xu | Variations in axial magma lens properties along the East Pacific Rise (9830 0 N-10800 0 N) from swath 3-D seismic imaging and 1-D waveform inversion[END_REF], we have reliable constraints on the depth of the axial melt lens (AML) beneath the hydrothermal field ($1,500 m, [START_REF] Detrick | Multi-channel seismic imaging of a crustal magma chamber along the East Pacific rise[END_REF] and on the relatively low thickness of layer 2A ($150 m, [START_REF] Christeson | Extrusive thickness variability at the East Pacific Rise, 98-108N: Constraints from seismic techniques[END_REF][START_REF] Harding | A multichannel seismic investigation of upper crustal structure at 98N on the East Pacific Rise: Implications for crustal accretion[END_REF][START_REF] Sohn | Fine-scale seismic structure of the shallow volcanic crust on the East Pacific Rise at 9850 0 N[END_REF][START_REF] Vera | Seismic imaging of oceanic layer 2A between 9830 0 N and 108N on the East Pacific Rise from two-ship wide-aperture profiles[END_REF]. We have phase lag data for a specific vent within the 9850 0 N field: L-vent. Cross-spectral analysis yields phase lag estimates at this vent of 2078 6 58 at M2 frequency and 1888 6 78 at K1 frequency [START_REF] Barreyre | Poroelastic response of mid-ocean ridge hydrothermal systems to ocean tidal loading: Implications for shallow permeability structure[END_REF].

Lucky Strike Hydrothermal Field (LSHF)

The Lucky Strike hydrothermal field is located above a volcanic center on the slow-spreading Mid-Atlantic Ridge (MAR, 2.2 cm yr 21 ; [START_REF] Cannat | Mid-Atlantic Ridge-Azores hotspot interactions: Along-axis migration of a hotspot-derived magmatic pulse 14 to 4 ma ago[END_REF], and provides the most extensive set of hydrothermal discharge time series data acquired to date. Seismic studies have also provided strong constraints for both the magma chamber depth ($3,400 m) and the thickness of the extrusive layer 2A ($600 m) [START_REF] Arnulf | Strong seismic heterogeneity in layer 2A near hydrothermal vents at the Mid-Atlantic Ridge[END_REF][START_REF] Crawford | Hydrothermal seismicity beneath the summit of Lucky Strike volcano, Mid-Atlantic Ridge[END_REF][START_REF] Seher | Seismic layer 2A variations in the Lucky Strike segment at the Mid-Atlantic Ridge from reflection measurements[END_REF][START_REF] Singh | Discovery of a magma chamber and faults beneath a Mid-Atlantic Ridge hydrothermal field[END_REF]. It has been shown that there exists a difference in phase lag between vents located on the west (1558 6 58 at M2 frequency-no estimate at K1 frequency) versus the east (1738 6 48 at M2 frequency and 1688 6 78 at K1 frequency) side of the LSHF [START_REF] Barreyre | Poroelastic response of mid-ocean ridge hydrothermal systems to ocean tidal loading: Implications for shallow permeability structure[END_REF]. The two sets of vents are separated by a solidified lava lake ($ 300 m in diameter) that forms an impermeable cap on the extrusive layer [START_REF] Barreyre | Structure, temporal evolution, and heat flux estimates from the Lucky Strike deep-sea hydrothermal field derived from seafloor image mosaics[END_REF][START_REF] Fouquet | Atlantic lava lakes and hot vents[END_REF][START_REF] Ondr Eas | Recent volcanic events and the distribution of hydrothermal venting at the Lucky Strike hydrothermal field[END_REF]. Note that we had to lower the coherency threshold to c 2 ð Þ ! 0:7 at the K1 frequency for LSHF-E in order to obtain a phase lag estimate. While less stringent than in all other sites of the compilation, this criterion still ensures a coherent and stable result.

Main Endeavour Field (MEF)

The Main Endeavour Field is located on the intermediate-spreading Juan de Fuca Ridge (JdFR, 6 cm yr 21 ; Riddihough, 1984) at 478N. Seismic studies provide estimates for both the magma chamber depth ($2,300 m) and the extrusive layer 2A thickness ($460 m) [START_REF] Van Ark | Seismic structure of the Endeavour segment, Juan de Fuca Ridge: Correlations with seismicity and hydrothermal activity[END_REF]. Intrafield variability is also evident at the MEF, where the Grotto site exhibits phase lags of 2178 6 2.58 at M2 frequency (and 2018 6 48 at K1), which exceeds the phase lag at the S&M site by $408 (178.58 6 1.58 at M2 frequency and 1608 6 58 at K1). These differences are much larger than the phase lag uncertainties and intersite variability, indicating that they arise from deterministic differences in the subsurface permeability structure.

Table 3 Average Phase Lag ( /) and Errors ( e / ) Estimates at M2 Semidiurnal Frequency and K1 Diurnal Frequency Estimated From Spectral Analysis Following the Methodology Described in [START_REF] Barreyre | Poroelastic response of mid-ocean ridge hydrothermal systems to ocean tidal loading: Implications for shallow permeability structure[END_REF] Hydrothermal fields Geochemistry, Geophysics, Geosystems 10.1002/2017GC007152

/M2 ð Þ e /M2 ð Þ /K1 ð Þ e /K1 ð Þ H 2A m ð Þ H 2B m ð Þ k 2A m 2 ð Þ k 2B m 2 ð Þ LSHF (MAR)

Results

It is important to note that the permeabilities estimated below and discussed throughout this paper are interpreted as corresponding specifically to discharge zone permeabilities, and therefore, cannot be interpreted as whole crustal permeabilities. In layer 2A the flow to a vent field is likely the result of flow within a limited volume and thus, may not be representative of layer 2A as a whole [START_REF] Cann | Modeling periodic megaplume emission by black smoker systems[END_REF]Fontaine et al., 2007), since the very presence of the vent field could result from an anomalously permeable underlying crust (e.g., [START_REF] Crone | Permeability structure of young ocean crust from poroelastically triggered earthquakes[END_REF].

4.1. Inverting for Discharge Zone Permeability Layering at Hydrothermal Fields Using the model described above, we computed the predicted phase lag of vent temperature relative to tidal pressure at the seafloor as a function of the permeabilities of layer 2A and 2B (Figure 3). We assumed that the respective thickness of the layers beneath each hydrothermal site was known from previous seismic studies summarized in the previous section, and in Table 3. We then used our phase lag estimates at both the M2 and K1 frequencies to constrain the model and reduce the permissible solution space to narrow contours (black contours in Figure 3) for each site.

Figure 3 highlights a key dichotomy previously reported by [START_REF] Barreyre | Poroelastic response of mid-ocean ridge hydrothermal systems to ocean tidal loading: Implications for shallow permeability structure[END_REF], between high layer 2A permeability systems (k 2A > $10 211 m 2 : LSHF-W and LSHF-E and MEF-S&M), and low layer 2A permeability systems (k 2A < $10 211 m 2 : EPR-Lvent and MEF-Grotto). Our two-layer model, however, reveals previously unknown effects due to coupling with a less permeable layer 2B (e.g., bending of the solution contour for higher layer 2B permeability). Overall, we find that the phase lag between exit-fluid temperature and tidal loading is relatively insensitive to the permeability of layer 2B, with the range of permeabilities that fit the phase lag data spanning several orders of magnitude (i.e., up to 3).

In order to further constrain our permeability estimates, we only retained combinations of (k 2A , k 2B ) that could jointly explain the observed phase lag at both the M2 and K1 tidal frequencies (except for LSHF-W, where robust phase lag estimates at the K1 frequency could not be obtained). Final constrained permeability ranges for layers 2A and 2B are shown in Figure 4 and summarized in Table 3. Our estimates are in broad agreement with previous studies, but provide further insight into the systematic variability of permeability layering at MOR hydrothermal sites (Figure 4b).

At the EPR, where the base of layer 2A is located at the relatively shallow depth of $155 mbsf [START_REF] Sohn | Fine-scale seismic structure of the shallow volcanic crust on the East Pacific Rise at 9850 0 N[END_REF], the phase lag estimates for L-vent (2078658 @ M2 and 1888678 @ K1) require a permeability for layer 2A ranging from $1.5 3 10 213 to 2.5 3 10 213 m 2 and from $10 215 to 2 3 10 214 m 2 for layer 2B.

By contrast, at the LSHF the smaller phase lags (1558 6 58 @ M2 for LSHF-W and 1738 6 48 @ M2 and 1688 6 78 @ K1 for LSHF-E) require a much higher effective permeability for layer 2A ranging from $1.5 3 10 210 to 10 29 m 2 for LSHF-W and $7 3 10 211 to 2 3 10 210 m 2 for LSHF-E. They also constrain layer 2B permeability between $10 215 and 1.5 3 10 212 m 2 for LSHF-W and between $10 215 and 5 3 10 213 m 2 for LSHF-E. Intermediate behaviors are observed at the MEF where the extrusive layer thickness is $450 m [START_REF] Van Ark | Seismic structure of the Endeavour segment, Juan de Fuca Ridge: Correlations with seismicity and hydrothermal activity[END_REF]. At Grotto vent, the phase lags of 2178 6 2.58 @ M2 and 2018 6 48 @ K1 require an effective permeability ranging from $6 3 10 213 to 8 3 10 213 m 2 for layer 2A and from $10 215 to 7 3 10 214 m 2 for layer 2B, whereas the smaller phase lags observed for S&M vent (178.58 6 1.58 @ M2 and 1608 6 58 @ K1) require an effective permeability of $1.5 3 10 210 m 2 for layer 2A and ranging from $10 215 to 1.5 3 10 213 m 2 for layer 2B.

Note that the estimated permeability ranges above are bounded by a minimum permeability of 10 215 m 2 for layer 2B. This cutoff permeability of 10 215 m 2 is an outcome of the coupling between layers in our model, and corresponds to the minimum permeability at which layer 2B exerts some influence on the observed phase lag. For any lower k 2B , layer 2B becomes too impermeable and the model behaves like a single-layer, homogenous model with an effectively impermeable boundary at the base of layer 2A (Figures 2a and2b). In section 5.3, we propose a lower bound on k 2B using heat flow constraints.

Effect of Layer Thickness on the Poroelastic Response at the Field Scale

In the previous section, we estimated the permeability structure for several hydrothermal fields that have well-constrained layer thicknesses from seismic studies. Here, we focus on variability within individual fields, and the potential contribution of changes in layer 2A/2B thickness to this variability. At the LSHF, the difference in phase lags estimated for the western (1558 6 58 @ M2 for LSHF-W) versus eastern Geochemistry, Geophysics, Geosystems (1738 6 48 @ M2 and 1688 6 78 @ K1 for LSHF-E) vents is consistent with a change in the extrusive layer thickness that has been imaged seismically (300 m for the west vents versus 600 m for the east vents, [START_REF] Arnulf | Strong seismic heterogeneity in layer 2A near hydrothermal vents at the Mid-Atlantic Ridge[END_REF]. Figure 5 shows the effect of layer 2A thickness on the M2 phase lag assuming a permeability of 10 214 m 2 in layer 2B and a range of permeabilities for layer 2A. Systematic phase lag differences are observed between vent sites hosted on the east versus west side of the LSHF, which can be explained entirely by a change in 2A thickness and requires no change in 2A permeability.

On the other hand, intermediate behaviors are observed at the MEF where the extrusive layer thickness is $450 m [START_REF] Van Ark | Seismic structure of the Endeavour segment, Juan de Fuca Ridge: Correlations with seismicity and hydrothermal activity[END_REF].

Since the seismic data indicate the extrusive layer thickness is the same for both of these vents [START_REF] Van Ark | Seismic structure of the Endeavour segment, Juan de Fuca Ridge: Correlations with seismicity and hydrothermal activity[END_REF], at the MEF the phase lag difference between the two sites appears to result from a variation in layer 2A permeability rather than layer 2A thickness. Magnetic data suggest that the S&M and Grotto vents, which are separated by a distance of $150 m, are fed by distinct fluid upwelling zones [START_REF] Tivey | Crustal magnetization reveals subsurface structure of Juan de Fuca Ridge hydrothermal vent fields[END_REF], and our results suggest these zones have different layer 2A permeabilities. properties, which precludes any tidally induced horizontal flow. Secondly, the use of a multilayer model with uniform permeability for each layer does not fully account for the progressive decrease in permeability with depth (e.g., cracks progressively closing with increasing depth, [START_REF] Carlson | The influence of porosity and crack morphology on seismic velocity and permeability in the upper oceanic crust[END_REF] nor heterogeneities. Thirdly, the model assumes that the hydrothermal fluid and the matrix are in thermal equilibrium at all depths. Fourthly, we assume that the steady state temperature gradient C results from adiabatic cooling alone (because of depressurization on ascent-i.e., without any heat exchange between the upwelling zone and adjacent rocks), because lateral conductive heat loss cannot be incorporated into a 1-D model, and vertical conductive heat loss is negligible in a thermal plume [START_REF] Phillips | Flow and reactions in permeable rocks[END_REF]. Finally, we assume that the properties of the interstitial fluid are those of pure (liquid) water. These last three assumptions preclude a discussion of how the temperature distribution inside a real convection cell might complicate the flow.

Discussion

The first two assumptions dealing with spatial homogeneity are common to most permeability estimates and models, even though the matrix hosting hydrothermal circulation at MORs is fractured, faulted, fissured, and brecciated to high degrees (e.g., [START_REF] Escart In | Interplay between faults and lava flows in construction of the upper oceanic crust: The East Pacific Rise crest 9825 0 -9858 0[END_REF][START_REF] Wright | Crustal fissuring on the crest of the southern East Pacific Rise at 17815'-40[END_REF], which can impart significant heterogeneity to the permeability structure. To this point, most models, including that employed here, address this issue by considering length scales much larger than the heterogeneity scale, which allows the system to be characterized by an effective permeability that implicitly averages heterogeneities. Incorporating heterogeneous permeability fields into flow models is complex, computationally expensive, and more importantly, ad hoc, lacking in situ reliable constraints (particularly at depth). On the other hand, incorporating multiple layers into poroelastic flow models is more straightforward. The third assumption regarding thermal equilibrium is necessary, because we are using exit-fluid temperature fluctuations (as opposed to velocity fluctuations, for example) to constrain permeability [START_REF] Jupp | A poroelastic model for the tidal modulation of seafloor hydrothermal systems[END_REF] and is justified by the fact that the thermal diffusion length scale over a tidal cycle ($16 cm @ M2 and $20 cm @ K1) is less than the spacing of flow paths (i.e., faults and cracks) (e.g., [START_REF] Bohnenstiehl | Faulting patterns near 19830 0 S on the East Pacific rise: Fault formation and growth at a superfast spreading center[END_REF][START_REF] Wright | Crustal fissuring on the crest of the southern East Pacific Rise at 17815'-40[END_REF]. The fourth assumption is valid in situations where there is little-to-no heat exchange between upwelling fluids and the surrounding matrix, which appears to be a valid assumption for the hydrothermal fields discussed in this study based on exit-fluid chemistry (e.g., [START_REF] Langmuir | Hydrothermal vents near a mantle hot spot: The lucky strike vent field at 378N on the Mid-Atlantic Ridge[END_REF][START_REF] Pester | Subseafloor phase equilibria in high-temperature hydrothermal fluids of the Lucky Strike Seamount (Mid-Atlantic Ridge, 37817 0 N)[END_REF]. The last assumption regarding fluid properties is justified as long as there are no phase changes in the tidally pumped flow regime because the fluid dynamic properties of hydrothermal fluids are similar to those of pure water.

Measured Range of Layers 2A and 2B Permeability and Comparison With Previous Work

Hydrothermal systems exhibiting phase lags between tidal loading and exit-fluid temperature greater than $2008 are generally well explained by a poroelastic model in which the shallow (extrusive) crust has a moderate permeability (k 2A $ 10 213 -10 212 m 2 ). By contrast, systems with phase lags 1808 require a higher layer 2A permeability (k 2A ! 10 210 m 2 ). The dichotomy observed in layer 2A arises from the sensitivity of phase lag to extrusive layer permeability (Figures 2c and2d), and suggests that phase lag constitutes a firstorder proxy for the permeability of the shallowest portion of the upflow zone underlying a given hydrothermal site. This is consistent with the results of [START_REF] Barreyre | Poroelastic response of mid-ocean ridge hydrothermal systems to ocean tidal loading: Implications for shallow permeability structure[END_REF].

High permeability in the shallow crust (!10 210 m 2 )-as we infer beneath Lucky Strike and the S&M site of the MEF-has previously been proposed for the Endeavour segment near the Raven field $1 km north of the MEF through a quantitative analysis of seafloor fissuring [START_REF] Hearn | Surficial permeability of the axial valley seafloor: Endeavour Segment, Juan de Fuca Ridge[END_REF]. Such high permeabilities are similar to those measured through borehole flow calculations in shallow subaerial Hawaiian basalts (10 211 -10 29 m 2 ) [START_REF] Ingebritsen | The hydrogeology of Kilauea volcano[END_REF], and upper igneous oceanic crustal units (Becker & Davis, 2004, and references therein). They are similar to in situ permeability measurements in shallow submarine wellbores (yielding values as large as 10 211 m 2 ) [START_REF] Fisher | Permeability within basaltic oceanic crust[END_REF], and in wellbores from sedimented ridge flanks of both the MAR and the JdFR ($10 210 m 2 ) [START_REF] Davis | Formation-scale hydraulic and mechanical properties of oceanic crust inferred from pore pressure response to periodic seafloor loading[END_REF]. Permeabilities inferred from pressure diffusion along borehole transects of CORKs (i.e., Circulation Obviation Retrofit Kits) yield similarly high values (10 210 -10 29 m 2 ) [START_REF] Davis | An episode of seafloor spreading and associated plate deformation inferred from crustal fluid pressure transients[END_REF], and so do large-scale numerical models of fluid flow and heat transport (10 211 -10 29 m 2 ) (e.g., [START_REF] Stein | Observations and models of lateral hydrothermal circulation on young ridge flank: Numerical evaluation of thermal and chemical constraints[END_REF], and ophiolite studies (10 212 -10 28 m 2 , [START_REF] Nehlig | Flow porosities, permeabilities and preliminary data on fluid inclusions and fossil thermal gradients in the crustal sequence of the Sumail ophiolite (Oman)[END_REF][START_REF] Van Everdingen | Fracture characteristics of the sheeted dike complex, Troodos ophiolite, Cyprus: Implications for permeability of oceanic crust[END_REF].

Interestingly, our model predicts lower 2A permeability beneath the Grotto site ($10 212 m 2 ), located between S&M and Raven. This suggests that layer 2A permeability can vary substantially from one discharge zone to the next, even zones separated by less than a km. Our model also predicts lower permeability ($10 213 m 2 ) in layer 2A beneath EPR's L-vent (9846 0 N), which is in agreement with the 10 213.2 m 2 estimated by [START_REF] Marjanovic | Seismic signatures of hydrothermal pathways along the East Pacific Rise between 9816' and 9856 0 N[END_REF] from seismic velocities and porosity/permeability relationships [START_REF] Carlson | The influence of porosity and crack morphology on seismic velocity and permeability in the upper oceanic crust[END_REF]. This contrasts with the high permeability values ($10 210 -10 29 m 2 ) proposed by [START_REF] Crone | Permeability structure of young ocean crust from poroelastically triggered earthquakes[END_REF] for the upflow zone beneath the 9850 0 N cluster of vents located $10 km north of L-vent. This estimate was obtained by assuming that the timing of microearthquakes directly reflects the diffusion of poroelastic stresses related to tidal cycles. This approach yielded a ''background'' permeability of the crust close to 10 212 m 2 , and locally as low as 10 214 m 2 . Significantly higher permeability ($10 210 -10 29 m 2 ) was found within narrow areas of hydrothermal upflow and downflow, which highlighted the high lateral heterogeneity of the crust on length scales < 3 km. This model, however, as the one used by [START_REF] Marjanovic | Seismic signatures of hydrothermal pathways along the East Pacific Rise between 9816' and 9856 0 N[END_REF] did not make a distinction between layers 2A and 2B permeability.

The more moderate permeability values in the 10 213 -10 212 m 2 range, which we propose for layer 2A beneath the L-vent (EPR) and Grotto (MEF) sites are consistent with estimates of 3 3 10 213 -6 3 10 212 m 2 derived by modeling flow rate perturbations generated by earthquake swarms near the MEF [START_REF] Crone | Flow rate perturbations in a black smoker hydrothermal vent in response to a midocean ridge earthquake swarm[END_REF]. They also fall within the range determined through drill-string packer experiments at the eastern flank of the JdFR (3 3 10 214 -2 310 211 m 2 ; [START_REF] Becker | Permeability of upper oceanic basement on the eastern flank of the Juan de Fuca Ridge determined with drill-string packer experiments[END_REF].

Our two-layer model places less stringent constraints on layer 2B permeability, but does provide reliable upper bounds that increase with decreasing spreading rate: $2 3 10 214 , $10 213 , and $10 212 m 2 at L-vent (EPR), MEF, and Lucky Strike, respectively (Figure 4). These lower values are closer to (but generally slightly lower than) the estimates of [START_REF] Lowell | Characteristics of magma-driven hydrothermal systems at oceanic spreading centers[END_REF] who used a simple single-pass convection model to infer the effective, depth-averaged permeability of discharge zones based on heat flow estimates at hydrothermal sites. This suggests that the effective permeability that sets the pace of heat extraction represents some average of layers 2A and 2B permeability that is primarily influenced by layer 2B [START_REF] Rosenberg | The relationship between flow and permeability field in seafloor hydrothermal systems[END_REF]. We further explore this notion in the following section. It should of course be noted that all the above permeability estimates, including ours, were derived by fitting a model to an observable (e.g., phase lag, heat flow, earthquake distribution). The accuracy of these estimates largely relies on the model applicability and accuracy.

Reconciling Permeability Estimates From Tidal Modulation With the Heat Output of Hydrothermal Fields

By definition, permeability sets the Darcy velocity of hydrothermal fluids subjected to a given pressure gradient (equation ( 11)). Permeability thus exerts a very strong modulation on the thermo-chemical fluxes associated with hydrothermal circulation [START_REF] Lowell | Hydrothermal processes at mid-ocean ridges: Results from scale analysis and single-pass models[END_REF], as well as on venting temperature [START_REF] Driesner | The interplay of permeability and fluid properties as a first order control of heat transport, venting temperatures and venting salinities at mid-ocean ridge hydrothermal systems[END_REF]. Here we explore the possibility of further constraining layers 2A and 2B permeability by constructing a model that jointly accounts for (1) the tidal response of hydrothermal systems, (2) their steady-state heat flow, and (3) the associated venting temperatures. We adapt the theoretical framework developed by [START_REF] Driesner | The interplay of permeability and fluid properties as a first order control of heat transport, venting temperatures and venting salinities at mid-ocean ridge hydrothermal systems[END_REF] to calculate the heat flow that would be transported by upwelling fluids in a stratified hydrothermal discharge zone of known k 2A and k 2B . In our conceptual model, buoyant upflow is fueled by an imposed heat flux resulting from the accretion of hot, crystallizing crust, and leads to venting of a total flux Q at the seafloor. In low-permeability layer 2B, the upflow is assumed to occur within a cylindrical zone of radius R 2B . In layer 2A, the upflow is likely faster and thus occupies a narrower area (radius R 2A ) [START_REF] Rosenberg | The relationship between flow and permeability field in seafloor hydrothermal systems[END_REF]. Under the assumption of steady-state heat transfer, the heat flow through layer 2A (q 2A , in W m 22 ) and 2B (q 2B ) must match the venting flux (Q 5 pR 2 2A q 2A 5 pR 2 2B q 2B ). We write T H and P H the temperature and pressure of fluids at the base of layer 2B, and assume that the properties of fluids upwelling throughout layer 2B are well described by that of seawater at (T H , P H ). Likewise, we assume that the temperature (T* T H ) and pressure (P* P H ) that characterize the interface between layers 2B and 2A provide a reasonable approximation of the properties of fluids throughout layer 2A, and that T* can be considered a reasonable proxy for the seafloor venting temperature.

The temperature difference driving the flow through layer 2A is that between the hot fluid at temperature T* and the ocean (T C , at seafloor pressure P C ). The heat flow through layer 2A thus writes:

q 2A 5 k 2A g q C 2q à ð Þ l à q à h à 2h C ð Þ ; ( 20
)
where l is the dynamic viscosity of the fluid, q is its density, and h is its enthalpy (in JÁkg 21 ). Likewise, the heat flow through layer 2B is

Geochemistry, Geophysics, Geosystems

10.1002/2017GC007152 q 2B 5 k 2B g q à 2q H ð Þ l H q H h H 2h C ð Þ : (21)
In equations ( 20) and ( 21), subscripts ''H'', ''*'', and ''C'' refer to pressure-temperature conditions (P H , T H ), (P*, T*), and (P C , T C ).

For simplicity, we assume that T H is known, and for a given Q solve equations ( 20) and ( 21) to obtain combinations of k 2A and k 2B that satisfy Q 5 pR 2 2A q 2A 5 pR 2 2B q 2B . Each (k 2A , k 2B ) pair is associated with a unique value of T* in the region of parameter space considered. This procedure is carried out using the tabulated thermodynamic properties of pure water and steam (IAPWS-IF97, [START_REF] Wagner | The IAPWS industrial formulation 1997 for the thermodynamic properties of water and steam[END_REF], which constitute a reasonable approximation for hydrothermal fluid properties. Viscosity is estimated based on the tables of International Association for the Properties of Water and Steam (2003; Revised release on the IAPS formulation 1985 for the viscosity of ordinary water substance). Seafloor conditions are taken as T C 5 48C and P C 5 230 bar. P* and P H are chosen as 244 and 299 bar, respectively, but their exact value has little influence on the result compared to the choice of T H . We choose T H 5 4508C but recognize that any value between 420 and 6008C would be equally acceptable as modeling studies have shown that this range of basal temperatures is compatible with observed seafloor venting temperatures (Fontaine & Wilcock, 2007;[START_REF] Fontaine | Permeability changes due to mineral diagenesis in fractured crust: Implications for hydrothermal circulation at mid-ocean ridges[END_REF]. Finally, to simplify our exploration of a vast parameter space, we consider two cases with R 2B 5 500 and 100 m, and assume R 2A 5 R 2B /10.

Using this relatively simple model, we constrain combinations of k 2A and k 2B that are consistent with three values of Q which we consider representative of EPR's L-vent (160 MW), MEF (450 MW), and Lucky Strike (600 MW) (Lowell et al., 2013 and references therein). Of course, these estimates have large error bars, potentially larger than 650%. Further, we must use the estimated heat output of the EPR 9850 0 N hydrothermal field as a proxy for heat flow in the neighboring L-vent area, as no estimate of the heat output of that region is currently available. The exercise attempted here should therefore be considered as a proof-ofconcept for a methodology that reconciles the permeability constraints of tidal modulation, steady-state heat transfer and venting temperature. It is however useful in this context to compare the characteristics of various hydrothermal sites in a quantitative fashion.

Figure 6 shows two sets of curves indicating which permeability combinations can explain the observed heat outputs under two distinct assumptions for the radius of the cylindrical upflow zone in layer 2B: R 2B 5 500 and 100 m. Focusing first on the R 2B 5 500 m case, which is consistent with the structure of upflow zones inferred from microseismicity studies at EPR 9850 0 N [START_REF] Marjanovic | Seismic signatures of hydrothermal pathways along the East Pacific Rise between 9816' and 9856 0 N[END_REF][START_REF] Tolstoy | Seismic identification of along-axis hydrothermal flow on the East Pacific Rise[END_REF], it is apparent that heat flow is primarily sensitive to the permeability of layer 2B, which is consistent with earlier studies (e.g., [START_REF] Rosenberg | The relationship between flow and permeability field in seafloor hydrothermal systems[END_REF]. Further, this permeability is likely greater than 10 215 m 2 as a lower value would not permit the observed heat fluxes. At L-vent, a permeability of 6 3 10 215 m 2 in layer 2B and 2 3 10 213 m 2 in layer 2A (red star in Figure 6) would enable a heat flow of 160 MW over a 50 m radius upwelling zone in layer 2A, while also producing the observed phase lag between tidal loading and vent temperature. Interestingly, a much higher permeability ($10 210 m 2 ) as suggested by [START_REF] Crone | Permeability structure of young ocean crust from poroelastically triggered earthquakes[END_REF] throughout the crust beneath EPR 9850 0 N would likely result in unrealistically high heat flow and potentially lower temperature discharge [START_REF] Driesner | The interplay of permeability and fluid properties as a first order control of heat transport, venting temperatures and venting salinities at mid-ocean ridge hydrothermal systems[END_REF]. More targeted studies of the fine permeability structure of EPR will be required to assess the applicability, sensitivity, and potential biases of each approach.

Similarly, k 2B 5 10 214 m 2 and k 2A 5 7 3 10 213 m 2 would be a plausible combination for the Grotto site of the MEF. Our model is in both cases compatible with a shallow upwelling temperature in excess of 3008C, which is consistent with high-temperature venting occurring at both sites. By contrast, models satisfying heat flow and tidal modulation constraints at S&M and Lucky Strike require much greater layer 2A permeability (on the order of 10 210 m 2 ), and layer 2B permeabilities similar to those inferred beneath L-vent and Grotto (between 10 214 and 10 215 m 2 ). However, such models would predict moderate-to-low venting temperatures (<1508C) instead of the high temperatures observed at these sites. Alternatively, one can assume a narrower upwelling zone throughout the crust (R 2B 5 100 m), in which case models can jointly explain the phase lag, heat output and venting temperatures in excess of 2008C with higher k 2B values ($2 3 10 213 m 2 , blue and green stars in Figure 6).

Large uncertainties on T H , Q, and the size of crustal-scale upflow zones preclude any definitive estimate of k 2A and k 2B at the hydrothermal sites considered. However, the methodology outlined here could be applied at a more local scale beneath a particular venting site, where permeability layering could be better constrained. In the context of our study, this method provides useful guidelines to evaluate the plausibility of our permeability estimates. The most important one is that layer 2B permeability is likely greater than Geochemistry, Geophysics, Geosystems 10.1002/2017GC007152 $10 215 m 2 . This lower bound common to all sites provides a reasonably narrow range of permeability for layer 2B beneath L-vent at the EPR: 10 215 -2 3 10 214 m 2 . Figure 6 shows that layer 2B permeability could fall within that same low range at all other sites and satisfy the heat flow constraint. In that scenario, the variability in permeability structure from site to site-and potentially across spreading rates-would be primarily confined to layer 2A, with moderate-to-low k 2A at the fast-spreading EPR, both high and low k 2A in neighboring upflow zones on the intermediate-spreading JdFR, and high k 2A at the slow-spreading MAR. However, as the model also shows, high 2A permeabilities are difficult to reconcile with high-temperature venting [START_REF] Driesner | The interplay of permeability and fluid properties as a first order control of heat transport, venting temperatures and venting salinities at mid-ocean ridge hydrothermal systems[END_REF], unless upwelling zones are relatively focused (radius 100 m). If that were the case, greater 2B permeability ($10 213 m 2 ) could enable greater venting temperatures at the intermediate and slow-spreading sites, and a greater Rayleigh number characterizing the convection system, which would be consistent with narrower upflow zones (Fontaine & Wilcock, 2007;[START_REF] Lowell | Hydrothermal processes at mid-ocean ridges: Results from scale analysis and single-pass models[END_REF]. In that alternate scenario, both layers 2A and 2B permeability would increase with decreasing spreading rate.

An increase in upflow zone permeability at slower spreading rates is broadly consistent with the tectonomagmatic conditions that prevail in such settings, especially with regards to layer 2A. Along slow-spreading MORs (e.g., at the LSHF), the upper crust and its extrusive portion (i.e., layer 2A) are relatively thick [START_REF] Cannat | How thick is the magmatic crust at slow-spreading oceanic ridges[END_REF][START_REF] Hooft | Crustal thickness and structure along three contrasting spreading segments of the Mid-Atlantic Ridge, 33.58-358N[END_REF][START_REF] Hussenoeder | Upper crustal seismic structure of the slow spreading Mid-Atlantic Ridge, 35 degrees N: Constraints on volcanic emplacement processes[END_REF][START_REF] Smith | Building the crust at the Mid-Atlantic Ridge[END_REF] and tectonized by faults and fissures that accommodate a large fraction (>0.5) of plate separation [START_REF] Behn | Magmatic and tectonic extension at mid-ocean ridges: 1. Controls on fault characteristics[END_REF][START_REF] Buck | Modes of faulting at mid-ocean ridges[END_REF]. A greater degree of faulting likely results in more widespread damage zones, which increases the matrix permeability. By contrast, at fast-spreading ridges (e.g., the EPR), magmatic processes accommodate a greater fraction of plate spreading [START_REF] Behn | Magmatic and tectonic extension at mid-ocean ridges: 1. Controls on fault characteristics[END_REF][START_REF] Buck | Modes of faulting at mid-ocean ridges[END_REF][START_REF] Cowie | Fault strain and seismic coupling on mid-ocean ridges[END_REF], resulting in a thinner, less tectonized upper crust and a layer 2A that is frequently repaved by lava flows (e.g., [START_REF] Fornari | Submarine lava flow emplacement at the East Pacific Rise 9850 0 N: Implications for uppermost ocean crust stratigraphy and hydrothermal fluid circulation[END_REF][START_REF] Fornari | The East Pacific rise between 98N and 108N: Twenty-five years of integrated, multidisciplinary oceanic spreading center studies[END_REF]. Robust magmatism at fast-spreading ridges thus hinders the formation of high-permeability upflow zones, potentially providing a simple explanation for the differences in phase lag and permeability structure observed between the LSHF and the EPR vent fields. This simple conceptual framework is illustrated in Figure 7.

It is of course difficult to infer general trends in permeability layering versus spreading rate with the limited data set at our disposal. This is particularly true given that our estimates for the JdFR suggest considerable variability in upflow zone permeability between neighboring upflow zones. Such short-wavelength (< 1 km) heterogeneity is unlikely to reflect changes in tectono-magmatic conditions, and may instead relate to hydrothermal alteration and mineralization processes active on a local scale. In any case, applying the methodology presented in this study to many other sites along the global MOR system will provide meaningful constraints on the natural variability of upflow zone permeability and its possible correlation with spreading rate.

Conclusions

We have developed an analytical multilayer poroelastic model and applied it to a comprehensive data set of in situ exit-fluid temperature time series data from sites where the phase lag between these data and ocean tidal loading can be robustly constrained, and the first-order crustal layering is well-constrained by seismic data. Applying our model to these data allows us to constrain the depth-dependent permeability structure of upflow zones across MOR hydrothermal fields. Our main conclusions are as follows: 1. Overall, we model the shallow crust beneath vents with phase lags ! 2008 as a low-permeability k 2A $ 10 213 210 212 m 2 À Á extrusive layer 2A, and the shallow crust beneath vents with phase lags 1808 as a relatively high-permeability k 2A $ ! 10 210 m 2 À Á extrusive layer 2A. 2. Our model and methodology places less stringent constraints on deeper crust (i.e., layer 2B) permeability, but does provide reliable upper bounds: $10 214 , $10 213 , and $10 212 m 2 at L-vent (EPR), MEF, and Lucky Strike, respectively, and plausible lower bound of $10 215 m 2 to satisfy heat flow constraints. 3. Our results are compatible with a scenario in which both layers 2A and 2B upflow zone permeability increase with decreasing spreading rate. This is consistent with the tectono-magmatic conditions that characterize different spreading rates (i.e., a greater degree of faulting at slow-spreading MORs increasing matrix permeability, and a greater degree of magmatism at fast-spreading MORs likely resulting in a less tectonized crust that is frequently repaved by lava flows.)

Figure 1 .

 1 Figure 1. (a) Conceptual illustration of the perturbation exerted on hydrothermal upflow by tidal loading. (b, c) Setup of the 1-D poroelastic model used to relate the phase lag between vent temperatures and tides to the permeability of the oceanic crust. (b) Single-layer model (Jupp & Schulz, 2004); (c) two-layer model (this study). AML refers to the axial melt lens providing the heat to drive hydrothermal circulation.

Figure 2 .

 2 Figure 2. Controls on the phase lag between vent temperatures and M2 tides. (a) In a single layer (2A) of thickness 155 m, a phase lag of 2078 requires a permeability of $10 213 m 2 (thick contours). (b) By contrast, in a two-layer model of total thickness 1,500 m, the same phase lag requires a specific combination of layers 2A and 2B permeabilities. The two-layer model asymptotes the single-layer model as k 2B decreases. (c) Typical phase lag map illustrating the four key regimes (R1-R4) described in section 2.3. (d) Sensitivity of 2 phase lag contours (2108 in red and 1608 in blue) to changes in the thickness of the layers. (e) Sensitivity of the same contours to changes in the assumed average pressure and temperature that set the thermodynamic properties of the fluid in both layers (see section 2.3).

  (1999); and[START_REF] Fontaine | Permeability changes due to mineral diagenesis in fractured crust: Implications for hydrothermal circulation at mid-ocean ridges[END_REF]. b Density of hydrothermal fluids calculated for pressure at 330 bar and hydrothermal fluid temperature at 3508C; from Holzbecher (1998);[START_REF] Rabinowicz | Thermal convection in a vertical permeable slot: Implications for hydrothermal circulation along mid-ocean ridges[END_REF][START_REF] Fontaine | Permeability changes due to mineral diagenesis in fractured crust: Implications for hydrothermal circulation at mid-ocean ridges[END_REF]. c Calculated for hydrothermal fluid temperature at 3508C; fromFontaine et al. (2001). d Sohn et al. (

  Phase lag angles are estimated for coherency c 2 ð Þ ! 0:85 at the M2 frequency for LSHF, MEF, and EPR and at the K1 frequency for MEF and EPR; and for coherency c 2 ð Þ ! 0:7 at the K1 frequency for LSHF. a Estimated for lower coherency c 2 ð Þ ! 0:7. b Arnulf et al. (2011). c Singh et al. (2006); Crawford et al. (2013). d Van Ark et al. (2007). e Sohn et al. (2004). f Detrick et al. (1987).

Figure 3 .

 3 Figure 3. Contours of measured phase lags at the Lucky Strike Hydrothermal field (LSHF, Western and Eastern vents), the East Pacific Rise (EPR, L-vent), and the Main Endeavour Field (MEF, S&M, and Grotto vents), plotted as a function of k 2A and k 2B assuming layers 2A and 2B geometries summarized in Table3. The left and right columns correspond to M2 and K1 tidal frequencies, respectively.

  5.1. Preamble: Limitations of the 1-D Poroelastic Modulation ModelOur multilayer poroelastic model makes several key assumptions that we briefly review. Firstly, the 1-D formulation does not account for horizontal gradients in fluid pressure, matrix displacement, or matrix

Figure 4 .

 4 Figure 4. (a) Combined constraints on layers 2A and 2B permeability for each site obtained by combining phase lag data at M2 and K1 frequencies, whenever available. (b) Synthesis of our results compared with permeability estimates from previous studies. k eff denotes an effective permeability that is not explicitly tied to either layers 2A nor 2B permeabilities.

Figure 5 .

 5 Figure 5. Effect of layer 2A thickness on the M2 phase lag, assuming a permeability of 10 214 m 2 in layer 2B and a range of permeabilities for layer 2A. Squares indicate measured values. See section 4.2 for details.

Figure 6 .

 6 Figure 6. Joint constraints on layers 2A and 2B permeability. Colored areas mark possible combinations of k 2A and k 2B compatible with the observed phase lag between tidal loading and vent temperatures. Colored lines indicate combinations of k 2A and k 2B compatible with a heat flow of 160, 450, or 600 MW in a cylindrical upflow zone of radius R 2B and R 2A , under the model assumptions presented in section 5.3. Two cases are shown-i.e., for R 2B 5 500 and 100 m. Lines are color-coded by the associated venting temperature T*. Colored stars indicate possible combinations of k 2A and k 2B that jointly satisfy constraints from tidal modulation and heat flow.

Figure 7 .

 7 Figure7. A conceptual model for changes in permeability layering and effective heat output in basalt-hosted systems across ridge spreading rates. Schematic cross sections represent typical upper oceanic crust stratification with corresponding permeabilities. In this interpretation, the variability in layer 2A upflow zone permeability across spreading rates primarily reflects the pervasiveness of tectonic processes in the shallow extrusive layer.

Table 2

 2 Baseline Fluid and Poroelastic Parameters

	Fluid and poroelastic parameters	Symbols	Values
	M 2 tidal angular frequency	x M2 (rad s 21 )	1.408 3 10 24
	K 1 tidal angular frequency	x K1 (rad s 21 )	7.2921 3 10 25
	Density of seawater Density of hydrothermal fluids Fluid viscosity	q 0 (kg m 23 ) q f (kg m 23 ) l (Pa s)	1,047.3 a 700.5 b 7:9310 25c
			Layer 2A	Layers 2B/2C
	Ratio of heat capacities	R	1.5	1.3
	P wave velocity	V p (m s 21 )	2,200 d	5,500 e
	S wave velocity	V s (m s 21 )	431	2,940
	Porosity	U	0.2 f	0.03 g
	Bulk density	q (kg m 23 )	2,500.1	2,882.5
	Grain bulk modulus	K g (GPa)	50 h	70 i
	Fluid bulk modulus	K f (GPa)	0.29 j	0.29 j
	Matrix (drained) bulk modulus	K m (GPa)	10.6	53.4
	Storage compressibility for 1-D loading	S 1 (Pa 21 )	7.46 3 10 210	1.05 3 10 210

a Density of seawater calculated for pressure at 330 bar and background seawater temperature at 38C; from Holzbecher (1998);

Rabinowicz et al.

10.1002/2017GC007152

Acknowledgments

This research was funded by the National Science Foundation grants (OCE-1536705 & OCE-1536943) and postdoctoral scholarships from Woods Hole Oceanographic Institution and Lamont-Doherty Earth Observatory.

Data used in this study can be found at Marine Geoscience Data System (MGDS-http://www.marine-geo.org/ index.php), PANGAEA (http://doi. pangaea.de/10.1594/PANGAEA. 820343), European Multidisciplinary Seafloor and water column Observatory (EMSO-http://www. emso-fr.org/), and Ocean Networks Canada (ONC-http://dmas.uvic.ca/ DataSearch). The authors are grateful to William Wilcock and Lars R€ upke for detailed and constructive reviews that significantly improved the manuscript and to Editor Joshua Feinberg. Finally, I dedicate this paper to B.T. and A.T.B.

Appendix A: Synthesis of Relevant Poroelastic Parameters

Here, we list the equations used to estimate the properties of the hydrothermal fluid and oceanic crust (elastic matrix):

1. Bulk density: q5/q f 1 12/ ð Þq s . 2. Shear modulus: G5qV 2 S .

3. Matrix drained bulk modulus: K5 G 4 3 2

4. Fluid bulk modulus:

, H is specific enthalpy and P is the steady state pore pressure (Pa).

5. Dimensionless poroelastic stress coefficient [START_REF] Detournay | Fundamentals of poroelasticity[END_REF]:

Þ . 6. Drained Poisson's ratio of the matrix frame: m5 3K22G 2 3K1G ð Þ . 7. Biot-Willis parameter-coefficient of effective stress-defined as a measure of the relative magnitude of the matrix and grain bulk moduli [START_REF] Nur | An exact effective stress law for elastic deformation of rock with fluids[END_REF]): a512 K Ks . 8. Skempton's coefficient: b5a a1/K 1 K f 2 1 Ks h i 21 .


Moussa Gomina
email: moussa.gomina@ensicaen.fr

Antoine Barbulée

Joël Bréard

Davy Duriatti

4 th ICNF -Smart Sustainable Solutions 1 MODELING DRAW PROCESS OF FLAX FIBER

In this study we experimentally determined the length distributions of the fiber bundles after each drawing. The detailed analysis of bundles lengths distributions before and after each drawing made it possible to understand the modes of rupture and to determine their kinetics of shortening and refining. Interestingly, the refining of the rovings is accompanied by a decrease in the rate of bark (set of very hydrophilic compounds), which will promote the wetting of the fibers by the polymer during the implementation of the composite and thus will improve significantly the hydric behavior of the material.

INTRODUCTION

In the composites industry, it is essential that the linear density of fibers (expressed in tex, i.e. the mass in grams of one kilometer of fiber) be perfectly controlled because the balance of the fabric depends on it. However, in the case of plant fibers, the linear density of the tape at the end of combing fluctuates according to the supply at the entrance. Draw is the industrial process of refining a roving to obtain a lower and more consistent count. A thinner roving gives better strength and flexibility to the fabric, allows easier infiltration of the resin and provides a more homogeneous structure to the composite, which improves its mechanical performance and life. During successive draws, the fiber bundles in the dry roving are stressed in tension, compression and in-plane shear, resulting in a degradation of the mechanical strength, thus inducing premature ruptures during weaving. It is therefore important to study the morphology of the fiber bundles as a function of draw. From a combed tape with count T1, an industrial draw device is used to obtain a finer roving with count T2 = p * T1 / q (p is the number of rovings with count T1 at the input and q is the ratio between the output speed and the input speed (Fig. 1)). At the n-th draw, it will be possible to produce a roving with count Tn much smaller than T1. To determine the length distributions of the fiber bundles in the roving according to the draw we used a testing setup purposely built in our laboratory for plant fibers. To describe the changes in the length and fineness of the fiber bundles, we have used the volume distributions instead of the frequency distributions used by some authors (Martin, 2014) because the mechanical contribution of a bundle depends on its length and cross section area.

RESULTS AND CONCLUSIONS

Figure 2a shows experimental values of average fiber bundles lengths and associated standard deviations, as a function of draw. The average length decreases with a loss of about 2 cm after each draw and the relative standard deviation of the dispersions is 44%. The analysis of the distributions of fiber bundles broken transversely (shortening) or longitudinally (refining) are shown in Figures 2b and2c, respectively. The combination of these results makes it possible to represent the changes in the aspect ratio of the fiber bundles (average length / mean diameter ratio) as a function of the drawing (Figure 2d). This study thus allows predicting the morphology of the fiber bundles during draw, which is an important input for the manufacture of reinforcements for composites. Future work will focus on the structural and mechanical characterization of the draw rovings.

Fig. 1

 1 Fig.1 Diagram of the drawing device

Fig. 2

 2 Fig.2 Evolution of average bundles lengths in rovings (a), details of transversely broken (b)or longitudinally broken (c) bundles and variation of aspect ratio as function of draw (d)

ACKNOWLEDGMENTS

The authors gratefully acknowledge the funding of PIAVE TRAMPLIN project by BPIFRANCE, France, under contract n° 130852.


Moussa Gomina
email: moussa.gomina@ensicaen.fr

Matthias Rouch

Annette Roy

Davy Duriatti

CORRELATION BETWEEN ACCELERATED AGING AND NATURAL AGING OF A BIO-SOURCED FLAX FIBER / UNSATURATED POLYESTER COMPOSITE

come

Correlation between accelerated aging and natural aging

of a bio-sourced flax fiber/unsaturated polyester composite

INTRODUCTION

Despite the high specific mechanical properties of plant fibers, their use for the manufacture of semi-structural composite parts is reluctant because of their thermal and hydric properties. Indeed water absorption of the fibers activates physical and chemical mechanisms identical to those usually reported for polymers (swelling, plasticization, leaching). It follows that composites incorporating plant fibers are much more sensitive to hygrothermal conditions than are materials incorporating synthetic fibers [START_REF] Azwa | A review on the degradability of polymeric composites based on natural fibres[END_REF]. This is manifested by anisotropic dimensional variations and degradation of mechanical properties [START_REF] Li | Hydrothermal ageing mechanisms of unidirectional flax fabric reinforced epoxy composites[END_REF].

The study focused on flax fiber rovings (Aretha variety), unsaturated isophthalic polyester (NORESTER 823) and the quasi-unidirectional biosourced composite combining these two constituents (V f = 32%; density 220 g/m 2 , 7% porosity). Accelerated ageing consisted of immersion in deionized water at 23°C or 70°C for a maximum of 70 days. For natural aging, the materials were exposed to natural conditions on the roof of the laboratory for 24 months.

RESULTS AND CONCLUSIONS

The volume of water absorbed by the composite immersed in water at 23°C as well as the induced swellings are represented as a function of time in Figure 1. The strong dimensional anisotropy observed is explained by the orthotropic hydro expansion behavior of the plant fibers and the elastic deformation recovery in wet environment related to the process of implementation of the composite [START_REF] Gamstedt | Moisture induced softening and swelling of natural cellulose fibres in composite applications[END_REF]. The analysis of the damage mechanisms of the composite and the results of tensile mechanical tests show a similarity between accelerated aging in water at 23°C and natural ageing even if the reduction of properties is less severe in the latter case (Fig. 2). On the other hand, the accelerated ageing in the air at 70°C induces very different damages from those recorded for a natural aging and the reduction of properties are much more important. These differences are explained by the hydrothermal behavior of flax fibers and more precisely relate to the nature of organic compounds on the surface of the fibers.

Fig. 1 Fig. 2

 12 Fig.1 Dimensional variations (green symbols) and volume of water absorbed (blue circle) by the composite during accelerated ageing in water at 23°C

ACKNOWLEDGMENTS

The authors gratefully acknowledge the funding of PIAVE TRAMPLIN project by BPIFRANCE, France, under contract n° 130852.