Jump to Search Jump to Navigation Jump to Content

Numerical approach

Numerical approach

For most of the simulations shown here, the 2D code I2VIS developed primarily by Taras Gerya is used. It combines a conservative finite difference method with non-diffusive marker-in-cell techniques. In the models visco-plastic rheology of the rocks is used, details of which can be found in “Rheology of the Earth" by Ranalli, 1995. The governing conservation equations and the constitutive relationships between stress and strain-rate (needed in the creeping flow regime) are solved on a staggered grid in Eulerian configuration. Conservation of mass is described by the Lagrangian continuity equation for a compressible fluid (Gerya and Yuen 2007; Sizova et al. 2010; Gerya et al. 2010). Physical properties are transported on Lagrangian markers, which move according to the velocity field interpolated from the fixed grid. For interpolation between markers and nodes, a distance-dependent bilinear averaging scheme (Gerya and Yuen, 2003) is used according to which markers located closer to a node have higher statistical weights.

Mechanical properties of the rocks used in the models are presented in Table 1.


The rheologies used in the models are visco-plastic. The viscous creep viscosities of rocks are represented as a function of temperature, pressure and stress in terms of deformation invariants by experimentally determined flow laws (Table 1, Ranalli 1995):


is the second invariant of the strain rate tensor. AD (pre-exponential factor), E (activation energy), V (activation volume), n (creep exponent) are experimentally determined flow law parameters and R is the Gas (gas) constant. Plasticity is implemented using the following yield criteria, which limits the creep viscosity, altogether yielding an effective visco-plastic rheology:

The local plastic strength of a rock depends on the mean stress, Ptotal = P (dynamic pressure), the ambient brittle/plastic strength c, which is the strength at P = 0, and the effective internal friction angle, j, which is calculated from the friction angle of dry rocks, jdry and the pore fluid pressure factor λfluid. This factor is interpreted as:

The pore fluid pressure Pfluid reduces the yield strength σyield of fluid-containing porous or fractured rocks. The weakening effect of ascending melt is implemented in a similar manner. During a melt extraction episode, the yield strength σyield of rocks located in the column between the area of the melt extraction and emplacement is decreased. In our experiments both λfluid, λmelt are constant: λfluid = 0.1 and λmelt = 0.001(see discussion in Gerya and Meilick 2011, Vogt, Gerya et al. 2012). The values of c, jdry, λfluid, λmelt depend on the mineral composition. In our numerical experiments 1018 and 1026 Pa×s are the lower and upper cut-off values for viscosity of all types of rocks.

Petrological Model

Stable mineral assemblages, as a function of pressure and temperature are computed based on thermodynamic data and Gibbs free energy minimization using the software Perplex (Connolly 2005) . Examples of the procedure and computational strategy may be found in Connolly (2005) and Kerrick and Connolly (2001). All properties (i.e. water content, effective density, isobaric heat capacity, thermal expansion, adiabatic heat and latent heat) were subsequently saved on a grid with a resolution of 5K and 25 MPa. At every time step these precomputed properties are updated for all the Lagrangian markers.


The topographic evolution accounts for the effects of erosion and sedimentation. The crust/air interface (surface) evolves according to a transport equation, which is solved at each time-step on the Eulerian grid (Gorczyk, Willner et al. 2007)

∂ yes / ∂ t = vy − vx ∂ yes / ∂ x − vs + ve,

where yes is the vertical position of the surface as a function of the horizontal distance x; vy and vx are the vertical and horizontal components of the material velocity vector at the surface (y is positive downward, y = 0 at the top of the box); vs and ve are the sedimentation and erosion rates, respectively and are defined as follows (Gerya and Meilick 2011):

vs = 0 mm/a               ve = 0.3 mm/a for y < 9 km

vs = 0.03 mm/a                      ve = 0 mm/a for y >10 km

Hydration and water migration

Dehydration reactions and associated water release are computed based on the physicochemical conditions and the assumption of thermodynamic equilibrium (Gorczyk, Gerya et al. 2007, Nikolaeva, Gerya et al. 2008, Sizova, Gerya et al. 2010, Gerya and Meilick 2011). Expelled water is stored in a newly generated water marker that moves independently. The velocity of the fluid is computed accordingly to pressure gradients as:

vx(water) and vy(water) are the fluid velocities in x and y direction; vx and vy indicate the local velocity of the solid mantle; A is a water percolation constant, vpercolation = 10 cm/yr, a presumed standard water percolation velocity (Peacock 1990, Gorczyk, Willner et al. 2007); g=9.81 m/s is the gravitational acceleration; rmantle=3300 kg/m3;  and rfluid=1000 kg/m3 is the density of the mantle and fluid.

When a given water marker encounters a lithology capable of absorbing water by hydration or partial melting, the moving water marker is absorbed.

Partial melting and melt extraction (if melt extraction present)

Because the water transport model does not permit complete hydration of the peridotitic mantle, the mantle solidus is intermediate between the wet and dry peridotite solidi. To account for this behavior we assume that the degree of both hydrous and dry melting is a linear function of pressure and temperature (i.e.: Gerya and Yuen 2003).

For a given pressure and rock composition the volumetric degree of melting M0 is:

M0 = 0           when T < Tsolidus,

M0 = (T - Tsolidus) / (Tliquidus - Tsolidus) when Tsolidus< T < Tliquidus,

M0 = 1           when T > Tliquidus,

Tsolidus and Tliquidus are solidus temperature and dry liquidus temperature (wet and dry solidi are used for the hydrated and dry mantle, respectively) at a given pressure and rock composition (Table 2).

The effective density, reff, of partially molten rocks is computed according to the relation:

where rsolid and rmelt  are the densities of solid and molten rocks at given P (MPa) and T (K).

The effect of latent heating due to equilibrium melting/crystallization is included implicitly by increasing the effective heat capacity (Cpeff) and the thermal expansion (αeff) of melting/crystallizing rocks:

Hl is the latent heat of melting.

The effective viscosity of partially molten rocks (M > 0.1) is calculated according to (Pinkerton 1992, Bittner and Schmeling 1995):

Here η0 is an empirical parameter depending on rock composition where η0= 1013 Pa×s for molten mafic rocks, and η0=5x1014 Pa×s for molten felsic rocks. The viscosities of the rocks that contain minor melt fractions (M < 0.1) are computed according to experimentally determined flow laws (Table 1), where the effective creep viscosities are a function of strain rate, pressure, and temperature (see Rheology section 2.1. for details).

Although melt might accumulate and form large melt reservoirs, it is more likely that melt collects in channels or dykes and leaves the melting zone before reaching high melt fractions (Schmeling, Monz et al. 1999, Schmeling, Babeyko et al. 2008). Hence, melt exceeding a predefined melt threshold  of (i.e.: Schmeling, Babeyko et al. 2008) Mmax = 4 % is extracted and only a non-extractable amount of melt Mmin = 2 % remains at the source (Nikolaeva, Gerya et al. 2008, Sizova, Gerya et al. 2010, Gerya and Meilick 2011). Markers track the amount of extracted melt during the evolution of an experiment. The total amount of melt, M, for every marker takes into account the amount of previously extracted melt and is calculated as:

M = M0 − Σn Mext

Σn Mext is the total melt fraction extracted during the previous n extraction episodes. Rocks are considered refractory when the extracted melt fraction is larger than the standard one, when Σn Mext > M0. If the total amount of melt exceeds Mmax, the melt fraction Mext = M − Mmin is extracted and Σn  Mext is updated.

Extracted melt is transmitted instantaneously to emplacement areas, as melt migrates faster than rocks deform (Elliott, Plank et al. 1997, Hawkesworth, Turner et al. 1997). Although the average and range of intrusive:extrusive (I:E) volume ratios for different petrotectonic settings may differ (ranging from 1:1 to 16:1), a ratio of 4:1 can be viewed as common to most magmatic systems considering the uncertainities (White, Crisp et al. 2006). Accordingly, intrusive rocks, which comprise 80% of melt, are emplaced in areas of highest possible intrusion emplacement rate, i.e., highest possible local crustal divergence rate. divcrust is evaluated above the melt extraction region by locally computing the ratio of the effective melt overpressure and the effective viscosity of the crust:

divcrust = [Pymelt –g rmelt(ymelt-y)-Py ]/ηy                       

Pymelt is the pressure at the extraction level ymelt , Py is the pressure at the current level y, g is the gravitational acceleration in the y-direction [m/s2], rmelt is the melt density, and ηy is the effective local crustal viscosity. Extracted melts are thus emplaced at the level y where the computed possible crustal divergence rate divcrust is the highest. Both the effects of matrix compaction in the melt extraction area and crustal divergence in the melt emplacement area are taken into account in the compressible continuity equation (i.e.: Gerya and Yuen 2007), which provides correct coupling between local and global flow fields.

References for numerical method

Bittner, D. and H. Schmeling (1995). "Numerical modeling of melting processes and induced diapirism in the lower crust." Geophysical Journal International 123: 59-70.

Connolly, J. A. D. (2005). "Computation of phase equilibria by linear programming: a tool for geodynamic modeling and an application to subduction zone decarbonation." Earth and Planetary Science Letters 236: 524-541.

Elliott, T., T. Plank, A. Zindler, W. White and B. Bourdon (1997). "Element transport from slab to volcanic front at the Mariana arc." Journal of Geophysical Research-Solid Earth 102(B7): 14991-15019.

Gerya, T. V. and F. I. Meilick (2011). "Geodynamic regimes of subduction under an active margin: effects of rheological weakening by fluids and melts." Journal of Metamorphic Geology 29(1): 7-31.

Gerya, T. V. and D. A. Yuen (2003). "Characteristics-based marker-in-cell method with conservative finite-differences schemes for modeling geological flows with strongly variable transport properties." Physics of the Earth and Planetary Interiors 140(4): 293-318.

Gerya, T. V. and D. A. Yuen (2007). "Robust characteristics method for modelling multiphase visco-elasto-plastic thermo-mechanical problems." Physics of the Earth and Planetary Interiors 163(1-4): 83-105.

Gorczyk, W., T. V. Gerya, J. A. D. Connolly and D. A. Yuen (2007). "Growth and mixing dynamics of mantle wedge plumes." Geology 35(7): 587-590.

Gorczyk, W., A. P. Willner, T. V. Gerya, J. A. D. Connolly and J. P. Burg (2007). "Physical controls of magmatic productivity at Pacific-type convergent margins: Numerical modelling." Physics of the Earth and Planetary Interiors 163(1-4): 209-232.

Hawkesworth, C. J., S. P. Turner, F. McDermott, D. W. Peate and P. vanCalsteren (1997). "U-Th isotopes in arc magmas: Implications for element transfer from the subducted crust." Science 276(5312): 551-555.

Kerrick, D. M. and J. A. D. Connolly (2001). "Metamorphic devolatilization of subducted oceanic metabasalts: implications for seismicity, arc magmatism and volatile recycling." Earth and Planetary Science Letters 189(1-2): 19-29.

Nikolaeva, K., T. V. Gerya and J. A. D. Connolly (2008). "Numerical modelling of crustal growth in intraoceanic volcanic arcs." Physics of the Earth and Planetary Interiors 171(1-4): 336-356.

Peacock, S. M. (1990). "Fluid processes in subduction zones." Science 248: 329-337.

Pinkerton, H., Stevenson, R.J. (1992). "Methods of determining the rheological properties of magmas at subliquidus temperatures." Journal of Volcanology and Geothermal Research 53: 47–66.

Ranalli, G. (1995). Rheology of the Earth. London, Chapman & Hall.

Schmeling, H., A. Y. Babeyko, A. Enns, C. Faccenna, F. Funiciello, T. Gerya, G. J. Golabek, S. Grigull, B. J. P. Kaus, G. Morra, S. M. Schmalholz and J. van Hunen (2008). "A benchmark comparison of spontaneous subduction models-Towards a free surface." Physics of the Earth and Planetary Interiors 171(1-4): 198-223.

Schmeling, H., R. Monz and D. C. Rubie (1999). "The influence of olivine metastability on the dynamics of subduction." Earth and Planetary Science Letters 165(1): 55-66.

Sizova, E., T. Gerya, M. Brown and L. L. Perchuk (2010). "Subduction styles in the Precambrian: Insight from numerical experiments." Lithos 116(3-4): 209-229.

Vogt, K., T. V. Gerya and A. Castro (2012). "Crustal growth at active continental margins: Numerical modeling." Physics of the Earth and Planetary Interiors 192–193(0): 1-20.

White, S. M., J. A. Crisp and F. J. Spera (2006). "Long-term volumetric eruption rates and magma budgets." Geochemistry Geophysics Geosystems 7.