# Calculations with melt transport#

The original formulation of the equations in
Basic equations 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 \(\phi\)), and also a change of the
mechanical part of the system. The latter is implemented using the approach of
Keller *et al.* [2013] and changes the Stokes system to

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 \(\mathbf{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-\phi) (p_s-p_f)\). \(K_D\) is the Darcy coefficient, which is defined as the quotient of the permeability and the fluid viscosity and \(\Gamma\) is the melting rate. \(\eta\) and \(\xi\) 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 \(\eta = (1-\phi) \eta_0 e^{-\alpha_\phi \phi}\) with \(\alpha_\phi \approx 25...30\) and \(\xi = \eta_0 \phi^{-n}\) with \(n \approx 1\).

To avoid the density gradients in Equation (30), 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 The Boussinesq approximation (BA)) 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_{\text{static}} \approx \rho_s \textbf{g}\). This finally allows us to write

where \(\beta_s\) is the compressibility of the solid. In the paper that describes the implementation [Dannberg and Heister, 2016], \(\kappa\) is used for the compressibility. We change the variable here to be consistent throughout the manual.

For the fluid pressure, choosing a good approximation depends on the model parameters and setup (see Dannberg and Heister [2016]). Hence, we make \(\nabla \rho_{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

The melt velocity is computed as

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
Keller *et al.* [2013] formulation. Hence, we do not consider the elastic
deformation terms that would appear on the right hand side of Equations
(29) and (31) and that
include the elastic and compaction stress evolution parameters \(\xi_{\tau}\) and
\(\xi_p\). Moreover, our viscosity parameters \(\eta\) and \(\xi\) 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 \(\eta_{\text{eff}}\) and \(\xi_{\text{eff}}\) combining a failure
criterion and shear and compaction visco-elasticities.

Moreover, melt transport requires an advection equation for the porosity field \(\phi\):

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

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

so we get

More details on the implementation can be found in Dannberg and Heister [2016]. A benchmark case demonstrating the propagation of solitary waves can be found in The solitary wave benchmark.

# Calculations with Darcy flow#

To calculate fluid transport, ASPECT can advect compositional fields with the fluid velocity computed using Darcy’s Law (derived in McKenzie 1984). Currently, this method approximates the fluid pressure gradient as the lithostatic pressure, so the expression used for the fluid velocity is:

Where \(\mathbf{u}_s\) is the solid velocity, \(K_D\) is the Darcy Coefficient, \(\phi\) is the porosity, \(\mathbf{g}\) is the gravity vector, \(\rho_s\) is the solid density, and \(\rho_f\) is the fluid density. The implementation of this method shares some similarities with the melt implementation, but is much simpler and as a consequence has several limitations. Currently, in the absence of solid motion, the fluid phase will only be advected parallel to the gravity vector based on buoyancy forces. An additional limitation is that solid compaction with varying porosity is not calculated with the current implementation, and so mass is not conserved as fluid leaves the solid matrix at a given point. For small porosities, this does not pose too much of an issue, but this approximation would break down at high porosities. The advantage of this method is that it is relatively cheap, and it can be used in scenarios where porosities are small and fluid pressures can be approximated by \(\rho_s \mathbf{g}\).