Calculations with melt transport

2.13 Calculations with melt transport

The original formulation of the equations in Section 2.1 describes the movement of solid mantle material. These computations also allow for taking into account how partially molten material changes the material properties and the energy balance through the release of latent heat. However, this will not consider melt extraction or any relative movement between melt and solid and there might be problems where the transport of melt is of interest. Thus, ASPECT allows for solving additional equations describing the behavior of silicate melt percolating through and interacting with a viscously deforming host rock. This requires the advection of a compositional field representing the volume fraction of melt present at any given time (the porosity $ϕ$ ), and also a change of the mechanical part of the system. The latter is implemented using the approach of [KMK13] and changes the Stokes system to

\begin{align} - \nabla \cdot [2 η (ε (u_{s}) - \frac{1}{3} (\nabla \cdot u_{s}) 1)] + \nabla p_{f} + \nabla p_{c} & = ρ g & in Ω, & (40) \\ \nabla \cdot u_{s} - \nabla \cdot K_{D} \nabla p_{f} - K_{D} \nabla p_{f} \cdot \frac{\nabla ρ_{f}}{ρ_{f}} & = - \nabla \cdot K_{D} ρ_{f} g \\ + Γ (\frac{1}{ρ_{f}} - \frac{1}{ρ_{s}}) & (41) \\ - \frac{ϕ}{ρ_{f}} u_{s} \cdot \nabla ρ_{f} - \frac{1 - ϕ}{ρ_{s}} u_{s} \cdot \nabla ρ_{s} \\ - K_{D} g \cdot \nabla ρ_{f} & in Ω, \\ \nabla \cdot u_{s} + \frac{p_{c}}{ξ} & = 0 . & (42) \end{align}

We use the indices $s$ to indicate properties of the solid and $f$ for the properties of the fluid. The equations are solved for the solid velocity $u_{s}$ , the fluid pressure $p_{f}$ , and an additional variable, the compaction pressure $p_{c}$ , which is related to the fluid and solid pressure through the relation $p_{c} = (1 - ϕ) (p_{s} - p_{f})$ . $K_{D}$ is the Darcy coefficient, which is defined as the quotient of the permeability and the fluid viscosity and $Γ$ is the melting rate. $η$ and $ξ$ are the shear and compaction viscosities and can depend on the porosity, temperature, pressure, strain rate and composition. However, there are various laws for these quantities and so they are implemented in the material model. Common formulations for the dependence on porosity are $η = (1 - ϕ) η_{0} e^{- α_{ϕ} ϕ}$ with $α_{ϕ} \approx 25 . . . 30$ and $ξ = η_{0} ϕ^{- n}$ with $n \approx 1$ .

To avoid the density gradients in Equation (41), which would have to be specified individually for each material model by the user, we can use the same method as for the mass conservation (described in Section 2.10.3) and assume the change in solid density is dominated by the change in static pressure, which can be written as $\nabla p_{s} \approx \nabla p_{static} \approx ρ_{s} g$ . This finally allows us to write

\frac{1}{ρ_{s}} \nabla ρ_{s} \approx \frac{1}{ρ_{s}} \frac{\partial ρ_{s}}{\partial p_{s}} \nabla p_{s} \approx \frac{1}{ρ_{s}} \frac{\partial ρ_{s}}{\partial p_{s}} \nabla p_{s} \approx \frac{1}{ρ_{s}} \frac{\partial ρ_{s}}{\partial p_{s}} ρ_{s} g \approx κ_{s} ρ_{s} g .

For the fluid pressure, choosing a good approximation depends on the model parameters and setup (see [DH16]). Hence, we make $\nabla ρ_{f}$ a model input parameter, which can be adapted based on the forces that are expected to be dominant in the model. We can then replace the second equation by

\begin{array}{l} \nabla \cdot u_{s} - \nabla \cdot K_{D} \nabla p_{f} - K_{D} \nabla p_{f} \cdot \frac{\nabla ρ_{f}}{ρ_{f}} & = - \nabla \cdot (K_{D} ρ_{f} g) \\ + Γ (\frac{1}{ρ_{f}} - \frac{1}{ρ_{s}}) \\ - \frac{ϕ}{ρ_{f}} u_{s} \cdot \nabla ρ_{f} - (u_{s} \cdot g) (1 - ϕ) κ_{s} ρ_{s} \\ - K_{D} g \cdot \nabla ρ_{f} . \end{array}

The melt velocity is computed as

u_{f} = u_{s} - \frac{K_{D}}{ϕ} (\nabla p_{f} - ρ_{f} g),

but is only used for postprocessing purposes and for computing the time step length.

Note: Here, we do not use the visco-elasto-plastic rheology of the [KMK13] formulation. Hence, we do not consider the elastic deformation terms that would appear on the right hand side of Equation (40) and Equation (42) and that include the elastic and compaction stress evolution parameters

ξ_{τ}

and

ξ_{p}

. Moreover, our viscosity parameters

η

and

ξ

only cover viscous deformation instead of combining visco-elasticity and plastic failure. This would require a modification of the rheologic law using effective shear and compaction viscosities

η_{eff}

and

ξ_{eff}

combining a failure criterion and shear and compaction visco-elasticities.

Moreover, melt transport requires an advection equation for the porosity field $ϕ$ :

\begin{align} ρ_{s} \frac{\partial (1 - ϕ)}{\partial t} + \nabla \cdot [ρ_{s} (1 - ϕ) u_{s}] & = - Γ & in Ω, i = 1 \dots C & (43) \end{align}

In order to solve this equation in the same way as the other advection equations, we replace the second term of the equation by:

\nabla \cdot [ρ_{s} (1 - ϕ) u_{s}] = (1 - ϕ) (ρ_{s} \nabla \cdot u_{s} + \nabla ρ_{s} \cdot u_{s}) - \nabla ϕ \cdot ρ_{s} u_{s}

Then we use the same method as described above and assume again that the change in density is dominated by the change in static pressure

\frac{1}{ρ_{s}} \nabla ρ_{s} \cdot u_{s} \approx κ_{s} ρ_{s} g \cdot u_{s}

so we get

\frac{\partial ϕ}{\partial t} + u_{s} \cdot \nabla ϕ = \frac{Γ}{ρ_{s}} + (1 - ϕ) (\nabla \cdot u_{s} + κ_{s} ρ_{s} g \cdot u_{s}) .

More details on the implementation can be found in [DH16]. A benchmark case demonstrating the propagation of solitary waves can be found in Section 5.4.14.

[next] [prev] [prev-tail] [front] [up]