Computational framework for complex flow and transport in heterogeneous porous media

We present a flexible scalable open-source computational framework, named SECUReFoam, based on the finite-volume library OpenFOAM(R), for flow and transport problems in highly heterogeneous geological media and other porous materials. The framework combines geostatistical pre- and post-processing tools with specialised Partial Differential Equations solvers. Random fields, for permeability and other physical properties, are generated by means of continuous or thresholded Gaussian random fields with various covariance/variogram functions. The generation process is based on an explicit spectral Fourier decomposition of the field which, although more computationally intensive than Fast Fourier Transform methods, allows a more flexible choice of statistical parameters and can be used for general geometries and grids. Flow and transport equations are solved for single-phase and variable density problems, with and without the Boussinesq approximation, and for a wide range of density, viscosity, and dispersion models, including dual-continuum (dual permeability or dual porosity) formulations. The mathematical models are here presented in details and the numerical strategies to deal with heterogeneities, equation coupling, and boundary conditions are discussed and benchmarked for the heterogeneous Henry and Horton-Rodgers-Lapwood problems, and other test cases. We show that our framework is capable of dealing with large permeability variances, viscous instabilities, and large-scale three-dimensional transport problems.


Introduction
Partial differential equations (PDEs) models in porous media are an essential part of many engineering, environmental and biological applications.They play a fundamental role in geothermal energy utilization [1][2][3], aquifer remediation [4], drug delivery in tissues [5] or bio-film formation [6,7], to name a few.Fluid flow, energy and solute transport in porous media are complex because of the variation in fluid properties and the heterogeneous nature of the porous medium.Fluid properties such as density and viscosity vary with the solute concentration and fluid temperature.These changes can lead to fluid instabilities greatly affecting the migration of dissolved substances, the mixing of fluids and chemical reactions [8][9][10][11].Heterogeneity is present across all scales.The interaction between the porous structure and the fluid flow manifest itself in the formation of preferential flow paths and stagnant regions.This in turn affects the migration and residence times of solutes, which often display an anomalous behavior [12].
Given the complexity of the phenomena related to porous media, numerical simulations provide the possibility to explore a multitude of configurations that can help analysing field and experimental data.However, producing accurate and reliable numerical data requires, considering the current state-of-the-art, advanced knowledge of numerical schemes and programming languages.Therefore, there is a need for robust and precise simulation tools capable of providing high quality results, but also flexible enough to deal with the enormous variety of porous media applications.Reproducibility of numerical results is often as important as accuracy.Reproducibility can only be achieved if data sets, results and simulation tools are accessible to the public.This can be achieved through the use of an open-source license and by building the simulation tools over well-known, well-maintained open-source libraries.
In this work we present the open-source computational framework SECUReFoam for the simulation of complex flow and transport in porous media.The framework is developed in C++ within the OpenFOAM ® library and includes new solvers, boundary conditions, pre-and post-processing utilities, and new numerical schemes.We detail the implementation of the geostatistical pre-processing library and solution strategites for saturated flow and transport solvers that account for variable fluid properties (e.g., density, viscosity, relative permeability).Heterogeneity can be accounted directly through the definition of random permeability fields and by means of double porosity models.The formers are based on multi-Gaussian and thresholded Gaussian random fields (also known as pluri-Gaussian) [13,14].Although this kind of fields can be generated with other available geostatistics toolboxes (GSLib [15], T-PROGS [16]), the integration with the flow and transport solvers provides the computational framework the capability to tackle a wide variety complex problems.
Various other open-source packages are available for porous media applications, such as OpenGeoSys [17], porousMultiphaseFoam [18], GeoChemFoam (focused on pore-scale) [19], DuMux [20], MRST [21].Our work differs from these packages in a number of ways and provides unique contributions to the community, such as: 1. fully integrated workflow from geostatistics to complex flows, including post-and co-processing; 2. focus on a clear and simple mathematical formulation and general-purpose PDE solvers with extra features developed into object-oriented external modules; 3. novel numerical methods for equation coupling and stable finite volume formulation; 4. pluri-Gaussian discontinuous fields, multi-scale features, and dualporosity models; 5. Out-of-the-box parallel scalability and C++ implementation.
The paper is organised as follows, Section 2 presents the mathematical formulation of flow and transport models in porous media at the Darcy scale.Section 3 describes the generation of thresholded Gaussian random fields.Then the numerical details of the implementation are discussed in Section 4 and illustrated with some examples in Section 5. Finally, we present some conclusions.

Flow and transport models 2.1 Darcy flow
In all our models we assume the validity of the Darcy's law for porous medium fully saturated by an incompressible fluid, which reads as follows: where is the acceleration of gravity vector.We assume that ρ and µ are non constant but dependent only of the solute concentration c [M L −3 ] although extensions to double diffusive models (where there is an additional dependency on temperature), compressible, and multiphase descriptions are currently being developed in the code.
Darcy's law can be rewritten in terms of a reduced pressure p ρgh = p−ρg•x removing the hydrostatic pressure contribution1 as follows:

Continuity equation
The continuity equation is the balance of mass of fluid per volume of porous medium.It can be formulated as where φ [1] is the medium porosity and there is a sink/source flow We can expand the first term of (3) as where is the storativity which takes into account the linearised porous matrix compressibility.

Transport equation
The transport equation gives the solute mass balance per volume of porous medium.We write the transport in terms of the solute concentration c [M L −3 ] as where c * is the concentration of the sink/source and D [L 2 T −1 ] is the hydrodynamic dispersion tensor with I being the identity matrix, D m [L 2 T −1 ] the diffusion coefficient and α L , α T [L] the longitudinal and transverse dispersion coefficients respectively.

Density and viscosity models
Fluid viscosity and density can be written as a function of the concentration of dissolved solutes as well as temperature.In the following we will focus on linear relations for ρ(c) and µ(c), namely: where f is either density or viscosity.A more extensive choice of non-linear models is available in the code (see section C).The same structure can be used for temperature-dependent viscosity and densities and extended to simultaneous dependence on concentration and temperature (double diffusive model).

Boussinesq approximation
It is sometimes convenient to rewrite the continuity equation as a condition for the divergence of the velocity field to make use of incompressible flow splitting algorithms.In this case, eq. ( 3) can be reformulated as: where it can be noticed that the right hand side contains the Lagrangian derivative of c.Therefore, by either assuming ρ ρ 1 or advection-dominated solute transport, we obtain the Boussinesq approximation, which, for nondeformable porous skeleton (i.e., S 0 = 0) and in absence of sources/sinks, reduces to the divergence free condition

Dual-porosity and multi-continuum models
The equations above are no longer a good approximation when the porous medium is highly-heterogeneous.For example, the pore space could be composed of large highly permeable pores (such as fractures) connected to a system of narrow low permeability (micro)-pores.In these cases, two separate continuity and momentum equations can be considered for two overlapping continua representing the porous spaces.Under the assumptions of eq. ( 9) and eq.(3), the system of equations for the two pressures is: where • denotes the variables in the second continuum (e.g.fractures).The transfer term between the two continua is in general an unclosed integral term depending on the connectivity between the domains, but it is commonly approximated as a linear transfer term where τ 0 is a linear transfer coefficient between matrix and fracture continua.
Another key difference is that our random field generators do not make use of the fast Fourier transform or other discrete transform.This can make the generation more computationally intensive but this is mostly overcome by the efficient C++ implementation and parallel scalability.In our tests, in fact, although the random field generation can be significantly expensive, it is nevertheless often negligible compared to the cost of the flow and transport solvers, and these steps can be fully integrated in a single run.

Continuous and pluri-Gaussian truncated fields
Gaussian random fields (GRF) have often been adopted in geostatistics to mimic the spatial distributions of geological properties due to their mathematical properties.Generally these fields well fit the purpose of geological descriptors as long as the spatial transition of the geological properties is smooth.The choice of correlation function can tune the smoothness of the function but always in a continuous manner.However, sediments' deposits rarely show a continuous structure.In these cases a continuous GRF can be post-processed with thresholding or binning to allow for abrupt transitions and discontinuities.This process is called single-truncation rule.However, if the sediment pattern is characterised by non-layered or non-stratified geometries where one category share its border with more than other two, the single-truncation model does no longer provide realistic results.To reproduce non-stratified geometries the thresholding needs to be done on a multivariate Gaussian random field.When two (or more) independent Gaussian random fields are thresholded according to a multi-dimensional truncation rule, the resulting field is called a pluri-Gaussian random field.The underlying idea is to simulate two or more GRFs on a domain and compare them through a series of inequalities which allows us to assign a unique value to each cell (Fig. 1).An example of a two-dimensional truncation rule with three (t 1 , t 2 , t 3 ) and two (s 1 , s 2 ) thresholds in the fist and second dimension, respectively, applied to two independent random realisations of GRF ( Z 1 (x) and Z 2 (x)) is depicted  If the thresholds are expressed in terms of percentiles of the GRF, the probability, i.e. the proportion of the facies i, j is Given I and J thresholds, respectively, the problem of choosing a truncation rule for N facies is underdetermined if N < (I + 1)(J + 1).Extra constraints might come from the connectivity or surface area of each facies.In the previous examples a total of twelve regions are identified by the Cartesian truncation diagram but these have been then merged into a total of four facies.This can be solved by constrained minimisation problem involving (I + 1)(J + 1) error functions.In other words, in this case, there are infinite threshold combinations which honour a given set of volume proportions and some constraints need to be defined beforehand to allow for a single solution.A more detailed analysis of this problem can be found elsewhere [26].In the computational framework we present here, the user needs to specify directly the thresholds rather than the facies volumes.

Numerical implementation
In this section we present a few more details about some of the most important solvers and tools implemented in SECUReFoam, namely: • setRandomField: a pre-processing utility to generate random fields for permeability, porosity and other properties; • rhoDarcyFoam: a variable-density flow and transport solver; • dualSimpleDarcyFoam: a dual-porosity Darcy solver; together with boundary conditions, sources and post-processing tools.Other solvers and utilities included the library but not described here include • simpleDarcyFoam: a simple single-phase Darcy solver; • multiRateScalarTransportFoam: a multi-rate mass transfer model [27,28]; • poroelasticFoam: a linear Biot poroelasticity solver; • dualRhoDarcyFoam: a dual-porosity variable density solver; • mesh importing utilities from gslib, ijk, and grdecl format ; • spatialPdf: an utility to create linear and log-spaced histograms from spatial data; • fieldMetrics: co-and post-processing utility to compute statistics of spatial data.This is a specialised and extended version of the functionObject structure of OpenFOAM ® .

Meshing and random field generator setRandomField
OpenFOAM ® is an unstructured finite-volume library that can deal with arbitrary cell types.Cartesian grids are treated internally as non-structured meshes.In the present work we have limited ourselves to Cartesian meshes generated by the OpenFOAM ® native blockMesh utility.However, included in SECUReFoam are utilities for importing corner-point reservoir models, based on the work done within the MRST project [21].Alternatively more complex geometries can be meshed combining from elementary volume objects or CAD files.Once the grid has been created, the material properties can be populated with the random field generator setRandomField.As we have seen in the previous section, this is based on the generation of GRFs.Several wellknown methods have been implemented in literature for generating GRFs [15,22,29,30].According to Mandelbrot and Van Ness (1968) [31], a GRF can be represented using a stochastic Fourier integral: where a are frequencies, dW (a) is a complex-valued white noise random measure and S(a) is the amplitude of the spectral measure.The latter can be found as the Fourier spectrum of correlation functions written in spherical coordinates.To allow for deformations along the x, y, z directions (i.e., different correlation lengths in each direction), a rescaling is implemented to the three-dimensional frequency space.This can be easily extended by including rotations, or more general transformation matrices, allowing therefore more flexibility in the orientation of the facies.The discrete form of ( 14) becomes Instead of relying on the FFT algorithm to compute (15), which would significantly reduce the computational cost of computing the Fourier integral, we have directly implemented the formula for a number of reasons.First of all, FFT can hardly be extended to non-Cartesian unstructured meshes.These are very common in geology and reservoir simulation.The other disadvantage of the FFT is that the random realisation becomes dependent on the spatial discretisation used.This makes (deterministic) convergence studies (for a given random field) hardly achievable.With (15) instead, for a given (pseudo-)random set of Gaussian variables W i , we can generate multiple random fields at different spatial resolution.Finally this approach allows us to relax the periodicity assumptions by including larger wave numbers to the sum, simulating therefore non-stationary random fields.
Once the GRFs are generated as above, they are either scaled to obtain the desired mean and variance, or thresholded, as explained in section 3.1, to obtain discontinuous fields.The options and parameters for the random field generation are described in section B.

Variable density solver rhoDarcyFoam
rhoDarcyFoam solves variable-density flow and transport problems with the models presented in section 2. Different formulation of the equations can be used, based on standard pressure (1) or the reduced pressure (2), with or without Boussinesq approximation (( 9) vs ( 8)).The transport model is based on (5).Multiple parameters can be defined as heterogeneous fields, including permeability, porosity, dispersion parameters, and storativity.Various dispersion, density, and viscosity models are included by means of a object-oriented modular structure and selectable through simple input file (following Open-FOAM ® standard, most of the physics settings have been included into the transportProperties dictionary).
The resulting system is a non-linear coupled PDE system.An under-relaxed Picard iteration loop is employed to couple the transport and Darcy equations.This is based on the pimple class in OpenFOAM ® , that allows the control of outer iterations to exit the loop when the residuals fall below a certain value or when the maximum number of iterations is reached.This also allows to remove the relaxation in the last iteration to prevent it from delaying the dynamic evolution of the fields.Using existing OpenFOAM ® capabilities, fully implicit (through Picard iterations) and explicit time stepping (by setting one single iteration) can be implemented with first and second-order backward integration.Adaptive time-stepping is included based on a maximum Courant number computed using the total fluxes.The whole formulation is based on the fluxes across faces, rather than the cell-averaged velocities.This ensures exact mass conservation.Optionally, each linear system can be solved multiple times to better incorporate explicit corrections due to non-orthogonal meshes or flux-limiter schemes, as is the case for most of the standard solvers shipped with OpenFOAM ® .
The permeability and dispersivity field can be specified as a symmetric matrix field on the faces (e.g.inverse of face transmissibility) or at the faces.While the former does not require any interpolation, the latter instead requires an interpolation to the faces.The standard interpolation schemes in OpenFOAM ® are overridden to enforce weighted harmonic interpolation for diffusive fluxes.Gradients at the faces are approximated with two-point or multi-point flux approximation (TPFA, MPFA).The latter is based on least-squares high-order gradient approximation, included in OpenFOAM ® .

Dual-porosity solver dualSimpleDarcyFoam
Dual-porosity solvers are available for simple constant density Darcy flow (dualSimpleDarcyFoam) as well as an extension of the variable density solver for dual-porosity problems (dualRhoDarcyFoam).We limit ourselves here to the description of the Darcy solver.The dual-porosity model (10) is solved in a segregated (iterative) manner.To increase the convergence of the inner iteration, we make use of a new splitting scheme [32] recently developed for coupled systems.These are based on an approximate Schur-complement that allows to increase the coupling between the equations.

Boundary conditions and other modules
Various new elements as separate modules to be used with various solvers and applications.This includes many new boundary conditions, source terms, and post-processing utilities.Most new boundary conditions have been implemented based on a new general-purpose Robin boundary condition (Robin) with linear [33] or exponential [34] reconstruction of the solution near the wall.This is then extended to include advective fluxes (RobinPhi) and then used as base class for more complex BCs.All the BCs are implemented using the standard OpenFOAM ® class structure for fvPatchFields.We describe here two particular types of BCs which will be extensively used in the numerical results.
Darcy-based boundary conditions.Specific boundary conditions have been implemented for single-phase and variable-density Darcy solvers.These include a darcyFixedVelocity condition to impose a pressure gradient to ensure a fixed inlet or outlet velocity is obtained.This automatically switches between different formulations of the pressure equation, anisotropic permeabilities, different solvers, and the presence of gravity.hydrostaticPressure is a condition that allows to specify the dynamic pressure while using a solver that includes the hydrostatic pressure in the total pressure.

Flux-based boundary conditions.
In practical applications it is often convenient to specify a boundary condition for the total fluxes at a boundary for flow and transport.To this aim, we developed the boundary conditions named darcyFixedVelocity and fixedTrans-portFlux.The former adapts automatically the specific formulation used by the solver (see section 4.2) to impose a given fixed velocity.The latter can be written in general, for a given field f , where, the first and second terms of the left hand side represent the diffusive and advective normal fluxes, respectively, and the left hand side contains the external fluxes imposed by the user.These are an advective flux with a given inlet concentration f a and the true fluid velocity u n , a discrete diffusive flux computed with a diffusion coefficient d and an external concentration f d (that accounts for the given half-cell distance from the boundary ∆), and an explicit flux F .

Numerical examples
In the following we present a series of numerical example to illustrate the capabilities of SECUReFoam.The heterogeneous and dispersive Henry problem (subsection 5.1) shows how the platform deals with variable density in heterogeneous porous media.The performance of the under-relaxed Picard is demonstrated by solving a strongly non-linear viscous fingering instability (subsection 5.2) and a highly heterogeneous Horton-Rogers-Lapwood problem (subsection 5.3), which also shows the use of truncated pluri-Gaussian fields.Finally, we present a quarter five-spot injection problem to demonstrate the dual-porosity formulation and the capability to deal with discontinuous highly heterogeneous fields.

Heterogeneous Henry problem
The Henry problem [35] is an abstraction of the seawater intrusion in a coastal aquifer.It has been extensively used to understand the interaction between the aquifer and the sea and as benchmark of variable density groundwater flow codes [36].In the Along with the original Henry problem (figure 4 left), we also consider a dispersive case (figure 4 right) and a heterogeneous case (figure 5) to illustrate the role of hydrodynamic dispersion and the heterogeneity of the permeability.The parameters for all cases can be found in Table 1.As shown in figure 4, hydrodynamic dispersion causes the interface between the seawater and the fresh water to become flatter.Dispersion also affects the movement of the  saltwater wedge, which travels further inland in the dispersive case than in the diffusive case.The same effect on the wedge's toe position is observed for the heterogeneous case.Heterogeneity, however, distorts more strongly the geometry of the saltwater-seawater interface.The low permeability zones near the sea (right) boundary modify the discharge of the freshwater and lead to the formation of high concentration zones along the boundary.The boundary conditions on that boundary are switching between Dirichlet and Neumann, therefore creating a seemingly oscillatory patterns which is however stable and purely due to the heterogeneities and the resulting horizontal velocity which switches between positive and negative.This also shows the limitation of this simple testcase in highly heterogeneous systems.

Viscous fingering
Viscous fingering is a flow instability that appears when a fluid displaces another one of different viscosity [8].The instability is caused by the difference in mobility between the fluids, which leads to the formation of fingering patterns.The patterns displays a complex tip-splitting, shielding and coalescence dynamics, which is affected by the medium heterogeneity [37][38][39] and chemical reactions [40][41][42].Figure 6 shows an example of miscible viscous fingering in a 3D cylindrical geometry.In this case viscosity is a function of concentration µ(c) = µ 0 e Rc with R = −3.A fluid with c = 1 is injected at constant rate from the left boundary displaces the resident fluid (c = 0).Pressure and a zero gradient for concentration are prescribed at the outlet.The system is characterised by the Péclet number P e = q 0 L/φD m = 10 3 , where q 0 is the velocity at the inlet, L the domain length, φ the porosity and D m the diffusion coefficient.We observe how the interface deforms and the fingering pattern appears (figure 6, left).At late times (figure 6, right), the fingers merge and the pattern is form by fewer thicker fingers.The tip splitting and shielding mechanisms are also reproduced in the simulation.These results have been obtained with 200 × 200 × 1000 Cartesian grid which is cut and adapted to the cylinder walls using the OpenFOAM ® native meshing tool snappyHexMesh.The adaptive time step is chosen to be proportional to the local mesh size and the inverse of the local velocity magnitude, with a proportionality constant of 0.1.The instability patterns are however only qualitatively independent on the mesh size and time step as no initial perturbation is imposed.Particularly important here (and even more for heterogeneous permeability fields) is the harmonic interpolation of the viscosity and diffusive coefficients to properly characterise the fingering.

Unstable flow in highly heterogeneous media
The Horton -Rogers -Lapwood (HRL) problem [43,44] is a heat transport problem in which a temperature difference is prescribed between the top and bottom boundaries of a rectangular domain.Fluid density decreases linearly with temperature so that an unstable density stratification is form that triggers a Rayleigh -Bénard instability.The system is characterised by the Rayleigh number Ra = kg∆ρH φµK th (17) where K th is the thermal conductivity, k permeability, µ viscosity, ∆ρ the density difference between the top and bottom boundaries, H the height of the domain and φ the medium porosity.For a square domain the system is stable for Ra < 4π 2 .For 4π 2 < Ra < 1300 the system becomes unstable and convection cells that occupy the whole domain form [45].For Ra > 1300 the convection regimes becomes chaotic and flow organises itself in columnar patterns [46].During the convection-dominated regimes, mixing and heat fluxes through boundaries are significantly increased with respect to the stable regime [47].
For illustration purposes we consider a 2 × 1 domain discretised using 1024 × 512 cells.Temperature 1 and 0 is prescribed at the bottom and top boundaries respectively and zero temperature gradient is prescribed at the lateral boundaries.Initially temperature varies smoothly with depth.The initial time step is set to 0.1 and the maximum Courant number to 0.5.Two heterogenous cases are solved.First, a log-normally distributed permeability field with σ 2 = 2 and correlation lengths l x = 0.2; l z = 0.05 (Figure 7 left).Second, a truncated permeability field with thresholds 2.7 • 10 −1 , 58.9, 1.83 • 10 −2 , and 3.99) chosen from the previous permeability distribution (Figure 8    It is important to notice, compared to the classical HRL problem with constant permeability, that the heterogeneities tend to stabilise the flow.In fact, both results in figs.7 and 8 are steady state, although for other statistically equivalent realisations it has been observed that small local fluctuations could still happen for variance of the log-permeability up to 2. The discontinuities in the permeability field creates, as expected, more defined structures and further stabilise the flow.

Dual-porosity and discontinuous permeabilities
Dual porosity models are a convenient representation of fractures porous media [48].In the following example we consider a quarter of a five spot geometry [49] in which a low permeability matrix is traversed by two fractures with high anisotropic permeability.Although this model is suitable to describe complex relatively homogeneous networks of fractures, for demonstration purposes, we focus here on a testcase where the fractures are localised in a cross-like structure.Therefore, the fractures permeability tensors have only one non-zero component corresponding to the orientation of the fracture.That is, only the K yy component is non-zero for the vertical fracture and the K xx component for the horizontal one.An injection takes place in the lower left corner and an extraction of the same magnitude on the upper right corner (Figure 9).The underlying matrix permeability is isotropic and constant and set to an intermediate value K mat = 10 −11 m 2 .The matrix porosity is chosen to be 0.3 for the matrix, dropping to 10 −5 in the regions where fractures are present, while the fracture porosity is close to one (0.99) in the fracture-dominated regions and close to zero (10 −5 s −1 ) elsewhere.The linear transfer coefficient τ is constant and equal to 10 −5 .A Cartesian mesh of 100 × 100 is used.The resulting system has a very strong coupling as the flow can only traverse the domain by communicating through the fracture system.This means that the iterative solution of the coupled system via the sequential solutions of matrix and fracture pressure would be very slow (up to thousands of iterations).Thanks to the Schur-based decoupling scheme implemented (the interested reader is referred to [32]), the convergence is achieved in only four iterations.Figure 10 show the resulting velocity field with the sharp transitions between matrix-dominated and fracture-dominated region.For this type of system, where fractures are very localised, the dual-porosity model assumptions breaks down, and similar results could have been obtained with a single highly heterogeneous permeability field.However, the aim here is to test the robustness of the method for solving highly heterogeneous dual-porosity permeability fields.For smoother transitions and larger-scale testcases, where one has a non-negligible fracture and matrix porosity (and permeability) everywhere, the dual-porosity model assumptions would be instead verified.Our geostatistical generation of the fields could be used to produce more realistic scenarios but with the additional difficulty of prescribing a relation or correlation between, not only porosity and permeability of the matrix (as commonly done), but also transfer coefficient and fracture porosity and permeability.

Conclusions
In this work we have presented the general-purpose open-source framework SECUReFoam for the computational modelling of various porous media flow and transport problems.An important element of this work is the combination of geostatistical tools with partial differential equations solvers.We focused here on a very general yet simple approach to generate discontinuous random fields (for porosity, permeability and other properties), namely the the truncated pluri-Gaussian simulation.We have limited our attention to unconditioned fields although it is possible to extend this to assimilate real-data by conditioning the random fields.The numerical implementation of continuous and truncated random fields is based on the explicit spectral representation evaluated in each mesh points rather than Discrete Fourier Transform, therefore independent from the space discretisation and making it suitable for the use with non-structured, non-orthogonal, and locally refined meshes, which are very common in geosciences.Future work will include the conditioning on real data (kriging) and the extension to arbitrary anisotropic correlation matrices including rotations, or more general transformation matrices, allowing therefore more flexibility in the orientation of the facies.
Several mathematical models are implemented, focusing on Darcy's flow, including dual-porosity media, and variable-density flow with different correlation models available for viscosity, density and other parameters.The computational framework also include more advanced physics, such as poromechanics, unsaturated flows, and phase-field methods but, due to the physical complexity of these models, these will be the focus of later works.All equations are solved sequentially within a fully implicit Picard-type iteration which deals with non-linearities and coupling between the equations.Operatorbased preconditioning and relaxation is used to improve the coupling between the equations.Future work will include the implementation of Newton and Krylov-based iterations for coupled non-linear operators.
We test the framework with well-known benchmark problems, namely the Henry problem for variable-density Darcy flow, the Horton-Rodgers-Lapwood testcase for Rayleigh-Benard instability, a three-dimensional pipe flow for viscous fingering, and a quarter five-spot problem for the dual-porosity model.We demonstrate how the computational framework is robust to highly heterogeneous media, instabilities and coupled problems.All solvers and testcases are available open-source online [50,51].

Fig. 1 :
Fig. 1: Pluri-Guassian simulation method: (a) two continuous Gaussian random fields are created; (b) at every position in space, the values of the continuous GRFs are used as coordinates to enter the truncation rule; (c) an heterogeneous non-Gaussian field is generated.

Fig. 2 :
Fig. 2: Qualitative truncation rule.Note that the areas of the different facies do not correspond to the their respective proportions in the simulations, not because of the qualitative nature of the image, but because the underlying continuous variables are not uniformly distributed but Gaussian.

Fig. 3 :
Fig. 3: Two-dimensional lognormal random field with (a) Gaussian correlation (b) exponential correlation and (c) Matérn correlation with ν = 1 (see appendices for the mathematical definition of these covariances)

Fig. 4 :
Fig. 4: Solution for the concentration for the original Henry (left) and dispersive Henry problem (right) at time 6000 s.As a result of hydrodynamic dispersion the freshwater-saltwater interface is flatter and advances more inland.

Fig. 5 :
Fig. 5: Map of log-permeabiity field (left) used to solve the heterogenous Henry problem permeability and the resulting concentration (right) at time 6000 s.

Fig. 6 :
Fig. 6: 3D viscous fingering at P e = 10 3 at dimensionless times 0.3 and 0.7.A high viscosity fluid (red) displaces a low viscosity fluid forming a fingering instability.As time passes the initially small fingers merge into larger ones.Shielding and tip splitting can be observed in the fingers of the panel on the right.
left).Flow and transport parameters are chosen so that Ra = 10 4 .Solutions for temperature at time 5 are shown in the right panels of figures 7 and 8.

Fig. 10 :
Fig. 10: Velocity field for the quater of a five-spot injection problem.Injected water moves preferentially through the matrix (white arrows) until it reaches the high permeability factures (yellow arrows) which focus the flow through the domain.Water leaves the system through the extraction point in the upper right corner.

Table 1 :
Henry problem the aquifer is represented by a 2 × 1 m rectangle.The right side of the domain is the boundary with the sea where hydrostatic pressure is prescribed using seawater density for the flow equation and an inlet/outlet boundary condition is used for the transport equation.On the left boundary a freshwater flow is simulated by prescribing flow Q in .The rest of the boundaries are impervious.Density depends linearly on the salt concentration ρ(c) = ρ 0 + βc and viscosity is constant.The boundary conditions for the concentration are homogeneous Neumann everywhere but the right (sea) boundary where a constant Dirichlet (c = 1) is imposed.However, to avoid the formation of a boundary layer when an outward flux develops, an automatic switching is adapted to impose Neumann condition when the velocity is pointing outwards.Henry problem parameters for the diffusive and dispersive cases.