Title: Permeability structure of young ocean crust from poroelastically triggered earthquakes
Abstract: Geophysical Research LettersVolume 38, Issue 5 Solid EarthFree Access Permeability structure of young ocean crust from poroelastically triggered earthquakes Timothy J. Crone, Timothy J. Crone [email protected] Lamont-Doherty Earth Observatory, Columbia University, Palisades, New York, USASearch for more papers by this authorMaya Tolstoy, Maya Tolstoy Lamont-Doherty Earth Observatory, Columbia University, Palisades, New York, USASearch for more papers by this authorDanielle F. Stroup, Danielle F. Stroup Lamont-Doherty Earth Observatory, Columbia University, Palisades, New York, USASearch for more papers by this author Timothy J. Crone, Timothy J. Crone [email protected] Lamont-Doherty Earth Observatory, Columbia University, Palisades, New York, USASearch for more papers by this authorMaya Tolstoy, Maya Tolstoy Lamont-Doherty Earth Observatory, Columbia University, Palisades, New York, USASearch for more papers by this authorDanielle F. Stroup, Danielle F. Stroup Lamont-Doherty Earth Observatory, Columbia University, Palisades, New York, USASearch for more papers by this author First published: 10 March 2011 https://doi.org/10.1029/2011GL046820Citations: 20Find it @CityULibAboutSectionsPDF ToolsRequest permissionExport citationAdd to favoritesTrack citation ShareShare Give accessShare full text accessShare full-text accessPlease review our Terms and Conditions of Use and check box below to share full-text version of article.I have read and accept the Wiley Online Library Terms and Conditions of UseShareable LinkUse the link below to share a full-text version of this article with your friends and colleagues. Learn more.Copy URL Abstract [1] Permeability is a primary control on fluid flow within mid-ocean ridge hydrothermal systems and strongly influences the transfer of energy and mass between the ocean and the lithosphere. Little is known about how this parameter might vary in zero-age crust even though such variations may determine the locations and areal extents of upflow and downflow zones. Typically, estimates of permeability in seafloor environments are given as a single value (or range of values) for entire systems. Here we model crustal stresses inferred from poroelastically triggered earthquake patterns to estimate the two-dimensional permeability structure within a hydrothermal system on the East Pacific Rise at 9°50′N. We show that permeability in young ocean crust may vary by several orders of magnitude over horizontal scales of hundreds of meters with values ranging from 10−13.4 to 10−9.4 m2. Such values are consistent with other estimates of permeability in ocean crust. These variations may prescribe the geometry of hydrothermal convection and should be considered in future models of these systems. [2] The permeability of the lithosphere strongly controls fluid flow within mid-ocean ridge hydrothermal systems [Fisher, 1998, 2004]. This property helps determine the intensity of hydrothermal convection and the geometry of subseafloor flow, and it helps control the transfer of heat and chemicals between the crust and the ocean. For these reasons, permeability plays a key role in mediating the activity of chemosynthesizing microorganisms that support the diverse deep sea ecosystems often found in mid-ocean ridge settings [Kelley et al., 2002]. [3] Despite its importance, permeability is the most poorly constrained of the hydrological parameters, with little known about its spatial variation on the scales of vent fields or hydrothermal convection cells [Fisher, 2004]. Direct measurements of permeability have not been made because drilling into zero-age crust is extremely difficult. Models of these systems, ranging from analytical analyses [Wilcock and McNabb, 1996; Wilcock and Fisher, 2004; Lowell and Germanovich, 2004] to sophisticated three-dimensional simulations [Coumou et al., 2008, 2009], have provided substantial insights into crustal permeability, but most models assume uniform or nearly-uniform permeability distributions and have not explored the effects of abrupt permeability variations. [4] Because mid-ocean ridges are active spreading margins with extensive tectonism and magmatism, and hydrothermal flow helps drive the precipitation and dissolution of minerals, heterogeneous permeability distributions are expected. Characterizing the spatial variations of crustal permeability on the scales relevant to hydrothermal circulation will be a critical step toward understanding the plumbing of these systems and their influence on the chemical, geological, and biological processes within the deep ocean and Earth's youngest lithosphere. In this paper we use a two-dimensional poroelastic model of tidally induced stress perturbations, constrained by observations of tidal triggering patterns [Stroup et al., 2009], to map the permeability structure within a hydrothermal system on the East Pacific Rise (EPR). [5] During a seven month period in 2003–2004 a microearthquake survey was conducted near a hydrothermal vent field on the EPR at 9° 50′N [Tolstoy et al., 2008]. Using relative relocation techniques [Waldhauser and Ellsworth, 2000] 6,050 earthquakes were located along this section of ridge with a hypocentral accuracy of 50 m or less [Tolstoy et al., 2008] (Figure 1). The majority of these events were located within 500 m of the ridge's spreading axis and extended from very near the seafloor to a depth of about 1500 m, or to just above the axial magma chamber [Kent et al., 1993]. As previously observed at mid-ocean ridges [Wilcock, 2001; Tolstoy et al., 2002], earthquakes in this catalog tended to occur near times of peak predicted extensional stress [Stroup et al., 2007]. However, a detailed spatial analysis of these events [Stroup et al., 2009] revealed a more complex triggering pattern. Averaged over time, these earthquakes tended to occur in a "wave-like" pattern along axis, with events near the southern and central part of the study area tending to occur as much as four hours before those in the intervening sections of crust (Figures 1 and 2). The events occurring before maximal predicted extension are situated within the inferred downflow and upflow zones of a hydrothermal convection cell [Tolstoy et al., 2008]. Figure 1Open in figure viewerPowerPoint Map showing the epicenters of the 6,050 double difference relocated earthquakes used in this study and the locations of the ocean bottom seismometers. Each dot marking each event epicenter has been colored according to the average phase (relative to predicted maximal extension) of all events within a 100-m radius in plan view. The colors thus represent a vertically averaged phase. Figure 2Open in figure viewerPowerPoint (a) Colored contours of the mean phase of the closest 100 earthquakes to each grid cell using a three-dimensional search radius [Stroup et al., 2009], and (b) the data projected onto the along-axis dimension. Areas with low earthquake density in Figure 2a are colored white, gray dots indicate the locations of the 6,050 microearthquakes projected onto this plane, and red stars mark the locations of high temperature hydrothermal vents when these data were collected. (Vent locations are available through the Ridge 2000 Data Portal at http://www.marine-geo.org/portals/ridge2000.) A downflow zone (marked by arrows) is inferred near 9°49.3′N where earthquakes occur very near the seafloor [Tolstoy et al., 2008], and an upflow zone is presumed to exist beneath the hydrothermal vents. The wave-like pattern suggests that pore pressure perturbations travel along axis on each tidal cycle leading to enhanced triggering in different parts of the crust at different phases of the tide. [6] In continental systems, earthquake swarms can migrate through the lithosphere ostensibly triggered by aqueous fluid pressure perturbations [e.g., Ingebritsen and Manning, 2010]. As a potential explanation for the pattern observed at the EPR, we hypothesize that pore pressure perturbations generated by ocean tidal loading and Earth tides propagate laterally along axis from the inferred upflow and downflow zones into the intervening section of crust. As this pressure transient propagates through the crust, it raises the potential for fault rupture by lowering the normal stresses on fault planes. The effects of pore pressure on the mechanics of faulting is well studied [Beeler et al., 2000; Scholz, 2002; Shapiro et al., 2003], and although the process is complex, increases in pore pressure generally favor induced seismicity regardless of fracture orientation [Scholz, 2002]. Pore pressure gradients that drive the diffusion of along-axis pressure perturbations can be generated by variations in crustal permeability that allow pressure changes at the seafloor to enter the lithosphere at different rates and travel to different depths over a tidal cycle. [7] To test this hypothesis and estimate the required permeability structure, we developed a two-dimensional poroelastic model to predict the relative phases of extensional stress maxima within a fluid saturated crust subjected to Earth tides and ocean tidal loading [Crone and Wilcock, 2005]. The model domain (Figure 3) represents a 3200-m along-axis section of crust from the seafloor to a depth of 1500 m. The model medium is assigned elastic properties typical of zero-age seismic layer 2B lithosphere and a porosity of 5% [Crone and Wilcock, 2005; Stroup et al., 2009]. The model is isothermal (0°C) with the fluid assigned a viscosity of pure water [Holzbecher, 1998] and a density and bulk modulus of 3.2 wt% NaCl-H2O solution [Anderko and Pitzer, 1993; Pitzer et al., 1984]. The model is forced with a time-dependent fluid load on the top boundary to simulate ocean tidal loading. It is forced with a time-dependent solid load on one side boundary to simulate Earth tides, as these forces generate strains that are predominantly horizontal near the seafloor. We allow the ocean tidal loading stress to lag the Earth tide stress by 158 degrees, the value predicted by global tide models for this area [Matsumoto et al., 2000; Stroup et al., 2007]. The sides of the domain are symmetry boundaries and the bottom is closed to flow because a melt lens located just below 1500 m likely defines the base of the hydrothermal system along this section of the ridge [Kent et al., 1993]. Because strong horizontal pressure gradients are not generated near the side boundaries, the choice of closed side boundaries does not significantly affect our results. Figure 3Open in figure viewerPowerPoint Model geometry and boundary conditions used in this study. The poroelastic model domain is closed to fluid flow on the bottom and side boundaries, and open to flow across the top boundary, which represents the seafloor. The right-hand side boundary is forced with a sinusoidal solid stress with a period of 12 hours to simulate the effects of the solid Earth tide. The top boundary is forced with a sinusoidal fluid stress function with a similar period (lagging the Earth tide by 158 degrees) to simulate the effects of ocean tides. No boundary sustains any shear stresses. The elastic properties of the medium are assigned values that are typical of seismic layer 2B [Crone and Wilcock, 2005; Stroup et al., 2009]. For the finest resolution model runs, the poroelastic medium is divided into 64 vertical slices (gray lines), each of which can be assigned a different permeability value. We use a genetic algorithm minimization scheme to find the permeability distribution that minimizes the difference between the modeled extensional stress phases (Figure 4) and the inferred stresses from earthquake triggering (Figure 2). [8] We used forward modeling and a genetic algorithm technique [Goldberg, 1989] to find the permeability distribution that minimized the difference (i.e. the sum of the residuals where model and data overlap) between the modeled extensional stresses and the extensional stresses inferred from the triggering pattern (Figure 2). The minimization process proceeded in four stages: it started with a course grid of eight 400-m wide vertical permeability zones, and progressed to a finer grid of 64 50-m wide zones by doubling the number and halving the width of zones after each stage. The permeability within each zone was allowed to range between 10−16–10−8 m2. Resulting populations from each stage were interpolated onto a finer grid and used to seed the populations of subsequent stages. This approach allowed us to compute solutions with large numbers of individuals in the early stages to adequately sample the solution domain, and focus on solution refinement in the later stages as the computational cost of obtaining solutions increased. Because the residual sum function was relatively flat near the minimum, we repeated this procedure 32 times to obtain a rough estimate of the uncertainty of our permeability solution. [9] Figure 4 shows the results from this modeling procedure. Here we plot the average of the 32 best permeability distributions, contours of the resulting phases of maximal extensional stress, and horizontal cross-sections of the phases projected onto the along-axis dimension for comparison with the data from Figure 2. Error bars on the permeability distribution indicate the standard deviation of the permeability at each node from the 32 minimization runs. About two-thirds of the values in the permeability distribution are near 10−12 m2. However two sections have significantly higher permeability values. One section located about 400–700 m along axis has values that are near 10−11.5 m2, and another located about 2100–2500 m has values near 10−10 m2. Another section near the right-hand side boundary has values that are lower, with a minimum of about 10−13.5 m2. The baseline permeability of about 10−12 m2 is similar to (and bracketed by) values that have been predicted for hydrothermal systems using analytical models constrained by heat flow data [Wilcock and Fisher, 2004; Lowell and Germanovich, 2004], models of flow rate perturbations generated by earthquake swarms [Crone et al., 2010], and measurements and models in ridge flank systems [Becker and Fisher, 2008; Fisher et al., 2008; Hutnak et al., 2006]. Permeabilities in off-axis basement measured over very large spatial scales have been estimated to be as high as 1.7 × 10−10 m2 [Davis et al., 2000], a value similar to the highest permeability predicted by our model. Figure 4Open in figure viewerPowerPoint (a) Colored contours of the relative phases of modeled extensional stress maxima, (b) several horizontal cross sections of the data from the lower two-thirds of the model domain projected onto the along-axis dimension (blue dots) along with the earthquake triggering phases from Figure 2b (gray dots), and (c) the along-axis permeability distribution that generated the best fit between the model and the data. The standard deviation of the residuals (where model and data overlap) is 15 degrees. [10] In our model, this pattern in the timing of maximal extensional stress is generated by the Earth and ocean tidal forces through an interplay between the loading stresses and the fluid pressures which are strongly influenced by the permeability structure. At the beginning of each tidal cycle, as the ocean tide rises and compresses the poroelastic medium, fluid flows into the model from the top boundary, with more fluid flowing into the domain where the permeability is elevated. At nearly the same time, Earth tide stresses are acting to bring the system into extension. Where fluid pressures are greatest (i.e. within the high permeability zones), maximal extension is reached first. As fluids flow into lower permeability crust, maximal extension is achieved later in the cycle. With the permeability distribution shown in Figure 4c, the resulting stress pattern is similar to the inferred stresses from the triggering pattern in Figure 2. [11] This model can only generate large differential phase lags to match the earthquake data when the ocean tidal loading and Earth tide stress functions (Figure 3) have similar magnitudes (i.e. within about 25% of one another). The global tide model we use [Matsumoto et al., 2000] provides estimates of both phase and magnitude of these forcing functions, and because the EPR at 9°50′N is near an amphidromic point where ocean tides are small, the magnitude of the Earth tide stresses are predicted to be nearly twice that of the ocean loading stresses. However, although we have high confidence in the tide model's predictions of phase, we have lower confidence in its predictions of amplitude. The tide model's predictions are based on a radially-symmetric reference Earth model, which is a reasonable approximation for many locations on and within the Earth. It is probably not a good approximation for the lithosphere at a fast spreading mid-ocean ridge with a relatively shallow axial magma chamber and the presence of off-axis magma bodies [Kent et al., 1993; Canales et al., 2008]. Our results suggest that near 9°50′N on the EPR, the relative magnitude of Earth tide stresses may be smaller than predicted by a radially-symmetric reference Earth model, which is consistent with tilt measurements conducted in other ridge environments [Tolstoy et al., 1998]. A comprehensive geodetic survey of this section of the ridge using tilt meters and bottom pressure recorders may help resolve this issue. [12] Our results illustrate how horizontal variations in permeability could contribute to the pattern of earthquake triggering observed at 9°50′N on the EPR. The predicted permeability variations are large and occur over relatively small distances, but they are within the range of modeled values from other studies as well as values measured in off-axis environments. Such large permeability variations may be generated by tectonic stresses associated with the fourth-order ridge crest discontinuity near 9°49.3′N [Haymon, 1996], and the geochemically and seismically defined discontinuity at 9°50.3′N [Tolstoy et al., 2008; Von Damm and Lilley, 2004]. That these high permeability zones appear to coincide with areas of inferred downflow and upflow suggests that the permeability structure may strongly affect the geometry of hydrothermal flow along this section of the ridge. Future modeling studies of hydrothermal processes at mid-ocean ridges should consider the potential effects of large horizontal permeability variations in the lithosphere. These structural variations may ultimately control heat transport, chemical processes, and biological productivity in the subseafloor. Acknowledgments [13] We thank Andrew Fisher and one anonymous reviewer for helpful reviews of this manuscript. This research was supported by the National Science Foundation under grants OCE-0928181, OCE-0327283, OCE-0649538, and OCE-0732569. References Anderko, A., and K. S. Pitzer (1993), Equation of state representation of phase equilibria and volumetric properties of the system NaCl-H2O above 573 K, Geochim. Cosmochim. Acta, 57, 1657– 1680. Becker, K., and A. T. Fisher (2008), Borehole packer tests at multiple depths resolve distinct hydrologic intervals in 3.5-Ma upper oceanic crust on the eastern flank of Juan de Fuca Ridge, J. Geophys. Res., 113, B07105, doi:10.1029/2007JB005446. Beeler, N. M., R. W. Simpson, S. H. Hickman, and D. A. Lockner (2000), Poro fluid pressure, apparent friction, and coulomb failure, J. Geophys. Res., 105, 25,533– 25,542. Canales, J., et al. (2008), Discovery of off-axis melt lenses at the RIDGE-2000 East Pacfic Rise integrated studies site, Eos Trans. AGU, 89(53), Fall Meet. Suppl., Abstract B21A-0319. Coumou, D., T. Driesner, and C. A. Heinrich (2008), The structure and dynamics of mid-ocean ridge hydrothermal systems, Science, 321, 1825– 1828. Coumou, D., T. Driesner, S. Geiger, A. Paluszny, and C. A. Heinrich (2009), High-resolution three-dimensional simulations of mid-ocean ridge hydrothermal systems, J. Geophys. Res., 114, B07104, doi:10.1029/2008JB006121. Crone, T. J., and W. S. D. Wilcock (2005), Modeling the effects of tidal loading on mid-ocean ridge hydrothermal systems, Geochem. Geophys. Geosyst., 6, Q07001, doi:10.1029/2004GC000905. Crone, T. J., W. S. D. Wilcock, and R. E. McDuff (2010), Flow rate perturbations in a black smoker hydrothermal vent in response to a mid-ocean ridge earthquake swarm, Geochem. Geophys. Geosyst., 11, Q03012, doi:10.1029/2009GC002926. Davis, E. E., K. Wang, K. Becker, and R. E. Thomson (2000), Formation-scale hydraulic and mechanical properties of oceanic crust inferred from pore pressure response to periodic seafloor loading, J. Geophys. Res., 105, 13,423– 13,435. Fisher, A. T. (1998), Permeability within basaltic oceanic crust, Rev. Geophys., 36(2), 143– 182. Fisher, A. T. (2004), Rates of flow and patterns of fluid circulation, in Hydrogeology of the Oceanic Lithosphere, edited by E. E. Davis, and H. Elderfield, pp. 339– 377, Cambridge Univ. Press, Cambridge, U. K. Fisher, A. T., E. E. Davis, and K. Becker (2008), Borehole-to-borehole hydrologic response across 2.4 km in the upper oceanic crust: Implications for crustal-scale properties, J. Geophys. Res., 113, B07106, doi:10.1029/2007JB005447. Goldberg, D. E. (1989), Genetic Algorithms in Search Optimization and Machine Learning, Addison-Wesley, Reading, Mass. Haymon, R. (1996), The response of ridge-crest hydrothermal systems to segmented, episodic magma supply, in Tectonic, Magmatic, Hydrothermal and Biological Segmentation of Mid-Ocean Ridges, edited by C. J. MacLeod et al., Geol. Soc. Spec. Publ., 118, 157– 168, doi:10.1144/GSL.SP.1996.118.01.09. Holzbecher, E. O. (1998), Modeling Density-Driven Flow in Porous Media, Springer, Berlin. Hutnak, M., A. T. Fisher, L. Zühlsdorff, V. Spiess, P. Stauffer, and C. W. Gable (2006), Hydrothermal recharge and discharge guided by basement outcrops on 0.7–3.6 Ma seafloor east of the Juan de Fuca Ridge: Observations and numerical models, Geochem. Geophys. Geosyst., 7, Q07O02, doi:10.1029/2006GC001242. Ingebritsen, S. E., and C. E. Manning (2010), Permeability of the continental crust: Dynamic variations inferred from seismicity and metamorphism, Geofluids, 10, 193– 205, doi:10.1111/j.1468-8123.2010.00278.x. Kelley, D. S., J. A. Baross, and J. R. Delaney (2002), Volcanoes, fluids, and life at mid-ocean ridge spreading centers, Annu. Rev. Earth Planet. Sci., 30, 385– 491, doi:10.1146/annurev.earth.30.091201.141331. Kent, G. M., A. J. Harding, and J. A. Orcutt (1993), Distribution of magma beneath the East Pacific Rise between the Clipperton Transform and the 9°17′n deval from forward modeling of common depth point data, J. Geophys. Res., 98, 13,971– 13,995. Lowell, R. P., and L. N. Germanovich (2004), Hydrothermal processes at mid-ocean ridges: Results from scale analysis and single-pass models, in Mid-Ocean Ridges: Hydrothermal Interactions Between the Lithosphere and Oceans, Geophys. Monogr. Ser., vol. 148, edited by C. R. German et al., pp. 219– 244, AGU, Washington, D. C. Matsumoto, K., T. Takanezawa, and M. Ooe (2000), Ocean tide models developed by assimilating TOPEX/POSEIDON altimeter data into hydrodynamical model: A global model and a regional model around Japan, J. Oceanogr., 56, 567– 581. Pitzer, K. S., J. C. Peiper, and R. H. Busey (1984), Thermodynamic properties of aqueous sodium chloride solutions, J. Phys. Chem. Ref. Data, 13, 1– 106. Scholz, C. H. (2002), The Mechanics of Earthquakes and Faulting, 2nd ed., Cambridge Univ. Press, Cambridge, U. K. Shapiro, S. A., R. Patzig, E. Rothert, and J. Rindschwentner (2003), Triggering of seismicity by pore-pressure perturbations: Permeability-related signatures of the phenomenon, Pure Appl. Geophys., 160, 1051– 1066. Stroup, D. F., D. R. Bohnenstiehl, M. Tolstoy, F. Waldhauser, and R. T. Weekly (2007), Pulse of the seafloor: Tidal triggering of microearthquakes at 9°50′N East Pacific Rise, Geophys. Res. Lett., 34, L15301, doi:10.1029/2007GL030088. Stroup, D. F., M. Tolstoy, T. J. Crone, A. Malinverno, D. R. Bohnenstiehl, and F. Waldhauser (2009), Systematic along-axis tidal triggering of microearthquakes observed at 9°50′N East Pacific Rise, Geophys. Res. Lett., 36, L18302, doi:10.1029/2009GL039493. Tolstoy, M., S. Constable, J. Orcutt, H. Staudigel, F. K. Wyatt, and G. Anderson (1998), Short and long baseline tiltmeter measurements on axial seamount, Juan de Fuca Ridge, Phys. Earth Planet. Inter., 108, 129– 141. Tolstoy, M., F. L. Vernon, J. A. Orcutt, and F. K. Wyatt (2002), Breathing of the seafloor: Tidal correlations of seismicity at Axial volcano, Geology, 30(6), 503– 506. Tolstoy, M., F. Waldhauser, D. R. Bohnenstiehl, R. T. Weekly, and W. Y. Kim (2008), Siesmic identification of along-axis hydrothermal flow on the East Pacific Rise, Nature, 451, 181– 184, doi:10.1038/nature06424. Von Damm, K. L., and M. D. Lilley (2004), Diffuse flow hydrothermal fluids from 9°50′N East Pacific Rise: Origin, evolution and biogeochemical controls, in The Subseafloor Biosphere at Mid-Ocean Ridges, Geophys. Monogr. Ser., vol. 144, edited by W. S. D. Wilcock et al., pp. 245– 268, AGU, Washington, D. C. Waldhauser, F., and W. L. Ellsworth (2000), A double-difference earthquake location algorithm: Method and application to the northern Hayward fault, California, Bull. Seismol. Soc. Am., 6, 1353– 1368. Wilcock, W. S. D. (2001), Tidal triggering of microearthquakes on the Juan de Fuca Ridge, Geophys. Res. Lett., 28(20), 3999– 4002. Wilcock, W. S. D., and A. T. Fisher (2004), Geophysical constraints on the sub-seafloor environment near mid-ocean ridges, in The Subseafloor Biosphere at Mid-Ocean Ridges, Geophys. Monogr. Ser., vol. 144, edited by W. S. D. Wilcock et al., pp. 51– 74, AGU, Washington, D. C. Wilcock, W. S. D., and A. McNabb (1996), Estimates of crustal permeability on the Endeavour segment of the Juan de Fuca mid-ocean ridge, Earth Planet. Sci. Lett., 138(1–4), 83– 91. Citing Literature Volume38, Issue516 March 2011 FiguresReferencesRelatedInformation