Fluid flow channeling and mass transport with discontinuous porosity distribution

The flow of fluids within porous rocks is an important process with numerous applications in Earth sciences. Modeling the compaction-driven fluid flow requires the solution of coupled nonlinear partial differential equations that account for the fluid flow and the solid deformation within the porous medium. Despite the nonlinear relation of porosity and permeability that is commonly encountered, natural data show evidence of channelized fluid flow in rocks that have an overall layered structure. Layers of different rock types have discontinuous hydraulic and mechanical properties. We present numerical results obtained by a novel space-time method, which can handle discontinuous initial porosity (and permeability) distributions efficiently. The space-time method enables straightforward coupling to models of mass transport for trace elements. Our results indicate that, under certain conditions, the discontinuity of the initial porosity influences the distribution of incompatible trace elements, leading to sharp concentration gradients and large degrees of elemental enrichment. Finally, our results indicate that the enrichment of trace elements depends not only on the channelization of the flow but also on the interaction of fluid-filled channels with layers of different porosity and permeability.

1 Introduction

The flow of fluids in the Earth's subsurface is important for many applications. Examples of such applications include, but are not limited to, the migration of magma (McKenzie, 1984; Barcilon and Richter, 1986), the flow of glaciers (Fowler, 1984), the integrity of subsurface reservoirs (Räss et al., 2018; Yarushina et al., 2022), and the efficiency of geothermal systems (Utkin and Afanasyev, 2021). A distinctive aspect of the fluid flow within the deep Earth is that rocks cannot be treated as purely elastic or rigid, requiring consideration of their bulk (volumetric) viscous deformation (McKenzie, 1984; Scott and Stevenson, 1986). In fact, recent experiments have confirmed that the viscous/viscoelastic behavior of rocks can be observed also at near-surface conditions (Sabitova et al., 2021). Thus, the volumetric deformation and the associated fluid flow need to be considered in a coupled fashion since (de)compaction can drive fluid flow and vice versa (Connolly and Podladchikov, 1998; Vasilyev et al., 1998). In the latter studies, porosity waves were observed numerically. Such waves reflect the propagation of porosity perturbations (and the associated volumetric deformation) in a wave-like fashion with minimal dissipation. The transport of fluid-filled porosity in a non-dissipative fashion has been at the focus of research by geoscientists since it has important implications for the geochemical anomalies that are observed near the surface of the Earth (Richter, 1986; Navon and Stolper, 1987; Jordan et al., 2018).

The shape of porosity waves has been shown to depend very sensitively on the nonlinear behavior of the bulk (volumetric) viscosity. For example, in cases where the compaction/decompaction behavior is associated with significant changes in the effective viscosity, porosity waves take a channel-like shape (in two or three dimensions) that is responsible for the focusing of the flow towards the Earth's surface (Räss et al., 2018; Connolly and Podladchikov, 2007; Räss et al., 2014, 2019; Yarushina and Podladchikov, 2015; Yarushina et al., 2015, 2020). The focusing of the flow produces “chimney-like” features that resemble geophysical observations (Räss et al., 2018; Yarushina et al., 2020). The occurrence of such features is very important in the quantification of fluid flow and the associated geochemical anomalies (Spiegelman and Kelemen, 2003).

An essential feature of geological formations is that rocks are typically found in layers (strata). The layers are often composed of rock types that have different physical properties, such as porosity and permeability. It is exactly this change in permeability that is responsible for the formation of geological reservoirs. For example, a typical underground reservoir must be composed of rocks of high porosity (and permeability) and must be covered by rocks of negligible porosity (and permeability) that act as a “seal” to the underlying rock units. This configuration typically requires the consideration of porosity (and permeability) jump discontinuities across the lithological boundaries. However, the methods used to solve the respective poro-viscoelastic equations numerically cannot handle a discontinuous initial porosity, and hence only approximate it by a continuous function with steep gradient. This approach leads to smoothing effects and does not preserve the discontinuous nature of solutions. Resolving the solution behavior next to a discontinuity is crucial in all the applications where the quantification of the fluid flow is needed and can thus be important for safety analyses in geoengineering applications (Yarushina et al., 2022).

Here, we consider a poro-viscoelastic model that generalizes the one introduced in Connolly and Podladchikov (1998) and Vasilyev et al. (1998) for the interaction of porosity and pressure. For modeling sharp transitions between materials, as caused, for example, by stacked rock layers, it is important to be able to treat porosities with jump discontinuities. These discontinuities turn out to be determined mainly by the initial condition, as it was shown in Bachmayr et al. (2023) based on results from Simpson et al. (2006) for smooth initial porosities. In addition, we utilize a newly-developed space-time method that has been shown to be more accurate in solving this particular problem in the presence of jump discontinuities (Bachmayr and Boisserée, 2025). Our approach can be used to benchmark methods that do not include discontinuities and quantify the error between the two approaches. An additional advantage of the space-time method is that it can be coupled to simple models of chemical-tracer transport (see, for example, Richter, 1986; Jordan et al., 2018) as a post-processing step, since the entire porosity-pressure history is saved and the chemical transport problem does not feedback into the porosity-pressure (hydro-mechanical) model. The results obtained from this coupling allow us, for the first time, to investigate the evolution of chemical anomalies in the presence of channelized fluid flow, and their interactions with porosity/permeability discontinuities. In particular, our results are relevant to the formation of ore deposits and to the transport of trace elements in the subsurface.

1.1 The governing equations

The model for poro-viscoelastic flow on which we focus in this work reads

$$\begin{array}{}\text{(1a)}& {\displaystyle }{\partial }_t\mathit{\varphi }=-(\mathrm{1}-\mathit{\varphi })\left({\mathit{\varphi }}^m{\displaystyle \frac{p}{{\mathit{\eta }}_{\mathrm{b}}\mathit{\sigma }\left(p\right)}}+Q{\partial }_{t}p\right),\mathit{\varphi }(\mathrm{0},\cdot )={\mathit{\varphi }}_{\mathrm{0}},\text{(1b)}& {\displaystyle }\begin{array}{rl}{\partial }_{t}p& ={\displaystyle \frac{\mathrm{1}}{Q}}\left({\mathrm{div}}_x\left({\displaystyle \frac{k_{\mathrm{b}}}{\mathit{\mu }{\mathit{\varphi }}_{\mathrm{b}}^n}}{\mathit{\varphi }}^n\left({\mathrm{\nabla }}_{x}p+(\mathrm{1}-\mathit{\varphi })\mathit{\delta }\mathit{\rho }g{\mathit{e}}_{\mathrm{d}}\right)\right)\right.\\ & \left.-{\mathit{\varphi }}^m{\displaystyle \frac{p}{{\mathit{\eta }}_{\mathrm{b}}\mathit{\sigma }\left(p\right)}}\right),p(\mathrm{0},\cdot )=p_{\mathrm{0}},\end{array}\end{array}$$

as previously described in Connolly and Podladchikov (1998) and Vasilyev et al. (1998). Here, ϕ denotes the porosity (void ratio), p is the effective pressure, σ accounts for decompaction weakening (Räss et al., 2018, 2019), Q is the compressibility (equal to K⁻¹, where K is the bulk modulus), and $\mathit{\delta }\mathit{\rho }={\mathit{\rho }}^{\mathrm{s}}-{\mathit{\rho }}^{\mathrm{f}}$ the density difference. Furthermore, ${\mathrm{\nabla }}_{x}f=({\partial }_{x_{\mathrm{1}}}f,\mathrm{\dots },{\partial }_{x_d}f{)}^{\top }$ and ${\mathrm{div}}_x\mathit{f}=\sum _{i=\mathrm{1}}^d{\partial }_{x_i}{f}_i$ for functions $f:{\mathbb{R}}^d\to \mathbb{R}$, $\mathbf{\text{f}}:{\mathbb{R}}^d\to {\mathbb{R}}^d$ as usual. For any function g(t,x), we denote $g(\mathrm{0},\cdot )$ as the function g at a fixed time t=0 with varying x. Finally, t is time and e_d is the vector indicating the gravity acceleration direction (all symbols and the respective units are given in Table 1). The problem is furthermore supplemented with initial porosity ${\mathit{\varphi }}_{\mathrm{0}}:\mathrm{\Omega }\to (\mathrm{0},\mathrm{1})$ and initial effective pressure $p_{\mathrm{0}}:\mathrm{\Omega }\to \mathbb{R}$.

Table 1 Variables and physical quantities.

An extension of the hydro-mechanical model from Eqs. (1a) and (1b) is to consider the transport of a chemical tracer (such as a trace element) as described in Jordan et al. (2018, Sect. 3). In particular, the trace-element transport equations are chosen since we consider that the abundance of trace elements does not affect the mechanical or the hydraulic properties of the rock. As a consequence, the trace-element transport problem depends on the hydromechanical problem, but the opposite is not true. This allows us to treat the chemical transport problem as a post-processing step after we have calculated the respective fluid velocities and the porosity distribution. The amount of tracer is quantified using the total concentration

$$\begin{array}{}\text{(2)}& \mathcal{C}=\mathit{\varphi }{\mathit{\rho }}^{\mathrm{f}}{\mathit{\chi }}^{\mathrm{f}}+(\mathrm{1}-\mathit{\varphi }){\mathit{\rho }}^{\mathrm{s}}{\mathit{\chi }}^{\mathrm{s}},\end{array}$$

and fulfills the transport equation

$$\begin{array}{}\text{(3)}& {\partial }_t\mathcal{C}+{\mathrm{div}}_x\left({\mathit{v}}^{\mathrm{e}}\mathcal{C}\right)=\mathrm{0}.\end{array}$$

Here v^e denotes the effective velocity and at the limit where v^s≈0 holds, is given by

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\mathit{v}}^{\mathrm{e}}={\displaystyle \frac{{\mathit{v}}^{\mathrm{f}}\mathit{\varphi }}{\mathit{\varphi }+(\mathrm{1}-\mathit{\varphi })K_D}},\\ & {\mathit{v}}^{\mathrm{f}}={\displaystyle \frac{\mathrm{1}}{\mathit{\varphi }}}{\displaystyle \frac{k_{\mathrm{b}}}{\mathit{\mu }{\mathit{\varphi }}_{\mathrm{b}}^n}}{\mathit{\varphi }}^n\left({\mathrm{\nabla }}_{x}p+(\mathrm{1}-\mathit{\varphi })\mathit{\delta }\mathit{\rho }g{\mathit{e}}_{\mathrm{d}}\right),\end{array}\end{array}$$

where $K_D=\frac{{\mathit{\rho }}^{\mathrm{s}}{\mathit{\chi }}^{\mathrm{s}}}{{\mathit{\rho }}^{\mathrm{f}}{\mathit{\chi }}^{\mathrm{f}}}$ describes the concentration ratio of the tracer which is assumed to be constant as already indicated in Table 1. Note that, in this case, K_D is a ratio of concentrations and not of mass fractions. Furthermore, Eq. (3) assumes that porosity is continuous and its derivation can be found in Appendix A. For cases where porosity is discontinuous, the jump condition must guarantee the conservation of mass at the discontinuity (see Appendix B for details).

1.2 Applicability of assumptions

The previous hydro-mechanical system from Eqs. (1a) and (1b) results from the simplification of the multiphase-, viscoelastic-Stokes' equations at the static limit. The static limit occurs when no far-field stresses are imposed at the boundaries, and the buoyancy stresses are relatively small within the model domain. This limit is justified in cases where the effective pressure is close to zero. In such cases, the shear stresses that rocks can support are very small and, in many applications, can be assumed to be negligible (Aharonov et al., 1997; Connolly and Podladchikov, 1998, 2007; Scott and Stevenson, 1984). Being close to the static limit implies that the solid velocity for the mechanical problem is taken at the limit where ${\mathit{v}}^{\mathrm{f}}\gg {\mathit{v}}^{\mathrm{s}}\approx \mathrm{0}$ (but generally div_xv^s≠0).

One particular process where the trace-element transport is important is when melt is ascending within the Earth's mantle (Richter, 1986). In regions such as in the mantle wedge or within a mantle plume, the temperature does not change significantly. In these geodynamic environments, the confining pressure is large and the melt-filled porosity of the mantle rock is very small, typically in the order of 0.001–0.01 (Sims et al., 1999). To model the trace-element equilibrium and the chemical interaction between solid and fluid, we use the partition coefficient K_D. The partition coefficient changes as a function of the mineralogy of the rock, its pressure and its temperature. However, for a given material, the variation of the partition coefficient with pressure is very weak and can be considered constant over several GPa of pressure (Taura et al., 1998).

Having v^f from Eqs. (1a) and (1b) allows the solution of Eq. (3). The specific form of Eq. (3) is valid at the limit where grain-scale chemical diffusion and hydro-dynamic dispersion are ignored. Previous studies indicate that, on the large scales considered here, these phenomena can be neglected (Richter, 1986; Stavropoulou et al., 1998).

1.3 Existing numerical methods

Various methods have been proposed to solve the hydro-mechanical problem from Eqs. (1a) and (1b) numerically, for example finite difference schemes with implicit time-stepping in Connolly and Podladchikov (1998) and adaptive wavelets in Vasilyev et al. (1998). In a number of recent works, pseudo-transient schemes based on explicit time stepping in a pseudo-time variable have been investigated. Due to their compact stencils, low communication overhead and simple implementation, such schemes are well suited for parallel computing on GPUs, so that very high grid resolutions can be achieved to compensate the low order of convergence, as shown for example in Räss et al. (2014, 2018, 2019), Utkin and Afanasyev (2021), Yarushina et al. (2020) and Reuber et al. (2020). Even though all of these schemes are observed to work well for smooth initial porosities ϕ₀, their convergence can be very slow in problems with non-smooth ϕ₀, in particular in the presence of discontinuities (Bachmayr and Boisserée, 2025). Examples of this behavior are also shown in Appendix C. In such cases, due to the smoothing that is implicit in the finite difference schemes, accurately resolving sharp localized features can require extremely large grids that are computationally inefficient.

1.4 Novel contributions

Our approach considers the utilization of a space-time method to solve the hydro-mechanical problem from Eqs. (1a) and (1b). This method was introduced for this particular problem in Bachmayr and Boisserée (2025), and has the advantage that the entire solutions of porosity and effective pressure fields can be stored in a space-time grid. In addition, this approach can handle discontinuities in the porosity ϕ without approximating it by a continuous function with steep gradient. As a result, smaller grids and less computational effort are needed compared to continuous schemes such as finite differences. Since the method generates efficient approximations of the entire time history of a solution to Eqs. (1a) and (1b) in a sparse format, its coupling to the problem of chemical-tracer transport (CT) given by Eq. (3) becomes straightforward. This is because the CT problem does not give feedback to the model from Eqs. (1a) and (1b), and thus solving it can be seen as post-processing step.

1.5 Outline

Since our results from the HM model in Eqs. (1a) and (1b) are uncoupled to the results of the CT problem in Eq. (3), we begin with a short description of the methods used to solve the HM model. In Sect. 2 we introduce the methods to obtain the numerical results both for the HM model in Sect. 3 and for the CT problem in Sect. 4. We finish with a discussion regarding the implication of our results for the porous fluid transport in natural systems.

2 Methods

2.1 Hydro-mechanical model (HM)

To solve Eqs. (1a) and (1b) we consider the space-time adaptive method which was introduced in Bachmayr and Boisserée (2025) based on a combination of a Picard iterations for Eq. (1a) and a particular adaptive least squares discretization of Eq. (1b) which itself is based on Führer and Karkulik (2021) and Gantner and Stevenson (2021, 2024). To make this more precise, we start by introducing the new variable

$$\begin{array}{}\text{(5)}& \mathit{\phi }=-\log (\mathrm{1}-\mathit{\varphi }),\end{array}$$

so that $\mathit{\varphi }=\mathrm{1}-e^{-\mathit{\phi }}$. The previous transformation allows the investigation of cases where the porosity is larger than the typical “small-porosity limit” (Vasilyev et al., 1998). The system in Eqs. (1a) and (1b) can then be written in the form

$$\begin{array}{}\text{(6a)}& {\displaystyle }\begin{array}{rl}& {\partial }_t\mathit{\phi }=-\left(\mathit{\beta }\left(\mathit{\phi }\right){\displaystyle \frac{p}{\mathit{\sigma }\left(p\right)}}+Q{\partial }_{t}p\right),\mathit{\phi }(\mathrm{0},\cdot )={\mathit{\phi }}_{\mathrm{0}},\end{array}\text{(6b)}& {\displaystyle }\begin{array}{rl}{\partial }_{t}p=& {\displaystyle \frac{\mathrm{1}}{Q}}\left({\mathrm{div}}_x\left(\mathit{\alpha }\left(\mathit{\phi }\right)\left({\mathrm{\nabla }}_{x}p+\mathit{\zeta }\left(\mathit{\phi }\right)\right)\right)-\mathit{\beta }\left(\mathit{\phi }\right){\displaystyle \frac{p}{\mathit{\sigma }\left(p\right)}}\right),\\ & p(\mathrm{0},\cdot )=p_{\mathrm{0}},\end{array}\end{array}$$

where

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& \mathit{\alpha }\left(\mathit{\phi }\right)={\displaystyle \frac{k_{\mathrm{b}}}{\mathit{\mu }{\mathit{\varphi }}_{\mathrm{b}}^n}}{\left(\mathrm{1}-e^{-\mathit{\phi }}\right)}^n,\mathit{\beta }\left(\mathit{\phi }\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\eta }}_{\mathrm{b}}}}{\left(\mathrm{1}-e^{-\mathit{\phi }}\right)}^m,\\ & \mathit{\zeta }\left(\mathit{\phi }\right)=e^{-\mathit{\phi }}\mathit{\delta }\mathit{\rho }g{\mathit{e}}_{\mathrm{d}}.\end{array}\end{array}$$

To solve Eq. (6b) for a fixed $\stackrel{\mathrm{‾}}{\mathit{\phi }}$ we consider a linearization, that is, we solve

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{\partial }_t{p}^{\left(k\right)}& ={\displaystyle \frac{\mathrm{1}}{Q}}({\mathrm{div}}_x\left(\mathit{\alpha }\left(\stackrel{\mathrm{‾}}{\mathit{\phi }}\right)\left({\mathrm{\nabla }}_x{p}^{\left(k\right)}+\mathit{\zeta }\left(\stackrel{\mathrm{‾}}{\mathit{\phi }}\right)\right)\right)\\ & -\mathit{\beta }\left(\stackrel{\mathrm{‾}}{\mathit{\phi }}\right){\displaystyle \frac{p^{\left(k\right)}}{\mathit{\sigma }\left(p^{(k-\mathrm{1})}\right)}}),p^{\left(k\right)}(\mathrm{0},\cdot )=p_{\mathrm{0}},\end{array}\end{array}$$

for p^(k) given the previous iterate p^(k−1). By defining

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& G\left[p^{(k-\mathrm{1})}\right]\left(p^{\left(k\right)},{\mathit{\psi }}^{\left(k\right)}\right)\\ & =\left(\begin{array}{c}\frac{\mathrm{1}}{Q}\left(\mathrm{div}\left(p^{\left(k\right)},{\mathit{\psi }}^{\left(k\right)}\right)+\mathit{\beta }\left(\stackrel{\mathrm{‾}}{\mathit{\phi }}\right)\frac{p^{\left(k\right)}}{\mathit{\sigma }\left(p^{(k-\mathrm{1})}\right)}\right)\\ {\mathit{\psi }}^{\left(k\right)}+\mathit{\alpha }\left(\stackrel{\mathrm{‾}}{\mathit{\phi }}\right){\mathrm{\nabla }}_x{p}^{\left(k\right)}\\ p^{\left(k\right)}(\mathrm{0},\cdot )\end{array}\right),\\ & R=\left(\begin{array}{c}\mathrm{0}\\ -\mathit{\alpha }\left(\stackrel{\mathrm{‾}}{\mathit{\phi }}\right)\mathit{\zeta }\left(\stackrel{\mathrm{‾}}{\mathit{\phi }}\right)\\ p_{\mathrm{0}}\end{array}\right),\end{array}\end{array}$$

with $\mathrm{div}(p,\mathit{\psi })={\partial }_{t}p+{\mathrm{div}}_x\mathit{\psi }$, we can reformulate Eq. (8) as first-order system

$$\begin{array}{}\text{(10)}& G\left[p^{(k-\mathrm{1})}\right]\left(p^{\left(k\right)},{\mathit{\psi }}^{\left(k\right)}\right)=R.\end{array}$$

Note that the second row of Eq. (9) can be rewritten as ${\mathit{\psi }}^{\left(k\right)}=-\mathit{\alpha }\left(\stackrel{\mathrm{‾}}{\mathit{\phi }}\right)({\mathrm{\nabla }}_x{p}^{\left(k\right)}+\mathit{\zeta }(\stackrel{\mathrm{‾}}{\mathit{\phi }}\left)\right)$. Plugging ψ^(k) into div(p^(k),ψ^(k)) in the first row of Eq. (9) yields Eq. (8).

Numerically we now use the approach presented in Führer and Karkulik (2021) and Gantner and Stevenson (2021, 2024), and thus calculate a least squares minimizer with respect to an appropriately chosen norm. Using the numerical approximation of $p\left[\stackrel{\mathrm{‾}}{\mathit{\phi }}\right]$ from Eq. (10) we solve Eq. (6a) by discretizing the iteration

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{\mathit{\phi }}^{(k+\mathrm{1})}(t,\cdot )& =Q\left(p_{\mathrm{0}}-p\left[{\mathit{\phi }}^{\left(k\right)}\right](t,\cdot )\right)\\ & -\underset{\mathrm{0}}{\overset{t}{\int }}\mathit{\beta }\left({\mathit{\phi }}^{\left(k\right)}\right){\displaystyle \frac{p\left[{\mathit{\phi }}^{\left(k\right)}\right](s,\cdot )}{\mathit{\sigma }\left(p\left[{\mathit{\phi }}^{\left(k\right)}\right](s,\cdot )\right)}}\mathrm{d}s.\end{array}\end{array}$$

which is based on integrating Eq. (6a) in time. For more details including proofs of convergence we refer to Bachmayr and Boisserée (2025, Sects. 3 and 4).

Figure 1 Example of a porosity channel at T=1.5 Myr (a) with the associated adaptive space-time grid (b); the color of each grid cell denotes its refinement level.

This scheme can generate space-time grids corresponding to spatially adapted time steps; an example of this can be found in Fig. 1. Furthermore, the method provides computable a-posteriori estimates of the error with respect to the exact solution of the coupled nonlinear system from Eqs. (6a) and (6b). Therefore, this method can be used to steer an adaptive grid refinement routine which yields efficient approximations of localized features of solutions, in particular in the presence of discontinuities. In addition, one obtains optimal convergence rates for ϕ and p independent of the presence of discontinuities in ϕ, as observed in Bachmayr and Boisserée (2025, Sect. 5.2).

2.2 Chemical-tracer transport model (CT)

To solve the chemical transport in Eq. (3), we follow its characteristics. Namely, we consider

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& {\partial }_{t}x\left(t\right)={\mathit{v}}^{\mathrm{e}}(t,x(t\left)\right),\\ & c\left(t\right)=\mathcal{C}(t,x(t\left)\right).\end{array}\end{array}$$

Then we calculate

$$\begin{array}{}\text{(13)}& \begin{array}{rl}{\partial }_{t}c\left(t\right)& ={\partial }_t\mathcal{C}(t,x(t\left)\right)+{\mathrm{\nabla }}_x\mathcal{C}(t,x(t\left)\right)\cdot {\partial }_{t}x\left(t\right)\\ & =-{\mathrm{div}}_x\left({\mathit{v}}^{\mathrm{e}}(t,x(t\left)\right)\mathcal{C}(t,x(t\left)\right)\right)\\ & +{\mathrm{\nabla }}_x\mathcal{C}(t,x(t\left)\right)\cdot {\mathit{v}}^{\mathrm{e}}(t,x(t\left)\right)\\ & =-{\mathrm{div}}_x\left({\mathit{v}}^{\mathrm{e}}(t,x(t\left)\right)\right)\mathcal{C}(t,x(t\left)\right)\\ & =-{\mathrm{div}}_x\left({\mathit{v}}^{\mathrm{e}}(t,x(t\left)\right)\right)c\left(t\right),\end{array}\end{array}$$

the previous yields a coupled system of ordinary differential equations (ODEs)

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& {\partial }_{t}x\left(t\right)={\mathit{v}}^{\mathrm{e}}(t,x(t\left)\right),\\ & {\partial }_{t}c\left(t\right)=-{\mathrm{div}}_x\left({\mathit{v}}^{\mathrm{e}}(t,x(t\left)\right)\right)c\left(t\right),\end{array}\end{array}$$

that is used to solve Eq. (3). The solution is provided along the characteristics given by v^e starting with some initial value x₀ and initial concentration c₀. We use an explicit Euler scheme to solve Eq. (14) for many different starting values x₀. Note that this approach is highly parallelizable, since we need to solve a high number of independent ODEs for each starting value. By exploiting this we usually achieve very low wall-clock times, even for many starting values corresponding to a high resolution. Note furthermore that this approach only conserves the quantity 𝒞 if v^e is continuous. For the discontinuous cases one needs to ensure continuity of the flux and we refer to Sect. B for more details.

2.3 Model parameters

The model parameters can be derived by non-dimensionalizing the physical models from Eqs. (1a), (1b) and (3) with values given in Table 1. Choosing the independent scales

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& x^{\mathrm{sc}}={\mathrm{10}}^{\mathrm{4}}\mathrm{m},\mathit{\delta }{\mathit{\rho }}^{\mathrm{sc}}{g}^{\mathrm{sc}}=\mathrm{5}\times {\mathrm{0}}^{\mathrm{3}}\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{s}}^{-\mathrm{2}},\\ & {\mathit{\eta }}_{\mathrm{b}}^{\mathrm{sc}}={\mathrm{10}}^{\mathrm{19}}\mathrm{Pa}\mathrm{s},\end{array}\end{array}$$

yields the dependent scales

$$\begin{array}{}\text{(16)}& \begin{array}{rl}& p^{\mathrm{sc}}=\mathit{\delta }{\mathit{\rho }}^{\mathrm{sc}}{g}^{\mathrm{sc}}{x}^{\mathrm{sc}}=\mathrm{5}\times {\mathrm{10}}^{\mathrm{7}}\mathrm{Pa},\\ & t^{\mathrm{sc}}={\displaystyle \frac{{\mathit{\eta }}_{\mathrm{b}}^{\mathrm{sc}}}{p^{\mathrm{sc}}}}=\mathrm{2}\times {\mathrm{10}}^{\mathrm{11}}\mathrm{s},\\ & {\displaystyle \frac{k_{\mathrm{b}}^{\mathrm{sc}}}{{\mathit{\mu }}^{\mathrm{sc}}}}={\displaystyle \frac{{\left(x^{\mathrm{sc}}\right)}^{\mathrm{2}}}{{\mathit{\eta }}_{\mathrm{b}}^{\mathrm{sc}}}}={\mathrm{10}}^{-\mathrm{11}}{\mathrm{m}}^{\mathrm{2}}{\mathrm{Pa}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}.\end{array}\end{array}$$

Hence we end up with the nondimensional parameters

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& \widehat{{\mathit{\eta }}_{\mathrm{b}}}=\mathrm{1},{\displaystyle \frac{\widehat{k_{\mathrm{b}}}}{\widehat{\mathit{\mu }}{\mathit{\varphi }}_{\mathrm{b}}^n}}=\mathrm{1000},\widehat{\mathit{\delta }\mathit{\rho }}\widehat{g}{\mathit{e}}_{\mathrm{d}}=\left(\begin{array}{c}\mathrm{0}\\ \mathrm{1}\end{array}\right),\\ & B\widehat{Q}=Q\mathit{\delta }{\mathit{\rho }}^{\mathrm{sc}}{g}^{\mathrm{sc}}{x}^{\mathrm{sc}}={\displaystyle \frac{\mathrm{1}}{\mathrm{60}}},\end{array}\end{array}$$

where B=Q^scδρ^scg^scx^sc denotes a non-dimensional number that is the ratio of buoyancy stress to bulk modulus. Note that this may look similar to the Deborah number defined in Connolly and Podladchikov (1998), however, the lengthscale x^sc is taken as an independent quantity in our approach, while in other studies it is derived from the compaction length (Connolly and Podladchikov, 1998; Vasilyev et al., 1998).

Due to the given length scale and time scale, we can directly translate physical times and domain sizes into our model parameters if we divide by x^sc or t^sc, respectively. For T=1.5 Myr, this corresponds to

$$\widehat{T}=T/t^{\mathrm{sc}}={\displaystyle \frac{\mathrm{1.5}\times {\mathrm{10}}^{\mathrm{6}}\times \mathrm{365.25}\times \mathrm{24}\times \mathrm{60}\times \mathrm{60}}{\mathrm{2}\times {\mathrm{10}}^{\mathrm{11}}}}=\mathrm{236.682}.$$

Note that, in the following, the “$\widehat{}$” symbols are omitted for convenience. Furthermore, we consider σ, as suggested in Räss et al. (2018, 2019), which is an expression of the form

$$\begin{array}{}\text{(18)}& \begin{array}{rl}\mathit{\sigma }\left(v\right)& =\mathrm{1}-{\displaystyle \frac{\mathrm{1}-c_{\mathrm{1}}}{\mathrm{2}}}\left(\mathrm{1}+\tanh \left(-{\displaystyle \frac{v}{c_{\mathrm{2}}}}\right)\right)\\ & ={\displaystyle \frac{c_{\mathrm{1}}+\exp \left(\mathrm{2}v/c_{\mathrm{2}}\right)}{\mathrm{1}+\exp \left(\mathrm{2}v/c_{\mathrm{2}}\right)}},v\in \mathbb{R},\end{array}\end{array}$$

and provides a phenomenological model for decompaction weakening. Here $c_{\mathrm{1}}\in (\mathrm{0},\mathrm{1}]$ and c₂>0, where 1+tanh can be regarded as a smooth approximation of a step function taking values in the interval (0, 2). In the most well-studied case c₁=1, as considered in Vasilyev et al. (1998), one observes the formation of porosity waves, whereas the case of c₁<1 with appropriate problem parameters and initial conditions, one can observe the formation of channels. With parameters from the stated ranges, σ as in Eq. (18) fulfills (Bachmayr et al., 2023, Assumptions 1), and hence we refer to Bachmayr et al. (2023, Sect. 4) for the well-posedness of Eqs. (1a) and (1b).

In this work we consider two cases, namely c₁, c₂=1 (no decompaction weakening) and c₁, c₂=0.002 (including decompaction weakening). The resulting functions read

$$\begin{array}{}\text{(19)}& {\mathit{\sigma }}^{\mathrm{a}}\left(v\right)=\mathrm{1},{\mathit{\sigma }}^{\mathrm{b}}\left(v\right)={\displaystyle \frac{\mathrm{0.002}+\exp \left(\mathrm{1000}v\right)}{\mathrm{1}+\exp \left(\mathrm{1000}v\right)}},v\in \mathbb{R}.\end{array}$$

3 Hydro-mechanical model results

In this part we show numerical results for solving the hydro-mechanical problem from Eqs. (1a) and (1b) using the method described in Sect. 2.1. We compare three different initial porosities which are shown in Fig. 2, reflecting our main focus on investigating the cases of jump discontinuities in the initial porosity distribution. Taking an initial homogeneous porosity as a reference case (Fig. 2a), we consider a drop (along x₂; Fig. 2b) versus an increase (along x₂; Fig. 2c) in the initial porosity distribution. For the effective pressure p, we assume homogeneous initial data p₀(x)=0. Note that all calculations in this section are carried out on larger grids, namely grids of size 30 km in x₂ direction, to avoid problems near the upper boundary. This, however, does not affect the computation time significantly since the additional grid cells located above 20 km are generally not refined because ϕ and p are almost constant there.

Figure 2 Three initial porosity distributions ${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{a}}$ (a), ${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{b}}$ (b) and ${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{c}}$ (c).

We start with the well-known scenario of the smooth initial porosity ${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{a}}$ as in Fig. 2a and compare the results without decompaction weakening (σ^a) and with decompaction weakening (σ^b) in Fig. 3. The resulting plots show the expected spreading of the fluid front in the case without decompaction weakening in Fig. 3a and c. In contrast, the fluid flow is focused in the presence of weakening as one can see in Fig. 3b and d. These results are used as reference and will not be discussed further since they confirm previous findings (Räss et al., 2018; Connolly and Podladchikov, 1998, 2007; Yarushina et al., 2015).

Figure 3 Porosity (a, b) and effective pressure (c, d) after T=1.5 Myr for a smooth initial condition (${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{a}}$) without decompaction weakening (σ^a) (a, c) and with decompaction weakening (σ^b) (b, d).

Figure 4 shows the results of the two initial conditions ${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{b}}$ and ${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{c}}$ from Fig. 2b and c that consider an initial porosity discontinuity.

Figure 4 Porosity (a, b) and effective pressure (c, d) without decompaction weakening (σ^a) after T=1.5 Myr for ${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{b}}$ (a, c) and ${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{c}}$ (b, d).

Both results consider the case without decompaction weakening. One can see very sharp transitions of ϕ at the locations of the initial discontinuities. Note that the discontinuity itself cannot move since the model from Eqs. (1a) and (1b) was derived under the assumption that v^s≈0 and porosity is a property of the solid. This agrees with the theoretical results shown in Bachmayr et al. (2023, Thm. 4.6). As it is especially visible for the porosity distribution, the sign of its transition (positive or negative) depends on whether the initial porosity of the upper layer was smaller (negative jump) or larger (positive jump) compared to the porosity of the underlying layer. Figure 5 shows a cross section of Fig. 4 that shows the discontinuities in ϕ more clearly. In contrast, the solution for p is continuous which aligns with the theory derived in Bachmayr et al. (2023, Sect. 4).

Figure 5 Cross section of Fig.  4 for x₁=5 km with porosity (a, b) and effective pressure (c, d).

Figure 6 shows the effect of decompaction weakening on the same initial discontinuous configurations (case σ^b). In the case of the negative jump (${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{b}}$), the channel has a slightly higher maximal porosity compared to the continuous case. Furthermore, Fig. 6a shows that there is a very steep increase in porosity at the place of the discontinuity. On the other hand, for the positive jump (${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{c}}$), the channel focuses significantly before spreading in the high-porosity/permeability zone as it is shown in Fig. 6b.

Figure 6 Porosity (a, b) and effective pressure (c, d) with decompaction weakening (σ^b) after T=1.5 Myr for ${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{b}}$ (a, c) and ${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{c}}$ (b, d).

This can be expected: we see a narrower channel in the domain that has smaller porosity, but the channel spreads quickly once the fluid enters the domain of high porosity and permeability. This shows that the fluid does not need to channelize as much as in the less porous domain in order to travel upwards.

4 Chemical-transport model results

An extension of the general model from Eqs. (1a) and (1b) is to consider the transport of a chemical tracer by solving Eq. (3). This is achieved by following its characteristics as described in Sect. 2.2. The natural range of partition coefficients can be very large (Irving, 1978) and their magnitudes depend on several parameters (Karato, 2016). However, as it was described in Sect. 1.2, it is reasonable to assume that the partition coefficient can be approximated as constant given a limited range of temperatures and constant mineralogical composition (that is implicitly assumed in our models). Without loss of generality, we examine the case where $K_D={\mathrm{10}}^{-\mathrm{3}}$ (as already indicated in Table 1) to consider incompatible elements. Incompatible elements are those that partition preferentially in the fluid. Solving for 𝒞 and using the prescribed values for ρ^s, ρ^f and K_D from Table 1 directly yields χ^s and χ^f as well.

In this part, we will only plot the normalized chemical tracer $\mathcal{C}/{\mathcal{C}}_{\mathrm{0}}$. This allows us to quantify the overall enrichment or depletion of a trace element with respect to the initial configuration. Note that since we use a characteristics-based approach for the advection of chemical elements, the chemical evolution is calculated only for the areas where the characteristics are initialized. As a result, the domain where the chemical evolution is calculated, changes in time depending on the effective velocity v^e. For simplicity, we consider constant initial data 𝒞₀(x)=1 since 𝒞 in Eq. (3) can be scaled arbitrarily without affecting the solution.

Figure 7 shows the normalized tracer compositions $\mathcal{C}/{\mathcal{C}}_{\mathrm{0}}$ connected to the solutions shown in Fig. 3 (with ${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{a}}$) for $K_D={\mathrm{10}}^{-\mathrm{3}}$.

Figure 7 $\mathcal{C}/{\mathcal{C}}_{\mathrm{0}}$ after T=1.5 Myr with an initially continuous porosity (${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{a}}$) without decompaction weakening (σ^a) (a) and with decompaction weakening (σ^b) (b).

We see the distribution of an incompatible element that prefers to stay with the fluid, and hence, it gets transported efficiently while draining the area of origin. The role of decompaction weakening becomes more apparent in the case of the channelization of the fluid flow, as shown in Fig. 7b. In that case, we observe a more pronounced enrichment in the region defined by the fluid-rich channel. It is important to note that this enrichment occurs in both the solid and the fluid, and it occurs at the expense of the trace element's distribution in the source region.

The solution of the CT problem having initial discontinuous porosity ϕ₀ is shown in Fig. 8.

Figure 8 $\mathcal{C}/{\mathcal{C}}_{\mathrm{0}}$ without decompaction weakening (σ^a) after T=1.5 Myr for ${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{b}}$ (a) and ${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{c}}$ (b).

This figure is calculated based on the HM model (porosity-pressure evolution) shown in Fig. 4 and assumes no decompaction weakening (σ^a). The results generally agree with the previous findings that show that the incompatible elements ($K_D={\mathrm{10}}^{-\mathrm{3}}$) travel further and enrich the upper layer. Furthermore, this enrichment seems to be traveling slightly faster in the region where the initial fluid content was higher (central region of the domain). This requires that the propagation velocity of the enrichment front is not constant and moves further from the location of the porosity discontinuity, which is located exactly at the middle of the domain (x₂=10 km). In contrast, for the case of the negative initial-porosity jump, the enrichment is negligible and is located in the area just above the discontinuity. A marked feature of the discontinuous models shown in Fig. 8 is that, exactly at the discontinuity, we observe a significant enrichment or depletion of 𝒞 depending whether we have a drop (${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{b}}$) or an increase (${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{c}}$) in the initial porosity.

Figure 9 $\mathcal{C}/{\mathcal{C}}_{\mathrm{0}}$ with decompaction weakening (σ^b) after T=1.5 Myr for ${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{b}}$ (a) and ${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{c}}$ (b).

Figure 10 Cross section of Fig. 9 for x₁=5 km with $\mathcal{C}/{\mathcal{C}}_{\mathrm{0}}$ on the y axis.

Finally, Figs. 9 and 10 show the resulting normalized tracer element $\mathcal{C}/{\mathcal{C}}_{\mathrm{0}}$ for the case of decompaction weakening (σ^b) and an initially discontinuous porosity ϕ₀. The associated HM model can be found in Fig. 6. The resulting cases show marked differences and can be summarized as follows. The case with negative jump discontinuity (${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{b}}$) shows a marked enrichment with respect to the incompatible element. In particular, there is a marked enrichment at the discontinuity (at x₂=10 km), and within the channel in general. Interestingly, for the case of positive jump discontinuity (${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{c}}$), the enrichment of the incompatible element is localized close to the discontinuity location (but is smaller at the discontinuity itself). This is explained by the fact that the fluid spreads beyond this point as it was shown in Fig. 6b and d.

5 Discussion and conclusions

We have presented results for the case of compaction-driven fluid flow in relation to fluid migration in the deep subsurface. Our method aims to resolve the effects of discontinuous porosity distributions as already discussed in Bachmayr and Boisserée (2025). The models confirm previous findings for the cases of homogeneous initial porosity (ϕ₀) distribution (Räss et al., 2018; Connolly and Podladchikov, 1998, 2007; Yarushina et al., 2015). However, for the cases when the ϕ₀ has jump discontinuities, our method predicts discontinuous solutions without artificial smoothing due to numerical diffusion. Such results are useful for cases where the mechanical variables, such as the effective and fluid pressure, need to be quantified in applications (Räss et al., 2018; Yarushina et al., 2022), and thus, our approach can be used to provide a reference case for numerical benchmarks.

An additional advantage of the space-time method is that the one-way coupling of the HM problem to the CT problem can be easily solved using the pre-calculated results of the HM problem. This allows for the investigation of the behavior of various trace elements and the overall mass transport in rock formations that have discontinuous porosity. Our results confirm previous data which suggest that incompatible elements are the most mobile and can travel together with the fluid (Richter, 1986). This selective enrichment in incompatible elements becomes more prominent in cases where the flow is channelized, leading to the formation of localized geochemical and mineralogical anomalies. Although channeling mechanisms have been discussed in previous works (Aharonov et al., 1997; Spiegelman and Kelemen, 2003; Schiemenz et al., 2011), the mechanism for the channeling in our case is different. In the aforementioned studies, the formation of channels was due to the selective dissolution of matrix minerals (Schiemenz et al., 2011; Spiegelman et al., 2001). In contrast, in our case the channeling is the result of decompaction weakening (Connolly and Podladchikov, 2007; Yarushina et al., 2015, 2020). In any case, it becomes apparent that, whatever the localization mechanism may be, the localization of the fluid amplify the enrichment of incompatible elements significantly. Furthermore, our new results also show the interaction of a fluid-filled channel with a jump discontinuity in the initial porosity. This example is very relevant for the case of fluid transport across heterogeneous layers. In particular, the results show a marked enrichment of the incompatible trace elements at the initial porosity discontinuity for the cases where the initial porosity exhibits a negative jump (i.e. porosity drops sharply at the transition). In the case of a positive jump, we observe a marked depletion at exactly the same location. The results indicate that both porosity and the incompatible-element enrichment, that are associated to the discontinuity, do not move over time and remain at the same location.

Our results indicate that the effects of the channeling of the flow together with the presence of initial discontinuities will produce a variety of element-enrichment patterns that can be investigated in future studies that focus on particular element behavior. These results can be very important in targeted mineral exploration and in the understanding of ore-formation processes.

Appendix A: Derivation of the chemical model

Starting with the conservation of mass in the two phases, we get

$$\begin{array}{}\text{(A1a)}& {\displaystyle }& {\displaystyle }{\partial }_t\left(\mathit{\varphi }{\mathit{\rho }}^{\mathrm{f}}{\mathit{\chi }}_i^{\mathrm{f}}\right)+{\mathrm{div}}_x\left(\mathit{\varphi }{\mathit{\rho }}^{\mathrm{f}}{\mathit{\chi }}_i^{\mathrm{f}}{\mathit{v}}_i^{\mathrm{f}}\right)={\mathrm{\Gamma }}_i^{\mathrm{f}},\text{(A1b)}& {\displaystyle }& {\displaystyle }{\partial }_t\left((\mathrm{1}-\mathit{\varphi }){\mathit{\rho }}^{\mathrm{s}}{\mathit{\chi }}_i^{\mathrm{s}}\right)+{\mathrm{div}}_x\left((\mathrm{1}-\mathit{\varphi }){\mathit{\rho }}^{\mathrm{s}}{\mathit{\chi }}_i^{\mathrm{s}}{\mathit{v}}_i^{\mathrm{s}}\right)={\mathrm{\Gamma }}_i^{\mathrm{s}},\end{array}$$

for i=1, …, n chemical elements, where ${\mathrm{\Gamma }}_i^{\mathrm{f}}$, ${\mathrm{\Gamma }}_i^{\mathrm{s}}$ denote reaction terms. Note that $\sum _{i=\mathrm{1}}^n{\mathit{\chi }}^{\mathrm{f}}=\sum _{i=\mathrm{1}}^n{\mathit{\chi }}^{\mathrm{s}}=\mathrm{1}$ and ${\mathrm{\Gamma }}_i^{\mathrm{f}}+{\mathrm{\Gamma }}_i^{\mathrm{s}}=\mathrm{0}$ for all i=1, …, n, since chemical elements can only be exchanged between the two phases. Defining the barycentric velocities

$$\begin{array}{}\text{(A2)}& {\mathit{v}}^{\mathrm{f}}=\sum _{i=\mathrm{1}}^n{\mathit{\chi }}_i^{\mathrm{f}}{\mathit{v}}_i^{\mathrm{f}},{\mathit{v}}^{\mathrm{s}}=\sum _{i=\mathrm{1}}^n{\mathit{\chi }}_i^{\mathrm{s}}{\mathit{v}}_i^{\mathrm{s}},\end{array}$$

we rewrite Eqs. (A1a) and (A1b)

$$\begin{array}{}\text{(A3a)}& {\displaystyle }\begin{array}{rl}{\partial }_t\left(\mathit{\varphi }{\mathit{\rho }}^{\mathrm{f}}{\mathit{\chi }}_i^{\mathrm{f}}\right)& +{\mathrm{div}}_x\left(\mathit{\varphi }{\mathit{\rho }}^{\mathrm{f}}{\mathit{\chi }}_i^{\mathrm{f}}{\mathit{v}}^{\mathrm{f}}\right)\\ & +{\mathrm{div}}_x\left(\mathit{\varphi }{\mathit{\rho }}^{\mathrm{f}}{\mathit{\chi }}_i^{\mathrm{f}}\left({\mathit{v}}_i^{\mathrm{f}}-{\mathit{v}}^{\mathrm{f}}\right)\right)={\mathrm{\Gamma }}_i^{\mathrm{f}},\end{array}\text{(A3b)}& {\displaystyle }\begin{array}{rl}{\partial }_t\left((\mathrm{1}-\mathit{\varphi }){\mathit{\rho }}^{\mathrm{s}}{\mathit{\chi }}_i^{\mathrm{s}}\right)& +{\mathrm{div}}_x\left((\mathrm{1}-\mathit{\varphi }){\mathit{\rho }}^{\mathrm{s}}{\mathit{\chi }}_i^{\mathrm{s}}{\mathit{v}}^{\mathrm{s}}\right)\\ & +{\mathrm{div}}_x\left((\mathrm{1}-\mathit{\varphi }){\mathit{\rho }}^{\mathrm{s}}{\mathit{\chi }}_i^{\mathrm{s}}\left({\mathit{v}}_i^{\mathrm{s}}-{\mathit{v}}^{\mathrm{s}}\right)\right)={\mathrm{\Gamma }}_i^{\mathrm{s}},\end{array}\end{array}$$

and use that for trace elements the diffusion fluxes obey the Fickean limit

$$\begin{array}{}\text{(A4a)}& {\displaystyle }& {\displaystyle }{\mathit{\rho }}^{\mathrm{f}}{\mathit{\chi }}_i^{\mathrm{f}}\left({\mathit{v}}_i^{\mathrm{f}}-{\mathit{v}}^{\mathrm{f}}\right)=-D_i^{\mathrm{f}}{\mathrm{\nabla }}_x\left({\mathit{\rho }}^{\mathrm{f}}{\mathit{\chi }}_i^{\mathrm{f}}\right),\text{(A4b)}& {\displaystyle }& {\displaystyle }{\mathit{\rho }}^{\mathrm{s}}{\mathit{\chi }}_i^{\mathrm{s}}\left({\mathit{v}}_i^{\mathrm{s}}-{\mathit{v}}^{\mathrm{s}}\right)=-D_i^{\mathrm{s}}{\mathrm{\nabla }}_x\left({\mathit{\rho }}^{\mathrm{s}}{\mathit{\chi }}_i^{\mathrm{s}}\right)\end{array}$$

for both the fluid and solid phase. For advection dominated problems, the diffusion coefficients D_i are very small and hence we can cancel the corresponding terms in Eqs. (A3a) and (A3b). Adding the resulting equations yields

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}{\partial }_t(\mathit{\varphi }{\mathit{\rho }}^{\mathrm{f}}{\mathit{\chi }}_i^{\mathrm{f}}& +(\mathrm{1}-\mathit{\varphi }){\mathit{\rho }}^{\mathrm{s}}{\mathit{\chi }}_i^{\mathrm{s}})\\ & +{\mathrm{div}}_x(\mathit{\varphi }{\mathit{\rho }}^{\mathrm{f}}{\mathit{\chi }}_i^{\mathrm{f}}{\mathit{v}}^{\mathrm{f}}+(\mathrm{1}-\mathit{\varphi }\left){\mathit{\rho }}^{\mathrm{s}}{\mathit{\chi }}_i^{\mathrm{s}}{\mathit{v}}^{\mathrm{s}}\right)\\ & ={\mathrm{\Gamma }}_i^{\mathrm{f}}+{\mathrm{\Gamma }}_i^{\mathrm{s}}=\mathrm{0}.\end{array}\end{array}$$

Next, we define the concentration ratio $K_D^i=\frac{{\mathit{\rho }}^{\mathrm{s}}{\mathit{\chi }}_i^{\mathrm{s}}}{{\mathit{\rho }}^{\mathrm{f}}{\mathit{\chi }}_i^{\mathrm{f}}}$ as well as the total concentration ${\mathcal{C}}_i=(\mathit{\varphi }+(\mathrm{1}-\mathit{\varphi }\left)K_D^i\right){\mathit{\rho }}^{\mathrm{f}}{\mathit{\chi }}_i^{\mathrm{f}}$ which allows us to rewrite Eq. (A5) further as

$$\begin{array}{}\text{(A6)}& {\partial }_t{\mathcal{C}}_i+{\mathrm{div}}_x\left({\mathcal{C}}_i{\mathit{v}}_i^{\mathrm{e}}\right)=\mathrm{0}\end{array}$$

where ${\mathit{v}}_i^{\mathrm{e}}=\frac{\mathit{\varphi }{\mathit{v}}^{\mathrm{f}}+(\mathrm{1}-\mathit{\varphi }){\mathit{v}}^{\mathrm{s}}}{\mathit{\varphi }+(\mathrm{1}-\mathit{\varphi })K_D^i}$. Note that this equation is equivalent to Eq. (3) where we omitted the index i for convenience and assumed that v^s≈0.

Appendix B: Handling discontinuous velocities

In order to solve Eq. (14) for a discontinuous velocity field v^e, we need to ensure mass balance at the discontinuity. This is normally done via the Rankine-Hugoniot jump condition (see, for example, Anderson, 1990, Sect. 4.3 or LeVeque, 2002, Sect. 11.8), which in the case of Eq. (3) leads to

$$\begin{array}{}\text{(B1)}& c_{+}{\mathit{v}}_{+}^{\mathrm{e}}\cdot \mathit{n}=c_{-}{\mathit{v}}_{-}^{\mathrm{e}}\cdot \mathit{n}\end{array}$$

where c₊, c₋, ${\mathit{v}}_{+}^{\mathrm{e}}$, ${\mathit{v}}_{-}^{\mathrm{e}}$ denote the values of c and v^e on both sides of the discontinuity and n is the normal vector with respect to the discontinuity. In the test cases shown in Figs. 8 and 9 we have n=e₂, which simplifies the numerical calculations.

In practice, when running the explicit Euler code to solve Eq. (14), we ensure that each time step does not advect the total concentration accross the discontinuity. Once the total concentration reaches the discontinuity, we recalculate the mass flux and use it to evaluate the concentration jump (without loss of generality we call it c₋), as follows

$$\begin{array}{}\text{(B2)}& c_{+}=c_{-}{\displaystyle \frac{{\mathit{v}}_{-}^{\mathrm{e}}\cdot \mathit{n}}{{\mathit{v}}_{+}^{\mathrm{e}}\cdot \mathit{n}}}.\end{array}$$

Appendix C: Comparison with finite difference code

Here we compare our space-time approach with a classical finite difference scheme. Note that for simplicity we chose a one-dimensional test case without decompaction weakening. Hence, this is similar to the tests shown in Figs. 3a, c and 4. In Fig. C1 one can see the error of the finite difference scheme at the terminal time; note the very different rates for the discontinuous and continuous test cases. This difference is more pronounced for the porosity in Fig. C1a since the correct solution is discontinuous whereas the corresponding effective pressure is still continuous with a kink at the location of the discontinuity. Hence, the convergence rates for the pressure in Fig. C1b are higher than for the porosity, even though they are still lower than in the case of a continuous initial porosity. Note also that the rates in the continuous case are the theoretically optimal ones. In comparison, the space-time approach does not suffer from slow rates in discontinuous cases as one can see in Fig. C2. In addition, the rates here are optimal as well. However, we would like to emphasize that a direct comparison is not possible. This is because the space-time approach is fundamentally different from the finite difference one. For the finite differences, we measure the L₂(Ω)-error at the terminal time, whereas for the space-time approach we considered a more complicated space-time error norm. In addition, the number of degrees of freedom (dofs) is not directly comparable since in Fig. C1 they correspond to the number of dofs at the terminal time (even though before there were many time steps involved) and in Fig. C2 they correspond to the total number of dofs for the entire space-time grid. Finally, we note that the norms measuring the errors of the porosity in Fig. C2a and of the effective pressure in Fig. C2b are of different kind. In summary, even though the convergence rates of the two methods are not directly comparable, the new space-time approach does not suffer from reduced convergence rates in the presence of discontinuous initial porosities.

Figure C1 Errors of a finite difference approximation of the one-dimensional hydromechanical model without decompaction weakening for porosity (a) and effective pressure (b).

Figure C2 Space-time errors of the one-dimensional hydromechanical model without decompaction weakening for porosity (a) and effective pressure (b).

Appendix D: Continuous approximation

In this section we want to compare the approximation of the hydromechanical and chemical model for the discontinuous initial function ${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{b}}$ and a continuous approximation of it. In Fig. D1, we plot a cross section of both the discontinuous initial function and its smooth approximation.

Figure D1 Cross section of initial porosity ${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{b}}$ (a) and its smooth approximation (b) for x₁=5 km.

In Fig. D2a and c, one can see the hydromechanical model solution to discontinuous case (as in Fig. 6a and c) whereas in Fig. D2b and d the results were obtained using the continuous (even though very steep) initial function shown in Fig. D1b.

Figure D2 Porosity (a, b) and effective pressure (c, d) with decompaction weakening after T=1.5 Myr for ${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{b}}$ (a, c) and its continuous approximation (b, d).

Here, one can see no major difference between the two approaches. This shows that the solution of the continuous approximation indeed approximates the discontinuous one if the continuous function is steep enough.

Figure D3 Cross section of Figure D2 for x₁=5 km with porosity (a, b) and effective pressure (c, d).

However, as it is shown in the cross section in Fig. D3, the continuous approach still misses most of the steep gradient of ϕ at the position of the discontinuity. This can be improved by using an even steeper approximation of ${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{b}}$, which, on the other side, increases the computational complexity and slows the computation times.

A very similar behavior can be observed when looking at the chemical enrichment patterns connected to the hydromechanical model results. The chemical enrichment patterns are shown in Fig. D4 and the corresponding hydromechanical models are shown in Fig. D2.

Figure D4 $\mathcal{C}/{\mathcal{C}}_{\mathrm{0}}$ with decompaction weakening (σ^b) after T=1.5 Myr for ${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{b}}$ (a) and its continuous approximation (b).

Here, the smooth approximation results in a slightly blurred version compared to the discontinuous problem. Only when looking at the cross section in Fig. D5, we see a considerable difference of the enrichment at the discontinuity (at x₂=10 km) which can, for example, have a significant impact on the creation of ore deposits at specific layers in the subsurface.

Figure D5 Cross section of Fig. D4 for x₁=5 km with $\mathcal{C}/{\mathcal{C}}_{\mathrm{0}}$ on the y axis. (a) corresponds to the discontinuous solution whereas (b) corresponds to the continuous approximation.