# Effects of fluid influx, fluid viscosity, and fluid density on fluid migration in the mantle wedge and their implications for hydrous melting

Nestor G. Cerpa, Ikuko Wada, Cian Wilson

► To cite this version:

Nestor G. Cerpa, Ikuko Wada, Cian Wilson. Effects of fluid influx, fluid viscosity, and fluid density on fluid migration in the mantle wedge and their implications for hydrous melting. Geosphere, 2018, 15 (1), pp.1-23. 10.1130/GES01660.1. hal-02130796

HAL Id: hal-02130796

https://hal.umontpellier.fr/hal-02130796

Submitted on 2 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

THE GEOLOGICAL SOCIETY
OF AMERICA®

This paper is published under the terms of the
CC-BY-NC license.

© 2018 The Authors

# Effects of fluid influx, fluid viscosity, and fluid density on fluid migration in the mantle wedge and their implications for hydrous melting

Nestor G. Cerpa<fnref n="1" />, Ikuko Wada<fnref n="2" />, and Cian R. Wilson<fnref n="3" />

## ABSTRACT

The migration pathways of hydrous fluids in the mantle wedge are influenced by the compaction of the porous mantle matrix, which depends on the matrix permeability, fluid viscosity, and fluid density. Experimental studies show that when fluids are interconnected, the permeability depends on mineral grain size and porosity, the latter of which depends on the amount of fluids introduced into the system (fluid influx). Here, we investigate the role of fluid influx, fluid viscosity, and fluid density in controlling fluid migration in the mantle wedge, using a 2-D numerical model accounting for the effects of grain-size variation and matrix compaction. Our models predict that fluid influx and fluid viscosity are key controls on fluid pathways, while fluid density plays a secondary role. Temperature dependence of fluid viscosity promotes downdip drag of fluids at the base of the forearc mantle toward the subarc region. High fluid influx at postarc depths promotes updip flow near the base of the mantle wedge, guiding the fluids arcward. The model that is applied to northern Cascadia predicts upward fluid migration focused beneath the arc but cannot explain high electrical conductivity observed slightly west of the upward fluid migration. We estimate the amount of hydrous melt that can be produced in the mantle wedge using calculated fluid distributions. Up to a few percent partial melting is predicted in a relatively small region in the core part of the subarc mantle wedge in most subduction settings, including northern Cascadia, and beneath the backarc in old-slab subduction zones.

## 1. INTRODUCTION

Generation of arc magmas in subduction zones occurs primarily because the addition of water from the dehydrating slab to the overlying mantle lowers the solidus of the mantle (Tatsumi et al., 1986; Grove et al., 2006). Experimental studies and numerical thermo-petrological models predict that the depths of slab dehydration depend on the thermal state of the slab and can vary widely

among different subduction zones (Schmidt and Poli, 1998; Hacker et al., 2003; van Keken et al., 2011), potentially promoting magma generation over a wide region. However, present-day arcs are located in a relatively narrow region that lies ~100 km above the top of the slab regardless of the thermal structure of the subduction zone (England et al., 2004; Syracuse and Abers, 2006; England and Katz, 2010), and it is unclear what controls the location of the arc. What are the pathways of hydrous fluids in the mantle wedge? How do fluids become focused beneath the arc?

Theoretical and experimental studies predict that fluids flow through interconnected pores at grain edges at the pressure and temperature conditions relevant to the ductile part of the mantle wedge (Bulau et al., 1979; von Bargen and Waff, 1986; Mibe et al., 1999), and the governing equations that describe the percolation of fluids through a viscously deformable porous matrix (Sleep 1974; McKenzie 1984; Scott and Stevenson, 1984; Fowler, 1985) have been widely used to model the migration of slab-derived fluids and mantle melts in subduction zones. In this physical model, fluid migration is controlled by solid flow and fluid flow through the solid matrix, the latter of which is driven by the buoyancy of the fluid and the pressure gradients that are induced by compaction and dilation of the viscous matrix over a characteristic length, called the compaction length (Spiegelman, 1993). Fluid flow also depends on fluid mobility, defined here as the ratio of the matrix permeability to fluid viscosity. The former is related to the geometry of the pores network and is proportional to the square of mineral grain size and to a power (generally between 2 and 3) of the fluid fraction in mantle rocks (von Bargen and Waff, 1986; Wark et al., 2003; Simpson et al., 2010; Miller et al., 2014). Based on this model of fluid transport, several hypotheses regarding how fluids are focused beneath the arc have been introduced.

In an earlier study, Spiegelman and McKenzie (1987) proposed the role of dynamic pressure as the main driving force to redirect fluids toward the subarc region. In their model, large dynamic pressure gradients induced by mantle corner flow drove fluids toward the tip of the corner flow. However, the magnitude of pressure gradients was overestimated by the use of a uniform

<footnote n="1">¹Géosciences Montpellier, UMR 5243, Université de Montpellier, CNRS, 34095 Montpellier, France</footnote>

<footnote n="2">²Department of Earth Sciences, University of Minnesota, Minneapolis, Minnesota 55455-0149, USA</footnote>

<footnote n="3">³Department of Terrestrial Magnetism, Carnegie Institute for Sciences, Washington, District of Columbia 20015-1305, USA</footnote>
mantle viscosity of 10²¹ Pa s, which largely exceeds current estimates of tem-
perature-dependent and strain-rate-dependent viscosity thought to better rep-
resent mantle rheology (e.g., Kelemen et al., 2003; van Keken et al., 2008).

Wilson et al. (2014) studied the role of matrix compaction in controlling fluid pathways in subductions zones, using a temperature-dependent solid rheology in their fluid migration models. Based on their results, they found two possible mechanisms that focus fluids that are liberated from the subducting slab at postarc depths toward the subarc region, emphasizing the role of compaction in controlling fluid pathways. The first mechanism occurs in the uppermost part of the subducting slab where fluids are produced by dehydration and form a layer of high permeability parallel to the dip of the slab. Resulting compaction pressure gradients further drive nearby fluids into this region, resulting in updip channelized fluid flow through the uppermost part of the slab. Thereby, some of the fluids produced at postarc depth migrate updip and are subsequently introduced into the mantle wedge at subarc depths. The second mechanism occurs in the mantle wedge. The fluids that do not migrate updip are directly released into the mantle wedge. During their ascent through the mantle wedge, their pathways are deflected toward the trench. The deflection is promoted by an upward increase in solid shear viscosity in the shallowest part of the mantle wedge, leading to the formation of decompacting channels that are subparallel to the isotherms. The results showed that the fluid pathways depend strongly on the compaction length, over which fluids feel the effects of matrix compaction, and that the two above mechanisms are only effective when compaction lengths are sufficiently high. Because compaction depends on fluid mobility, the authors suggested that spatial variation of grain size is another possible factor that contributes to focusing of fluids beneath the arc because it affects the mantle permeability.

Using simplified models of porous flow, in which the compaction of the
matrix was neglected, several studies address the migration of hydrous fluids
in the mantle wedge (Iwamori, 1998, 2000; Arcay et al., 2005; Cagnioncle et al.,
2007; Hebert and Montesi, 2013; Angiboust et al., 2012; Wada and Behn 2015).
Cagnioncle et al. (2007) illustrated the influence of grain size on fluid flow in
the mantle wedge. Assuming a uniform grain size, the authors showed that
relatively small grain size could lead to trapping of fluids at the base of the
mantle wedge and their downdip transport by mantle outflow. Later, Wada et
al. (2011) estimated the grain-size distribution in the mantle wedge using the
"paleowatmetter" model of Austin and Evans (2007). Grain size was shown
to increase from a few µm near the tip of the corner flow to a few cm in the
hottest part of the mantle wedge. The calculated grain-size distribution was
adapted by Wada and Behn (2015) in a simplified porous flow model to test
its influence on fluid migration. The results showed that small grain size along
the base of the mantle wedge at shallow depth causes entrapment of fluids in
the downgoing mantle, providing a mechanism for downdip transport of fluids
from forearc depths to the subarc region. In their model, however, the effect of
matrix compaction was not incorporated.

In our previous study, we developed numerical models of fluid migration
in the mantle wedge accounting for both grain-size distribution and mantle

compaction (Cerpa et al., 2017). The modeling results showed that spatial vari-
ations of grain size and solid shear viscosity both serve to focus fluids beneath
the arc, and that their contributions as fluid-focusing mechanisms vary spa-
tially in the mantle wedge. Grain-size variation plays a key role at the base of
the mantle wedge near the tip of the corner flow and focuses fluids released
beneath the forearc toward the subarc region for all modeled subduction
zones as long as fluids are introduced beneath the forearc. The pathways of
fluid released at postarc depths in cold-slab subduction zones are deflected
toward the arc at shallow depths due to gradients of compaction pressure and
advection by mantle inflow. In the study, the pattern of fluid release from the
slab was calculated based on the thermal state of the slab, but only a frac-
tion of the total fluid released from the slab was introduced into the mantle
wedge such that upward fluid migration occurs in the model. Fluid influx is
an important parameter that determines the fluid fraction (i.e., porosity) at the
base of the mantle wedge and thus partly controls the permeability in this
region and subsequent evolution of fluid migration pathways in the overlying
mantle. Fluid viscosity and fluid density vary by several orders of magnitude
with temperature and composition (mainly the water content) at relatively high
pressure (e.g., Audétat and Keppler, 2004) and can also affect fluid migration,
yet their effects were not quantified in the previous study.

The present work focuses on quantifying the effects of fluid influx and
fluid properties on the pattern of fluid migration through the application of
the fluid migration model developed by Cerpa et al. (2017). In the following,
we first briefly summarize the results of theoretical and laboratory studies,
from which our model parameters are derived, and the modeling approach.
We then present a series of modeling results that illustrate the effects of fluid
influx, fluid viscosity, and fluid density. In the last three simulations presented
in this study, the subduction parameters are chosen to represent the North-
ern Cascades subduction zone, where the distribution of fluids in the mantle
wedge has been inferred from magnetotelluric (MT) observations. Finally, we
discuss our results with an emphasis on the implications for hydrous melting
in the mantle wedge.

2. PERMEABILITY AND FLUID PROPERTIES IN THE MANTLE WEDGE

## 2.1. Permeability of Mantle Rocks

Fluid migration by porous flow occurs if the fluids form an interconnected
network of channels or pores within the hosting solid rock. Indeed, an inter-
granular fluid phase can be distributed along grain edges forming the chan-
nels or at the grain corners forming nodes of fluid, depending on which con-
figuration minimizes the free energy at grain boundaries (Bulau et al., 1979).
The dihedral angle θ, that is, the angle between two fluid-solid interfaces at
a corner of a fluid-filled pore, controls the topology and interconnectivity of
the fluid phase in texturally equilibrated rocks. When θ < 60°, channels are
interconnected, and when θ > 60°, the channels are closed off, and fluid is trapped in isolated pockets at grain corners, preventing fluid migration (Bulau et al., 1979; Watson, 1982). However, von Bargen and Waff (1986) computed dihedral angles for single-phase isotropic tetrakaidekahedral crystals of constant size and found that the transition from open to closed channels ("pinch-off" criterion) only applies at relatively low fluid fractions (<<0.01), and at relatively high porosities (>>0.01), grain boundaries may become wetted (Waff and Faul, 1992; Miller et al., 2014), particularly when the bulk undergoes high shear stress conditions (Hier-Majumder and Kohlstedt, 2006). At the base of the mantle wedge, a relatively high fluid fraction and a high shear stress condition are expected, likely promoting interconnectivity.

Laboratory measurements of the dihedral angle in synthetic rock samples hosting CO₂-H₂O fluids at 1000 °C and 1 GPa were carried out by Watson and Brenan (1987). Their results showed high dihedral angles (above 70°), suggesting that hydrous fluids could not form an interconnected network in the mantle at relatively low porosities. However, in a subsequent study, Mibe et al. (1998) measured the dihedral angle in aqueous fluid-forsterite systems at 1000 °C and 3-5 GPa and reported a mean dihedral angle of 40°. They emphasized that θ may drop below 60° at 2 GPa as a result of changes in pressure-dependent olivine solubility that affects the composition of hydrous fluids. Therefore, at the conditions of the flowing part of the mantle wedge, fluids likely occupy interconnected pores and migrate by porous flow (Mibe et al., 1999).

At textural equilibrium, the relationship between permeability, fluid fraction, and grain-size follows a power law of the form of

$$K(b, \phi) = \frac{\phi^n b^2}{C}, \qquad (1)$$

where b is the grain size diameter or channel spacing, φ is fluid volume fraction or porosity, n is a porosity exponent, and C is a geometrical factor (McKenzie, 1984; von Bargen and Waff, 1986; Turcotte and Schubert, 2014). In this study, we use C = 270 (Wark et al., 2003).

In the creeping regime, grain size can evolve through grain growth and dynamic recrystallization (Karato, 1984). Grain growth occurs by grain boundary migration, minimizing the interfacial free energy between grains (Atkinson, 1988; Evans et al., 2001). Under applied stresses, dynamic recrystallization occurs by reorganizations of dislocations to form new subgrain boundaries and/or by migration of existing grain boundaries to form new grains.

Austin and Evans (2007) derived a model for the evolution of grain size called the "paleowattmeter" model, in which the rate of mean grain-size reduction through dynamic recrystallization is related to the mechanical work done to the grains. In the model, dynamic recrystallization and grain growth are assumed independent of each other, and when their rates are balanced, the mean grain size reaches its steady state (Behn et al., 2009). Wada et al. (2011) coupled the "paleowattmeter" model to thermomechanical models of subduction zones to estimate the steady-state grain-size distribution in the mantle wedge. They found that the grain size increases from a few µm near the tip of the corner flow to a few cm in the hottest part of the flowing mantle wedge

and that this downdip increase in grain size by approximately two orders of magnitude occurs regardless of the subduction parameters. The model likely overestimates the rate of downdip increase in grain size with depth at the base of the mantle wedge where the effect of advection, which is neglected in the grain-size calculation, becomes relatively important (Wada et al., 2011), and therefore the fine grain size may persist to even greater depths than the model predicts. Furthermore, the "paleowattmeter" model is based on experiments performed with single-phase (olivine) systems and neglects the effects of a secondary phase on grain growth. The presence of a secondary phase, such as pyroxene, likely limits the growth of olivine grains to a few cm (Evans et al., 2001, and references therein). Nonetheless, Wada et al. (2011) tested the effects of imposing an upper limit of their calculated grain size to a few cm and found little effect of this value on the grain-size distribution in the regions where the grain size is smaller than a few cm.

The effects of grain-size distribution on fluid migration have been studied by Wada and Behn (2015) and Cerpa et al. (2017), who showed that fluids become trapped in small grain size at shallow depths and become dragged downdip by solid flow until grain size becomes large enough to allow upward fluids migration, potentially redirecting shallow fluids toward the subarc region. While these studies focused on the effect of grain size on permeability, the role of fluid influx at the base of the mantle wedge was not investigated. The imposed fluid influx controls the fluid fraction and thus permeability at the base of the mantle wedge. Thus, in this study, we focus on the role of fluid fraction at the base of the mantle wedge in controlling the evolution of fluid migration.

The parameter n from Equation (1) that relates permeability to fluid fraction has been estimated experimentally and numerically for partially molten rocks. There has been little experimental work done on the quantification of n for a hydrous-fluid-bearing rock. However, the theoretical models used for the prediction of n are based primarily on geometrical considerations, and therefore the results have been applied to models of water-rich fluids migration. A similar reasoning applies to the geometrical factor C.

Estimations of n range from ~2 for idealized isotropic systems of uniform grain size (Frank, 1968; Maaløe and Scheie, 1982; von Bargen and Waff, 1986; Simpson et al., 2010) to ~3 for more complex systems (anisotropic, non-uniform grain sizes) (Wark et al., 2003; Connolly et al., 2009). More recently, Miller et al., (2014) conducted laboratory estimations of permeability in synthetic partially molten rocks using X-ray microtomography and have inferred an exponent of n = 2.6 from their data set. These authors observed two topological regimes that could be represented by two end members: n = 2, when melt resides principally along channels at relatively low porosities and n = 3, when melt tends to form films at grain boundaries at relatively high porosities (see also Rudge, 2018). As discussed by the authors, their inferred value of n may represent the coexistence of both regimes in different sub-volumes of their samples. In our study, we are primarily interested in melt migration from the bottom of the mantle to the region of flux melting where a relatively high fluid fraction is expected. Furthermore, the mantle wedge corner flow is a high-deformation environment, which is likely to promote wetting of grain
boundaries (Hier-Majumder and Kohlstedt, 2006). Hence, in our fluid-flow models, we assume $n = 3$.

The fluid influx at the base of the mantle wedge controls the overall porosity in the system. To assess the effects of fluid fraction on permeability and on fluid-flow pathways, we systematically vary the fraction of the thermopetrologically predicted slab-derived fluids that are allowed to flux into the mantle wedge in our model.

## 2.2. Fluid Viscosity and Density

A series of dehydration reactions occurs in the downgoing slab, producing a hydrous fluid that can migrate into the mantle wedge. The hydrous fluid likely dissolves silicate minerals as it percolates the subducted lithosphere. The solubility of silicates in hydrous fluids increases with pressure and temperature (Manning, 2004, and references therein), and thus the distribution of dissolved silicates in the fluids that migrate into the mantle wedge depends on the thermal state of the slab and history of the fluid. Fluids that migrate into the hot part of the mantle wedge trigger melting of the mantle (e.g., Grove et al., 2006). The resulting silicate-rich melts contain some amount of water as inferred from the H₂O contents found in arc lavas (Wallace, 2005, and references therein), but the exact H₂O contents in arc magmas are still debated (Plank et al., 2013). Hence, from the slab to the overlying crust, the composition of the migrating fluids can vary significantly, resulting in spatial variations in the physical properties of the fluids. Of those, viscosity and density affect the mobility and the buoyancy of fluids, respectively, and are critical parameters to fluid migration.

Fluid viscosity in the mantle wedge is expected to vary by several orders of magnitude, spanning from values similar to pure water 10⁻⁴ Pa s near the base of the mantle wedge to those of highly silicic melts 10¹¹ Pa s as they reach the overlying crust (Audétat and Keppler, 2004; Giordano et al., 2008; Hack and Thompson, 2011). Audétat and Keppler (2004) measured the viscosities of H₂O-rich fluids (>20 wt%) with varying water and dissolved silica contents at high temperatures (600–950 °C) and high pressures (1–2 GPa). They found that fluid viscosity is insensitive to pressure variation but is strongly dependent on temperature and composition and proposed an empirical relationship for the viscosity (μ) of fluids as a function of dissolved silica concentration ($c_{sil}$) (in wt%) and temperature (T) (in K) in subduction zones:

$$ \log(\mu) = -4.2 + \left(0.008 c_{sil} \left(\frac{10000}{T} - 2.0\right)\right). \quad (2) $$

The authors suggest that their model reproduces the viscosities of fluids that contain up 80 wt% of dissolved silicate in experiments. Their model predicts that fluid viscosity varies from 10⁻⁴ Pa s (e.g., viscosity of pure water) for $c_{sil} = 0$ wt% to 10 Pa s for $c_{sil} = 80$ wt% at T = 700°C (Fig. 1). In this study, we test the role of fluid viscosity in controlling fluid migration pathways by applying a range of constant fluid viscosities in the range of 0.001–10 Pa s in our fluid

migration models, which corresponds approximately to a temperature range of 700–1300 °C and dissolved silica concentrations of >20 wt%.

In contrast to fluid viscosity, fluid density is not expected to vary by more than a factor of two in the mantle wedge. We will illustrate the effects of density through three models with different density values.

## 3. MODELING APPROACH

We use the modeling approach of Cerpa et al. (2017), which incorporates the effects of grain-size distribution and mantle compaction on fluid migration. Given the relatively low porosities that are calculated in our models, it is reasonable to assume that the effects of the fluid phase on the solid-state mantle flow are negligibly small (Wilson et al., 2014). Based on this assumption, we neglect the effect of compaction on the divergence of solid velocity field (e.g., Dymkova and Gerya, 2013), the effects of porosity on the buoyancy of the mantle, and the weakening of the solid due to the presence of a free-fluid phase (Hirth and Kohlstedt, 2003). The advection of heat by the fluid phase (Rees-Jones et al., 2018) is also neglected in the temperature calculations.

<figure><img src="image_4.png" /><figcaption>Figure 1. Fluid viscosity as a function of silica content in the fluid for various temperatures (black solid lines) calculated using the Arrhenian model of Audétat and Keppler (2004). Gray dashed lines indicate the fluid viscosity values used in models of fluid migration (labels).</figcaption></figure>
The governing equations are the same as in Cerpa et al. (2017) and are briefly summarized in this section. Parameters with a prime are non-dimensional, while other parameters are dimensional. The subscript 0 indicates that the parameter is a reference parameter used for the non-dimensionalization of the equations.

## 3.1. Governing Equations

### 3.1.1. Solid Flow

We assume that the solid-state mantle flow has reached a steady state. The governing equations for the solid flow are:

$$ \nabla \cdot (2\eta\dot{\epsilon}') - \nabla p' = 0 \qquad (3) $$

$$ \nabla \cdot v_s' = 0 \qquad (4) $$

$$ v_s' \cdot \nabla T' - \nabla \cdot (\kappa' \nabla T') = Q', \qquad (5) $$

where $v_s$ is mantle (solid) velocity vector, $p$ is dynamic pressure, $T$ is temperature, $\kappa$ is thermal diffusivity, $Q$ is radiogenic heat source, and $\dot{\epsilon}$ is strain rate tensor defined as

$$ \dot{\epsilon}' = \frac{1}{2}(\nabla v_s' + \nabla v_s'') \qquad (6) $$

The mantle (solid) shear viscosity $\eta$ accounts for both diffusion creep and dislocation creep (see Appendix A). The above governing equations are exactly the same as those in Cerpa et al. (2017) except that we include radiogenic heat in the heat transfer equation.

As in Cerpa et al. (2017), the thermomechanical model is coupled to the paleowattmeter model (Austin and Evans, 2007) to calculate the grain-size distribution in the mantle wedge. Our model assumes that the timescale of grain-size evolution is small compared to the timescale of the fluid migration, and thus we neglect advection in the calculation of the steady-state grain size (Behn et al., 2009; Wada et al., 2011). The parameters used in the calculations are given in Appendix A.

### 3.1.2. Fluid Flow

In the fluid-flow calculations, we neglect mass exchanges between phases, and we do not account for the effects of the dynamic pressure due to mantle shear. The non-dimensional governing equations for the fluid flow are thus

$$ \frac{\partial \phi'}{\partial t'} + v_s' \cdot \nabla \phi' - \frac{h_0^2 v_0}{\delta_0^2 v_{s_0}} \frac{\mathcal{P}'}{\zeta'} = 0 \qquad (7) $$

$$ \frac{h_0^2 \mathcal{P}'}{\delta_0^2 \zeta'} - \nabla \cdot \left( \frac{\tilde{K}}{\mu'} \nabla \mathcal{P}' \right) - \nabla \cdot \left( \frac{\tilde{K}}{\mu'} \hat{g} \right) = 0, \qquad (8) $$

where $\phi$ is porosity or fluid fraction, $\mathcal{P}$ is compaction pressure, $\zeta$ and $K$ are bulk viscosity and permeability of the solid, respectively, $\mu$ is fluid viscosity, $\hat{g}$ is the unit vector in the direction of gravity, $v_{s_0}$ and $v_{s_0}$ are fluid and solid reference velocities, respectively. $h_0$ is a reference length scale and $\delta_0$ is a reference compaction length. The relationships among the reference parameters that are important in our study are given in Appendix B along with their values (see also Cerpa et al., 2017). In particular, in our models the reference compaction length is defined as

$$ \delta_0^2 = \eta_0 \left( \frac{\Psi_0^2 K_0}{(\Delta \rho g)^2 \mu_0} \right)^{\frac{1}{3}}, \qquad (9) $$

where $\Delta \rho g$ is fluid buoyancy, $\eta_0$, $K_0$ and $\mu_0$ are reference solid viscosity, permeability, and fluid viscosity, respectively. $\Psi_0$ is a reference flux used for the non-dimensionalization of the fluid influx boundary condition as described below. As described by Wilson et al. (2014), compaction lengths govern the mode of fluid flow. Low compaction lengths favor a fluid transport by solid advection, while high compaction lengths promote a buoyancy-driven flow and enhance the effects of compaction pressure. In our models, the compaction length varies spatially and temporally, and the non-dimensional compaction length is defined as

$$ \delta^{\prime 2} = \sqrt{\frac{K}{\mu'} \zeta'}, \qquad (10) $$

where permeability ($K'$) and bulk solid viscosity ($\zeta'$) are

$$ K' = b'^2 \phi'^3 \quad \text{and} \qquad (11) $$

$$ \zeta' = \eta' \phi'^{-2}. \qquad (12) $$

We regularize the above two parameters to avoid numerical singularities (Wilson et al., 2014; Cerpa et al., 2017) and express regularized values as $K'$ and $\zeta'$.

The fluid velocity, $v_{s,}$, is the vectorial sum of the solid velocity and the fluid velocity relative to the mantle matrix based on Darcy's law:

$$ v_{s,}' = \frac{v_{s_0}}{v_{s_0}} v_{s,}' - \frac{1}{\phi_{s,}'} \frac{K'}{\mu'} (\nabla \mathcal{P}' + \hat{g}) = 0. \qquad (13) $$

## 3.2. Model Setup

The model setup for solid-state mantle flow is similar to that used in Cerpa et al. (2017), except that the Cascadia model accounts for internal radiogenic heat in the crust (Wada and Wang, 2009).
<figure><img src="image_5.png" /></figure>

In the solid-state mantle flow model, the subduction velocity ($v_{sub}$) is prescribed at the top and the bottom of the subducting slab. We assume the slab and the overlying mantle are decoupled at shallow depths and become fully coupled at a 75 km depth (hereafter referred to as the maximum decoupling depth [MDD]) to be consistent with previous studies (e.g., Wada and Wang, 2009). The rigid overriding crust is fixed. The geotherm applied on the trench-side vertical boundary of the model is calculated using the GDH1 plate cooling model, which assumes a 95-km-thick slab (Stein and Stein, 1992). The geotherm imposed on the shallow part of the backarc-side boundary is calculated using a one-dimensional conductive heat model with an average surface heat flux of 80 µW/m⁻². As stated above, internal heat production is included in Cascadia models (Table 1). On the deeper part of the backarc side boundary, down to the approximate depth of the transition between mantle inflow and mantle outflow, we impose a geotherm calculated by a potential temperature of 1350 °C and an adiabatic gradient of 0.3 °C km⁻<fnref n="1" />. The surface temperature ($T_{sur}$) is 0 °C, and slab bottom temperature ($T_{bot}$) is 1450 °C.

Following Wada et al. (2012) and Cerpa et al. (2017), we estimate the pattern of fluid release from the slab by first computing the temperature field and then computing the water produced by dehydration reactions, using tabulated phase diagrams from the thermodynamic calculation code Perple_X (Connolly, 2009). We assume a hydrated layer at the uppermost part of the subducting plate arriving at the trench consisting of a 0.5-km-thick volcanic layer, a 1.5-km-thick layer of dikes, a 5-km-thick gabbros layer, and a 4-km-thick peridotitic uppermost mantle layer with initial water contents of 2.6, 1.8, 0.8, and 2.0 wt%, respectively. The uppermost part of the slab is divided into 1-km-wide columns, each consisting of 500-m-tall cells. Assuming vertical fluid flow within the slab, the integration over each column of the water release per cell provides an estimation of the volumetric flux $Q_{H_2O}$ of fluids released at the top of the slab as a function of depth (z). Because our study

focuses on fluid migration in the creeping mantle wedge, we exclude the "cold nose," that is, the cold corner tip between the subducting slab and the overlying crust that does not deform, from our calculations, and we estimate fluid release from the slab only below the MDD. In the calculation, the effects of updip migration (Wilson et al., 2014) and rehydration reactions (Wada et al., 2012) within the slab are neglected. The fluid influx to be applied in the model ($Q(z)$) is a product of the calculated $Q_{H_2O}(z)$ and a tuning parameter $\omega$, which controls the fraction of slab-derived fluids that are introduced into the mantle wedge to account for the fact that all fluids released from the slab do not migrate instantaneously into the mantle wedge. In Cerpa et al. (2017), we applied a few different fluid influx patterns $Q_{H_2O}(z)$ to test their effects on fluid pathways in the mantle wedge. The tuning parameter value that controls the magnitude of fluid influx was fixed for all the models and was chosen as the minimum value that allowed upward migration in all the models. Here, we investigate the influence of the tuning parameter $\omega$.

In the fluid-flow calculations, Equations (7) and (8) are not solved in the slab nor in the overriding crust. The non-dimensional volumetric flux $Q'(z) = Q(z)/\Psi_0$ is imposed as a boundary condition at the base of the mantle wedge. The top of the fluid domain through which fluids are allowed to escape is set at a 45 km depth. Above this depth, fluid migration mechanisms involving mechanical failure of rocks may play a role (Havlin et al., 2013; Keller et al., 2013) but are neglected in our models. As described in Cerpa et al. (2017), in the fluid calculations, we use the mantle shear viscosity from the solid calculation to which we impose a cap value of $10^{21}$ Pa s. Similarly, we impose a lower bound of 10 µm and an upper bound of 2 cm to the grain size in the fluid-flow models.

All calculations are performed using the software package TerraFERMA (Transparent Finite Element Rapid Model Assembler) (Wilson et al., 2017).

## 4. EFFECTS OF FLUID INFLUX AND FLUID PROPERTIES ON FLUID MIGRATION

We first present models with a simple subduction zone geometry with a constant slab dip of 30° and a 30-km-thick overriding crust. For the calculation of the thermal and solid-state mantle flow fields, we impose a subduction rate of 4 cm/yr and a slab age of 50 Ma. In the model (S1), full coupling between the slab and the overlying mantle downdip of the MDD induces a corner flow that replenishes the mantle wedge with the hot mantle from the backarc region (Fig. 2A; data used to create figures in Supplemental Material<fnref n="1" />). The model-predicted distributions of mantle wedge flow, shear viscosity, and grain size have been described in Cerpa et al. (2017). The key observations (Figs. 2B and 2C) that are important to fluid migration include (1) low shear viscosity (~$10^{18}$ Pa s) along the line separating the shallow incoming mantle flow toward the trench from the deeper downgoing mantle flow due to relatively high shear stresses and temperature; (2) moderate viscosity ($10^{19}$–$10^{20}$ Pa s) in the inflow region where the strain rate is relatively low; (3) increasing viscosity (>$10^{20}$ Pa s) toward the surface

<table><caption>TABLE 1. MAIN DIMENSIONAL MODEL PARAMETERS USED IN SOLID FLOW CALCULATIONS</caption><thead><tr><th>Models</th><th>Slab age (Ma)</th><th><i>v</i><sub>sub</sub><br>(cm/yr)</th><th>Radiogenic<br>heat<sup>a</sup></th></tr></thead><tbody><tr><td>S1</td><td>50</td><td>4.0</td><td>no</td></tr><tr><td>SC1</td><td>7</td><td>3.9</td><td>yes</td></tr></tbody></table>

aWhen radiogenic heat is applied, we use the following parameters:

Radiogenic heat (µW/m²):

1.3 (upper crust <15 km depth)

0.4 (lower crust)

0.02 (mantle and slab)

Density (kg/m³):

2750 (upper and lower crust)

3300 (mantle and slab)

Specific heat capacity (J/kg/K): 1250

<footnote n="1">¹Supplemental Material. Data used to generate figures. Please visit https://doi.org/10.1130/GES01660.S1 or access the full-text article on www.gsapubs.org to view the supplemental data. The code TerraFERMA is open-source and available at http://terraferma.github.io.</footnote>
<figure><img src="image_6.png" /></figure>

Figure 2. Solid-state mantle flow in model S1. The color scales represent the temperature (A); bounded mantle shear viscosity structure as imposed in the fluid-flow calculations (cap value of 10²¹ Pa s) (B); and bounded grain size as imposed in the fluid-flow calculations (cap value of 10²¹ Pa s) (B); and bounded grain size as imposed in the fluid-flow calculations (C). Thin white arrows are velocity fields. The location of the maximum decoupling depth (MDD) is given by a thick white arrow. Red triangle represents the idealized location of the arc defined as the point at the surface that lies at a 100 km depth above the top of the slab.

above the mantle inflow region because of an upward decrease in temperature; and (4) downdip and arcward increase in grain size from a few µm near the tip of the corner flow to a few cm in the hottest part of the mantle wedge.

The thermal structure of model S1 leads to release of fluids from the slab over two distinct depth ranges (Fig. 3). The first peak in the fluid release between 75 and 100 km depths is associated with the dehydration of the upper

crust while the second peak between 135 and 150 km depths is associated with the dehydration of the subducting mantle.

We use the above solid shear viscosity structure, grain-size distribution, and pattern of fluid release in the fluid-flow calculations with a simple slab geometry (models F1–11) as described in Section 3. In the calculations, we apply uniform fluid properties unless otherwise stated.

<figure><img src="image_7.png" /></figure>

Figure 3. (A) Maximum water content within the slab. (B) Volumetric flux of fluid release from the slab (bottom x axis) and cumulative water loss (top x axis) with depth (y axis) estimated at the top of the subducting slab.
4.1. Effects of Fluid Influx

We use a model with a fluid viscosity of 0.1 Pa s and a tuning parameter of 0.5% as a reference model (model F1). The model is identical to model 50E in Cerpa et al. (2017), and therefore we only summarize the main features. Fluids introduced during the first peak dehydration are initially dragged down-

dip (Fig. 4A) because of the small grain size near the MDD and the presence
of a layer of relatively high solid shear viscosity that develops just above
the slab. As grain size increases downdip, the mantle permeability and the
compaction length increase, and fluids start to migrate upwards at subarc
depths (Figs. 4A and 4B). Fluids released by the second peak dehydration
migrate largely upward because of the relatively high grain size in the deep

<figure><img src="image_8.png" /></figure>

Figure 4. (A–D) Porosity distribution (color) and fluid velocity below 60 km depth (white arrows) at time 2.0, 4.0, 6.0, and 10.0, respectively in model F1. (E) Porosity and fluid velocity time-averaged fields over $t \in [0; t]$ with $t$ the final time given in Table 2. Yellow and magenta lines are isocontours of grain size and solid shear viscosity, respectively. Red triangle indicates the idealized horizontal position of the volcanic front.
mantle wedge (Figs. 4A–4C). In the mantle inflow region, fluid pathways are deflected trenchward under the effects of the compaction pressure gradients and the solid advection (Fig. 4C). In the shallowest part of the mantle wedge where solid flow vanishes, fluids tend to pond due to the vertical increase in solid shear viscosity and decrease in grain size. In our fluid-flow models, fluids initially migrate through a dry mantle. As we are rather interested in mature fluid migration pathways (Figs. 4C and 4D), we use the time-averaged porosity field (Fig. 4E) to represent the “long-term” behavior of the fluids.

We now consider a model with a higher tuning parameter of 5% (model F2), i.e., a larger flux of fluids introduced at the base of the mantle wedge. All other parameters remain identical to model F1. In this model, fluids from the first peak dehydration are first dragged downdip as in model F1 and then start to migrate upward at subarc depths (Fig. 5A). However, the distance over which these fluids are dragged downdip is shorter compared to that in model F1. Some fluids from the second peak dehydration migrate immediately upward as in model F1, but there are also some fluids that first migrate updip along the base of the mantle wedge (Fig. 5A). Above the second peak dehydration, solid viscosity is relatively high, and compaction pressure gradients guide fluids to migrate updip through high permeability layers oriented parallel to the dip of the slab. Those fluids start to move upwards at a ~130 km depth, but during their ascent, some of them become entrained and dragged downdip by the outflowing mantle, merging back with the rising fluids from the second peak dehydration (Fig. 5B). This complex pattern equilibrates progressively, eventually leading to the formation of two main vertical fluid pathways separated by a horizontal distance of ~30 km from each other (Figs. 5C and 5D). In model F2, upward migration of fluids is less affected by mantle inflow and variations in solid shear viscosity in flowing mantle wedge than in model F1 owing to higher permeability due to increased fluid influx. In the cold shallowest mantle wedge, fluids tend to pond as in model F1. The time-averaged porosity field (Fig. 5E) shows wider fluid pathways than in model F1. In the following, only time-averaged porosity is discussed.

We test the effects of a high tuning parameter (50%) in model F3 (Fig. 6), while all other parameters remain identical to models F1 and F2. In this model, the distance over which fluids from the first peak dehydration are dragged is even shorter compared to that in models F1 and F2. Conversely, the distance of updip migration of some fluids from the second peak dehydration is longer than that in model F2, resulting in two vertical fluid pathways about a 50 km distance apart. In model F3, solid advection plays a modest role in controlling fluid migration.

## 4.2. Effects of Fluid Viscosity

### 4.2.1. Uniform Fluid Viscosity

Relative to the reference model F1, we decrease and increase fluid viscosity by two orders of magnitude to 0.001 in models F4 and 10 Pa s in F5, respectively (Fig. 6). All other parameters are identical to model F1. The fluid viscosities in models F4, F1, and F5 correspond to dissolved silicates concentrations ($c_{sil}$) of 20.5, 54.6, and 88.8 wt%, respectively, at a temperature of 800 °C based on the empirical law for aqueous fluid viscosity of Audôtat and Keppler (2004) (Fig. 1).

Decreasing fluid viscosity promotes greater upward migration from the first peak dehydration and updip migration of fluids from the second peak dehydration. Conversely, increasing fluid viscosity by two orders of magnitude leads to the entrapment of fluids at the base of the mantle wedge and their downdip drag by the mantle flow.

Although the absolute values of fluid fraction in model F4 differs from that in model F2, the two models show identical fluid migration pathways, indicating that decreasing fluid viscosity by two orders of magnitude has the same effect as increasing fluid influx by one order of magnitude. In other words, models with the same $\omega^2/\mu$ produce the same fluid pathways although the distributions of fluid fraction are different due to the difference in the amount of fluids introduced into the system. This scaling relationship in controlling compaction-driven flow is consistent with the relationship described in Equation (9) for the reference compaction length, fluid flux, and fluid viscosity; the compaction length is scaled with $\Psi_0^2/\mu_0$, where $\Psi_0$ plays an equivalent role to $\omega$ as the latter controls fluid influx.

We include two additional models F6 and F7 to demonstrate this scaling relationship. In these models, $\omega$ and/or $\mu$ take different values from those in models F1–F5, but $\omega^2/\mu$ in F6 is identical to that of model F1, and that in F7 is identical to those in model F2 and F4, and models that share the same value show identical fluid migration paths. The tuning parameter $\omega$ controls fluid influx and thus the porosity and permeability at the base of the mantle wedge, affecting fluid mobility and compaction length. As a consequence, for a given grain-size distribution, the impact of matrix compaction is mainly controlled by $\omega^2/\mu$.

Models with $\omega^2/\mu < 2.5 \times 10^{-6} Pa^{-1} s^{-1}$ (model F5) lead to entrapment of fluids at the base of the mantle wedge and little upward fluid migration. In models with $\omega^2/\mu < 2.5 \times 10^{-4} Pa^{-1} s^{-1}$, downdip drag is still prominent for fluids from the first peak dehydration, but upward fluid migration develops. Models with $\omega^2/\mu < 2.5 \times 10^{-2} Pa^{-1} s^{-1}$ show a lesser effect of solid advection and eventually promote some updip fluid migration of fluids from the second peak dehydration at the base of the mantle wedge.

### 4.2.2. Temperature-Dependent Fluid Viscosity

In nature, fluid viscosity is expected to vary spatially because of its dependence on fluid composition, pressure, and temperature. Variation in the composition of hydrous fluids that are traveling through the mantle wedge depends largely on silicate solubility, which increases with both the temperature and the pressure (e.g., Newton and Manning, 2002). As fluids travel to shallower depths from the base of the mantle wedge to its hot core, pressure decreases but temperature increases. Given the competing effects of
<figure><img src="image_9.png" /><figcaption>Figure 5. (A–D) Porosity distribution (color) and fluid velocity below 60 km depth (white arrows) at time 0.8, 2.2, 4.0, and 6.0, respectively, in model F2. (E) Porosity and fluid velocity time-averaged fields over $t \in [0; t]$ with $t$ the final time given in Table 2. Please refer to Figure 4 for other legend details.</figcaption></figure>
<figure><img src="image_10.png" /></figure>

Figure 6. Diagram showing porosity (color) and fluid velocity below 60 km depth (white arrows) for models F1–F7 as a function of volumetric flux tuning parameter (ω) and fluid viscosity (μ). Models are averaged over $t \in [0; t_f]$ with $t_f$ the final time given in Table 2. Please refer to Figure 4 for other legend details.

temperature and pressure on silica solubility, the change in fluid composition may be relatively small. The direct effect of pressure on fluid viscosity is expected to be also relatively small, whereas temperature has a large impact on fluid viscosity as fluids travel between the relatively cold base and the hot core of the mantle wedge (Audôtat and Keppler, 2004). Here, we investigate the effects of temperature-dependent fluid viscosity on fluid migration.

Two end-member models are studied: model F8 with relatively low silicate content $c_{sil} = 30$ wt% and model F9 with relatively high silicate content $c_{sil} = 70$ wt% (Fig. 7). In both models, $ω = 5\%$. Near the “cold nose,” fluid viscosity is expected to be relatively high, while it is expected to be relatively low in the flowing part of the mantle wedge (Figs. 7A–7C). The relatively cold temperature in the shallow mantle wedge beneath the overlying crust and immediately above the 50 Ma slab also results in relatively high fluid viscosity. Fluids are thus expected to be less mobile at the base of the mantle wedge than in its hottest part.

In model F8 with $c_{sil} = 30$ wt%, fluid viscosity at the base of the mantle wedge is 0.01–0.1 Pa s, similar to that in model F2, and it becomes lower (~10⁻³ Pa s) in the hot part of the mantle wedge, resulting in a higher fluid mobility in this region compared to model F2. Therefore, model F8 exhibits some downdip drag of fluids from the first peak dehydration at the base of the mantle wedge over a distance similar to that in model F2. The two vertical fluid pathways

beneath the backarc are ~40 km apart and are spaced closer than those in model F2. This indicates that fluids from the second peak dehydration travel farther updip compared to model F2. This is consistent with the higher mobility of fluids expected in the model with temperature-dependent viscosity as they migrate away from the base of the mantle wedge. Fluid pathways at shallow depths are narrower in model F8 than in model F2 because buoyancy-driven upward fluid flow becomes more dominant than fluid migration due to mantle inflow and compaction-driven flow.

In model F9 with $c_{sil} = 70$ wt%, fluid viscosity in the hottest part of the mantle wedge is predicted to be comparable (~0.1 Pa s) to that of model F2 but higher (10–100 Pa s) at the base of the mantle wedge. As a consequence, fluids from the first peak dehydration become dragged downdip farther than that in model F2, and updip migration above the second peak dehydration does not occur. At the base of the mantle wedge, the effects of solid advection and spatial variations in solid shear viscosity are similar to that in model F2.

Models F8 and F9 show that a temperature-dependent viscosity promotes lateral transport of fluids away from the trench at the base of the mantle wedge. This results in more focusing of fluids that are released beneath the forearc region toward the subarc region. On the contrary, a temperature-dependent viscosity reduces trenchward deflection and subsequent spread of
<figure><img src="image_11.png" /></figure>

Figure 7. (A) Predicted fluid viscosity as a function of temperature for (A) model F8 and (C) model F9. (B) Time-averaged porosity field for model F8 and (D) model F9. Porosity is averaged over $t \in [0; t_f]$ with $t_f$ the final time given in Table 2. Please refer to Figure 4 for other legend details.

upward fluid migration pathways by mantle inflow at shallow depths. However, as discussed by Cerpa et al. (2017), fluid migration pathways in and above the hottest part of the mantle wedge may differ from our predictions since melting may produce more viscous and less dense fluids.

## 4.3. Effects of Fluid Density

To test the effects of fluid buoyancy, we reduce the density contrast between the fluid and the solid from the reference value of 2000 kg m⁻³ to 1000 kg m⁻³ in model F10 and to 200 kg m⁻³ in model F11 (Fig. 8). Because $\Delta\rho$ is used for the non-dimensionalization of governing equations, we modify the reference flux $\Psi_0$ and thus the tuning parameter $\omega$ accordingly (see Table 2) so that the reference compaction length and boundary conditions in models F10 and F11 are identical to those in model F2 (see Appendix B).

Decreasing the density contrast by a factor of two (model F10) compared to model F2 generates fluid pathways that are similar to those obtained in models F1 and F6. Decreasing density contrast even further, by one order of magnitude (model F11), does have a noticeable effect on fluid pathways, which is comparable to the effect of increasing fluid viscosity by two orders of magnitude (model F5).

However, such small density contrast between the fluids and the mantle rock is unlikely, particularly below the core part of the mantle wedge. Therefore, the effect of fluid density variation on fluid migration pathways from the base to the core part of the mantle wedge may be small compared to the effects of variations in fluid viscosity.

## 5. FLUID PATHWAYS IN THE NORTHERN CASCADIA SUBDUCTION ZONE

The Cascades subduction zone (39°N–53°N) is formed by the subduction of the Juan de Fuca plate beneath the North American plate and is often considered as an end-member warm-slab subduction zone because of the relatively young age (ca. <20 Ma) of the subducting plate (Kirby et al., 1991; Peacock and Wang, 1999; Hacker et al., 2003). Thermal models for Cascadia predict that the subducting Juan de Fuca slab releases most of its water before it reaches the subarc depths (Wada and Wang, 2009; Syracuse et al., 2010; van Keken et al., 2011). However, geochemical and geophysical observations indicate the addition of slab-derived fluids to the subarc mantle, leading to the generation of arc magmas beneath the Cascades.
<figure><img src="image_12.png" /><figcaption>Figure 8. Porosity field time-averaged over $t \in [0; t_f]$ for models F10 and F11 with $t_f$ the final time given in Table 2. Please refer to Figure 4 for other legend details.</figcaption></figure>

Based on geochemical analyses of primitive basalts, Schmidt et al. (2008) divided the Cascades margin into four different segments (North, Columbia, Central, and South segments). In particular, the composition of the most primitive calk-alkaline basalts in the Columbia segment between Mount Rainier (46.8°N, 121°W) and Mount Jefferson exhibits high $^{87}$Sr/$^{86}$Sr and Ba/Ce ratios compared to other basalts from the Cascades and may reflect the effects of flux melting (Schmidt et al., 2008, and references therein). Although less pronounced in the Columbia segment compared to the southern part of the Cascades, relatively high ratios of fluid-mobile to fluid-immobile trace elements indicate the addition of slab-derived fluids. There is also some evidence of dry melting in the backarc region ~46°N (Leeman et al., 2005).

<table><caption>TABLE 2. MAIN DIMENSIONAL MODEL PARAMETERS USED IN THE FLUID-FLOW CALCULATIONS</caption><thead><tr><th>Models</th><th>ω<br>(%)</th><th>μ<br>(Pa s)</th><th>Δρ<br>(kg m<sup>-3</sup>)</th><th>Ψ<sub>0</sub><br>(m s<sup>-1</sup>)</th><th>t<sub>f</sub><br>(Myr)</th></tr></thead><tbody><tr><td>F1</td><td>0.5</td><td>0.1</td><td>2000</td><td>10<sup>-12</sup></td><td>12</td></tr><tr><td>F2</td><td>5</td><td>0.1</td><td>2000</td><td>10<sup>-12</sup></td><td>8</td></tr><tr><td>F3</td><td>50</td><td>0.1</td><td>2000</td><td>10<sup>-12</sup></td><td>6</td></tr><tr><td>F4</td><td>0.5</td><td>0.001</td><td>2000</td><td>10<sup>-12</sup></td><td>8</td></tr><tr><td>F5</td><td>0.5</td><td>10</td><td>2000</td><td>10<sup>-12</sup></td><td>8</td></tr><tr><td>F6</td><td>5</td><td>10</td><td>2000</td><td>10<sup>-12</sup></td><td>12</td></tr><tr><td>F7</td><td>50</td><td>10</td><td>2000</td><td>10<sup>-12</sup></td><td>8</td></tr><tr><td>F8</td><td>5</td><td>var</td><td>2000</td><td>10<sup>-12</sup></td><td>8</td></tr><tr><td>F9</td><td>5</td><td>var</td><td>2000</td><td>10<sup>-12</sup></td><td>14</td></tr><tr><td>F10</td><td>0.5</td><td>0.1</td><td>200</td><td>0.1 × 10<sup>-12</sup></td><td>20</td></tr><tr><td>F11</td><td>2.5</td><td>0.1</td><td>1000</td><td>0.5 × 10<sup>-12</sup></td><td>18</td></tr><tr><td>FC11</td><td>5</td><td>0.1</td><td>2000</td><td>10<sup>-12</sup></td><td>10</td></tr><tr><td>FC12</td><td>5</td><td>10</td><td>2000</td><td>10<sup>-12</sup></td><td>11</td></tr><tr><td>FC13</td><td>50</td><td>10</td><td>2000</td><td>10<sup>-12</sup></td><td>10</td></tr></tbody></table>

Note: var refers to the use of the empirical law of Audôtat and Keppler (2004).

Magnetotellurics (MT) has been used to infer the distribution of fluids in subduction zones, given the sensitivity of electrical conductivity (resistivity) to the presence of fluids and melts (Ni et al., 2011; Pommier and Garnero, 2014). The electrical resistivity model for the region beneath Mount Rainier of the central Cascades shows a high-conductivity anomaly originating just above the slab at ~80 km depth and extending both upward and downdip with decreasing values in those directions (McGary et al., 2014). Those anomalies were interpreted as indicators of the presence of slab-derived fluids and melts that are generated by flux melting. Another electrical resistivity model for the same region also indicates a deeper high-conductivity anomaly above the slab situated beneath the backarc region at 150–200 km depths (Wannamaker et al., 2014), which is interpreted as a region of decompression melting of upwelling mantle, consistent with geochemical studies in this region (e.g., Leeman et al., 2005). Both MT models of McGary et al. (2014) and Wannamaker et al. (2014) predict the presence of fluids in a region that extends from the base of the mantle wedge to the overlying crust and is deflected toward the trench at shallow depths. This feature was suggested by McGary et al. (2014) to represent the pathways of buoyant diapirs deflected by mantle inflow. The MT results, however, do not provide information on the direction of fluid migration. Quantitative models of fluid flow in the mantle wedge will therefore help to improve the interpretation of MT observations and other geophysical observations, such as seismic velocity and attenuation structures. Here, as a first step, we develop fluid migration models for Cascadia along an E-W profile at 47°N (Fig. 9) to complement the existing MT observations and test the effects of ω²/μ on fluid migration and melt generation.
<figure><img src="image_13.png" /></figure>

Figure 9. (A) Map of Cascadia subduction zone. Gray lines are iso-contours of slab surface depth with corresponding values in km (from McCrory et al., 2012). Red triangles are major arc volcanoes. Blue line is a W-E slab profile at latitude 47°N extracted for constructing slab geometry. (B) Slab surface in cross-sectional area. Solid line is profile from McCrory et al. (2012), and dashed line is an extrapolation of slab surface to 170 km depth. (C) Resistivity model of McGary et al. (2014) along the profile line shown in (A) (modified from their figure 2b). The areas of primary fluid migration in the mantle wedge as inferred from low-resistivity anomalies are labeled A to D. Red triangle and red circles indicate the projected locations of Mount Rainier and seismicity within 20 km distance from the profile line, respectively.
For the model construction, we adopt the slab geometry model for Cascadia from McCrory et al. (2012). The model extends only to a 95 km depth, and we extend the slab to the bottom of our model, assuming a constant slab dip of 30.3°. In the solid flow calculation, we apply a subduction rate of 3.9 cm/yr, which corresponds to a profile-parallel component of the subduction velocity of the Juan de Fuca plate (4.5 cm/yr) (Wada and Wang, 2009). The subducting plate age at the trench is 7 Ma (Syracuse et al., 2010), and the thickness of the

overriding crust is 40 km. Other model parameters for the solid flow calculations are identical to those listed in Table 1.

Compared to the reference model S1 with a 50 Ma slab, the slab in the Cascadia model is much warmer. For instance, the 600 °C isotherm is advected to ~70 km depth in the Cascadia model (Fig. 10A), while in model S1, the same isotherm extends beyond the bottom of the model domain at 160 km depth. Overall, the grain-size distribution is similar to that in model S1, but grain size

<figure><img src="image_14.png" /></figure>

**Figure 10.** (A–C) Solid-state mantle flow for Cascadia model SC1 with temperature field (A), bounded mantle shear viscosity structure as imposed in the fluid-flow calculations (B), and bounded grain size as imposed in the fluid-flow calculations (C). Thin white arrows are velocity fields. Thick black arrow indicates the maximum decoupling depth (MDD). A red triangle indicates an idealized horizontal position of Mount Rainier from Figure 4A in (McCrory et al., 2012). (D) Estimated volumetric flux $Q_{H,0}$ at the top of the slab and cumulativater $M_{H,0}$ loss with depth. The gray area corresponds to depths above the MDD that are neglected in our study.
is slightly larger particularly at the base of the mantle wedge, owing to the warmer slab thermal state (Fig. 10B). Similarly, the mantle wedge viscosity structure in model SC1 (Fig. 10C) is similar to that in model S1, but the viscosity is slightly lower in model SC1.

The volumetric flux of fluids from the top of the slab in model SC1 is calculated using the same approach described above (Fig. 10D). Because of the high temperatures within the slab, dehydration of the subducting crust occurs beneath the "cold-nose," but a peak dehydration of the subducting mantle occurs between 75 and 95 km depths, consistent with earlier petrological modeling results for Cascadia (Wada and Wang, 2009; van Keken et al., 2011).

We present three fluid migration models for Cascadia (models FC11-FC13; Fig. 11). In model FC11 and FC13, $ω$ is 5 and 50%, respectively, and μ is 0.1 and 10 Pa s, respectively, yielding $ω^2/μ$ of $2.5 \times 10^{-2}$ Pa$^{-1}$ s$^{-1}$. In model F12, $ω$ is 5, and μ is 10 Pa s, yielding $ω^2/μ$ of $2.5 \times 10^{-4}$ Pa$^{-1}$ s$^{-1}$. As discussed above, the exact values of $ω$ and μ are not critical to the fluid migration pathways, and the value of $ω^2/μ$ depicts the fluid migration paths. In all three models, a constant density contrast of 2000 kg m$^{-3}$ is assumed. In models FC11 and FC13, fluids immediately flow upward at subarc depths with little initial downdip transport along the base of the mantle wedge. During their ascent, fluids are little affected by solid advection and spatial variations in solid shear viscosity in the inflow region. In model FC12, fluids are transported downdip and then migrate

upward at postarc depths (~90-100 km depths), and some fluids are dragged farther downdip. The deflection of fluids toward the trench by the incoming mantle flow is more pronounced in model FC12. All three models predict some focusing of fluids near the subarc region.

## 6. DISCUSSION

### 6.1. Physical Factors Influencing Fluid Pathways in the Mantle Wedge

Wilson et al. (2014) investigated the influence of mantle compaction on fluid migration pathways in subduction zones. In their model, grain size was assumed constant, and mantle permeability varied only with fluid fraction. Thereby, the compaction length was controlled by the reference fluid mobility (i.e., the ratio of matrix permeability to fluid viscosity $K_0/μ_0$). They showed that, with relatively high fluid mobility ($K_0/μ_0 \sim 10^{-7}$ m$^3$ s$^{-1}$) and thus a high compaction length ($\delta_0 \sim 20$ km), fluid migration paths are strongly influenced by the solid rheology, while with relatively low fluid mobility ($K_0/μ_0 \sim 10^{-11}$ m$^3$ s$^{-1}$) and thus a low compaction length ($\delta_0 \sim 2$ km), the influence of solid advection becomes important.

<figure><img src="image_15.png" /></figure>

Figure 11. Diagram showing porosity (color) and fluid velocity below 60 km depth (white arrows) for models FC11-FC13 as a function of tuning parameter ($ω$) and fluid viscosity (μ). Models are averaged over $t \in [0; t]$ with $t$, the final time given in Table 2. Please refer to Figure 6 for other legend details.
In our recent work (Cerpa et al., 2017), we have investigated the combined effects of grain-size distribution and mantle wedge viscosity structure on fluid-flow pathways. In these models, the reference mobility is ~10⁻⁹ m³ s kg⁻¹, which is equivalent to moderate mobility in Wilson et al. (2014). However, the spatial variation in grain size caused the fluid mobility and thus the compaction length to vary spatially. Fluid mobility is relatively low near the MDD where grain size is small and increases with the distance away from the MDD. As a result, fluids introduced beneath the forearc are initially trapped at the base of the mantle wedge and dragged downdip by mantle flow. The downdip increase in permeability increases fluid mobility and thus decreases the effectiveness of solid advection, allowing the upward migration of fluids at subarc depths as shown in our model F1 (Fig. 4). In other regions of the mantle wedge, where permeability is relatively high, fluids are largely unaffected by solid advection except in regions where variations in solid shear viscosity promote lateral fluid migration through its effects on compaction pressure gradients and solid advection. In the latter regions, some fluids migrate toward the trench and may contribute to the focusing of some fluids released at postarc depths toward the subarc region. Overall, spatial variations in both grain size and mantle viscosity help to focus fluids beneath the arc, but their contribution to the focusing varies spatially.

In the present work, we have investigated the effects of fluid influx in terms of ω and fluid viscosity on fluid migration pathways. The results illustrate the scaling relationship between ω and μ such that the fluid migration pathways are identical in models with the same ω²/μ value, although the distribution of fluid fraction varies with ω. A low ω²/μ value (2.25 × 10⁻⁶ Pa⁻¹s⁻¹) leads to the entrapment of fluids at the base of the mantle wedge and their downdip transport by the mantle flow, inconsistent with the predictions of flux melting occurring beneath the arcs. A very high ω²/μ value (2.5) leads to mostly upward fluid migration. At intermediate to high ω²/μ values (2.5 × 10⁻⁴–2.5 × 10⁻² Pa⁻¹s⁻¹), fluids released beneath the forearc are transported downdip to the subarc depths by solid advection before migrating upward, and some fluids released during a second peak dehydration at postarc depths are laterally advected toward the arc in the mantle inflow region. Increasing ω²/μ (>2.5 × 10⁻² Pa⁻¹s⁻¹) reduces the lateral transport of fluids by solid advection and the widening of the pathways of rising fluids through the mantle wedge beneath the backarc.

Wilson et al. (2014) observed the formation of updip fluid-flow channels in the subducting slab in models with high fluid mobility and thus a high compaction length. Our models predict, instead, that some fluids released at postarc depths can flow updip near the base of the mantle wedge if ω²/μ is relatively high, providing a new mechanism for focusing fluids beneath the subarc region. Geophysical observations in NE Japan suggest the existence of a dipping low-velocity layer near the base of the mantle wedge (Kawakatsu and Watada, 2007. This layer may indicate the presence of free fluids, but whether fluids are migrating downdip or updip would be difficult to resolve. Our models show that the direction of fluid migration at the base of the mantle wedge at postarc depths depends on ω²/μ. Alternatively, it has been proposed that the

low-velocity layer just above the slab represents the downdip drag of hydrous phases that are stable up to 150 km depth (Iwamori, 1998; Grove et al., 2009; Horiuchi and Iwamori, 2016), although thermomechanical models with realistic rheologies for the mantle wedge predict that the temperatures are too high for most hydrous phases to remain stable beyond a depth of 70–80 km (e.g., van Keken et al., 2011; Wada et al., 2012).

A common simplification made in fluid migration models for subduction zone is the use of a constant fluid viscosity as in our models F1 to F7. In nature, however, fluid viscosity may change spatially because of its dependence on pressure, temperature, and composition of the fluids. As a first step in incorporating fluid viscosity variation, we applied temperature-dependent fluid viscosity in models F8 and F9. The results show that the relatively cold temperature just above the slab and the subsequent low mobility may contribute to the drag of fluid released beneath the forearc toward the subarc region. However, the fluid migration paths in these models do not differ significantly from those in the models with a constant fluid viscosity. Production of silicic melts in the fluxed region may produce important changes in fluid properties but is neglected in our models. Hence, fluid migration pathways in the shallow mantle wedge may differ from the predictions by our models.

Although a relatively large effect was observed when the fluid density was decreased by one order of magnitude, it is unlikely to decrease by more than a factor of two in reality (Hack and Thompson, 2011). We can thus conclude that in a realistic range of fluid density in the mantle wedge, its effect on fluid migration is small compared to the effects of fluid viscosity and fluid influx on fluid migration pathways. Nonetheless, the results suggest that the combination of high fluid density and fluid viscosity due to high silicate content in the fluxed region of the mantle can have a large impact on fluid flow in the shallow part of the mantle wedge, enhancing the effect of solid advection. Hence, as we have suggested in Cerpa et al. (2017), the focusing of fluids toward the subarc region may be greater than what is predicted by our models.

## 6.2. Fluid Distribution and Hydrous Melting

The supply of water from the subducting slab to the hot mantle wedge is thought to trigger hydrous melting (Sisson and Grove, 1993; Grove et al., 2006), and the degree of melting likely increases with the amount of water added to the system (Stolper and Newman, 1994; Hirschmann et al., 1999; Katz et al., 2003; Kelley et al., 2010). Although, our models do not incorporate melting in the calculations, we use the parameterized melting function of Kelley et al. (2010), following the approach of Wilson et al. (2014), to calculate the fraction of hydrous melts that can be produced in the mantle wedge based on the results of selected fluid-flow models that generate some upward fluid migration and time-averaged porosities ≥10⁻⁴ (models F1, F2, F3, F6, and F7). The degree of melting is calculated by using fluid content, temperature, and pressure predicted by our fluid migration models. The parameterized melting function is given in terms of temperature, pressure, and water content. We use
the temperature field from model S1 (Fig. 2). A lithostatic pressure based on solid density of 3300 kg m⁻³ (Δρ = 2000 kg m⁻³) is assumed for the calculation of pressure. Water concentration is calculated using the time-averaged fluid fraction of the fluid-flow models. We exclude from the calculations the shallowest part of the mantle wedge where fluid accumulation due to ponding is likely overestimated by our models without inclusion of the effects of brittle deformation on fluid extraction (Cerpa et al., 2017).

The fluid-flow models predict the formation of two or three main pathways through the hottest part of the mantle wedge. The first main fluid pathway forms at relatively shallow depth (<100 km) beneath the arc (models F1 and F6) or beneath the forearc at a distance up to ~25 km from the arc, if fluid mobility is high (model F3). A second main fluid pathway is observed in all models beneath the backarc region and generally forms above the location of the second peak dehydration situated between 135 and 150 km depths. A third main pathway eventually appears beneath the backarc region as a consequence of the updip flow of some fluids from the second peak dehydration, if fluid mobility is high (models F2, F3, and F7).

Our calculations predict that hydrous melting occurs beneath the arc and the backarc (Fig. 12). Beneath the arc, melt production is observed in a relatively narrow region (width <20 km) between 70 and 85 km depths near the updip sloped isocontours of solid shear viscosity in the inflow region of the mantle wedge. Beneath the backarc, melt production develops in a wider region than in the subarc as the fluid-flow calculations predict wider fluid pathways in the former region. Overall, melt fraction is relatively low (0.1%–1%), reaching 1%–4% only in models F5 and F7. At a given value of ω²/μ (identical fluid pathways), the highest melt fractions are observed for the models with the highest ω.

The degree of melting in the mantle wedge inferred from geochemical data on arc lavas and volcanics ranges from 1% to 20% (Portnyagin et al., 2007; Kelley et al., 2010; Watt et al., 2013). Our models generally predict a lower degree of melting in the mantle wedge (≤1%). This can be explained by underestimation of the amount of slab-derived fluids in the hottest part of the mantle wedge or the temperature or both. Increasing tuning parameter increases the predicted melt fraction but requires relatively high fluid viscosity to cause focusing of fluids beneath the arc (model F7). Alternatively, mantle potential temperature may be higher than that assumed in models F1–F11 ($T_m$ = 1350 °C) and may be closer to the dry solidus temperature of mantle peridotite. Moreover, experimental studies on melting of natural peridotite samples in the presence of water indicate that low-degree melts may be produced at the base of the mantle wedge at relatively low temperatures of ~800–900 °C below an 80 km depth (Grove et al., 2006; Till et al., 2012). The latent heat released by such melting would be advected by the rising melts, locally increasing the temperature in the core part of the mantle wedge and promoting further melting (Rees-Jones et al., 2018).

Our models with relatively old subducting slabs also predict hydrous melting beneath the backarc. In nature, volcanoes do form behind the volcanic front above relatively old subducting slabs, such as in Northeast Japan (Tatsumi et al., 1983), Marianas (Kelley et al. 2010), and Kamchatka (e.g., Tatsumi et al., 1994), indicating the presence of magma as predicted by our models As discussed in Cerpa et al. (2017), some of the melts that are produced beneath the backarc may also be redirected arcward by the incoming mantle flow, particularly when the fluid viscosity and density are relatively high (i.e., low fluid mobility and buoyancy). Gradients of compaction pressure induced by melting alone may also facilitate lateral

<figure><img src="image_16.png" /><figcaption>Figure 12. Diagram showing hydrous melt fraction (color) for models F1–F7 as a function of volumetric flux tuning parameter (ω) and fluid viscosity (μ) slab and fluid viscosity. Melt fraction is calculated below 60 km depth using the time-averaged fluid porosity (Fig. 6) and the parametric function of Kelley et al. (2010).</figcaption></figure>
transport of melts as proposed for melt focusing at mid-ocean ridges (Turner et al., 2017).

### 6.3. Northern Cascadia

Our fluid-flow models for Cascadia (models FC11–FC13) (Fig. 11) display fluid pathways comparable to those inferred from MT models. Our solid-state flow model for Cascadia predicts the release of fluids into the flowing mantle wedge between 75 and 95 km depths, consistent with the highest electrical conductivity observed just above the slab between 80 and 100 km depths by (McGary et al., 2014). The high-conductivity anomaly at the base of the mantle wedge extends downdip up to a 120 km depth, and there appears to be little fluid present beyond the 120 km depth, if any. Similarly, our models with the highest fluid mobility (FC11 and FC13) exhibit modest fluid fraction beyond 110 km depth. However, the MT model of Wannamaker et al. (2014) for the same Cascadia profile indicates another slab-parallel high-conductivity anomaly at the base of the mantle wedge between a 150 km and a 200 km depth. These authors interpreted this anomaly as a region of decompression melting of upwelling mantle. Our models (models FC11–FC13) predict that some fluids do

become dragged downdip by mantle flow to greater depth, potentially contributing to the deeper high-conductivity anomaly.

In all three Cascadia models, despite two different fluid mobilities (FC11 and FC13 versus FC12), upward migration of fluids occurs in the hottest part of the mantle wedge near the subarc region as observed by McGary et al. (2014) and Wannamaker et al. (2014). Their MT models indicate some lateral deflection of upward fluid pathways that may be explained by the combined effect of mantle inflow and vertical variations in mantle shear viscosity as shown in our models. However, our Cascadia models display upward fluid pathways east of Mount Rainier, while the MT models indicate upward pathways west of Mount Rainier. This discrepancy points to two possible explanations. One is that the depth of fluid influx is shallower in the actual system, potentially due to higher-than-expected slab temperature and thus shallower fluid release or the updip fluid migration within the subducting slab that leads to fluid influx at shallower depths than in the model. The other possible explanation is that the permeability structure at the base of the mantle wedge updip of the influx region is higher in the actual system, for example, due to shear-induced permeability (e.g., Précigout et al., 2017) or reaction-induced fracturing (Okamoto and Shimizu, 2015), allowing updip migration along the base of the mantle wedge.

<figure><img src="image_17.png" /><figcaption>Figure 13. Diagram showing hydrous melt fraction (color) for models FC11–FC13 below 45 km depth. Please see Figure 12 for other legend details.</figcaption></figure>
We use the time-averaged porosity distribution of our Cascadia models to predict the regions of hydrous melting (Fig. 13). The model with the lowest fluid viscosity (model FC11) predicts negligible melting. The two models with higher fluid viscosity (models FC12 and FC13), predict 0.1%–1% degree of melting in the subarc region. This relatively low degree of melting predicted by these models is consistent with the relatively high ratios of incompatible to compatible elements found in Cascadia arc lavas (e.g., Ce/Yb; Schmidt et al., 2008). In hot-subduction zones such as Cascadia, relatively high fluid viscosity is expected because of the temperature dependence of silicates solubility in fluids. Further, fluids released from the dehydration of the subducting mantle in the subarc region may trigger slab melting as they travel through the overlying crust as proposed to explain geochemical data in Southern Cascadia (Walowski et al., 2015), increasing the viscosity of the fluids before they migrate into the mantle wedge. Thus, as indicated by models FC12 and FC13 (Fig. 13), with relatively high fluid viscosity, a modest fraction of slab-derived fluids from the slab can explain 0.1%–1% melting beneath Mount Rainier, and melting occurs in a relatively small region compared to the distribution of fluids in the mantle wedge (Fig. 11).

## 7. CONCLUSIONS

In our previous study (Cerpa et al. [2017]), we showed that the spatial variations in grain size (permeability) and mantle shear viscosity play a key role on fluid migration in the mantle wedge. In the present study, we assess the influence of the magnitude of fluid influx through the parameter ω and fluid properties in our fluid migration models. In these models, fluids are generally introduced from the base of the mantle wedge beneath the forearc and the backarc. Our results suggest that fluid influx and fluid viscosity (μ) have a large effect on fluid migration pathways, while density contrast plays a secondary role. In our models, mantle compaction is mainly controlled by ω²/μ for a given grain-size distribution. Low values of ω²/μ produce entrapment of fluids at the base of the mantle wedge. Increasing ω²/μ allows upward migration of fluids. In particular, the downdip drag of fluid released beneath the forearc decreases with increasing ω²/μ. At high values, the models predict updip flow of fluids at the base of the mantle wedge beneath the backarc. Further, temperature-dependent fluid viscosity enhances the downdip drag of fluids by mantle flow. Based on the averaged porosities predicted by our models, we estimate relatively low degree of melting (0.1%–5%) occurring beneath the arc and beneath the backarc. Finally, our model for Northern Cascadia predicts fluid distribution in the mantle wedge that is consistent with electromagnetic imaging of the mantle wedge.

## ACKNOWLEDGMENTS

We thank T. Gerya and an anonymous reviewer for comments that helped to improve the manuscript. I.W. acknowledges financial support from the University of Minnesota in the form of

start-up funds, which also supported N.G.C. C.W. was supported by NSF grant OCE-1358091. The code TerraFERMA is open-source and available at http://terraferma.github.io.

# APPENDIX A. PARAMETERS IN SOLID-FLOW CALCULATIONS

## Solid Shear Viscosity

The solid-state mantle flow model accounts for diffusion creep and dislocation creep, which are assumed to be independent. The solid shear viscosity η is thus expressed as

$$ \frac{1}{\eta} = \left( \frac{1}{\eta_{\text{diff}}} + \frac{1}{\eta_{\text{disl}}} + \frac{1}{\eta_{\text{max}}} \right), \quad (14) $$

where ηmax is a maximum shear viscosity set to 10²⁴ Pa s, and ηdiff and ηdisl are diffusion and dislocation shear viscosity, respectively, and are defined as

$$ \eta_r = A b^{\alpha} C_{\text{oh}}^{\beta} \exp\left(\frac{E}{RT}\right) \left(\frac{1}{\nu_{\text{eff}}}\right)^{\frac{1}{\gamma}} \left(\frac{\epsilon_{\text{eff}}}{\epsilon_{\text{0}}}\right)^{\gamma}, \quad (15) $$

where the subscript i denotes creep mechanism, b is grain size, and $\epsilon_{eff}$ is the second invariant of strain rate. The rheological parameters λ, Coh, E, and R are pre-exponential factor, water content, activation energy, and gas constant, respectively. Parameters α, β, and γ are grain size exponent, water content exponent, and stress exponent, respectively. For these parameters, we use values that are determined experimentally for wet olivine (Hirth and Kohlstedt, 2003) (Table A1).

<table><caption>TABLE A1. RHEOLOGICAL PARAMETER VALUES FOR WET OLIVINE (HIRTH AND KOHLSTEDT, 2003)</caption><thead><tr><th></th><th>A<br>(Pa<sup>3</sup> s<sup>-1</sup>)</th><th>C<br>(H/10<sup>6</sup> Si)</th><th>α</th><th>β</th><th>γ</th><th>E<br>(J mol<sup>-1</sup>)</th></tr></thead><tbody><tr><td>Creep</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Diffusion</td><td>0.333</td><td>1000</td><td>3</td><td>1</td><td>1</td><td>3.35 × 10<sup>5</sup></td></tr><tr><td>Dislocation</td><td>3.0 × 10<sup>-20</sup></td><td>1000</td><td>0</td><td>1.2</td><td>3.5</td><td>4.80 × 10<sup>5</sup></td></tr></tbody></table>

## Grain Size

We calculate the steady-state grain size as (Behn et al., 2009)

$$ b = \left| \frac{G c \lambda \exp\left(\frac{E}{RT}\right)}{2 p_g \chi \sigma: \dot{\epsilon}_{\text{eff}}} \right|^{\frac{1}{\gamma}}, \quad (16) $$

where σ and $\dot{\epsilon}_{disl}$ are stress tensor and strain rate tensor associated with dislocation creep, respectively, G is grain growth constant, λ is average specific grain boundary energy, E_g is activation energy for grain growth, C is a geometric constant, p_g is grain growth exponent, and χ represents the fraction of work done by dislocation creep associated with changing grain boundary area. The value of the parameters used for the calculations are provided Table A2.

<table><caption>TABLE A2. PARAMETER VALUES FOR STEADY-STATE GRAIN-SIZE CALCULATION (BEHN ET AL., 2009)</caption><thead><tr><th>G<br>[m<sup>3</sup> s<sup>-1</sup>]</th><th>c</th><th>λ<br>[K m<sup>-2</sup>]</th><th>E<sub>g</sub><br>[J mol<sup>-1</sup>]</th><th>P<sub>g</sub></th><th>χ</th></tr></thead><tbody><tr><td>4.5 × 10<sup>-4</sup></td><td>3</td><td>1</td><td>3.5 × 10<sup>5</sup></td><td>3</td><td>0.1</td></tr></tbody></table>
# APPENDIX B. SCALING RELATIONSHIPS

The reference compaction length is given by Spiegelman (1993):

$$ \delta_0 = \sqrt{\frac{K_0 \phi_0^{n-m} n_0}{\mu_0}}, \quad (17) $$

where we have defined the reference porosity as:

$$ \phi_0 = \left| \frac{\psi_0 \mu_0}{\Delta p g K_0} \right|^{\frac{1}{n}}, \quad (18) $$

with $\psi_0$ a reference flux used for the non-dimensionalization of the boundary condition at the base of the mantle wedge following:

$$ Q'(z') = \frac{\omega Q_{n0}}{\psi_0}. \quad (19) $$

Note that the reference flux and density contrast are also used to define the reference fluid velocity:

$$ v_0 = \left| \frac{\mu_0 \Delta p g}{K_0} \right|^{\frac{1}{n}} \frac{n-1}{\psi_0^n}. \quad (20) $$

The reference permeability ($K_0 = 3.7 \times 10^{-5} m^2$), fluid viscosity ($\mu_0 = 1$ Pa s), and solid viscosity ($n_0 = 10^{19}$ Pa s) are fixed in all models.

## REFERENCES CITED

Angiboust, S., Wolf, S., Burov, E., Agard, P., and Yamato, P., 2012, Effect of fluid circulation on subduction interface tectonic processes: Insights from thermo-mechanical numerical modelling: Earth and Planetary Science Letters, v. 357, p. 238–248, https://doi.org/10.1016/j.epsl.2012.09.012.

Arcay, D., Tric, E., and Doin, M.P., 2005, Numerical simulations of subduction zones: Effect of slab dehydration on the mantle wedge dynamics: Physics of the Earth and Planetary Interiors, v. 149, no. 1–2, p. 133–153, https://doi.org/10.1016/j.pepi.2004.08.020.

Atkinson, H.V., 1988, Overview no. 65: Theories of normal grain growth in pure single phase systems: Acta Metallurgica, v. 36, no. 3, p. 469–491, https://doi.org/10.1016/0001-6160(88)90079-X.

Audétat, A., and Keppler, H., 2004, Viscosity of fluids in subduction zones: Science, v. 303, no. 5657, p. 513–516, https://doi.org/10.1126/science.1092282.

Austin, N.J., and Evans, B., 2007, Paleowattmeters: A scaling relation for dynamically recrystallized grain size: Geology, v. 35, no. 4, p. 343–346, https://doi.org/10.1130/G23244A.1.

Behn, M.D., Hirth, G., and Elsenbeck, J.R., 2009, Implications of grain size evolution on the seismic structure of the oceanic upper mantle: Earth and Planetary Science Letters, v. 282, no. 1, p. 178–189, https://doi.org/10.1016/j.epsl.2009.03.014.

Bulau, J.R., Wa, H.S., and Tyburczy, J.A., 1979, Mechanical and thermodynamic constraints on fluid distribution in partial melts: Journal of Geophysical Research. Solid Earth, v. 84, B11, p. 6102–6108, https://doi.org/10.1029/JB084IB11p06102.

Cagnioncle, A.-M., Parmentier, E.M., and Elkins-Tanton, L.T., 2007, Effect of solid flow above a subducting slab on water distribution and melting at convergent plate boundaries: Journal of Geophysical Research. Solid Earth, v. 112, no. B9, https://doi.org/10.1029/2007JB004934.

Cerpa, N.G., Wada, I., and Wilson, C.R., 2017, Fluid migration in the mantle wedge: Influence of mineral grain size and mantle compaction: Journal of Geophysical Research. Solid Earth, https://doi.org/10.1002/2017JB014046.

Connolly, J.A.D., 2009, The geodynamic equation of state: What and how: Geochemistry, Geophysics, Geosystems, v. 10, no. 10, https://doi.org/10.1029/2009GC002540.

Connolly, J.A.D., Schmidt, M.W., Solferino, G., and Bagdassarov, N., 2009, Permeability of asthenospheric mantle and melt extraction rates at mid-ocean ridges: Nature, v. 462, no. 7270, p. 209–212, https://doi.org/10.1038/nature08517.

Dymkova, D., and Gerya, T., 2013, Porous fluid flow enables oceanic subduction initiation on Earth: Geophysical Research Letters, v. 40, no. 21, p. 5671–5676, https://doi.org/10.1002/2013GL057798.

England, P., Engdahl, R., and Thatcher, W., 2004, Systematic variation in the depths of slabs beneath arc volcanoes: Geophysical Journal International, v. 156, no. 2, p. 377–408, https://doi.org/10.1111/j.1365-246X.2003.02132.x.

England, P.C., and Katz, R.F., 2010, Melting above the anhydrous solidus controls the location of volcanic arcs: Nature, v. 467, no. 7316, p. 700–703, https://doi.org/10.1038/nature09417.

Evans, B., Renner, J., and Hirth, G., 2001, A few remarks on the kinetics of static grain growth in rocks: International Journal of Earth Sciences, v. 90, no. 1, p. 88–103, https://doi.org/10.1007/s005310000150.

Fowler A. C., 1985, A mathematical model of magma transport in the asthenosphere: Geophysical & Astrophysical Fluid Dynamics, v. 33, no. 1–4, p. 63–96. https://doi.org/10.1080/03091928508245423.

Frank, F.C., 1968, Two-component flow model for convection in the Earth's upper mantle: Nature, v. 220, p. 350–352, https://doi.org/10.1038/220350a0.

Giordano, D., Russell, J.K., and Dingwell, D.B., 2008, Viscosity of magmatic liquids: A model: Earth and Planetary Science Letters, v. 271, no. 1, p. 123–134, https://doi.org/10.1016/j.epsl.2008.03.038.

Grove, T.L., Chatterjee, N., Parman, S.W., and Médard, E., 2006, The influence of H₂O on mantle wedge melting: Earth and Planetary Science Letters, v. 249, no. 1, p. 74–89, https://doi.org/10.1016/j.epsl.2006.06.043.

Grove, T.L., Till, C.B., Lev, E., Chatterjee, N., and Médard, E., 2009, Kinematic variables and water transport control the formation and location of arc volcanoes: Nature, v. 459, no. 7247, p. 694–697, https://doi.org/10.1038/nature08044.

Hack, A.C., and Thompson, A.B., 2011, Density and viscosity of hydrous magmas and related fluids and their role in subduction zone processes: Journal of Petrology, v. 52, no. 7–8, p. 1333, https://doi.org/10.1093/petrology/egg048.

Hacker, B.R., Abers, G.A., and Peacock, S.M., 2003, Subduction factory 1. Theoretical mineralogy, densities, seismic wave speeds, and H₂O contents: Journal of Geophysical Research. Solid Earth, v. 108, no. B1, https://doi.org/10.1029/2001JB001127.

Havlin, C., Parmentier, E.M., and Hirth, G., 2013, Dike propagation driven by melt accumulation at the lithosphere-asthenosphere boundary: Earth and Planetary Science Letters, v. 376, p. 20–28, https://doi.org/10.1016/j.epsl.2013.06.010.

Hebert, L.B., and Montési, L.G., 2013, Hydration adjacent to a deeply subducting slab: The roles of nominally anhydrous minerals and migrating fluids: Journal of Geophysical Research: Solid Earth, v. 118, https://doi.org/10.1002/2013JB010497.

Hier-Majumder, S., and Kohlstedt, D., 2006, Role of dynamic grain boundary wetting circulation beneath volcanic arcs: Geophysical Research Letters, v. 33, no. 8, https://doi.org/10.1029/2006GL025716.

Hirschmann, M.M., Asimow, P.D., Ghiorso, M.S., and Stolper, E.M., 1999, Calculation of peridotite partial melting from thermodynamic models of minerals and melts. III. Controls on isobaric melt production and the effect of water on melt production: Journal of Petrology, v. 40, no. 5, p. 831–851, https://doi.org/10.1093/petroj/40.5.831.

Hirth, G., and Kohlstedt, D., 2003, Rheology of the upper mantle and the mantle wedge: A view from the experimentalists, in Eiler, J., ed., Inside the Subduction Factory: American Geophysical Union, Geophysical Monograph Series, v. 138, p. 83–105, https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/138GM06.

Horiuchi, S.S., and Iwamori, H., 2016, A consistent model for fluid distribution, viscosity distribution, and flow-thermal structure in subduction zone: Journal of Geophysical Research. Solid Earth, v. 121, no. 5, p. 3238–3260, https://doi.org/10.1002/2015JB012384.

Iwamori, H., 1998, Transportation of H₂O and melting in subduction zones: Earth and Planetary Science Letters, v. 160, no. 1–2, p. 65–80, https://doi.org/10.1016/S0012-821X(98)00080-6.

Iwamori, H., 2000, Deep subduction of H₂O and deflection of volcanic chain towards backarc near triple junction due to lower temperature: Earth and Planetary Science Letters, v. 181, p. 41–46, https://doi.org/10.1016/S0012-821X(00)00180-1.

Karato, S.-I., 1984, Grain-size distribution and rheology of the upper mantle: Tectonophysics, v. 104, no. 1, p. 155–176, https://doi.org/10.1016/0040-1951(84)90108-2.

Katz, R.F., Spiegelman, M., and Langmuir, C.H., 2003, A new parameterization of hydrous mantle melting: Geochemistry, Geophysics, Geosystems, v. 4, no. 9, https://doi.org/10.1029/2002GC000433.
Kawakatsu, H., and Watada, S., 2007, Seismic evidence for deep-water transportation in the mantle: Science, v. 316, no. 5830, p. 1468–1471, https://doi.org/10.1126/science.1140855.

Kelemen, P.B., Rilling, J.L., Parmentier, E.M., Mehl, L., and Hacker, B.R., 2003, Thermal structure due to solid-state flow in the mantle wedge beneath arcs, in Eiler, J., ed., Inside the Subduction Factory: American Geophysical Union, Geophysical Monograph Series, v. 138, p. 293–311, https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/138GM13.

Keller, T., May, D.A., and Kaus, B.J.P., 2013, Numerical modelling of magma dynamics coupled to tectonic deformation of lithosphere and crust: Geophysical Journal International, v. 195, no. 3, p. 1406–1442, https://doi.org/10.1093/gji/ggt306.

Kelley, K.A., Plank, T., Newman, S., Stolper, E., Grove, T., Parman, S., and Hauri, E., 2010, Mantle melting as a function of water content beneath the Mariana arc: Journal of Petrology, v. 51, no. 8, p. 1711, https://doi.org/10.1093/petrology/egg036.

Kirby, S.H., Durham, W.B., and Stern, L.A., 1991, Mantle phase changes and deep-earthquake faulting in subducting lithosphere: Science, v. 252, no. 5003, p. 216, https://doi.org/10.1126/science.252.5003.216.

Leeman, W.P., Lewis, J.F., Evarts, R.C., Conrey, R.M., and Streck, M.J., 2005, Petrologic constraints on the thermal structure of the cascades arc: Journal of Volcanology and Geothermal Research, v. 140, no. 1, p. 67–105, https://doi.org/10.1016/j.jvolgeores.2004.07.016.

Maaløe, S., and Scheie, Å., 1982, The permeability controlled accumulation of primary magma: Contributions to Mineralogy and Petrology, v. 81, no. 4, p. 350–357, https://doi.org/10.1007/BF00371689.

Manning, C.E., 2004, The chemistry of subduction-zone fluids: Earth and Planetary Science Letters, v. 223, no. 1, p. 1–16, https://doi.org/10.1016/j.epsl.2004.04.030.

McCrorry, P.A., Blair, J.L., Waldhauser, F., and Oppenheimer, D.H., 2012, Juan de Fuca slab geometry and its relation to Wadati-Benioff zone seismicity: Journal of Geophysical Research. Solid Earth, v. 117, no. B9, https://doi.org/10.1029/2012JB009407.

McGary, R.S., Evans, R.L., Wannamaker, P.E., Elsenbeck, J., and Rondenay, S., 2014, Pathway from subducting slab to surface for melt and fluids beneath Mount Rainier: Nature, v. 511, no. 7509, p. 338, https://doi.org/10.1038/nature13493.

McKenzie, D., 1984, The generation and compaction of partially molten rock: Journal of Petrology, v. 25, no. 3, p. 713, https://doi.org/10.1093/petrology/25.3.713.

Mibe, K., Fujii, T., and Yasuda, A., 1998, Connectivity of aqueous fluid in the Earth's upper mantle: Geophysical Research Letters, v. 25, no. 8, p. 1233–1236, https://doi.org/10.1029/98GL00872.

Mibe, K., Fujii, T., and Yasuda, A., 1999, Control of the location of the volcanic front in island arcs by aqueous fluid connectivity in the mantle wedge: Nature, v. 401, no. 6750, p. 259–262, https://doi.org/10.1038/45762.

Miller, K.J., Zhu, W.-I., Mont'esi, L.G.J., and Gaetani, G.A., 2014, Experimental quantification of permeability of partially molten mantle rock: Earth and Planetary Science Letters, v. 388, p. 273–282, https://doi.org/10.1016/j.epsl.2013.12.003.

Newton, R.C., and Manning, C.E., 2002, Solubility of enstatite + forsterite in H₂O at deep crust/upper mantle conditions: 4 to 15 kbar and 700 to 900°C: Geochemica et Cosmochimica Acta, v. 66, no. 23, p. 4165–4176, https://doi.org/10.1016/S0016-7037(02)00998-5.

Ni, H., Keppler, H., and Behrens, H., 2011, Electrical conductivity of hydrous basaltic melts: Implications for partial melting in the upper mantle: Contributions to Mineralogy and Petrology, v. 162, no. 3, p. 637–650, https://doi.org/10.1007/s00410-011-0617-4.

Okamoto, A., and Shimizu, H., 2015, Contrasting fracture patterns induced by volume-increasing and-decreasing reactions: Implications for the progress of metamorphic reactions: Earth and Planetary Science Letters, v. 417, p. 9–18, https://doi.org/10.1016/j.epsl.2015.02.015.

Peacock, S.M., and Wang, K., 1999, Seismic consequences of warm versus cool subduction metamorphism: Examples from southwest and northeast Japan: Science, v. 286, no. 5441, p. 937–939, https://doi.org/10.1126/science.286.5441.937.

Plank, T., Kelley, K.A., Zimmer, M.M., Hauri, E.H., and Wallace, P.J., 2013, Why do mafic arc magmas contain ~4 wt% water on average?: Earth and Planetary Science Letters, v. 364, p. 168–179, https://doi.org/10.1016/j.epsl.2012.11.044.

Pommier, A., and Garnero, E.J., 2014, Petrology-based modeling of mantle melt electrical conductivity and joint interpretation of electromagnetic and seismic results: Journal of Geophysical Research. Solid Earth, v. 119, no. 5, p. 4001–4016, https://doi.org/10.1002/2013JB010449.

Portnyagin, M., Hoernle, K., Plechov, P., Mironov, N., and Khubunaya, S., 2007, Constraints on mantle melting and composition and nature of slab components in volcanic arcs from volatiles (H₂O, S, Cl, F) and trace elements in melt inclusions from the Kamchatka Arc: Earth and Planetary Science Letters, v. 255, no. 1, p. 53–69, https://doi.org/10.1016/j.epsl.2006.12.005.

Précigout, J., Prigent, C., Palasse, L. and Pochon, A., 2017, Water pumping in mantle shear zones: Nature Communications, v. 8, article no. 15736, https://doi.org/10.1038/ncomms15736.

Rees Jones, D.W., Katz, R.F., Tian, M., and Rudge, J.F., 2018, Thermal impact of magmatism in subduction zones: Earth and Planetary Science Letters, v. 481, p. 73–79, https://doi.org/10.1016/j.epsl.2017.10.015.

Rudge, J.F., 2018, Textural equilibrium melt geometries around tetrakaidecahedral grains: Proceedings of the Royal Society A, v. 474, https://doi.org/10.1098/rspa.2017.0639.

Schmidt, M.E., Grunder, A.L., and Rowe, M.C., 2008, Segmentation of the cascade arc as indicated by Sr and Nd isotopic variation among diverse primitive basalts: Earth and Planetary Science Letters, v. 266, no. 1, p. 166–181, https://doi.org/10.1016/j.epsl.2007.11.013.

Schmidt, M.W., and Poli, S., 1998, Experimentally based water budgets for dehydrating slabs and consequences for arc magma generation: Earth and Planetary Science Letters, v. 163, no. 1, p. 361–379, https://doi.org/10.1016/S0012-821X(98)00142-3.

Scott, D.R., and Stevenson, D.J., 1984, Magma solitons: Geophysical Research Letters, v. 11, no. 11, p. 1161–1164, https://doi.org/10.1029/20GL011i011p0161.

Simpson, G., Spiegelman, M., and Weinstein, M.I., 2010, A multiscale model of partial melts: 2. Numerical results: Journal of Geophysical Research. Solid Earth, v. 115, B4, https://doi.org/10.1029/2009JB006376.

Sisson, T.W., and Grove, T.L., 1993, Experimental investigations of the role of H₂O in calc-alkaline differentiation and subduction zone magmatism: Contributions to Mineralogy and Petrology, v. 113, no. 2, p. 143–166, https://doi.org/10.1007/BF00283225.

Sleep, N.H., 1974, Segregation of magma from a mostly crystalline mush: Geological Society of America Bulletin, v. 85, no. 8, p. 1225–1232, https://doi.org/10.1130/0016-7606(1974)85<1225:SOMFAM>2.0:CO2.

Spiegelman, M., 1993, Flow in deformable porous media. Part 1 Simple analysis: Journal of Fluid Mechanics, v. 247, p. 17–38, https://doi.org/10.1017/S00221102903000369.

Spiegelman, M., and McKenzie, D., 1987, Simple 2-D models for melt extraction at mid-ocean ridges and island arcs: Earth and Planetary Science Letters, v. 83, no. 1, p. 137–152, https://doi.org/10.1016/0012-821X(87)90057-4.

Stein, C.A., and Stein, S., 1992, A model for the global variation in oceanic depth and heat flow with lithospheric age: Nature, v. 359, no. 6391, p. 123–129, https://doi.org/10.1038/359123a0.

Stolper, E., and Newman, S., 1994, The role of water in the petrogenesis of Mariana trough magmas: Earth and Planetary Science Letters, v. 121, no. 3–4, p. 293–325, https://doi.org/10.1016/0012-821X(94)90074-4.

Syracuse, E.M., and Abers, G.A., 2006, Global compilation of variations in slab depth beneath arc volcanoes and implications: Geochemistry, Geophysics, Geosystems, v. 7, no. 5, https://doi.org/10.1029/2005GC001045.

Syracuse, E.M., van Keken, P.E., and Abers, G.A., 2010, The global range of subduction zone thermal models: Physics of the Earth and Planetary Interiors, v. 183, no. 1, p. 73–90, https://doi.org/10.1016/j.pepi.2010.02.004.

Tatsumi, Y., Sakuyama, M., Fukuyama, H., and Kushiro, I., 1983, Generation of arc basalt magmas and thermal structure of the mantle wedge in subduction zones: Journal of Geophysical Research. Solid Earth, v. 88, no. B7, p. 5815–5825, https://doi.org/10.1029/JB088iB07p05815.

Tatsumi, Y., Hamilton, D.L., and Nesbitt, R.W., 1986, Chemical characteristics of fluid phase released from a subducted lithosphere and origin of arc magmas: Evidence from high-pressure experiments and natural rocks: Journal of Volcanology and Geothermal Research, v. 29, no. 1–4, p. 293–309, https://doi.org/10.1016/0377-0273(86)90049-1.

Tatsumi, Y., Furukawa, Y., Kogiso, T., Yamanaka, Y., Yokoyama, T., and Fedotov, S.A., 1994, A third volcanic chain in Kamchatka: Thermal anomaly at transform/convergence plate boundary: Geophysical Research Letters, v. 21, no. 7, p. 537–540, https://doi.org/10.1029/94GL00214.

Till, C.B., Grove, T.L., and Withers, A.C., 2012, The beginnings of hydrous mantle wedge melting: Contributions to Mineralogy and Petrology, v. 163, no. 4, p. 669–688, https://doi.org/10.1007/s00410-011-0692-6.

Turcotte, D.L., and Schubert, G., 2014, Geodynamics: Cambridge, UK, Cambridge University Press, https://doi.org/10.1017/CBO9780511843877.

Turner, A.J., Katz, R.F., Behn, M.D., and Keller, T., 2017, Magmatic focusing to mid-ocean ridges: The role of grain-size variability and non-Newtonian viscosity: Geochemistry, Geophysics, Geosystems, v. 18, no. 12, p. 4342–4355, https://doi.org/10.1002/2017GC007048.

van Keken, P.E., Currie, C., King, S.D., Behn, M.D., Cagnioncle, A., He, J., Katz, R.F., Lin, S.-C., Parmentier, E.M., Spiegelman, M., and Kelin, W., 2008, A community benchmark for subduction zone modeling: Physics of the Earth and Planetary Interiors, v. 171, no. 1, p. 187–197, https://doi.org/10.1016/j.pepi.2008.04.015.
van Keken, P.E., Hacker, B.R., Syracuse, E.M., and Abers, G.A., 2011, Subduction factory: 4. Depth-dependent flux of H₂O from subducting slabs worldwide: Journal of Geophysical Research. Solid Earth, v. 116, no. B1, https://doi.org/10.1029/2010JB007922.

von Bargen, N., and Waff, H.S., 1986, Permeabilities, interfacial areas and curvatures of partially molten systems: Results of numerical computations of equilibrium microstructures: Journal of Geophysical Research. Solid Earth, v. 91, no. B9, p. 9261-9276, https://doi.org/10.1029/JB091iB09p09261.

Wada, I., and Behn, M.D., 2015, Focusing of upward fluid migration beneath volcanic arcs: Effect of mineral grain size variation in the mantle wedge: Geochemistry, Geophysics, Geosystems, v. 16, no. 11, p. 3905-3923, https://doi.org/10.1002/2015GC005950.

Wada, I., and Wang, K., 2009, Common depth of slab-mantle decoupling: Reconciling diversity and uniformity of subduction zones: Geochemistry, Geophysics, Geosystems, v. 10, no. 10, https://doi.org/10.1029/2009GC002570.

Wada, I., Behn, M.D., and He, J., 2011, Grain-size distribution in the mantle wedge of subduction zones: Journal of Geophysical Research. Solid Earth, v. 116, B10, https://doi.org/10.1029/2011JB008294.

Wada, I., Behn, M.D., and Shaw, A.M., 2012, Effects of heterogeneous hydration in the incoming plate, slab rehydration, and mantle wedge hydration on slab-derived H₂O flux in subduction zones: Earth and Planetary Science Letters, v. 353, p. 60-71, https://doi.org/10.1016/j.epsl.2012.07.025.

Waff, H.S., and Faul, U.H., 1992, Effects of crystalline anisotropy on fluid distribution in ultramafic partial melts: Journal of Geophysical Research. Solid Earth, v. 97, no. B6, p. 9003-9014, https://doi.org/10.1029/92JB00066.

Wallace, P., 2005, Volatiles in subduction zone magmas: Concentrations and fluxes based on melt inclusion and volcanic gas data: Journal of Volcanology and Geothermal Research, v. 140, no. 1, p. 217-240, https://doi.org/10.1016/j.jvolgeores.2004.07.023.

Walowski, K.J., Wallace, P.J., Hauri, E.H., Wada, I., and Clynne, M.A., 2015, Slab melting beneath the Cascade Arc driven by dehydration of altered oceanic peridotite: Nature Geoscience, v. 8, no. 5, p. 404-408, https://doi.org/10.1038/ngeo2417.

Wannamaker, P.E., Evans, R.L., Bedrosian, P.A., Unsworth, M.J., Maris, V., and McGary, R.S., 2014, Segmentation of plate coupling, fate of subduction fluids, and modes of arc magmatism in Cascadia, inferred from magnetotelluric resistivity: Geochemistry, Geophysics, Geosystems, v. 15, no. 11, p. 4230-4253, https://doi.org/10.1002/2014GC005509.

Wark, D.A., Williams, C.A., Watson, E.B., and Price, J.D., 2003, Reassessment of pore shapes in microstructurally equilibrated rocks, with implications for permeability of the upper mantle: Journal of Geophysical Research. Solid Earth, v. 108, B1, 10.1029/2001JB001575.

Watson, E.B., 1982, Melt infiltration and magma evolution: Geology, v. 10, no. 5, p. 236-240, https://doi.org/10.1130/0091-7613(1982)10<236:MIAME>2.0.CO;2.

Watson, E.B., and Brenan, J.M., 1987, Fluids in the lithosphere, 1. Experimentally-determined wetting characteristics of CO₂-H₂O fluids and their implications for fluid transport, host-rock physical properties, and fluid inclusion formation: Earth and Planetary Science Letters, v. 85, no. 4, p. 497-515, https://doi.org/10.1016/0012-821X(87)90144-0.

Watt, S.F.L., Pyle, D.M., Mather, T.A., and Naranjo, J.A., 2013, Arc magma compositions controlled by linked thermal and chemical gradients above the subducting slab: Geophysical Research Letters, v. 40, no. 11, p. 2550-2556, https://doi.org/10.1002/grl.50513.

Wilson, C.R., Spiegelman, M., van Keken, P.E., and Hacker, B.R., 2014, Fluid flow in subduction zones: The role of solid rheology and compaction pressure: Earth and Planetary Science Letters, v. 401, p. 261-274, https://doi.org/10.1016/j.epsl.2014.05.052.

Wilson, C.R., Spiegelman, M., and van Keken, P.E., 2017, Terraferma: The transparent finite element rapid model assembler for multiphysics problems in earth sciences: Geochemistry, Geophysics, Geosystems, https://doi.org/10.1002/2016GC006702.