A DEIM driven reduced basis method for the diffuse Stokes/Darcy model coupled at parametric phase-field interfaces

In this article, we develop a reduced basis method for efficiently solving the coupled Stokes/Darcy equations with parametric internal geometry. To accommodate possible changes in topology, we define the Stokes and Darcy domains implicitly via a phase-field indicator function. In our reduced order model, we approximate the parameter-dependent phase-field function with a discrete empirical interpolation method (DEIM) that enables affine decomposition of the associated linear and bilinear forms. In addition, we introduce a modification of DEIM that leads to non-negativity preserving approximations, thus guaranteeing positive-semidefiniteness of the system matrix. We also present a strategy for determining the required number of DEIM modes for a given number of reduced basis functions. We couple reduced basis functions on neighboring patches to enable the efficient simulation of large-scale problems that consist of repetitive subdomains. We apply our reduced basis framework to efficiently solve the inverse problem of characterizing the subsurface damage state of a complete in-situ leach mining site.


Introduction
Fluid flow through porous media is typically modeled with the Darcy equations.When there are large cracks and voids in the porous medium, then a homogenization of the material into a single permeability tensor is no longer appropriate.The creeping flow in those domains can more accurately be modeled with the Stokes equations.The coupling relations between the Darcy and Stokes flow regimes were first studied by Beavers and Joseph [1], and later supplemented by Saffman [2].The resulting Stokes/Darcy equations coupled via Beavers-Joseph-Saffman interface conditions play an important role in many disciplines, for instance modeling groundwater flow [3], petroleum and karst reservoirs [4,5], in-situ leach mining [6]), perfusion of blood through tissue [7,8], filtration devices [9,10], and chemical reactors [11].The wide range of applications has lead to significant interest in methods that computationally approximate the coupled Stokes/Darcy problem, including their numerical analysis and algorithmic treatment [12][13][14][15][16][17][18][19][20].Many of these applications are characterized by incomplete or uncertain data in terms of the geometry and topology of the two flow domains (e.g., uncertain subsurface soil characteristics, limited resolution of CT or MRI scans).At the same time, they generally involve a wide range of length scales (e.g., from small rock cavities to a complete mining site, from small capillaries to a full organ).
In this article, our objective is to enable the efficient simulation of such multiscale systems with parametrically defined free-flow Stokes domains at the lowest scale.To achieve this objective, we make use of model order reduction by means of reduced basis methods [21][22][23] on small repetitive subdomains.A reduced basis method replaces a computationally expensive high-dimensional finite element discretization of a parametrized partial differential equation (PDE) by a small low-dimensional set of basis functions that have high approximation power with respect to the solution manifold of the parametrized problem [24,25].We develop a reduced basis method for the coupled Stokes/Darcy model on a variable internal geometry, and model the large-scale system consisting of many repetitive subdomains by coupling together many such reduced basis functions.Similar approaches to coupling reduced basis functions on repetitive subdomains have been used in [26][27][28][29][30][31].
In a reduced basis context, the reduced basis functions themselves are defined on a high-dimensional finite element approximation space.Consequently, the high-resolution finite element mesh must remain fixed.This requirement conflicts with our aim of flexibly handling the geometry of the internal Stokes and Darcy domains.Usually, parametric geometry is handled by mapping back to a reference domain [22,23], but such an approach does not permit changes of topology as this would degenerate the Jacobian of the mapping.We therefore require a method that is able to model a topologically flexible Stokes domain that may merge and disperse within the Darcy domain.Diffuse interface methods [32][33][34][35][36], also known as diffuse domain or phase-field methods, offer such a flexible framework for solving coupled boundary value problems on non-boundary fitted meshes.The internal geometry is implicitly represented by a phase-field function, which smoothly transitions from zero to one.The phase-field indicator function and its gradient can be leveraged to replace integrals on subdomains or interfaces by weighted volumetric integrals on the complete domain.The resulting phase-field formulation is equivalent to the sharp-boundary interface problem when the width of the diffuse interface (controlled by a characteristic length-scale parameter) limits to zero.Phase-field geometry representations have been widely applied, for instance in growth modeling of tissues and crystals [37,38], for tracking the evolution of crack patterns [39][40][41][42], for enforcing boundary conditions in imaging data based analysis [43,44], for modeling multi-phase flow [45][46][47], for variational image processing and segmentation [48,49], or for modeling phase transition and segregation processes [50][51][52][53].
A crucial aspect of a reduced basis method is the affine decomposition of the linear and bilinear forms.Affine decomposition enables the precomputation of reduced-basis stiffness matrices such that the final reduced order model is completely independent of the size of the original high-fidelity model.However, the (bi)linear forms resulting from the diffuse representation of the internal geometry do not satisfy this criterion.This is a common issue for many relevant parametrized PDEs, such as those that feature nonlinearities [54][55][56][57], complex material laws [58][59][60], or in general, spatially varying model coefficients [61,62].For those cases, the discrete empirical interpolation method (DEIM) may be used [61,63].With DEIM, all non-affine parameter dependent fields are replaced by a low-dimensional interpolation on DEIM modes.For our diffuse representation, this method solves the problem of non-affine parameter dependence, but the interpolation of a domain indicator function is likely to produce Gibbs-type oscillations at regions with high gradients (i.e., the diffuse interfaces).Oscillations necessarily imply negative values and values larger than one, which in turn could produce nonphysical domain representations and unstable system matrices.In this paper, we mitigate this issue by introducing a non-negativity preserving version of DEIM.
Our article is structured as follows: in Section 2, we derive the diffusely coupled Stokes/Darcy equations.We show that all three Beavers-Joseph-Saffman conditions can be treated naturally within a diffuse interface framework.In Section 3, we propose a non-negativity preserving variation of the discrete empirical interpolation method for the dimensional reduction of the phase-field geometry representation and show its effectiveness for three benchmark problems.In Section 4, we use the same three benchmark problems to study the relation between the required number of phase-field DEIM modes and the number of reduced basis functions in the reduced order model.In Section 5, we apply our methodology to efficiently estimate the subsurface flow characteristics of an in-situ leach (ISL) mining site that consists of a large number of repetitive hexagonal units.In Section 6, we summarize our work and draw conclusions.

Diffusely coupled Stokes/Darcy equations on parametrically defined domains
We consider an incompressible fluid moving through a partially porous medium at velocities that are sufficiently small to neglect the convective components in the material derivative.Steady state equilibrium of the fluid is then described by: with σ the Cauchy stress tensor, u the fluid velocity, f the body force, and Ω the d-dimensional spatial domain.

Weak formulation of the coupled Stokes/Darcy equations
Our porous medium in Ω contains voids and cracks where the flow is unobstructed.We denote the union of these (potentially disconnected) subdomains Ω S , and in these subdomains we use the constitutive relation of an incompressible Newtonian fluid: ) with µ the viscosity and p the pressure field.Hence, in Ω S , the governing equations are the Stokes equations.
In the remaining domain, Ω D = Ω \ ΩS , the interaction between the porous medium and the fluid produces a flow resistance that is assumed to be sufficiently high such that the viscous effects in the Newtonian fluid may be neglected.The closing relations become: where κ is the second order permeability tensor.These assumptions lead to the Darcy equations in Ω D .The interface that couples the Stokes and Darcy domains is denoted Γ = ΩS ∩ ΩD .On Γ , the unit normal vectors n S and n D point out of the Stokes and Darcy domains, respectively.Across the interface, the solution fields are coupled due to the required balance of mass and balance of linear momentum.Additionally, the porous material introduces a shearing resistance to the fluid on the Stokes domain.These considerations lead to the so-called Beavers-Joseph-Saffman coupling conditions: where subscripts S and D refer to the solution fields in the Stokes and Darcy domains, respectively.The tensor T = (I − n S ⊗ n S ) is the tangential projector and α is the shear resistance parameter.
On the boundary of Ω, denoted ∂Ω, we permit the following boundary conditions: ) ) The first two conditions represent no-inflow with free-slip and the last condition is a pressure condition.These specific boundary conditions are chosen because they are valid essential and natural conditions for both the Stokes and the Darcy equations.The inflow condition is chosen homogeneous for simplicity, but it may just as well be set to a non-zero value.
We obtain a weak formulation by multiplying Eq. (2.1) by test functions v and q, integrating over the domain Ω, substituting the constitutive relations Eq. (2.2) and Eq. ( 2.3) on their respective domains, performing integration by parts, and by substituting the coupling and boundary conditions of Eqs.(2.4) and (2.5).This leads to the statement: where the function space H(Ω S , Ω D ) is defined as: This space ensures sufficient regularity on the vector valued functions in the Stokes and Darcy domains.

Diffuse interface representation
Next, we consider a parametrically defined domain Ω S (β) (and thus also where β is a point in the parameter space P. The parameter space P is finite dimensional and encodes in some sense the set of potential geometries of Ω S .For example, parameters in P could denote the number, position, size and/or orientation of voids and cracks through the domain Ω.To handle this extensive geometric flexibility, we introduce a diffuse representation of the interface geometry in the weak statement.The strength of such an implicit interface representation is that the computational mesh does not have to fit the internal interface.Let φ(β) ∈ C 1 (Ω) be a "phase-field" indicator function with values between 1 and 0. The function φ(β) tends to 1 in the Stokes domain and to 0 in the Darcy domain, and it monotonically decreases from 1 to 0 along straight lines that cross a thin region around the Stokes/Darcy interface.The thickness of the diffuse interface is characterized by the parameter δ, as illustrated in Fig. 2.1.Based on this indicator function, volume integrals on Ω S or Ω D and surface integrals on Γ can be approximated by volume integrals on Ω according to: ) In the following, we indicate dependencies on β only when relevant.The normal vector n S may be approximated as −∇φ/||∇φ||, such that: (2.9) Finally, we assume that the solutions u are such that ||T u D || ||T u S ||, whereby T u S on Γ can be approximated as twice the solution in the diffuse interface.
By making use of these approximations in Eq. (2.6), the new weak statement becomes: where we introduce the new second order tensor field α, defined as: and the function space H 0 (φ, ∂Ω u ), defined as: (2.12)

Relation with the Brinkman model
Diffuse geometry representations are often plagued by non-physical behavior in the diffuse interface region.
The only way to address this issue is to reduce the length-scale parameter δ and thereby reduce the interface width.As a result, significant adaptive mesh refinement is required to ensure that the induced error is below a tolerance threshold, which can easily lead to a prohibitive increase in computational cost.We would like to point out, however, that the weak formulation (2.10) corresponds to model equations with physical relevance even for a non-vanishing interface width.
Assuming that the obtained solution pair (u, p) is sufficiently smooth, we can perform integration by parts to arrive at the following statements: which hold for all v ∈ H(φ) and all q ∈ L 2 (Ω).The corresponding strong form equations are: with μ = φµ and κ−1 = (1−φ)µκ −1 +α.These are the Brinkman equations [64], which are used frequently for describing a sufficiently fast moving fluid in porous media [4,65,8].The field φ, which was originally a domain indicator, is now a material parameter, and the third Beavers-Joseph-Saffman coupling condition (2.4) has emerged as an additional orthotropic contribution to the material porosity.

Finite element approximation
A finite element formulation based on the weak formulation (2.10) no longer requires a mesh that fits the interface.Instead, we are able to use the same mesh for computing the solution for all domain configurations in the parameter space P.This is a crucial requirement for obtaining a reduced basis representation of the parametric problem later on.The mixed nature of the equations, however, warrants careful selection of the finite element spaces.This is particularly challenging for the coupled Stokes/Darcy equations: the well-posedness requirements, i.e., the Ladyzhenskaya-Babuška-Brezzi (LBB) conditions, in the Stokes limit are different from those in the Darcy limit [66].Since our domain is parametric, but our mesh and approximation spaces are fixed, we seek a pressure/velocity pair of elements that is stable for both a pure Stokes problem and a pure Darcy problem.So, we require approximation spaces V h × Q h that are subspaces of the relevant function spaces for all domain configurations: and satisfy the LBB conditions in the limiting cases of both the Stokes and the Darcy equations.One such pair is the combination of linear nodal elements for the construction of Q h and so-called MINI elements for the construction of V h [67,65].The MINI element is a simplicial element and consist of linear interpolation functions enriched with a bubble function [68].We use this combination of elements throughout the remainder of this article.For all simulations, we make use of the FEniCS finite element library [69] for computing the large stiffness matrices and the PETSc library for manipulating those matrices [70].
3 Non-negativity preserving discrete empirical interpolation of the phase-field To ensure efficient formation of the stiffness matrix in the reduced order model, it is imperative that the bilinear and linear forms under consideration depend on the parameter β in an affine sense.That is, we require: where B(•, • ; β) and L(• ; β) are the bilinear and linear forms corresponding to Eq. (2.10).The (bi)linear forms B i (•, •) and L i (•) are parameter independent, and the parameter dependency of B(•, • ; β) and L(• ; β) is captured through multiplication with the scalar-valued functions θ B i (β) and θ L i (β).

The discrete empirical interpolation method
In our case, the complex parameter dependency of the phase-field φ(β) does not naturally permit a decomposition of the bilinear form as shown in Eq. (3.1a).To remedy this deficiency, the affine decomposition can be approximated with the discrete empirical interpolation method (DEIM).DEIM produces the following low-dimensional representation of a field: where fi (x) are the parameter independent functions in space (the DEIM interpolation modes) and θ f i (β) are the corresponding parameter dependent weighting functions.In summary, the DEIM algorithm for obtaining such a decomposition is a two-step procedure: 1.The manifold M f = {f (x; β) : β ∈ P} is explored and the modes that best represent the complete manifold in some norm are identified.In practice, a discrete representation of the manifold is obtained as a snapshot matrix of solution vectors of the projection of the functions f (x; β) onto a finite element space for a sampling of the parameter space.The modes fi (x) follow from a principal orthogonal decomposition or singular value decomposition of the snapshot matrix.2. The functions θ f i (β) are designed such that Eq. (3.2) induces an interpolation of f (x) by the modes fi (x) at a selection of points.These points are iteratively determined: the n-th point is the location where the interpolation of the n-th mode by the first n − 1 modes produces the largest error.When the interpolation locations xj are defined for j = 1, • • • , N f , the weighting values θ f i (β) are computed by sampling f (x j ; β) and by solving the interpolation system of equations: The reader is referred to [61,62] for more details.

A non-negativity preserving alteration
The following properties of the phase-field φ and tensor field α guarantee that the stiffness matrix is positive semi-definite: These conditions are naturally satisfied by the true phase-field.To guarantee stability of the reduced basis method, they must also be satisfied by the DEIM approximation of these fields after the reconstruction strategy of Eq. (3.2).This is challenging due to the nature of the possible fields φ(β): they are mostly constants of 1 or 0 except at thin regions where they have steep gradients.The non-local support of the DEIM interpolation functions would produce Gibbs-like oscillations in the neighborhood of the interface.
To remedy this issue, we propose to rewrite the fields (3.4) in the following form: We will now approximate the new fields ξ, ζ and t i with Eq. (3.2) and the DEIM algorithm.The affine representation of the original fields then follows by expanding the summation multiplication: where, in Eq. (3.6c), d is taken 2 for simplicity of notation.
Since the final approximations are constructed from squares, these series expansions will satisfy the requirements of Eq. (3.4).The downside of this approach is the increase of the cost for the summations involved.

Benchmark problems
In the following, we investigate the relation between the quality of the approximation and the required number of modes for three benchmark problems.For each of these problems, we plot the first four DEIM modes φk and tk , and we show reconstructions of the true fields φ and t for representative parameter points.We also qualitatively compare our non-negativity preserving DEIM reconstruction to a standard DEIM approximation of that same field.

Circular domain with an inwards spiraling channel
The first benchmark problem consists of a circular domain of porous material with a circular void in the center.The phase-field, representing Stokes flow, is spiraling inwards.The parameter space P = [180 • , 540 • ] is one-dimensional and denotes the number of full rotations of the spiral.The material parameters are  , and the boundary pressures t n are set to 1000Pa and 0Pa at the exterior and interior circles respectively.The problem set-up is illustrated in Fig. 3.1a, and Figure 3.1b shows the high-fidelity mesh that we use for all subsequent computations.First, we obtain the DEIM approximation modes of the fields ξ(β), ζ(β) and t(β).We use a equidistantly spaced discrete approximation of the parameter space with 1001 samples for the parameter sampling: For each β ∈ P h we construct discrete approximations ξ h (β), ζ h (β) and t h (β) as interpolants of the true fields on the high-fidelity mesh.We collect the solution vectors corresponding to the interpolated fields in a matrix, and performing a singular value decomposition on this snapshot matrix.The first four left-singular-vectors of the fields ξ(β) are shown in Fig. 3.2.
As an example, we use these interpolation modes to reconstruct the field φ at the parameter point β = 360 • .The original field φ(360 • ) is shown in Figure 3.3a.The other subfigures show the reconstructions with the non-negativity preserving DEIM method for different numbers of modes.As may be observed from the color bar, the reconstruction of φ is non-negative everywhere, as guaranteed by our reconstruction strategy of Eq. (3.6).When fewer than ten modes are used for the reconstruction then there is still a significant mismatch between the approximated and the true phase-field.From ten onward, the reconstruction clearly resembles the true field, and the difference between the approximated field and the true field is minor when the reconstruction makes use of twenty modes.To demonstrate the importance of such a non-negativity preserving reconstruction, we employ the standard DEIM method for the reconstruction of the same field φ(360 • ) with the same interpolation modes.The result for five and ten interpolation modes are shown in Fig. 3.4.As the color bar indicates, there are banded regions around the interface in the approximated phase-field that are negative.Such undershoots would result in a negative diffusive term, potentially rendering the weak formulation (2.10) unstable.
Next we focus on the approximation of the field t(β).The first four interpolation modes are shown in Fig. 3.5, and an example reconstruction for the case β = 360 • with various numbers of modes is shown in Fig. 3.6.We observe a similar trend as for the approximation of φ: ten interpolation modes are required to yield an approximation that resembles the true field, and for twenty interpolation modes the approximation and the true field are qualitatively difficult to distinguish.
The singular values of the singular value decomposition of the different snapshot matrices can be leveraged to quantify the approximation power of a certain number of modes of either one of the fields ψ, ξ and t with respect to the entire (discrete) parameter set.This notion is made rigorous by the Eckart-Young-Mirsky theorem, which relates the singular values to the optimal relative error (in the Frobenius     norm) that can be achieved when representing the entire snapshot matrix with the first n interpolation modes [71].Based on this theorem, we consider the following quality measure: where σ i are the singular values.Figure 3.7 shows (n) for the different fields ξ, ψ and t.The overlap of the lines ξ and ψ for the majority of the for graph indicates that the approximation of ξ and ψ may be expected to perform equally well.This is to be expected since they are closely related.The graphs also shows that t is the most challenging field to represent accurately with a finite-dimensional linear approximation.Still, with a total of ∼30 modes, the measure for the relative error drops below 1%, which is more than sufficient for applications where interface geometries are described by a phase-field.The 20-mode reconstruction of φ, shown in Fig. 3.3d, corresponds to an ξ -value of approximately 0.01, and the 20-mode reconstruction of t, shown in Fig. 3.6d, corresponds to an t -value of 0.03.The difference between (n) for the fields ξ, ψ and t implies that a computationally optimal reduced basis method will require different numbers of interpolation modes to represent the different phase-field quantities.We explore this concept in more depth in Section 4.

Rotating channel
Next, we consider an angled straight channel of Stokes flow that is embedded in a square domain of porous material.The parameter is the angle of rotation of the channel, such that the parameter space  angles of the rotating channel, the phase-field is largely uncorrelated.Widely varying entries in the φ(β)-manifold means that the Kolmogorov n-width of the manifold is large and that linear dimensional reduction requires many modes to produce accurate approximations.The purpose of this benchmark is to investigate the impact of this problem.We are interested in the required number of DEIM modes such that we are still able to obtain acceptable reconstructions of the internal geometry.We again sample the parameter space with 1001 equidistantly spaced samples, producing as a discrete approximation of the full parameter space.Figure 3.9 shows the DEIM modes of the field ξ.Examples the reconstruction of the field φ(30 • ) are shown in Fig. 3.10 for the nonnegativity preserving DEIM approximation and in Fig. 3.11 for the standard DEIM approximation.We again observe large patches with negative phase-field values in Fig. 3.11, and no regions with negative           To quantify these qualitative observations, we plot the error measure for each of the three fields in Fig. 3.14.The slow decay of the -values in this graph confirm that the DEIM representation of the geometry for this benchmark is less effective than it was for the first benchmark problem.For example, the graph indicates that ∼70 modes are required to reach an -values of ∼1% for the t-field, whereas ∼30 modes were required to reach the same threshold for benchmark problem 1 (as shown in Fig. 3.3d).The acceptable reconstructions of Figs.3.10d and 3.13d again correspond to ξ ≈ 0.01 and t ≈ 0.03, implying that these are appropriate target values.method, which is based on a reinterpretation of the geometric variability as a mapping from a reference domain [22,23].
The parameter space P of this benchmark problem is three-dimensional.We construct the sampling set P h as a Cartesian product of 21 uniformly distributed samples along each axis, yielding a set of 9261 samples:  The approximations of φ and t still show defects even when 40 interpolation modes are used.Naturally, this is due to the higher dimensionality of the parameter space.We can determine roughly how many modes would be required to obtain adequate approximations by investigating the error measure (n).This measure is plotted in Fig. 3.21 for each of the three fields.When compared to Figs. 3.7 and 3.14, we indeed observe a significantly larger overall .Where benchmarks 1 and 2 require 30 and 70 modes respectively to reach an t of 1%, this threshold is not reached for benchmark 3 within the first 100 modes.
For benchmarks 1 and 2 -values of 1% and 3% are sufficient to obtain an acceptable approximation of φ and t respectively.For the current benchmark problem, these -values are reached for ∼65 and ∼80 modes respectively.We consider this an excessive amount.In the next section, we investigate the importance of the accuracy of the DEIM reconstructions of φ(β), ψ(β) and t(β) on the fidelity of the reduced basis formulation.4 Reduced basis method of the coupled equations on the parametric domain From Section 2, we have the following parameter dependent finite element formulation: where the bilinear and linear forms are those of Eq. (2.10).The mixed finite element space V h ×Q h is that described in Section 2.4 and could be kept independent of β.This approximation space is assumed to be sufficiently refined such that the finite element solutions (u h (β), p h (β)) can for all intents and purposes be considered the true solution (u(β), p(β)).The meshes illustrated in Figs.3.1b, 3.8b and 3.15b satisfy this requirement for the previous three benchmark problems.We now wish to find a low-dimensional subspace RB r ⊂ V h 0 × Q h (with r = dim(RB r )) that has high approximation power with respect to the full solution manifold M (u h ,p h ) = {(u h (β), p h (β)) : β ∈ P}.Once this low-dimensional subspace has been found, it may be used in place of V h 0 × Q h in Eq. (4.1) to produce a reduced order model for approximating any solution pair (u h , p h ) ≈ (u, p) in M (u h ,p h ) for any β ∈ P at a very low computational cost.The procedure for determining RB r again relies on a singular value decomposition of a snapshot matrix, which is in this case constructed from solution vectors of high-fidelity finite element solutions.We perform this computation for all three benchmark problems.The snapshot matrices are filled with solutions to Eq.The singular values corresponding to the reduced basis functions again relate to the approximation power of the cumulative basis compared to the full (discrete) solution manifold.This is quantified by the measure from Eq. (3.7), which is plotted for all three benchmark problems in Fig. 4.4.We observe the same trends as we did for the DEIM geometry representations in Section 3: benchmark 1 performs best, followed by benchmark 2, and benchmark 3 is worst.Again, this can be motivated from the complexity of the problems, and the related Kolmogorov n-width of the solution manifold.The higher irregularity of benchmark 2 produces a higher Kolmogorov n-width than the better-behaved benchmark problem 1, but the higher dimensionality of benchmark 3 has a more significant impact.Nevertheless, all three benchmark problems achieve an -value below 5% for less than 20 solution modes.In order to efficiently use these functions as basis functions in the reduced basis formulation of Eq. (4.1), we require an affine decomposition of the (bi)linear forms into parameter independent bilinear forms multiplied by parameter dependent weights.This permits a precomputation of the corresponding stiffness matrices during an offline phase, such that a full integration of the reduced basis functions can be avoided during the online phase.Based on the DEIM representation of all parameter dependent fields from Section 3, an affine decomposition of the bilinear form follows as: with: With the functions from Figs. 4.1 to 4.3 as approximation bases for the solution space, our final reduced order model reads: where u represents the solution tuple that combines the velocity and the pressure solution.The matrix representation of this problem becomes: Find ûRB ∈ R r s.t.: where all K ∈ R r×r , F ∈ R r , and ûRB is the vector of coefficients for each of the reduced basis functions.This reduced order model is operated by taking the following steps: 1.In the offline phase (i.e., before operation), the matrices corresponding to the bilinear forms B φ i , B ψ i and B α i are precomputed.2. Also in the offline phase, the interpolation matrix from Eq.  3. In the online phase (i.e., during operation), the weight vectors θ ξ (β), θ ζ (β) and θ t (β) are computed for the given parameter β. 4. From those weights, we compute the weights θ φ i (β), θ ψ i (β) and θ α i (β) for to the non-negativity preserving approximation of Eq. (3.6). 5.The summation from Eq. (4.4) is carried out and the resulting r × r system of equation is be solved.Subsequently, the obtained low-order approximation can be visualized in a postprocessing step by computing the weighted sum of the reduced basis functions and the computed coefficients.
Figures 4.5 to 4.7 show example results of the complete reduced order model for all three benchmark problems for various numbers of reduced basis functions (i.e., various r).In all cases, 50 modes were used for the DEIM approximation, which we consider the maximal number for maintaining computational efficiency.For all three benchmarks, 20 solution modes are sufficient to produce reduced basis approximations that are sufficiently accurate in the "eye-ball norm" for these particular parameter points.Only upon careful inspection of the velocity magnitudes are we able to distinguish the approximations from the high-fidelity computations.Interestingly, benchmark 1 and benchmark 3 produce seemingly acceptable results already for ten solution modes, whereas benchmark 2 clearly does not.The remaining open issue is the relation between the number of DEIM modes for each of the three fields, the number of solution modes in the reduced basis, and the fidelity of the reduced order model throughout the parameter space.Recall that (n) plotted in Figs.3.7, 3.14, 3.21 and 4.4 represents the lowest possible relative error that may be achieved when approximating all solutions in the snapshot matrix with the first n modes.This does not guarantee that the DEIM approximation and the reduced order model yield solutions close to that optimum.Additionally, the measure concerns an average over the parameter space.This says little about the maximum relative error, which is arguably the more important measure.Still, these computed values may guide the design of the reduced order model in terms of number of modes for the various reduced order approximations.
As a first indicator of the quality of the reduced order model, and its dependence on the DEIM approximation of the phase-field geometry representation, we compute the maximum L 2 -error of the velocity field of the reduced order model over the complete (discrete) parameter space as a function of r for various numbers of DEIM interpolation modes.The results are plotted in Figs.4.8a, 4.9a and 4.10a for each of the three benchmark problems.The general trend in each of the figures is the same: for small r values, the maximum L 2 -error is independent of the number of modes for the DEIM approximations.As r increases above a certain threshold, the L 2 -error starts to drop.This drop plateaus at different values for the different numbers of DEIM modes.For the first two benchmark problems, an acceptable maximum relative error of ∼3% is achieved for fewer than 25 solution modes.For benchmark 3, however, the maximum error does not drop below 10%.Apparently, 50 DEIM modes are insufficient to accurately represent the internal geometry.Indeed, the -values for the phase-field quantities of benchmark 3, previously graphed in Fig. 3.21, decayed significantly slower than those of benchmarks 1 and 2.
In Section 3, we show that the DEIM approximation of the three phase-field related quantities performs differently well.At the same time, the error in the reconstruction of the three fields may have different impact on the final solution error.These observations imply that a different number of modes may be required for the DEIM approximation of the three fields.A more targeted and optimized choice of the number of DEIM interpolation modes could reduce the computational expense of operating the reduced order model without adversely affecting the quality of the final solution.We propose to equate the -value for the particular choice of number of solution modes (i.e., those in Fig. 4.4) to the average of the -values for the DEIM interpolation modes (i.e., those in Figs.3.7, 3.14 and 3.21).That is, we use the minimal values n ξ , n ζ and n t for the DEIM reconstruction of the respective fields such that: which we solve by iteratively incrementing n ξ , n ζ or n t depending on which reduces the right-hand-side of Eq. (4.6) most.We limit the number of DEIM modes for each field to 50 to avoid excessive computational expense.
With the choice of number of DEIM interpolation modes for ξ, ζ and t based on Eq. (4.6), we recompute the L 2 -errors throughout the entire parameter space.The results are graphed for each of the three benchmark problems in Figs.4.8b, 4.9b and 4.10b.The graphs markers at intermediate points in the graphs state the predicted number of DEIM interpolation modes as (n ξ , n ζ , n t )-tuples.The figures also include a copy of the blue 50-modes lines from Figs. 4.8a, 4.9a and 4.10a as references.We observe that the maximum L 2 -errors of the optimized approach closely follows the ("overkill") 50-modes line while substantially reducing the required number of DEIM interpolation modes.The optimized mode number lines break shortly before the 50-modes lines plateau, confirming that more than 50 DEIM modes are required for at least one of the fields to maintain a drop of error.These results indicate that Eq. (4.6) provides an effective estimation for the required number of DEIM modes for all fields such that the DEIM approximations do not dominate the source of error of the reduced order model.

Example application: in-situ leach mining
In-situ leaching (ISL), also known as solution mining, is a mining process for recovering underground minerals such as copper and uranium.In 2019, ISL accounted for 57% of the total uranium mining worldwide [72].During ISL, a solution (typically a mixture of native groundwater, a complex agent and an oxidant) is pumped through an injection well down to an ore deposit.At this depth, the solution flows towards the production well, while dissolving minerals from the ore.The pregnant solution is then pumped to the surface at the production well, where the dissolved uranium is later extracted from the solution.The injection and production wells are distributed in regular patterns.The most frequently used pattern types are the 3-spot, 5-spot, and 7-spot patterns, shown in Figure 5.1.The choice of pattern depends on factors such as subsurface permeability, deposit type, ore grade and installation cost.The distance between the injection and productions wells can be as high as 50m−60m for the 3-spot pattern, and as low as 15 -30m for the 5-spot and 7-spot patterns.As a result, the 5-spot and 7-spot patterns have a higher installation cost but increased uranium recovery rate and operation flexibility.
Challenges that are faced while constructing and operating an ISL mining site include (i) ensuring that nearby groundwater is not contaminated by only operating sections of the mine and continuously taking measurements at monitoring wells surrounding the mine, and (ii) managing the permeability of the ore deposit with hydraulic fracturing or controlled explosives.In both cases, a predictive simulation tool with fast response time could facilitate decision making processes, may serve to optimize long-term operation of the mining site and enables more targeted interventions.We make use of the reduced order model described in the previous sections to develop such a simulation tool.We illustrate its operation capability by solving the inverse problem of predicting the damage state of the ore deposit based on inflow and outflow measurements at the injection and production wells.

The high-fidelity model
In the following, we examine the Honeymoon mining project illustrated in Fig. 5.2a.This mine is located in South Australia and was active from 2011 to 2013, after which operations were halted due to the low price point of uranium [73].The mine features a 7-spot pattern.A distinctive property of this pattern is its tiling into hexagonal segments, illustrated in Fig. 5.2b.The complete mining site consists of 225 such hexagons.We exploit this regular pattern in our computational model for the subsurface flow throughout the entire mining site: we develop a reduced basis method that can represent the flow in a single hexagon, and reuse this model in each tile.For each hexagon, the distance between the single production well to each of the six injection wells is 15m.The diameter of the injection and production wells is 0.5m.The uranium-bearing sand lies 100 -150m underground, under sheets of gravel, clay, and sand.We treat it as coarse-grained pyritic sand, with a permeability of 100 millidarcy (10 −13 m 2 ).Since the solution consists largely of groundwater, we assume a fluid viscosity of 1mPa•s.The boundary conditions that we assume for each hexagon are a variable positive pressures at the injection wells, a zero pressure at the production well and symmetry conditions on the remaining boundaries.These conditions correspond to the ones from Section 2 (i.e., Eq. (2.5)).
The rate at which uranium can be leached depends not only on the permeability of the sandstone but also on the crack patterns running through it.These cracks create pathways in the ore deposit for the solution to penetrate.Once cracks become too large they reduce the effectiveness of the uranium recovery.To include the crack patterns in our model, we construct an evolving damage field that runs from each of the six injection wells to the production well.The damage field depends only on a single parameter, the damage parameter D, and represents the severity of the cracked state for each sextant of the hexagon.The damage parameter varies between 0 and 1, signifying an undamaged state or a fully developed crack from injection well to production well.Figures 5.3 and 5.4 show the evolving damage field and the corresponding velocity profile for D between 0 and 1 in one sextant.For D = 0 there is no damage and the resulting velocity profile follows from pure Darcy flow.While D increases it creates a mixed Stokes/Darcy domain, which can be interpreted as a homogenization of smaller cracks in the sandstone.The resulting velocities in this created channel are therefore higher than in the pure Darcy domain.This difference increases with higher values of D. At D = 0.5 a thinner channel on top of the existing one is developing and represents the growth of a single straight crack.At D = 1 this thin crack is fully developed, represented by a domain of Stokes flow.

The reduced order model in one hexagon
Computing finite element approximations of the high-fidelity model on each of the 225 hexagons would yield an excessive amount of degrees of freedom for the complete system.Solving the resulting discrete system would be too time consuming for optimization purposes or for solving inverse problems.Moreover, for each new damage field in each different hexagon, the system would have to be reassembled, adding to the severe computational expense.Hence, we create a reduced basis method with a very limited number of degrees of freedom within each hexagon.To ensure that our reduced order model is capable of accurately representing the possible solution states up to a tolerance level, we perform the procedure discussed at the end of Section 4 to determine the required number of DEIM and RB modes.
In each hexagon, there are six parameters running from 0 to 1 that determine the overall damage state.Additionally, there are six pressure boundary conditions, which we vary between 0 and 10 5 .The resulting parameter space is thus twelve-dimensional: P = [0, 1] 6 × [0, 10 5 ] 6 .To explore the solution manifold, we combine a tensor product sampling grid with a random sampling set, giving P h = P h U P h R .We construct P h U as {0.25, 0.5, 0.75, 1} 6 × {10 5 } 6 such that it contains 4096 points representing different internal geometries (damaged states).To also include the variable pressure, we construct P h R from another 6000 points that are sampled with a uniform random distribution along every axis in P. For each parameter point in P h , we compute the phase-field related fields ξ, ζ, t, as well as the solutions (u h , p h ) on the highfidelity mesh.Figures 5.5a and 5.5b show the measures from Eq. (3.7) as a function of the number of modes for the phase-field quantities and the solution fields respectively.Based on these results we choose  25 solution modes per hexagon and, using the guided choice discussed above, (21,30,48) for the DEIM modes respectively.With this set-up, we obtain the following reduced order model in each individual hexagon: or, in matrix representation: where, of course, the parameter dependence of K is affine due to the DEIM approximation of the geometry representation.

Connecting reduced order models
Due to the symmetry conditions between the hexagons, there is no mass-flow through any of the boundary facets, nor is there a shear stress acting on those facets.These homogeneous conditions are satisfied by each of the solution snapshots in the snapshot matrix, and thus by each of the reduced basis functions.All coupling between the hexagonal cells then occurs through shared connectivity with the injection well at each of the corners.The available measurement data at each of these injection wells is the total amount of fluid that is pumped into the underground system from the surface.Since the distribution of the mass-flow among the neighbouring hexagons is unknown, this produces a coupled system of equations.
To introduce this connectivity, we add a Lagrange multiplier constraint that enforces a prescribed total inflow for each of the injection wells.For injection well 2 in the example of Fig. 5.6a, the variational statement corresponding to the inflow condition reads: where q 2 is the scalar test function that corresponding to the constraint.The inflow conditions of Eq. ( 5.3) can be recognized as the transpose of the pressure boundary conditions of Eq. (5.1a).Indeed, the Lagrange multipliers corresponding to the constraint of q 2 is the pressure p 2 , which now becomes a degree of freedom.
The full system of equations is then assembled as follows.The 25 degrees of freedom for the reduced basis approximation within each hexagon form a 25 × 25 block on the diagonal.For each of the injection wells, another row and column are added to the system.The column and the row are populated with three F jvectors (respectively their transposes) corresponding to the three neighboring hexagons.This procedure is illustrated in Fig. 5.6.For the complete Honeymoon mining-site shown in Fig. 5.2, this results in a system of equations of "merely" 6,077 degrees of freedom.Without the reduced basis representation of the flow in each hexagon, the full system would exceed 15 million degrees of freedom.For any given new subsurface damage state in any particular hexagon, the corresponding matrix K(D 1 , • • • , D 6 ) has to be recomputed (which is a cheap operation due to its affine dependence on the local damage parameters), and then the block on the system matrix can simply be overwritten.As a result, we are able to recompute a solution state for a completely new damage parameter in around 0.1s on an Intel i9-9900k @5Ghz, computing in serial.The majority of the compute time is spent on assembling the 225 diagonal blocks.Since these blocks are completely independent of one another, this computation could alternatively be performed in parallel very effectively.For our purposes, though, the 0.1s computation time suffices.

Solving an inverse problem: calibration of the subsurface state
The subsurface state affects the flow profile, for which measurement data is available in terms of mass throughputs at the injection and production wells.Determining the subsurface state based on the measurement data thus involves solving an inverse problem.Based on Sections 5.1 to 5.3, we have a low-order computational model for the subsurface flow through the entire Honeymoon mining site of Fig. 5.2.To produce an example problem, we artificially generate a distribution of inflow and outflow measurements.We do so by assuming smoothly varying inflow pressure values for the injection sites, ranging between 8.5 • 10 4 Pa and 1.0 • 10 5 Pa.We also assume a smoothly varying damage field.We overlay this damage pattern with random noise of about 1% and add some larger damage values to represent local cracks.We then compute the outflow field in /s of the described problem by multiplying the resulting throughput in each production well by a height of 5m.This is the assumed intake height of the production wells inside the uranium deposits.Finally, we add random noise to these computed values to ensures that our model will not trivially be able to obtain the exact solution.The resulting incompatibility between inflow and outflow also occurs in actual mining sites due to leakage of the solute into different ground layers and due to extraction of groundwater.The final pattern of outflow measurements is shown in Fig. 5.7.Given only this measurement data, we now solve an optimization problem for the damage parameter D in each sextant of each hexagon.The objective is to minimize the error between the given measurements of Fig. 5.7 and the resulting outflow at each production well of the coupled computation.As an initial state, the damage parameter is homogeneously set to 0. We employ a gradient-descent method to iteratively adjust each damage parameter, until the model predicts the correct outflows for the given inflow.The gradient is computed numerically with a central difference finite difference scheme, which results in two function evaluations per sextant for changes in the one-dimensional damage parameter.Such an approach is only feasible due to the reduced-order model.To track the achieved accuracy of the scheme, we determine the error between the computed outflow and the measured outflow according to: where u i,out is the measured outflow for production well i, and u h i,out is the computed outflow.The convergence of this relative error with respect to the number of iterations is shown in Fig. 5 graph indicates, we reach a root mean square relative error below 1.5% in under 10 iterations.This residual error can be attributed to the artificially added effects of leakage and groundwater extraction.
Having obtained a prediction for the damaged state, we can plot the resulting velocity field in Fig. 5.9, illustrated as a distribution of velocity magnitudes.These velocity magnitudes confirm that the crack patterns are more developed in the areas where higher mass flows are measured compared to the areas where low mass flows are measured.In particular, in the bottom center a row of hexagons developed cracks with damage parameters around 0.95.This area is of particular interest, because it is connected to the outside of the well field and could lead to the uranium-rich solution escaping into the groundwater.Such results should prompt the installment of additional monitoring wells in the vicinity of the area.In the top center of the mine, where low mass-flow (under 1.5 /s) for one hexagon is measured, the subsurface soil is highly impermeable with damage parameters around 0.05.Such results may instigate hydraulic fracturing measures at these areas.

Conclusions
In this article, we have developed a reduced basis method capable of handling the coupled Stokes/Darcy equations on variable internal geometry.We based our model on a diffuse interface representation of the Stokes and Darcy domains, which permits changes in topology of the parameter dependent internal geometry.We showed that all three of the Beavers-Joseph-Saffman coupling conditions can be treated straightforwardly in a diffuse geometry setting, and we highlighted the equivalency with a particular type of Brinkman model.We used the discrete empirical interpolation method (DEIM) to handle the non-affinity of the phasefield quantities that arise in the coupled equations.To ensure non-negativity of the domain indicator required by physics, we introduced a non-negativity preserving version of DEIM.We showed that in total, we require three such DEIM representations of phase-field related quantities that arise from the diffuse geometry representation and the Beavers-Joseph-Saffman coupling conditions.We studied the performance of the non-negativity preserving DEIM reconstruction for three benchmark problems.The benchmarks differed in complexity, which could also be observed in the convergence graphs for , a relative error measure depending on the singular values for the DEIM interpolation modes.
The same three benchmarks were used to investigate the performance of a complete reduced basis formulation.The solution modes that act as reduced basis functions were obtained from the singular value decomposition of a snapshot matrix.We investigated the dependence of the number of solution modes and the number of DEIM interpolation modes for each of the three fields on the maximum L 2 -error of the reduced order model throughout the parameter space, and we proposed a strategy for determining a required number of DEIM modes for each of the phase-field quantities for a given number of solution modes based on precomputed -values.Making use of this strategy, we obtained the optimum maximum relative L 2 -error for a given number of solution modes while significantly reducing the number of DEIM interpolation modes compared to a naive approach.
We then illustrated the application of our reduced order model for a large-scale computations via the analysis of an in-situ leach mining site.A reduced order model for the subsurface flow throughout the entire mining site was constructed from reduced basis approximations on subdomains and by subsequently coupling the different subdomains with Lagrange multipliers.We showcased the capability of the resulting model for predicting subsurface crack patterns for an existing in-situ leach mining site based on measured inflow and outflow data at injection and production wells.A complete model evaluation required approximately 0.1 seconds on a desktop machine, and involved solving a system of equation of 6,077 degrees of freedom.This constitutes a significant cost reduction compared to the 15 million degrees of freedom that are required to solve the full problem without the reduced order model at virtually the same fidelity level.

( a )
True field.(b) Reconstruction with five modes.(c) Reconstruction with ten modes.(d) Reconstruction with twenty modes.
(a) Reconstruction with five modes.(b) Reconstruction with ten modes.
(a) True field.(b) Reconstruction with five modes.(c) Reconstruction with ten modes.(d) Reconstruction with twenty modes.

1 0. 2 (
, and t n = 1000 at the left boundary and t n = 0 at the right boundary.The problem is illustrated in Fig.3.8a, and the corresponding high-fidelity mesh is illustrated in Fig.3.8b.We may anticipate that the DEIM approximation is very ineffective for this type of problem: for different parameter points β, i.e., different 1 a) Schematic drawing.(b)High fidelity mesh.

( a )
True field.(b) Reconstruction with ten modes.(c) Reconstruction with twenty modes.(d) Reconstruction with forty modes.

( a )
Reconstruction with ten modes.(b) Reconstruction with twenty modes.

( a )
True field.(b) Reconstruction with ten modes.(c) Reconstruction with twenty modes.(d) Reconstruction with forty modes.

( a )
True field.(b) Reconstruction with ten modes.(c) Reconstruction with twenty modes.(d) Reconstruction with forty modes.

( a )
Reconstruction with ten modes.(b) Reconstruction with twenty modes.

( a )
True field.(b) Reconstruction with ten modes.(c) Reconstruction with twenty modes.(d) Reconstruction with fifty modes.
(4.1) on the meshes of Figs.3.1b, 3.8b and 3.15b for the same parameter samplings described in Sections 3.3.1 to 3.3.3.Figures 4.1 to 4.3 show the first four resulting reduced basis functions for each of the three benchmark problems.The figures show pressure fields overlapped by velocity fields.This is indicative of the nature of these functions: since we perform a singular value decomposition of a snapshot matrix of velocity-pressure pairs, each of the left singular vectors represents a combination of a velocity and a pressure field.

Fig. 4 . 4 :
Fig. 4.4: Convergence of for the solution fields with the number of modes.

( a )
True solution.(b) Approximation with five modes.(c) Approximation with ten modes.(d) Approximation with twenty modes.
Naive choice of number of DEIM modes.
Optimal number of DEIM modes (b) Guided choice of number of DEIM modes.

Fig. 4 . 8 :
Fig. 4.8: Convergence of the relative maximum L 2 -error over the complete parameter space as a function of the number of reduced basis modes for different number of DEIM modes for benchmark 1.
Naive choice of number of DEIM modes.
modes Optimal number of DEIM modes (b) Guided choice of number of DEIM modes.

Fig. 4 . 9 :
Fig.4.9:Convergence of the relative maximum L 2 -error over the complete parameter space as a function of the number of reduced basis modes for different number of DEIM modes for benchmark 2.
Naive choice of number of DEIM modes.
Guided choice of number of DEIM modes.

Fig. 4 . 10 :
Fig. 4.10: Convergence of the relative maximum L 2 -error over the complete parameter space as a function of the number of reduced basis modes for different number of DEIM modes for benchmark 3.

(a) 3 Fig. 5 . 1 :
Fig. 5.1: Patterns of injection and production wells used in practice for ISL mining.
Solution fields on one hexagon.

Fig. 5 . 5 :
Fig. 5.5: Convergence of as a function of the number of modes.

Fig. 5 . 6 :
Fig. 5.6: Example of three adjacent hexagons with corresponding coupling matrix (green bars indicating the position of the coupling vectors F j ).

Fig. 5 . 7 :
Fig. 5.7: Assumed measurements of outflow at production wells over entire well field.

Fig. 5 . 8 :
Fig. 5.8: Error between measurements and computed outflows over the number of iterations.