model setup
Our model, manifests a ocean-continent subduction zone cross-section (Fig. 1) from the bottom of upper mantle to lithospheric, within an area of 820 km in depth and 4000 km in length. The model is based on 2D thermal-mechanic code (I2ELVIS)18,42, processing by finite differences and marker-in-cell method. The 2D coupled thermal-mechanic model, simulates forced subduction of an oceanic plate beneath a continental plate. The model is divided into non-uniform rectangular grids with 2041×481 nodes, with more than 10 million randomly distributed markers. The subduction zone (from x = 1498 km to x = 3002 km) has a high resolution (1 km×1 km) both in horizontal and vertical, while the rest region adopts a relatively lower resolution, which horizontally varying from 1.0421 km to 5.0046 km according to the distance from the subduction zone. The initial set of materials is shown in Fig. 1, and detailed material properties are listed in Table.1 and Table.2. A wet low brittle/plastic strength rheological olivine weak shear zone, whose sin(φ) = 0.1 (where φ is the effective internal friction angle), is built between the overriding and subducting slab to prescribe the initialization of subduction. The accretionary sediments are settled between the oceanic and continent crust and overlying the weak shear zone (Fig. 1).
The bottom boundary of the model is vertically permeable, while the rest three boundaries are free slip surface in mechanics43,44. A boundary push (imposed on both subducting and overriding plates) at a pre-existing weak zone initiates the subduction, until the oceanic slab sinks into mantle approximately 100–150 km, when the subduction system can convert into a self-sustaining subduction. The initial thermal field in the subducting plate is set according to an oceanic geotherm45 for a 40 Ma lithospheric cooling age20. The initial temperature of the overriding continental plate linearly increases from 0 ℃ at the surface to 1367 ℃ at the lithosphere-asthenosphere boundary, below which an adiabatic temperature gradient of 0.5 ℃/km is used.
Table.1
Material parameters used in the model.
| Material | ρ0, kg/m3 | c0, MPa | sin(φdry) | QL, kJ/kg | Hr, µW/ m3 | Cp, J/kg·K | α, 1/K | β, 1/MPa | Flow law |
| Sediments, volcanics from sediments | 2600 (solid) 2400 (molten) | 10 | 0.15 | 300 | 2 | 1000 | 3×10− 5 | 1×10− 5 | Wet quartzite |
| Upper continental crust | 2700 (solid) 2400 (molten) | 10 | 0.15 | 300 | 1 | 1000 | 3×10− 5 | 1×10− 5 | Wet quartzite |
| Lower continental crust | 2700 (solid) 2400 (molten) | 10 | 0.15 | 380 | 1 | 1000 | 3×10− 5 | 1×10− 5 | Wet quartzite |
| Upper oceanic crust (basalt) | 3000 (solid) 2900 (molten) | 10 | 0.1 | 380 | 0.25 | 1000 | 3×10− 5 | 1×10− 5 | Wet quartzite |
| Lower oceanic crust (babbro) | 3000 (solid) 2900 (molten) | 10 | 0.6 | 380 | 0.25 | 1000 | 3×10− 5 | 1×10− 5 | Plagioclase An75 |
| Volcanics from wet molten and subducted basalts and gabbro | 3000 (solid) 2900 (molten) | 10 | 0.15 | 380 | 0.25 | 1000 | 3×10− 5 | 1×10− 5 | Wet quartzite |
| Lithoshpere-asthenosphere dry mantle | 3300 (solid) 2900 (molten) | 10 | 0.6 | 400 | 0.022 | 1000 | 3×10− 5 | 1×10− 5 | Dry olivine |
| Reference a | 1,2 | 3 | 3 | 1,2 | 1 | 1,2 | 1,2 | 1,2 | 3 |
a 1. DonaldL.Turcotte & JerrySchubert45; 2. Bittner. & Schmeling46; 3. Ranalli47.
Table.2
Temperature and pressure dependent parameters used in the model.
| Material | κ, W/m·K (at Tκ, P) | Tsolidus, K (at P) | Tliquidus, K (at P) |
| Sediments, volcanics from sediments | [0.64 + 807/(T + 77)]×exp(0.00004P) | 889 + 17900/(P + 54) + 20200/(P + 54)2 at P < 1200 MPa. 831 + 0.06P at P >1200 MPa | 1262 + 0.09P |
| Upper continental crust | [0.64 + 807/(T + 77)]×exp(0.00004P) | 889 + 17900/(P + 54) + 20200/(P + 54)2 at P < 1200 MPa. 831 + 0.06P at P >1200 MPa | 1262 + 0.09P |
| Lower continental crust | [0.64 + 807/(T + 77)]×exp(0.00004P) | 973-70400/(P + 354) + 77800000/(P + 354)2 at P < 1600 MPa. 935 + 0.0035P + 0.0000062P2 at P >1600 MPa | 1423 + 0.105P |
| Upper oceanic crust (basalt) | [1.18 + 474/(T + 77)]×exp(0.00004P) | 973-70400/(P + 354) + 77800000/(P + 354)2 at P < 1600 MPa. 935 + 0.0035P + 0.0000062P2 at P >1600 MPa | 1423 + 0.105P |
| Lower oceanic crust (gabbro) | [1.18 + 474/(T + 77)]×exp(0.00004P) | 973-70400/(P + 354) + 77800000/(P + 354)2 at P < 1600 MPa. 935 + 0.0035P + 0.0000062P2 at P >1600 MPa | 1423 + 0.105P |
| Volcanics from wet molten and subducted basalts and gabbro | [1.18 + 474/(T + 77)]×exp(0.00004P) | 973-70400/(P + 354) + 77800000/(P + 354)2 at P < 1600 MPa. 935 + 0.0035P + 0.0000062P2 at P >1600 MPa | 1423 + 0.105P |
| Lithoshpere-asthenosphere dry mantle | [0.73 + 1293/(T + 77)]×exp(0.00004P) | 273.15 + 1085.17 + 13.29P-0.051P2-43\({\text{X}}_{{\text{H}}_{2}\text{O}}^{0.75}\)b | 273.15 + 1475 + 8P-0.32 P2-43\({\text{X}}_{{\text{H}}_{2}\text{O}}^{0.75}\)b |
| Reference a | 1,2 | 3,4,5,6,7,8 | 3 |
a 1. Clauser. & Huenges48, 2. Hofmeister49, 3. Schmidt & Poli23, 4. Hess50, 5. Hirschmann51, 6. Johannes52, 7. Poli & Schmidt53, 8. Katz et al.54.
b \({\text{X}}_{{\text{H}}_{2}\text{O}}^{0.75}\) is the weight percent of water Katz et al.54.
hydration and dehydration process
Our model initially sets a homogeneously hydrated oceanic crust with a length of 7 km, which is composed of 2 km of hydrothermally altered hydrated basalts (4 wt.% H2O) and 5 km of gabbros (1.4 wt.% H2O). These water contents are used to replace fluids, which is mainly from the slab hydration process at the mid-ocean ridge. The main shallow hydration process in our model appears at the slab bending zone. Prior to subduction the oceanic plate begins to bend downward, and thus the near-surface rocks are placed in tension leading to generation of normal faults. The fluid deeply percolates into the bending oceanic plate along these normal faults, which favors pumping the fluid downward and further reactivation27 (Fig. 2). Another shallow hydration of oceanic plate is that seawater percolate into the porous and permeable basaltic crust. The maximum serpentinization depth at the trench is set down to sub-Moho (10–20 km). A finite brittle deformation threshold (e.g. ɛserp=0.05 in our experiments) is set to control the mantle serpentinization. When the deformation of lithospheric mantle exceeds the prescribed threshold, the material transforms to serpentinized rock (2 wt.% H2O) with 15% of serpentinization (Fig. 2). These initially maximum serpentinization depth dataset is based on the average value estimated from tomographic researches29,55,56 at the trench. At intermediate-deep depths, plastic behavior is not active. Serpentinization of the slab interior occurs intrinsically in our model controlling by a result of reaction between fluids in excess and the dry rocks through which they migrate.
Both porous compaction and dehydration reactions are attributed to expelled water from subducting oceanic plate. The free fluids in hydrous basalts and gabbros keep releasing from 4 wt.% H2O/1.4 wt.% H2O to zero at shallow (75 km)57:
$${\text{X}}_{{\text{H}}_{2}\text{O}}\left(\text{w}\text{t}.\text{\%}\right)=(1-0.013\bullet \varDelta \text{y})\bullet {\text{X}}_{{\text{H}}_{2}\text{O}}^{0}$$
,
where \({\text{X}}_{{\text{H}}_{2}\text{O}}^{0}\) is the weight percent of water at the surface, \(\varDelta \text{y}\) is the depth in km (0–75 km). We calculate the dehydration reactions/water release on the basis of physicochemical conditions and the assumption of thermodynamic equilibrium20,57,58. The hydration and dehydration processes are modeled in the form of marks. A dehydration curve is used to control the dehydration reaction, which is defined as23:
\(\text{T}=478+0.18\times \text{P}-0.000031\times {\text{P}}^{2}\) , at P<2.1 GPa
\(\text{T}=740-0.0018\times \text{P}-0.0000039\times {\text{P}}^{2}\) , at P>2.1 GPa
where T is the dehydration temperature in ℃, P is the current pressure in MPa. The released fluid is stored in a newly generated Lagrangian marker (named water marker) and propagates upward into the mantle wedge42,59. The migration velocity of water marker is calculated as44:
$${\text{v}}_{\text{x}}^{\text{w}\text{a}\text{t}\text{e}\text{r}}={\text{v}}_{\text{x}}^{\text{m}\text{a}\text{n}\text{t}\text{l}\text{e}}, {\text{v}}_{\text{y}}^{\text{w}\text{a}\text{t}\text{e}\text{r}}={\text{v}}_{\text{y}}^{\text{m}\text{a}\text{n}\text{t}\text{l}\text{e}}-{\text{v}}_{\text{y}}^{\text{p}\text{e}\text{r}\text{c}\text{o}\text{l}\text{a}\text{t}\text{i}\text{o}\text{n}},$$
where \({\text{v}}_{\text{x}}^{\text{m}\text{a}\text{n}\text{t}\text{l}\text{e}}\) and \({\text{v}}_{\text{y}}^{\text{m}\text{a}\text{n}\text{t}\text{l}\text{e}}\) are the local velocity of mantle in x and y direction, respectively. \({\text{v}}_{\text{y}}^{\text{p}\text{e}\text{r}\text{c}\text{o}\text{l}\text{a}\text{t}\text{i}\text{o}\text{n}}=10 \text{c}\text{m}/\text{a}\) describes the relative velocity of the upward percolation of the water. The water marker propagates upward independently, until it encounters mantle rocks, which can consumed the water by hydration reaction or partial melting. The seismic data60 implies that the majority mantle wedges contain average water contents less than ~ 4 wt.%, while saturated mantle rocks achieve ~ 8 wt.% H2O61. According to this fact, We assume an incomplete hydration mantle wedge of 2 wt.% H2O as a consequence of channelizaion of slab-derived fluids57,62,63. The stable hydrous phase in our experiments is calculated by P-T conditions of materials (see in Table. 2) and equilibrium water content (≤ 2 wt.% H2O).
Magma-related process
Like aforementioned hydration/dehydration process, our experiments use Lagrangian Markers to track the magma-related process during model evolution (i.e. partial melting, melt extraction, percolation and accretion). A simplified method64 is adopted to implement magma-related process in our 2D coupled thermal-mechanic models.
Katz54 presented a parameterization for melt fraction as a function of pressure, temperature, water content and modal cpx, based on recent thermodynamic and experimental investigations. This melt fraction is applied in our models as:
$${\text{M}}_{\text{c}\text{p}\text{x}}={\left(\frac{\text{T}-{\text{T}}_{\text{s}\text{o}\text{l}\text{i}\text{d}\text{u}\text{s}}}{{\text{T}}_{\text{l}\text{i}\text{q}\text{u}\text{i}\text{d}\text{u}\text{s}}^{\text{l}\text{h}\text{e}\text{r}\text{z}}-{\text{T}}_{\text{s}\text{o}\text{l}\text{i}\text{d}\text{u}\text{s}}}\right)}^{{{\beta }}_{\text{c}\text{p}\text{x}}} , \text{w}\text{h}\text{e}\text{n} {\text{T}}_{\text{s}\text{o}\text{l}\text{i}\text{d}\text{u}\text{s}}<\text{T}<{\text{T}}_{\text{c}\text{p}\text{x}-\text{o}\text{u}\text{t}}$$
$${\text{M}}_{\text{o}\text{p}\text{x}}={\text{M}}_{\text{c}\text{p}\text{x}-\text{o}\text{u}\text{t}}+\left(1-{\text{M}}_{\text{c}\text{p}\text{x}-\text{o}\text{u}\text{t}}\right){\left(\frac{\text{T}-{\text{T}}_{\text{c}\text{p}\text{x}-\text{o}\text{u}\text{t}}}{{\text{T}}_{\text{l}\text{i}\text{q}\text{u}\text{i}\text{d}\text{u}\text{s}}-{\text{T}}_{\text{c}\text{p}\text{x}-\text{o}\text{u}\text{t}}}\right)}^{{{\beta }}_{\text{o}\text{p}\text{x}}} , \text{w}\text{h}\text{e}\text{n} {\text{T}}_{\text{c}\text{p}\text{x}-\text{o}\text{u}\text{t}}<\text{T}<{\text{T}}_{\text{l}\text{i}\text{q}\text{u}\text{i}\text{d}\text{u}\text{s}}$$
where \({\text{M}}_{\text{c}\text{p}\text{x}}\) is the degree of melting prior to the exhaustion of cpx, \({\text{M}}_{\text{o}\text{p}\text{x}}\) is the degree of melting for\({\text{T}}_{\text{c}\text{p}\text{x}-\text{o}\text{u}\text{t}}<\text{T}<{\text{T}}_{\text{l}\text{i}\text{q}\text{u}\text{i}\text{d}\text{u}\text{s}}\), and the melting reaction consumes mostly opx after the exhaustion of cpx. T is the temperature in Kelvin, \({\text{T}}_{\text{s}\text{o}\text{l}\text{i}\text{d}\text{u}\text{s}}\) and \({\text{T}}_{\text{l}\text{i}\text{q}\text{u}\text{i}\text{d}\text{u}\text{s}}\) are the solidus and liquidus of mantle, respectively (shown in Table. 2). \({\text{T}}_{\text{l}\text{i}\text{q}\text{u}\text{i}\text{d}\text{u}\text{s}}^{\text{l}\text{h}\text{e}\text{r}\text{z}}\) (lherzolite liquidus) is applied to create a kinked melting function54. \({{\beta }}_{\text{c}\text{p}\text{x}}={{\beta }}_{\text{o}\text{p}\text{x}}=1.5\) are the equation exponents derived from "best fit" assemblage of the experimental solidus from Hirschmann51. \({\text{M}}_{\text{c}\text{p}\text{x}-\text{o}\text{u}\text{t}}\)=0.15/(0.5 + 0.08P) (P is the pressure in GPa) is the total cpx melting degree in a closed (batch) system. \({\text{T}}_{\text{c}\text{p}\text{x}-\text{o}\text{u}\text{t}}={\text{M}}_{\text{c}\text{p}\text{x}-\text{o}\text{u}\text{t}}^{{{\beta }}_{\text{c}\text{p}\text{x}}^{-1}}\left({\text{T}}_{\text{l}\text{i}\text{q}\text{u}\text{i}\text{d}\text{u}\text{s}}^{\text{l}\text{h}\text{e}\text{r}\text{z}}-{\text{T}}_{\text{s}\text{o}\text{l}\text{i}\text{d}\text{u}\text{s}}\right)+{\text{T}}_{\text{s}\text{o}\text{l}\text{i}\text{d}\text{u}\text{s}}\) is the temperature (in Kelvin) when cpx is exhausted through isobarically melting.
Partial melting can be observed in a wide area in our model, which is similar to the melt pooling region65. The melt is extracted and stored at the shallowest part of the partial melting zone, and this process is traced by Lagrangian markers in our models. The melt transportation is implemented by converting the markers (with rock type of molten mantle/asthenosphere) to a new marker type (molten basalt) in the shallowest part of the melting pool. The total amount of the extracted melt matches the volume of the converted markers. The total amount of melt for every marker is calculated at every modeling time step by the following equation:
$$\text{M}={\text{M}}_{0}-\sum _{\text{i}=1}^{\text{n}}{\text{M}}_{\text{i}}^{\text{e}\text{x}\text{t}}$$
,
where\({ \text{M}}_{0}\) is standard melt fraction (\({ \text{M}}_{0}={\text{M}}_{\text{c}\text{p}\text{x}}, \text{w}\text{h}\text{e}\text{n} {\text{T}}_{\text{s}\text{o}\text{l}\text{i}\text{d}\text{u}\text{s}}<\text{T}<{\text{T}}_{\text{c}\text{p}\text{x}-\text{o}\text{u}\text{t}}\), \({ \text{M}}_{0}={\text{M}}_{\text{o}\text{p}\text{x}}, \text{w}\text{h}\text{e}\text{n} {\text{T}}_{\text{c}\text{p}\text{x}-\text{o}\text{u}\text{t}}<\text{T}<{\text{T}}_{\text{l}\text{i}\text{q}\text{u}\text{i}\text{d}\text{u}\text{s}}\)), n indicates the time steps of previous extraction, \({\text{M}}_{\text{i}}^{\text{e}\text{x}\text{t}}\) represents the amount of extracted melt at i time step. The rock is considered refractory/non-molten, when the total extracted melt fraction (\(\sum _{\text{i}=1}^{\text{n}}{\text{M}}_{\text{i}}^{\text{e}\text{x}\text{t}})\) is larger than the standard one (\({\text{M}}_{0}\)). If the total amount of melt (\(\text{M}\)) exceeds a certain threshold, the melt in marker is extracted and \(\sum _{\text{i}=1}^{\text{n}}{\text{M}}_{\text{i}}^{\text{e}\text{x}\text{t}}\) is updated. The extracted melt is transmitted instantaneously to emplacement area (e.g. the shallowest part of melting pool and forms a magma chamber), since the extracted melt leaves the melting zone much faster than rock deform66,67. It is assumed that 80% of extracted melt emplace lower depth or underlying the continental plate, forming plutons in the continental crust in areas of highest possible intrusion emplacement. The remaining 20% can propagate upwards to the surface above the extraction zone, influencing the surface topography evolution (e.g. forming volcanoes). High possible local crustal divergence rate is used to predict the location of extracted melt intrusion emplacement, which is the ratio of the effective melt overpressure to the effective viscosity of the crust57:
$${\text{d}\text{i}\text{v}}_{\text{c}\text{r}\text{u}\text{s}\text{t}}=\left[{\text{P}}_{\text{m}\text{e}\text{l}\text{t}}-\text{g}{\rho }\left({\text{y}}_{\text{m}\text{e}\text{l}\text{t}}-\text{y}\right)-\text{P}\right]/{\eta }$$
,
where \({\text{P}}_{\text{m}\text{e}\text{l}\text{t}}\) and \(\text{P}\) are the pressure at the extraction depth (\({\text{y}}_{\text{m}\text{e}\text{l}\text{t}}\)) and at the current depth (\(\text{y}\)), respectively, \(\text{g}\) is the gravitational acceleration, \({\rho }\) is the melt density, and \({\eta }\) is the current effective local crustal viscosity. Extracted melts are emplaced at the depth where the computed possible local crustal divergence rate is highest. In order to obtain a correct coupling between local and global flow field, effective of matrix compaction in the melt extraction zone and crustal divergence in the melt emplacement area must be take into account. We deal these process with a compressible continuity Eq. 59.
Other effective values of physical parameters of materials during magma-related process are calculated by the equation modified based on the melt fraction. Latent heating effect due to melting/crystallization equilibrium is taken into account by increasing the effective heat capacity and thermal expansion in the energy conservation equation:
$${\text{X}}_{\text{e}\text{f}\text{f}}={\text{X}}_{\text{m}\text{o}\text{l}\text{t}\text{e}\text{n}}{\text{M}}_{\text{c}\text{p}\text{x}/\text{o}\text{p}\text{x}}+{\text{X}}_{\text{s}\text{o}\text{l}\text{i}\text{d}}(1-{\text{M}}_{\text{c}\text{p}\text{x}/\text{o}\text{p}\text{x}})$$
$${\text{C}}_{\text{p},\text{e}\text{f}\text{f}}={\text{C}}_{\text{p}}+{\text{Q}}_{\text{L}}{\left(\frac{\partial {\text{M}}_{\text{c}\text{p}\text{x}/\text{o}\text{p}\text{x}}}{\partial \text{T}}\right)}_{\text{P}}$$
$${{\alpha }}_{\text{e}\text{f}\text{f}}={\alpha }+{\rho }\frac{{\text{Q}}_{\text{L}}}{\text{T}}{\left(\frac{\partial {\text{M}}_{\text{c}\text{p}\text{x}/\text{o}\text{p}\text{x}}}{\partial \text{P}}\right)}_{\text{T}}$$
where \({\text{X}}_{\text{e}\text{f}\text{f}}\) represents the effective value of a physical parameter, \({\text{X}}_{\text{m}\text{o}\text{l}\text{t}\text{e}\text{n}}\) and \({\text{X}}_{\text{s}\text{o}\text{l}\text{i}\text{d}}\) are the physical parameter values of material in molten and solid phase, respectively. \({\text{Q}}_{\text{L}}\)is the latent heat, and \({\text{M}}_{\text{c}\text{p}\text{x}/\text{o}\text{p}\text{x}}\) is the melt fraction (\({\text{M}}_{\text{c}\text{p}\text{x}/\text{o}\text{p}\text{x}}={\text{M}}_{\text{c}\text{p}\text{x}}, \text{w}\text{h}\text{e}\text{n} {\text{T}}_{\text{s}\text{o}\text{l}\text{i}\text{d}\text{u}\text{s}}<\text{T}<{\text{T}}_{\text{c}\text{p}\text{x}-\text{o}\text{u}\text{t}}\) ,\({\text{M}}_{\text{c}\text{p}\text{x}/\text{o}\text{p}\text{x}}={\text{M}}_{\text{o}\text{p}\text{x}}, \text{w}\text{h}\text{e}\text{n} {\text{T}}_{\text{c}\text{p}\text{x}-\text{o}\text{u}\text{t}}<\text{T}<{\text{T}}_{\text{l}\text{i}\text{q}\text{u}\text{i}\text{d}\text{u}\text{s}}\)). All material parameters are shown in Tables. 1 and 2.