the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Modelling chemical advection during magma ascent

### Nicolas Riel

### Pierre Lanari

Modelling magma transport requires robust numerical schemes for chemical advection. Current numerical schemes vary in their ability to be mass conservative, computationally efficient, and accurate. This study compares four of the most commonly used numerical schemes for advection: an upwind scheme, a weighted essentially non-oscillatory (WENO-5) scheme, a semi-Lagrangian (SL) scheme, and a marker-in-cell (MIC) method. The behaviour of these schemes is assessed using the passive advection of two different magmatic compositions. This is coupled in 2D with the temporal evolution of a melt anomaly that generates porosity waves. All algorithms, except the upwind scheme, are able to predict the melt composition with reasonable accuracy, but none of them is fully mass conservative. However, the WENO-5 scheme shows the best mass conservation. In terms of total running time and when multithreaded, the upwind, SL, and WENO-5 schemes show similar performance, while the MIC scheme is the slowest due to reseeding and removal of markers. The WENO-5 scheme has a reasonable total run time, has the best mass conservation, is easily parallelisable, and is therefore best suited for this problem.

- Article
(4955 KB) - Full-text XML
- BibTeX
- EndNote

Mechanisms of magma ascent and emplacement within the lithosphere and upper asthenosphere remain largely unconstrained (e.g. Connolly and Podladchikov, 2007b; Katz et al., 2022). Studies have attempted to address this problem using techniques ranging from geophysical measurements of the present-day lithosphere to geochemical analysis of the rock record. However, geophysical studies are hampered by indirect measurements, and natural samples in geochemical studies represent only the end-product of the melting processes (Brown, 2013; Clemens et al., 2022; Johnson et al., 2021). Comparatively, numerical modelling allows the investigation of these processes at a range of scales in space and time (e.g. Keller, 2013; Katz and Weatherley, 2012).

To numerically model such open systems, it is necessary to be able to describe the chemical and physical processes responsible for magma ascent in a rock. At low melt fractions and in the absence of externally applied stress, the physical processes are based on the continuum formulation of two-phase flow. It takes into account the concurrent mechanisms of rock matrix compaction and buoyancy of partial melt in an interconnected porous network (e.g. Scott and Stevenson, 1984; McKenzie, 1984). This formulation is based on mass and momentum conservation and an appropriate set of constitutive relationships. In addition, conservation of energy needs to be ensured to link mechanical to chemical processes (e.g. Katz, 2008). Chemical processes, such as phase reactions, can be considered using thermodynamics and/or kinetics and relate the equilibration of the melt with the hosting rock (e.g. Omlin et al., 2017; Bessat et al., 2022). They contribute to the transport dynamics by changing rock properties, such as density, viscosity, porosity, and permeability (Jha et al., 1994; Aharonov et al., 1995b; Keller and Katz, 2016). However, the amount of melt interacting with the rock is also modulated by transport mechanisms (Kelemen et al., 1997; Spiegelman and Kenyon, 1992; Aharonov et al., 1995a). Therefore, the development of integrated models that successfully describe the complex interaction between reaction and transport is key to understanding melting and melt extraction at all scales.

Numerous numerical studies have investigated reactive melt transport. It has been shown that melts that partially crystallise or dissolve the host rock could be a viable mechanism for channelling flow and creating heterogeneities in the mantle in the context of mid-ocean ridges (Aharonov et al., 1997; Spiegelman et al., 2001) and sub-arc magmatism (Bouilhol et al., 2011). Concerning lower-crust melting, this approach has mainly been used to understand the processes of chemical differentiation and the compositional range of magmas in mafic systems (e.g. Jackson et al., 2005; Solano et al., 2012; Riel et al., 2019).

One challenge of reactive melt transport modelling is the advection of the melt composition through its ascent. This part, which is mathematically well understood, being described by a mass balance equation, is numerically challenging (e.g. LeVeque, 1992). This is mainly due to the fact that most numerical models are based on a Eulerian frame of reference, where the discretised space is fixed in space and in time. In contrast, transport is by essence better defined from a Lagrangian perspective, where the observer follows the particles of fluid as they move. In addition, two-phase flow models are at least 2D problems due to the formation of channels (e.g. Barcilon and Lovera, 1989; Connolly and Podladchikov, 2007b) and to the fact that mass cannot be transported efficiently in 1D in the melt (Jordan et al., 2018). This brings a limitation to the resolution of the models and hence requires accurate advection schemes.

This study compares four numerical schemes applied to the problem of the advection of magmatic composition: an upwind scheme, a weighted essentially non-oscillatory (WENO) scheme, a semi-Lagrangian (SL) scheme, and a marker-in-cell (MIC) method. This selection provides a representation of the different approaches to solving advection problems that are commonly used in a wide range of applications. We assess the performance of each scheme in terms of accuracy, mass conservation, and computational time. A 2D model coupling chemical advection with a two-phase flow solver is then used to evaluate which algorithm is best suited for this problem.

Chemical transport in two-phase flow systems is described by the four mass conservation equations of the system (e.g. Aharonov et al., 1997). The first two equations describe the conservation of the total mass of the solid and the liquid:

where f and s represent the fluid and solid phases, *t* is the time (in s), *ϕ* is the fluid-filled porosity, *ρ* is the density of the respective phase (kg s^{−3}), and ** v** is the velocity of the respective phase (m s

^{−1}). The last two equations express the conservation of each chemical component within the solid and fluid phases:

where *C*_{e} is the mass fraction of the chemical component *e* in the respective phase, ${D}_{e}^{\mathrm{s}}$ is the solid diffusion coefficient of the chemical component *e* (in m^{2} s^{−1}), and ${\mathit{D}}_{e}^{\mathrm{f}}$ is the hydrodynamic dispersion tensor of the chemical component *e* in the fluid (m^{2} s^{−1}). These four equations assume no mass transfer due to reactions between the solid and the liquid phases.

## Simplifications

In this study, the advection of the chemical components transported by the liquid phase is considered, and the diffusion term in Eqs. (3) and (4) is neglected. Since *ρ*_{s} is assumed to be constant and the host rock has a fixed composition, Eq. (3) is omitted.

Subtracting Eq. (2) in Eq. (4), and dividing by *ρ*_{f} and *ϕ* yields

Equation (5) is formally equivalent to Eq. (4) without the dispersion term. Moreover, Eq. (4) is written in conservative form, whereas Eq. (5) is expressed in Lagrangian or non-conservative form. Equation (5) removes the time dependence on *ϕ* and is linear. It is a common form used in the reactive transport modelling community (e.g. Carrera et al., 2022).

An expression for *v*_{f} can be derived by coupling Eqs. (1) and (2) to the momentum conservation equations (e.g. McKenzie, 1984; Bercovici et al., 2001). These are usually solved before Eq. (5); a description of the system used in this study is provided below in Sect. 5.1.

Solving an advection equation using a linear Eulerian scheme leads to high numerical diffusion for first-order schemes, such as the upwind scheme (Courant et al., 1952), and to oscillations on sharp gradients for higher-order schemes (LeVeque, 2002). The latter effect is described by Godunov's theorem (Godunov and Bohachevsky, 1959). This theorem states that linear Eulerian schemes with an order of accuracy greater than 1 cannot preserve the monotonicity of the solution for sharp gradients, discontinuities, or shocks. This has led to extensive developments in the design of high-order Eulerian non-linear schemes that can achieve high accuracy without bringing oscillations. Examples of such developments are the essentially non-oscillatory (ENO) methods (Harten et al., 1987) that later led to WENO schemes (Liu et al., 1994). These schemes are based on the idea of using a non-linear adaptive procedure to automatically choose the locally smoothest stencil, and early examples of applications include the modelling of shocks appearing in acoustics (e.g. Grasso and Pirozzoli, 2000) or solving the Hamilton–Jacobi equations (e.g. Jiang and Peng, 2000).

Another approach is to use schemes closer to the Lagrangian perspective, such as the MIC (or alternatively named marker-and-cell) method (e.g. Harlow et al., 1955; Gerya and Yuen, 2003a). It consists of tracking individual particles on a Lagrangian frame and reinterpolating them when needed on a Eulerian stationary grid. This approach has the advantage of producing little numerical diffusion, being unconditionally stable, and has been extensively used in geodynamic models (e.g. Gerya, 2019; van Keken et al., 1997; Duretz et al., 2011).

Finally, there are intermediate methods, such as semi-Lagrangian methods, trying to take advantages from both Eulerian and Lagrangian schemes (Robert, 1981; McDonald, 1984). These schemes look at different particles at each time step, considering only particles whose final trajectories correspond to the position of grid nodes. This has the advantage of only considering a number of particles equal to the resolution of the Eulerian grid and is computationally efficient. They are also unconditionally stable but have issues with mass conservation (Chandrasekar, 2022). They were first developed for atmospheric modelling (e.g. Robert, 1981; Staniforth and Côté, 1991) and later successfully used in the plasma modelling community (e.g. Sonnendrücker et al., 1999).

To solve for Eq. (5) in the context of two-phase flow, we implemented and tested four different advection schemes that are representative of the approaches described above: an upwind scheme, a WENO scheme, an SL scheme, and a MIC method.

## 3.1 Upwind scheme

The upwind scheme is among the simplest algorithms for solving an advection equation on a Eulerian grid (e.g. LeVeque, 1992). It is explicit and first order in space and in time. It consists of using a spatially biased stencil that depends on the direction of the flow (Fig. 1).

### 3.1.1 Spatial discretisation

Using a first-order spatially biased stencil, Eq. (5) can be approximated for one chemical element and in 1D as

where *i* is a spatial index in the *x* direction, Δ*x* is the constant grid spacing, and *v*_{f,i} is the *x* component of the fluid velocity.

### 3.1.2 Temporal discretisation

Combined with the first-order forward Euler method, we retrieve the classical upwind scheme from Eq. (6):

where Δ*t* is the time step.

It can also be rewritten in a more compact form:

where

This scheme is well-known to produce a lot of numerical diffusion and is bounded by the following Courant–Friedrichs–Lewy (CFL) condition for *p* dimensions:

where *Co*_{max} is the maximum Courant (*Co*) number (e.g. Hirsch, 2007). In addition, it is not mass conservative for non-constant *v*_{f}, especially for divergent flow.

## 3.2 Weighted essentially non-oscillatory scheme

Weighted essentially non-oscillatory schemes were developed by Liu et al. (1994). The reader can refer to Shu (2009) for a comprehensive review of the development of WENO schemes and Pawar and San (2019) for implementations using Julia.

They are high-order schemes able to resolve sharp gradients, produce little numerical diffusion, but also follow the same CFL condition as the upwind scheme. The key idea behind them is to use a non-linear adaptive procedure to automatically choose the locally smoothest stencil. This allows WENO schemes to dispose of oscillations when advecting sharp gradients.

We use a fifth-order-in-space finite-difference approach for non-conservative problems, referenced as WENO-5 hereafter.

### 3.2.1 Spatial discretisation

Equation (5) can be discretised in space using the WENO-5 scheme similarly to the upwind scheme, in 1D, for one chemical element and for a single grid point such as

where

Here, ${C}_{i-\frac{\mathrm{1}}{\mathrm{2}}}^{L}$ and ${C}_{i+\frac{\mathrm{1}}{\mathrm{2}}}^{R}$ are omitted to avoid redundancy. They can be obtained by shifting the index by −1 and 1, respectively.

The non-linear weights *w* are defined as

The values of the optimal weights ${d}_{k}^{L}$ and ${d}_{k}^{R}$ are given in Table 1. *ϵ* represents the machine epsilon, the relative approximation error due to rounding in floating-point arithmetic, and is used to avoid division by zero.

Smoothness indicators *β* are equal to

The WENO-5 scheme requires, in 1D, a stencil of five points biased towards the left for positive velocities and five points biased towards the right for negative velocities as shown in Fig. 2. This commonly requires two ghost points on each side of the model to apply the boundary conditions. To extend this scheme to 2D, two new terms can be added to Eq. (9) for the new positive and negative component of *v*_{f}. The expressions of *C* at half points of the new index can be derived using the same formulae as in 1D for the new direction.

### 3.2.2 Temporal discretisation

Weighted essentially non-oscillatory schemes are not stable using the standard forward Euler time integration method (Wang and Spiteri, 2007). The most commonly used discretisation is the third-order strong stability preserving (SSP) explicit Runge–Kutta method (e.g. Jiang and Shu, 1996; Ghosh and Baeder, 2012). Strong stability preserving schemes are used to fully capture discontinuous solutions and are therefore well-suited for solving hyperbolic partial differential equations (Gottlieb et al., 2001).

The third-order SSP Runge–Kutta for Eq. (5) for one chemical element can be written as

with *L* being the spatial discretisation operator:

With this formulation, the WENO-5 scheme is fifth order in space and third order in time.

## 3.3 Semi-Lagrangian schemes

Semi-Lagrangian schemes take a different approach than classical Eulerian methods and are related to tracer-based advection schemes. Semi-Lagrangian schemes aim to use the best of Lagrangian and Eulerian methods by solving the problem for particles whose trajectories pass through a fixed grid at the end of each time step rather than recording the full history of individual particles. They are therefore unconditionally stable. Two steps are usually required to implement SL schemes: trajectory tracing and interpolation back to the grid. In this study, a backward-in-time SL scheme is used for the trajectory tracing and the quasi-monotone scheme developed by Bermejo and Staniforth (1992) for the interpolation.

### 3.3.1 Trajectory tracing

The advantage of backward-in-time SL schemes is that the interpolant is defined from the Eulerian grid. In the case of a regular grid, this reduces the complexity of the implementation and the numerical cost of the interpolation function, since the interpolant is deﬁned on a regular grid. From a particle point of view, the goal is to retrieve the position of the particle at time *t*_{n} for which the position corresponds to a grid point at time *t*_{n+1}. Using Eq. (5) for one chemical element and in 1D, the following ordinary differential equation has to be solved:

Knowing $x\left({t}_{n+\mathrm{1}}\right)={x}_{i}$, where *i* is a grid point, *x*(*t*_{n})=*x*_{d}, where *d* is a departure point that needs to be found. In most practical cases, the velocity field varies greatly in time and in space between each time step, especially for porosity waves, so it is not easy to determine *x*_{d}. A common approach to accurately determine *x*_{d} is to use a linear multistep method such as the implicit mid-point scheme (Robert, 1981):

where *v*_{f} at time *t*_{n+½} is obtained by taking the mean between the velocity at *n* and *n*+1. The assumption behind the mid-point rule is that the velocity remains constant at the mid-point value during each time step. This ensures that each trajectory is linear, with the mid-point being the average of the positions of its endpoints (Fig. 3). This method is a second-order accurate trajectory method in both space and time.

Equation (11) must be solved implicitly because *x*_{d} is present on both sides of the equation and therefore requires iterations. It can be achieved for *r* iterations in this form:

A minimum of three iterations while using linear interpolation has been shown to be sufficient in most cases (e.g. McDonald, 1984).

### 3.3.2 Interpolation

In most cases, *x*_{d} does not correspond to a grid node (see Fig. 3). In this case, interpolation is required to retrieve the value of the unknown at *x*_{d}:

where ℒ is an interpolation operator and ${x}_{{i}_{k}}$ represents the nodes of the cell containing *x*_{d}.

Commonly, cubic interpolants are used as they offer a good compromise between performance and accuracy (e.g. Chandrasekar, 2022) and require in 1D four grid points ${x}_{{i}_{k}}$ per particle. Cubic B-splines are used in this study. Godunov's theorem still applies to linear SL schemes, and since cubic interpolation is third order in space, it introduces oscillations and overshoots for high gradients. To overcome this limitation, quasi-monotone (QM) SL schemes were developed by Bermejo and Staniforth (1992). The term QM means that the scalar field values cannot exceed the range of the previous time step but can still develop wiggles inside that range. Quasi-monotonicity is equivalent to the notion of being essentially non-oscillatory (Bermejo, 2001). A disadvantage of this method is an increased numerical diffusion, especially for a high *Co* number. A maximum *Co* number of 1.5 is generally used (e.g. Smith, 2000).

To implement quasi-monotone semi-Lagrangian (QMSL) schemes, let us define *C*^{−} and *C*^{+} as the minimum and maximum scalar values of the nodes of the cell containing *x*_{d} and *C*^{H} as the high-order non-monotone interpolant, here a cubic spline. Then, a local clipping can be applied at the end of each time step:

where *C*^{M} is the quasi-monotone interpolant. Equation (13) can be rewritten in a more compact way:

Formally, this formulation is equivalent to a linear combination between a high-order interpolant and a first-order (monotone) interpolant (Bermejo, 2001).

## 3.4 Marker-in-cell schemes

Marker-in-cell schemes share the same ambition as SL schemes, such as being unconditionally stable, but are closer to Lagrangian schemes. They record the complete history of individual particles, called markers and interpolate their values on a fixed grid. This approach has the advantage of greatly reducing numerical diffusion and making MIC schemes unconditionally stable. In addition to trajectory tracing and interpolation, the MIC schemes require markers to be generated within the domain of the model.

### 3.4.1 Initial marker generation and reseeding of particles

The number of markers per cell required can vary depending on the complexity of the problem, here 5 markers per cell dimension and effectively 25 in 2D. This is generally sufficient to achieve good accuracy (e.g. Gerya, 2019). The initial value in each marker can then be directly derived from the initial conditions or obtained by linear interpolation from the initial conditions of the Eulerian grid.

For highly divergent flows or sometimes strongly stretching flows, it is necessary to regenerate or remove markers during the simulation. For highly divergent flows, this is because particles will accumulate in zones with negative divergence values and create a gap in zones with positive divergence values. For highly stretching flow, increases or gaps in the density of the markers can be induced by preferential flow in a particular direction. For reseeding, a non-conservative strategy similar to Keller et al. (2013) is used. If the marker density per cell is less than 25 % of the initial density, new markers are generated and assigned the value of the nearest marker. The old markers are discarded after this step. For marker accumulation, the marker density cannot exceed twice the initial density. If it does, a quarter of the markers are discarded at random.

### 3.4.2 Trajectory tracing

The goal of trajectory tracing for MIC schemes is to determine the position of each marker at the next time step. The same equation as Eq. (10) is solved. However, compared to backward SL where the final position is known, the unknown in this case is the position of the arrival point at *t*^{n+1}. Also, the common way to solve this equation for MIC schemes is not using a linear multistep method such as the implicit mid-point scheme but rather Runge–Kutta schemes (Gerya, 2019). Using a second-order Runge–Kutta scheme, it consists of four steps. Interpolating *v*_{f} at *t*_{n} at the departure point *x*_{d} of the markers, finding the position of the markers at *t*_{n+½} using *v*_{f} at *t*_{n}, reinterpolating the velocity at this new position, and using this new velocity to compute the arrival point *x*_{a} of the markers at *t*_{n+1} from *t*_{n}.

Since classical interpolants do not retain the physical properties of the velocity field, such as its divergence, a simple bilinear interpolation may lead to unphysical clustering of markers on the timescale of a numerical model. To address this issue, Pusok et al. (2017) explored different interpolants and showed the advantages of using the LinP interpolation scheme (Gerya, 2019). The LinP interpolation scheme is an empirical relationship that combines two linear interpolants defined at the sides and at the centre of each cell. It is defined as

with *A* a constant commonly equal to $\mathrm{2}/\mathrm{3}$, ℒ a linear interpolant, and *x*_{side} and *x*_{center} the position of the sides and centre of the cell containing the marker.

Using this definition, we can rewrite the four steps of the second-order Runge–Kutta scheme in mathematical notation with four equations:

with *x*_{h} the intermediate position of the marker. Solving these four equations successively to obtain a value for *x*_{a} with this method is second order in space but only first order in time, as only the velocity at *t*_{n} is used.

### 3.4.3 Interpolation

After calculating the position of the markers, it is necessary to interpolate back on the Eulerian grid and/or to update the values of the markers from the Eulerian grid depending on the problem being solved. The step of updating the markers is not described in detail here, as it is not used in this study, but involves a simple interpolant when a regular grid is used. This step is more complex for updating the advected field on the Eulerian grid because the markers are not uniformly distributed for a non-trivial velocity field. Therefore, contrary to SL schemes, interpolation is performed on an unstructured grid as it is based on the position of the markers. In most cases, linear interpolants are used because they prevent oscillations, and the marker densities are high enough to prevent numerical diffusion. In this study, a weighted-distance-averaging linear interpolant is used (Gerya, 2019):

where *x*_{m} is the position of a marker *m*, Δ*x*_{m} is the distance between a marker *m* and the grid point *i*, and *w* is the weight of a marker *m*. All markers found in the cells surrounding grid point *i* are used for interpolation, as in Gerya and Yuen (2003b). The relationship between the markers and the grid in 2D is summarised in Fig. 4. One disadvantage of this interpolant is that it is prone to race conditions in shared memory systems, as it involves two sums, which require the use of atomic operations to parallelise the implementation.

To test the four advection schemes, two different numerical tests are performed in 2D: the pure rotation of a cylinder and the advection through a more complex velocity field mimicking a convection cell. In both cases, the domain is a square of size 1.0×1.0 with a constant spacing of $\mathrm{\Delta}x=\mathrm{\Delta}y=\mathrm{0.005}$ for a resolution of 201×201 nodes. The radius of the cylinder is 24Δ*x* with a mass fraction of 1.0 and is centred at coordinates (0.25,0.25). The initial conditions for both tests are shown in Fig. 5.

For the first test, the time increment is Δ*t*=400 with $\mathit{\omega}=\mathit{\pi}\times {\mathrm{10}}^{-\mathrm{5}}$, so it takes 500 time steps to make a full revolution. The velocity is defined as $\mathit{v}=(-\mathit{\omega}(y-\mathrm{0.5}),\mathit{\omega}(x-\mathrm{0.5}\left)\right)$, so the rotation is anti-clockwise, and the centre of it is at coordinates (0.5,0.5). The *Co* numbers inside the cylinder range between 0.45 and 0.8. The test is stopped after two revolutions. For the second test, ** v** is defined as $\mathit{v}=(-\mathrm{2}\mathit{\pi}\times \mathrm{sin}(\mathit{\pi}x)\times \mathrm{cos}(\mathit{\pi}y),\mathrm{2}\mathit{\pi}\times \mathrm{cos}(\mathit{\pi}x)\times \mathrm{sin}(\mathit{\pi}y\left)\right)$. The time increment Δ

*t*is fixed by constraining the

*Co*number to be 0.7 for a total time of 0.8 and 1016 time steps. At half of the total time, the opposite sign of

**is taken as the new value of**

*v***, such that the analytical solution of the problem corresponds to the initial conditions.**

*v*To compare and quantify the results of the different schemes, the following quantities were monitored: the mass conservation (*M*), the total error (*E*_{tot}), the maximum value of the final mass fraction (max(*C*_{f})), and the minimum elapsed computational time of one time step after 10000 runs, with one thread and with 24 threads.

The mass conservation is defined as

where *k* is a 2D grid point index, *K* is the total number of grid points, ${C}_{\mathrm{f}}^{k}$ is the final mass fraction at index *k*, and ${C}_{\mathrm{0}}^{k}$ is the initial mass fraction at index *k*.

The total error in the scheme is defined as the mean square error:

The results for both tests are reported in Table 2 and in Figs. 6 and 7 for the four different schemes. Both tests show the strong numerical diffusion of the upwind scheme and its high *E*_{tot} due to its first order in time and space. The WENO-5 scheme shows no oscillation and a good accuracy, being fifth order in space and third in time in both tests, with a small mass loss in the second test. The QMSL is not mass conservative for both problems. It shows a relatively good accuracy, being third order in space and second in time; is monotone; but shows deformation of the original cylinder at the end of the second test. Finally, for both tests, the MIC scheme is mass conservative and monotone and shows the best accuracy with almost no numerical diffusion. These simple tests highlight the properties of each scheme but use a velocity field that is divergence-free and without sharp variations. Coupling with a two-phase flow system is therefore necessary to assert which scheme is the more suitable in this case.

Solving Eq. (5) for concrete cases implies having an expression for *v*_{f} at each time step. In this section, Eq. (5) is coupled to a transport model based on two-phase flow formalism. This transport model is used to model magma ascent in a porous solid phase. The main mechanism of transport is decompaction weakening, buoyancy, and failure and combines the formulations of Connolly and Podladchikov (2007a) and Vasilyev et al. (1998). It considers a compressible viscoelastic matrix with incompressible solid grains and an incompressible fluid phase, and it neglects the effect of shear stresses on fluid flow and compaction.

## 5.1 Two-phase flow formulation

In the case of a laminar fluid flow, conservation of momentum for the fluid can be expressed using Darcy's law:

with *P*_{f} the fluid pressure (in Pa), *k* the permeability (m^{2}), a function of the filled porosity or melt fraction *ϕ*, *μ*_{f} the fluid viscosity (Pa s), and ** g** the gravity vector (m s

^{−2}).

The relation between permeability and the filled porosity is assumed to follow the Kozeny–Carman law (Carman, 1939; Costa, 2006):

where *a* is a proportionality constant.

The effective pressure *P*_{e} is defined as the difference between lithostatic pressure and fluid pressure:

with *P*_{lith} the lithostatic pressure or the vertical load (in Pa). Substituting Eq. (16) in Eq. (15) and assuming constant rock density, we obtain

Considering the solid phase as a Maxwell body, we introduce rheology as the sum of viscous and poroelastic deformation:

where *ζ* is the volume viscosity of the rock (in Pa s), *b* a constant, and *β*_{ϕ} the pore compressibility modulus (Pa^{−1}). The terms on the right-hand side represent viscous and poroelastic deformation. Equation (18) is valid on the basis that shear stress is neglected.

The volume viscosity *ζ* is defined as a function of *ϕ* and *P*_{e}:

with *μ*_{s} the shear viscosity of the rock (in Pa s), *m* a constant, and *R* the decompaction weakening factor defined as the inverse of the *R* factor in Connolly and Podladchikov (2007a). *H*(*P*_{e}) is originally defined as the Heaviside function but is here approximated by a hyperbolic tangent function as similarly done by Räss et al. (2018).

We approximate here *β*_{ϕ} as the inverse of *G*, the shear modulus of the rock (in Pa):

This is valid for cylindrical pores, as described by Yarushina and Podladchikov (2015).

Summing up the right-hand sides of Eqs. (1) and (2) describing mass conservation and neglecting the change in densities, we obtain the total volumetric flux of material. Applying the divergence operator, we can derive

We can then substitute Eqs. (17) and (18) in Eq. (19) to obtain

In addition, developing Eq. (1) with the assumption that *ϕ* is much smaller than unity and substituting with Eq. (18) yields

Equation (20) can be seen as the mass conservation equation of the system, relating the flux densities of the solid and fluid phases. Equation (21) relates the evolution of porosity with the deformation of the solid phase. Solving these two coupled equations for *P*_{e} and *ϕ* allows the calculation of *v*_{s} and *v*_{f} from Eqs. (18) and (17) at each time step, making the link with Eq. (5).

## 5.2 Non-dimensionalisation and numerical approach

To mitigate numerical errors, a dimensionless scaling of the system is applied. The scaling variables are defined in Table 3. Using the scaling variables with Eqs. (17), (18), (20), and (21) and rearranging, we obtain the dimensionless system of equations:

where *φ*, *p*, *t*_{c}, *u*_{s}, and *u*_{f} are the dimensionless porosity, the dimensionless effective pressure, the dimensionless time, the dimensionless solid velocity, and the dimensionless fluid velocity, respectively. The Deborah number *De* is formally the ratio of the relaxation time to the observation time (Reiner, 1964), and here it characterises the ratio between viscous and elastic deformation. In the limit of small porosities, *u*_{s} can be neglected and only Eqs. (22), (23), and (25) are here solved.

Equations (22) and (23) are strongly coupled and highly stiff due to the non-linearity of the system and require an efficient numerical solver. DifferentialEquations.jl (Rackauckas and Nie, 2017), a robust ordinary differential equation (ODE) solver package, was used. This package has the advantage of simplicity, both in concept and in coding, and allows arbitrary orders of accuracy in time to be easily tested using different ODE solvers.

Equations (22) and (23) are first discretised in space using finite differences on a uniform Cartesian grid in 2D and then integrated in time using the trapezoidal rule with the second-order backward difference formula (TR-BDF2) scheme, an implicit scheme suitable for highly stiff problems (Bank et al., 1985) using DifferentialEquations.jl. It uses adaptive time-stepping and the Newton–Raphson method as a non-linear solver, using forward automatic differentiation to compute the Jacobian matrix (Revels et al., 2016). Knowing *φ* and *p*, Eq. (25) is then solved to compute *u*_{f} at each time step. The boundary conditions are periodic in all directions for all models. The system is then dimensionalised back.

## 5.3 Application to magmatic system

To assess the behaviour of the four advection schemes coupled with a two-phase flow system, we model the ascent of a magmatic anomaly. The spatial domain is a 2D regular grid of 450 by 900 m, and the total physical time is 1.5 Myr. The initial melt fraction distribution is defined using the following 2D Gaussian function:

with *ϕ*_{0} the background porosity defined as 0.1 %, *ϕ*_{max} the maximum porosity defined as 5 %, and *x*_{0} and *z*_{0} the centre of the anomaly. The standard deviation *σ* of the Gaussian is 30 m. All physical parameters and corresponding scaling variables used are reported in Table 3. The evolution of porosity is shown in Fig. 8. All models were performed on a single computer with an Intel Xeon Gold 6128 processor and 128 GB of RAM using Julia version 1.10.2. All models were computed on a CPU with multithreading using 24 threads.

The melt fraction is associated with two different arbitrary chemical compositions: a basaltic composition for the background melt fraction and an andesitic composition for the anomaly, corresponding to a circle with a radius of 60 m. The aim is not to model a realistic magmatic system but to investigate how the advection schemes can numerically affect the predictions of the model. The two compositions are reported in Table 4. No feedback between the melt compositions and the physical properties of the melt was considered to prevent the advection schemes from influencing the two-phase flow. In real settings, the effect of melt composition on melt viscosity and density is not negligible for these conditions. The maximum time step allowed for the two-phase flow is constrained by the *Co* number associated with the melt velocity. Its maximum value allowed for the upwind and the WENO-5 schemes is 0.7, but values of 0.7 and 1.5 for the QMSL and the MIC schemes were both used to take advantage of the extended stability of these schemes. The results for the evolution of the silica content in the melt are shown in Fig. 9 for the *Co* number of 0.7 and for the four algorithms at a resolution of 500×1000.

^{a} Recalculated from Giordano and Dingwell (2003). ^{b} Recalculated from Neuville et al. (1993).

As there is no analytical solution to this particular problem, it is not possible to directly calculate the numerical error in the different advection schemes. Nevertheless, we can compute the mass conservation of the advected quantities. The total mass of the melt composition is conserved, as it is re-normalised to 100 % at each time step. This is a constant-sum constraint and is characteristic of compositional data (Aitchison, 1982). However, it is not necessarily the case for each individual oxide. In that light, similar to Eq. (14), we monitor the mass conservation of each individual oxide *M*_{ox} in the melt at each time step:

where *ϕ*^{k} and ${\mathit{\varphi}}_{\mathrm{0}}^{k}$ are the current and initial porosity at index *k*, and ${C}_{\text{ox}}^{k}$ and ${C}_{{\text{ox}}_{\mathrm{0}}}^{k}$ are the current and initial composition of the oxide of interest in the melt at index *k*.

The melt fraction *ϕ* is conserved through the models, as Eqs. (22) and (23) are solved using a conservative discretisation. Therefore, *M*_{ox} only monitors the effects of the advection schemes for the oxide of interest. To quantify how the mass conservation evolves for each individual oxide and for each advection algorithm, the same model was performed at five different resolutions: 100×200, 200×400, 300×600, 400×800, and 500×1000. The values of the mass conservation of silica content ${M}_{{\text{SiO}}_{\text{2}}}$ for each resolution are shown in Fig. 10 for all the models. The values of the mass conservation of each oxide for all the algorithms are shown in Fig. 11. The total running time of each model is reported in Fig. 12.

The numerical models produced allow a better understanding of the process of passive chemical transport in magma within porosity waves and the impact of each advection scheme on the magma composition over time. All models confirm two distinct composition domains at the top of the porosity waves at the end of the simulations (Fig. 9). It is effectively a mixing of the compositions from the initial background porosity and from the anomaly. This is because melt is incorporated by the waves as they rise. This is attributed to the fact that the velocity of the porosity waves is higher than the melt velocity and has also been reported in previous studies (e.g. Jordan et al., 2018).

Comparing the results of the four algorithms, it is clear that the upwind scheme has the highest amount of numerical diffusion, which increases chemical mixing for non-physical reasons. The WENO-5 and QMSL exhibit similar results in terms of numerical diffusion, while the MIC shows the lowest amount with almost purely advective behaviour (Fig. 9). This is consistent with the numerical tests (Figs. 6 and 7). In terms of mass conservation, the oxide content is not conserved in all four schemes (Figs. 10 and 11). The WENO-5 gives the best results, with a mass conservation ranging from 98.87 % to 100.51 % for the lowest resolution to 99.85 % to 100.06 % at the highest resolution for all the oxides. The MIC performs similarly at high resolution, ranging from 99.77 % to 100.11 % for all the oxides at the highest resolution and with a significant increase in mass conservation from the 200×400 resolution. An improvement in the mass conservation for a *Co* number of 1.5 compared to the value of 0.7 is also noticeable. The QMSL shows slightly lower mass conservation for a *Co* of 0.7, ranging from 99.52 % to 100.22 % at the highest resolution and 97.22 % to 101.26 % at the lowest. In contrast to the MIC, there is a significant decrease for the mass conservation for a *Co* number of 1.5, ranging from 97.22 % to 102.81 % at low resolution to up to 97.77 % to 101.01 % at high resolution. The upwind scheme shows the worst values for mass conservation, ranging from 92.14 % to 103.55 % at low resolution to values of 96.97 % to 101.37 % at the highest resolution (Figs. 10 and 11). The better mass conservation of the MIC for a higher *Co* number can be explained by less reseeding and removal of markers, as the approach used is not mass conservative. On the other hand, the poorer performance concerning the QMSL at higher *Co* number can be interpreted as showing the decrease in accuracy of the trajectory tracing with increasing time step. The differences in mass conservation observed in the different oxides through all the models show that the initial conditions play a role and that higher values in the anomaly at the beginning of the model leads to a loss of mass (e.g. SiO_{2} or Al_{2}O_{3}), whereas the opposite leads to an excess of mass (e.g. CaO or K_{2}O). Also, the greater the relative difference between the composition of the oxide from the anomaly and the background melt fraction, the greater the mass conservation loss or gain (e.g. H_{2}O or FeO). However, it is observed that the mass conservation of all the oxides appears to converge towards 1 with increasing resolution for all methods.

In terms of performance, all schemes, except the MIC, show a run time of the same order of magnitude for a *Co* number of 0.7 at all resolutions. This is explained by the multithreading approach, which allows high performance, even for more computationally expensive algorithms due to parallelism. The high computational cost of the MIC is mainly attributed to the reseeding and removal of the markers due to the highly divergent velocity field. This part was not fully parallelised due to race conditions caused by the removal and addition of memory at run time. Also, the MIC and QMSL perform better for a *Co* number of 1.5 compared to 0.7 (Fig. 12). This is explained by a larger adaptive time step used by the two-phase flow solver due to the extended stability domain, which means that fewer time steps are required to solve the system. All the calculations were performed on a single CPU, and the code could be further optimised, especially for the MIC. However, this result provides an idea of the cost of each method when parallelised and highlights the complexity of fully parallelising a MIC algorithm while dealing with a significant amount of reseeding and marker removal.

The upwind scheme is considered inadequate for this problem due to its high numerical diffusion, and its lack of mass conservation for highly divergent velocity fields. The MIC scheme shows very good results in terms of accuracy with the least amount of numerical diffusion and has no stability condition. It also demonstrates better mass conservation with a higher time step. However, it is expensive in terms of computation and memory, as it needs to keep track of the markers. As the velocity field *v*_{f} is strongly divergent, it requires frequently regenerating and deleting markers, which adds complexity to the implementation and additional numerical cost. As a result, we consider this scheme to be too costly for this particular problem but recognise its robust qualities for other geodynamic problems where diffusion is not acceptable, such as thermomechanical deformation or mantle convection (e.g. Duretz et al., 2011; Ueda et al., 2015; Trim et al., 2020). On the other hand, the QMSL scheme shows very good performance with its extended stability field and good accuracy but has very poor mass conservation for a high *Co*, which contradicts the purpose of this scheme. A potential approach to improve the mass conservation would be to improve the trajectory tracing step, either by using higher-order multistep methods (e.g. Filbet and Prouveur, 2016) or by using Runge–Kutta schemes, similar to the one used for the MIC. On the basis of these results, the WENO-5 advection scheme appears to be the most appropriate for this problem. Mass conservation is a critical property for studying mass balance and mass transport problems associated with magma transport at different scales on Earth, and this algorithm obtains the best results. It also has good accuracy and reasonable performance and is easy to extend to higher dimensions and to parallelise.

In this study, a series of tests were carried out to determine which advection scheme is the most suitable for modelling the chemical transport of magma. Four of the most commonly used algorithms in the literature were compared: the upwind, WENO-5, MIC, and QMSL schemes. To test them, we combined a 2D two-phase flow model, which describes the evolution of the melt fraction of magma over time, with the chemical advection of its composition.

All algorithms, with the exception of the upwind scheme, are able to predict the melt composition with reasonable accuracy. However, mass conservation of each individual oxide in the melt is not fully achieved for any of the schemes. The MIC, while showing the least amount of numerical diffusion, requires a very large amount of reseeding and removal of markers due to the strongly divergent melt velocity field. This procedure is costly and requires reallocating memory at run time, complicating the implementation. The QMSL scheme has the worst mass conservation of the three algorithms, especially at high *Co*. This could potentially be improved by refining the trajectory tracing step to make it a more valuable alternative. The WENO-5 scheme shows the best results in mass conservation, even at low resolutions, is explicit, is easy to implement, and extends in 3D, although it is constrained by the CFL condition. On the basis of these results, the WENO-5 scheme is the most suitable to use for transporting magma composition during magma ascent. This is also applicable to problems using similar formulations, such as chemical advection in aqueous fluids, and makes WENO-5 a suitable scheme for modelling reactive transport under crustal or mantle conditions.

The code used in this study allowing reproducibility of the data is available on GitHub (https://github.com/neoscalc/ChemicalAdvectionPorosityWave.jl, last access: 2 July 2024) and at a permanent DOI repository (Zenodo): https://doi.org/10.5281/zenodo.8411354 (Dominguez, 2024a). The code is written in the Julia programming language. Refer to the repository's README for additional information. The code is distributed under the GPL-3.0 license. The data produced and used during this study are available at a permanent DOI repository (Zenodo): https://doi.org/10.5281/zenodo.13306073 (Dominguez, 2024b).

HD and PL conceptualised the project. PL acquired the funding. HD conducted the study with input from NR. HD wrote the code and the original manuscript. NR and PL revised and edited the manuscript. All authors contributed to the discussions and interpretation of the results.

The contact author has declared that none of the authors has any competing interests.

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.

Hugo Dominguez thanks Thibault Duretz and Ludovic Räss for discussions on the MIC implementation and the two-phase flow, respectively. The authors are grateful to Marcin Dabrowski and Albert de Montserrat Navarro for their thorough and constructive review of the original draft and Mauro Cacace for his editorial handling.

Funding was provided by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant no. 850530).

This paper was edited by Mauro Cacace and reviewed by Marcin Dabrowski and Albert de Montserrat Navarro.

Aharonov, E., Whitehead, J. A., Kelemen, P. B., and Spiegelman, M.: Channeling Instability of Upwelling Melt in the Mantle, J. Geophys. Res.-Sol. Ea., 100, 20433–20450, https://doi.org/10.1029/95JB01307, 1995a. a

Aharonov, E., Whitehead, J. A., Kelemen, P. B., and Spiegelman, M.: Channeling Instability of Upwelling Melt in the Mantle, J. Geophys. Res.-Sol. Ea., 100, 20433–20450, https://doi.org/10.1029/95JB01307, 1995b. a

Aharonov, E., Spiegelman, M., and Kelemen, P.: Three-Dimensional Flow and Reaction in Porous Media: Implications for the Earth's Mantle and Sedimentary Basins, J. Geophys. Res.-Sol. Ea., 102, 14821–14833, https://doi.org/10.1029/97JB00996, 1997. a, b

Aitchison, J.: The Statistical Analysis of Compositional Data, J. Roy. Stat. Soc. B, 44, 139–160, https://doi.org/10.1111/j.2517-6161.1982.tb01195.x, 1982. a

Bank, R., Coughran, W., Fichtner, W., Grosse, E., Rose, D., and Smith, R.: Transient Simulation of Silicon Devices and Circuits, IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 4, 436–451, https://doi.org/10.1109/TCAD.1985.1270142, 1985. a

Barcilon, V. and Lovera, O. M.: Solitary Waves in Magma Dynamics, J. Fluid Mech., 204, 121, https://doi.org/10.1017/S0022112089001680, 1989. a

Bercovici, D., Ricard, Y., and Schubert, G.: A Two-Phase Model for Compaction and Damage: 1. General Theory, J. Geophys. Res.-Sol. Ea., 106, 8887–8906, https://doi.org/10.1029/2000JB900430, 2001. a

Bermejo, R.: Analysis of a Class of Quasi-Monotone and Conservative Semi-Lagrangian Advection Schemes, Numerische Mathematik, 87, 597–623, https://doi.org/10.1007/PL00005425, 2001. a, b

Bermejo, R. and Staniforth, A.: The Conversion of Semi-Lagrangian Advection Schemes to Quasi-Monotone Schemes, Mon. Weather Rev., 120, 2622–2632, https://doi.org/10.1175/1520-0493(1992)120<2622:TCOSLA>2.0.CO;2, 1992. a, b

Bessat, A., Pilet, S., Podladchikov, Y. Y., and Schmalholz, S. M.: Melt Migration and Chemical Differentiation by Reactive Porosity Waves, Geochem. Geophys. Geosyst., 23, e2021GC009963, https://doi.org/10.1029/2021GC009963, 2022. a

Bouilhol, P., Connolly, J. A., and Burg, J.-P.: Geological Evidence and Modeling of Melt Migration by Porosity Waves in the Sub-Arc Mantle of Kohistan (Pakistan), Geology, 39, 1091–1094, https://doi.org/10.1130/G32219.1, 2011. a

Brown, M.: Granite: From Genesis to Emplacement, Geol. Soc. Am. Bull., 125, 1079–1113, https://doi.org/10.1130/B30877.1, 2013. a

Carman, P. C.: Permeability of Saturated Sands, Soils and Clays, J. Agric. Sci., 29, 262–273, https://doi.org/10.1017/S0021859600051789, 1939. a

Carrera, J., Saaltink, M. W., Soler-Sagarra, J., Wang, J., and Valhondo, C.: Reactive Transport: A Review of Basic Concepts with Emphasis on Biochemical Processes, Energies, 15, 925, https://doi.org/10.3390/en15030925, 2022. a

Chandrasekar, A.: Numerical Methods for Atmospheric and Oceanic Sciences, Cambridge University Press, 1st edn., ISBN 978-1-00-911923-8 978-1-00-910056-4, https://doi.org/10.1017/9781009119238, 2022. a, b

Clemens, J. D., Bryan, S. E., Mayne, M. J., Stevens, G., and Petford, N.: How Are Silicic Volcanic and Plutonic Systems Related? Part 1: A Review of Geological and Geophysical Observations, and Insights from Igneous Rock Chemistry, Earth-Sci. Rev., 235, 104249, https://doi.org/10.1016/j.earscirev.2022.104249, 2022. a

Connolly, J. A. D. and Podladchikov, Y. Y.: Decompaction Weakening and Channeling Instability in Ductile Porous Media: Implications for Asthenospheric Melt Segregation, J. Geophys. Res., 112, B10205, https://doi.org/10.1029/2005JB004213, 2007a. a, b

Connolly, J. A. D. and Podladchikov, Y. Y.: Decompaction Weakening and Channeling Instability in Ductile Porous Media: Implications for Asthenospheric Melt Segregation, J. Geophys. Res., 112, B10205, https://doi.org/10.1029/2005JB004213, 2007b. a, b

Connolly, J. A. D. and Podladchikov, Yu. Yu.: Compaction-Driven Fluid Flow in Viscoelastic Rock, Geodinam. Acta, 11, 55–84, https://doi.org/10.1080/09853111.1998.11105311, 1998. a

Costa, A.: Permeability-Porosity Relationship: A Reexamination of the Kozeny-Carman Equation Based on a Fractal Pore-Space Geometry Assumption, Geophys. Res. Lett., 33, L02318, https://doi.org/10.1029/2005GL025134, 2006. a

Courant, R., Isaacson, E., and Rees, M.: On the Solution of Nonlinear Hyperbolic Differential Equations by Finite Differences, Commun, Pure Appl. Math., 5, 243–255, https://doi.org/10.1002/cpa.3160050303, 1952. a

Dominguez, H.: neoscalc/ChemicalAdvectionPorosityWave.jl: ChemicalAdvectionPorosityWave.jl 1.0.0 (v1.0.0), Zenodo [code], https://doi.org/10.5281/zenodo.8411354, 2024a. a

Dominguez, H.: Data produced and used in the publication “Modelling chemical advection during magma ascent” by Dominguez et al., 2024, Zenodo [data set], https://doi.org/10.5281/zenodo.13306073, 2024b. a

Duretz, T., May, D. A., Gerya, T. V., and Tackley, P. J.: Discretization Errors and Free Surface Stabilization in the Finite Difference and Marker-in-Cell Method for Applied Geodynamics: A Numerical Study, Geochem. Geophys. Geosyst., 12, 7, https://doi.org/10.1029/2011GC003567, 2011. a, b

Filbet, F. and Prouveur, C.: High Order Time Discretization for Backward Semi-Lagrangian Methods, J. Comput. Appl. Math., 303, 171–188, https://doi.org/10.1016/j.cam.2016.01.024, 2016. a

Gerya, T. V.: Introduction to Numerical Geodynamic Modelling, Cambridge University Press, 2nd edn., ISBN 978-1-316-53424-3 978-1-107-14314-2, https://doi.org/10.1017/9781316534243, 2019. a, b, c, d, e

Gerya, T. V. and Yuen, D. A.: Characteristics-Based Marker-in-Cell Method with Conservative Finite-Differences Schemes for Modeling Geological Flows with Strongly Variable Transport Properties, Phys. Earth Planet. In., 140, 293–318, https://doi.org/10.1016/j.pepi.2003.09.006, 2003a. a

Gerya, T. V. and Yuen, D. A.: Characteristics-Based Marker-in-Cell Method with Conservative Finite-Differences Schemes for Modeling Geological Flows with Strongly Variable Transport Properties, Phys. Earth Planet. In., 140, 293–318, https://doi.org/10.1016/j.pepi.2003.09.006, 2003b. a

Ghosh, D. and Baeder, J. D.: Compact Reconstruction Schemes with Weighted ENO Limiting for Hyperbolic Conservation Laws, SIAM J. Sci. Comput., 34, A1678–A1706, https://doi.org/10.1137/110857659, 2012. a

Giordano, D. and Dingwell, D. B.: Non-Arrhenian Multicomponent Melt Viscosity: A Model, Earth Planet. Sc. Lett., 208, 337–349, https://doi.org/10.1016/S0012-821X(03)00042-6, 2003. a

Godunov, S. K. and Bohachevsky, I.: Finite Difference Method for Numerical Computation of Discontinuous Solutions of the Equations of Fluid Dynamics, Matematičeskij sbornik, 47, 271–306, 1959. a

Gottlieb, S., Shu, C.-W., and Tadmor, E.: Strong Stability-Preserving High-Order Time Discretization Methods, SIAM Review, 43, 89–112, https://doi.org/10.1137/S003614450036757X, 2001. a

Grasso, F. and Pirozzoli, S.: Shock-Wave-Vortex Interactions: Shock and Vortex Deformations, and Sound Production, Theor. Comput. Fluid Dynam., 13, 421–456, https://doi.org/10.1007/s001620050121, 2000. a

Harlow, F. H., Evans, M., and Richtmyer, R. D.: A Machine Calculation Method for Hydrodynamic Problems, Los Alamos Scientific Laboratory of the University of California, 1955. a

Harten, A., Engquist, B., Osher, S., and Chakravarthy, S. R.: Uniformly High Order Accurate Essentially Non-oscillatory Schemes, III, in: Upwind and High-Resolution Schemes, edited by: Hussaini, M. Y., van Leer, B., and Van Rosendale, J., 218–290, Springer, Berlin, Heidelberg, ISBN 978-3-642-60543-7, https://doi.org/10.1007/978-3-642-60543-7_12, 1987. a

Hirsch, C.: Numerical Computation of Internal and External Flows: Fundamentals of Computational Fluid Dynamics, Elsevier/Butterworth-Heinemann, Oxford ; Burlington, MA, 2nd edn., ISBN 978-0-7506-6594-0, 2007. a

Jackson, M., Gallagher, K., Petford, N., and Cheadle, M.: Towards a Coupled Physical and Chemical Model for Tonalite–Trondhjemite–Granodiorite Magma Formation, Lithos, 79, 43–60, https://doi.org/10.1016/j.lithos.2004.05.004, 2005. a

Jha, K., Parmentier, E. M., and Phipps Morgan, J.: The Role of Mantle-Depletion and Melt-Retention Buoyancy in Spreading-Center Segmentation, Earth Planet. Sc. Lett., 125, 221–234, https://doi.org/10.1016/0012-821X(94)90217-8, 1994. a

Jiang, G.-S. and Peng, D.: Weighted ENO Schemes for Hamilton–Jacobi Equations, SIAM J. Sci. Comput., 21, 2126–2143, https://doi.org/10.1137/S106482759732455X, 2000. a

Jiang, G.-S. and Shu, C.-W.: Efficient Implementation of Weighted ENO Schemes, J. Comput. Phys., 126, 202–228, https://doi.org/10.1006/jcph.1996.0130, 1996. a

Johnson, T., Yakymchuk, C., and Brown, M.: Crustal Melting and Suprasolidus Phase Equilibria: From First Principles to the State-of-the-Art, Earth-Sci. Rev., 221, 103778, https://doi.org/10.1016/j.earscirev.2021.103778, 2021. a

Jordan, J. S., Hesse, M. A., and Rudge, J. F.: On Mass Transport in Porosity Waves, Earth Planet. Sc. Lett., 485, 65–78, https://doi.org/10.1016/j.epsl.2017.12.024, 2018. a, b

Katz, R. F.: Magma Dynamics with the Enthalpy Method: Benchmark Solutions and Magmatic Focusing at Mid-ocean Ridges, J. Petrol., 49, 2099–2121, https://doi.org/10.1093/petrology/egn058, 2008. a

Katz, R. F. and Weatherley, S. M.: Consequences of Mantle Heterogeneity for Melt Extraction at Mid-Ocean Ridges, Earth Planet. Sc. Lett., 335–336, 226–237, https://doi.org/10.1016/j.epsl.2012.04.042, 2012. a

Katz, R. F., Jones, D. W. R., Rudge, J. F., and Keller, T.: Physics of Melt Extraction from the Mantle: Speed and Style, Annu. Rev. Earth Planet. Sci., 50, 507–540, https://doi.org/10.1146/annurev-earth-032320-083704, 2022. a

Kelemen, P. B., Hirth, G., Shimizu, N., Spiegelman, M., and Dick, H. J.: A Review of Melt Migration Processes in the Adiabatically Upwelling Mantle beneath Oceanic Spreading Ridges, Philos. T. Roy. Soc. Lond. A, 355, 283–318, 1997. a

Keller, T.: Numerical Modeling of Magma Ascent and Emplacement in the Continental Lithosphere and Crust, PhD thesis, ETH Zurich, https://doi.org/10.3929/ETHZ-A-010192539, 2013. a

Keller, T. and Katz, R. F.: The Role of Volatiles in Reactive Melt Transport in the Asthenosphere, J. Petrol., 57, 1073–1108, https://doi.org/10.1093/petrology/egw030, 2016. a

Keller, T., May, D. A., and Kaus, B. J. P.: Numerical Modelling of Magma Dynamics Coupled to Tectonic Deformation of Lithosphere and Crust, Geophys. J. Int., 195, 1406–1442, https://doi.org/10.1093/gji/ggt306, 2013. a

LeVeque, R. J.: Numerical Methods for Conservation Laws, Lectures in Mathematics ETH Zürich, Birkhäuser Verlag, Basel, Boston, 2nd edn., ISBN 978-3-7643-2723-1 978-0-8176-2723-2, 1992. a, b

LeVeque, R. J.: Finite Volume Methods for Hyperbolic Problems, Cambridge Texts in Applied Mathematics, Cambridge University Press, Cambridge, ISBN 978-0-521-00924-9, https://doi.org/10.1017/CBO9780511791253, 2002. a

Liu, X.-D., Osher, S., and Chan, T.: Weighted Essentially Non-oscillatory Schemes, J. Comput. Phys., 115, 200–212, https://doi.org/10.1006/jcph.1994.1187, 1994. a, b

McDonald, A.: Accuracy of Multiply-Upstream, Semi-Lagrangian Advective Schemes, Mon. Weather Rev., 112, 1267–1275, https://doi.org/10.1175/1520-0493(1984)112<1267:AOMUSL>2.0.CO;2, 1984. a, b

McKenzie, D.: The Generation and Compaction of Partially Molten Rock, J. Petrol., 25, 713–765, https://doi.org/10.1093/petrology/25.3.713, 1984. a, b, c

Neuville, D. R., Courtial, P., Dingwell, D. B., and Richet, P.: Thermodynamic and Rheological Properties of Rhyolite and Andesite Melts, Contrib. Mineral. Petrol., 113, 572–581, https://doi.org/10.1007/BF00698324, 1993. a

Omlin, S., Malvoisin, B., and Podladchikov, Y. Y.: Pore Fluid Extraction by Reactive Solitary Waves in 3-D: Reactive Porosity Waves, Geophys. Res. Lett., 44, 9267–9275, https://doi.org/10.1002/2017GL074293, 2017. a

Pawar, S. and San, O.: CFD Julia: A Learning Module Structuring an Introductory Course on Computational Fluid Dynamics, Fluids, 4, 159, https://doi.org/10.3390/fluids4030159, 2019. a

Pusok, A. E., Kaus, B. J. P., and Popov, A. A.: On the Quality of Velocity Interpolation Schemes for Marker-in-Cell Method and Staggered Grids, Pure Appl. Geophys., 174, 1071–1089, https://doi.org/10.1007/s00024-016-1431-8, 2017. a

Rackauckas, C. and Nie, Q.: DifferentialEquations.Jl – A Performant and Feature-Rich Ecosystem for Solving Differential Equations in Julia, J. Open Res. Softw., 5, 15, https://doi.org/10.5334/jors.151, 2017. a

Räss, L., Simon, N. S. C., and Podladchikov, Y. Y.: Spontaneous Formation of Fluid Escape Pipes from Subsurface Reservoirs, Sci. Rep., 8, 11116, https://doi.org/10.1038/s41598-018-29485-5, 2018. a

Reiner, M.: The Deborah Number, Physics Today, 17, 62–62, https://doi.org/10.1063/1.3051374, 1964. a

Revels, J., Lubin, M., and Papamarkou, T.: Forward-Mode Automatic Differentiation in Julia, arXiv:1607.07892 [cs.MS], 2016. a

Riel, N., Bouilhol, P., van Hunen, J., Cornet, J., Magni, V., Grigorova, V., and Velic, M.: Interaction between Mantle-Derived Magma and Lower Arc Crust: Quantitative Reactive Melt Flow Modelling Using STyx, Geological Society, London, Special Publications, 478, 65–87, https://doi.org/10.1144/SP478.6, 2019. a

Robert, A.: A Stable Numerical Integration Scheme for the Primitive Meteorological Equations, Atmosphere-Ocean, 19, 35–46, https://doi.org/10.1080/07055900.1981.9649098, 1981. a, b, c

Scott, D. R. and Stevenson, D. J.: Magma Solitons, Geophys. Res. Lett., 11, 1161–1164, https://doi.org/10.1029/GL011i011p01161, 1984. a

Shu, C.-W.: High Order Weighted Essentially Nonoscillatory Schemes for Convection Dominated Problems, SIAM Review, 51, 82–126, https://doi.org/10.1137/070679065, 2009. a

Smith, C. J.: The Semi-Lagrangian Method in Atmospheric Modelling, PhD thesis, 2000. a

Solano, J. M. S., Jackson, M. D., Sparks, R. S. J., Blundy, J. D., and Annen, C.: Melt Segregation in Deep Crustal Hot Zones: A Mechanism for Chemical Differentiation, Crustal Assimilation and the Formation of Evolved Magmas, J. Petrol., 53, 1999–2026, https://doi.org/10.1093/petrology/egs041, 2012. a

Sonnendrücker, E., Roche, J., Bertrand, P., and Ghizzo, A.: The Semi-Lagrangian Method for the Numerical Resolution of the Vlasov Equation, J. Comput. Phys., 149, 201–220, https://doi.org/10.1006/jcph.1998.6148, 1999. a

Spiegelman, M. and Kenyon, P.: The Requirements for Chemical Disequilibrium during Magma Migration, Earth Planet. Sc. Lett., 109, 611–620, 1992. a

Spiegelman, M., Kelemen, P. B., and Aharonov, E.: Causes and Consequences of Flow Organization during Melt Transport: The Reaction Infiltration Instability in Compactible Media, J. Geophys. Res.-Sol. Ea., 106, 2061–2077, https://doi.org/10.1029/2000JB900240, 2001. a

Staniforth, A. and Côté, J.: Semi-Lagrangian Integration Schemes for Atmospheric Models – A Review, Mon. Weather Rev., 119, 2206–2223, https://doi.org/10.1175/1520-0493(1991)119<2206:SLISFA>2.0.CO;2, 1991. a

Trim, S. J., Lowman, J. P., and Butler, S. L.: Improving Mass Conservation With the Tracer Ratio Method: Application to Thermochemical Mantle Flows, Geochem. Geophys. Geosyst., 21, e2019GC008799, https://doi.org/10.1029/2019GC008799, 2020. a

Ueda, K., Willett, S., Gerya, T., and Ruh, J.: Geomorphological–Thermo-Mechanical Modeling: Application to Orogenic Wedge Dynamics, Tectonophysics, 659, 12–30, https://doi.org/10.1016/j.tecto.2015.08.001, 2015. a

van Keken, P. E., King, S. D., Schmeling, H., Christensen, U. R., Neumeister, D., and Doin, M.-P.: A Comparison of Methods for the Modeling of Thermochemical Convection, J. Geophys. Res.-Sol. Ea., 102, 22477–22495, https://doi.org/10.1029/97JB01353, 1997. a

Vasilyev, O. V., Podladchikov, Y. Y., and Yuen, D. A.: Modeling of Compaction Driven Flow in Poro-Viscoelastic Medium Using Adaptive Wavelet Collocation Method, Geophys. Res. Lett., 25, 3239–3242, https://doi.org/10.1029/98GL52358, 1998. a

Wang, R. and Spiteri, R. J.: Linear Instability of the Fifth-Order WENO Method, SIAM J. Numer. Anal., 45, 1871–1901, https://doi.org/10.1137/050637868, 2007. a

Yarushina, V. M. and Podladchikov, Y. Y.: (De)Compaction of Porous Viscoelastoplastic Media: Model Formulation, J. Geophys. Res.-Sol. Ea., 120, 4146–4170, https://doi.org/10.1002/2014JB011258, 2015. a