PyDDC: An Eulerian-Lagrangian simulator for density driven convection of $\mathrm{CO_{2}}$–-brine systems in saturated porous media

As shown in Fig. 1 the entire domain, $\mathrm{\Omega}$, is discretized into a finite number of cells, $\mathrm{N_{x}}$ and $\mathrm{N_{y}}$, in both x and y directions. Considering non-uniformity of the grid, a generic representation of cell centered co-ordinate, $x_{c}$, $y_{c}$, is given as:

with $x_{i}\mathrm{{}^{\prime}s}$ and $y_{j}\mathrm{{}^{\prime}s}$ representing the cell boundaries and $\Delta x_{i}=x_{i+1}-{x_{i}}$ and $\Delta y_{j}=y_{j+1}-{y_{j}}$ representing the cell widths.

The GSTools library developed by Müller et al. (2022) is used here to generate spatial random permeability fields with a specified covariance matrix, mean permeability, log-variance and axial correlation lengths. This package is employed into our flow module to obtain the Darcy flux and pressure field for a heterogeneous domain with a given set of user defined meta-parameters. The porosity field is determined from the permeability field via the Kozeny-Carman relation (Pape et al., 2000) $k(x)=\frac{\phi(x)^{3}}{c(1-\phi(x))^{2}S^{2}}$, where $\phi(x)$ is the porosity field, $c$ is the Kozeny constant and $S$ is a parameter depending on the grain size of the media. The realization of the random porosity-permeability ($k-\phi$) field along with their respective distributions are shown in Fig. 3.

For the FVM implementation, the discretized equations subjected to necessary boundary conditions are formulated as:

As the Darcy flux is entirely dictated by the pressure gradient and density term, we first solve for the pressure Poisson equation and use this computed pressure field to solve for the flux. Substituting equation 2 in equation 1 and integrating over the domain, we get:

where the above transformation from volume integral to surface integral is attained with the help of Green-Gauss theorem. This equation can be discretized as summing over all the facet elements belonging to each cell:

The boundary conditions on the pressure and the fluxes are represented in equation 6, where the subscript $D$ refers to Dirichlet boundary and subscript $N$ denotes Neumann boundary. The computational domain along with the different boundary conditions are highlighted in Fig. 3. Currently we only support imposition of one type of boundary condition, which is considered to be the most natural case for this simulation. Free slip boundary conditions are imposed on the Darcy flux at the top and bottom boundaries, and the pressure is taken to be hydro-static at the left and right boundaries. We consider the most common form of pressure gradient in an aquifer, represented as $P=\mathrm{P_{t}}+\rho gh$, where $\rho gh$ corresponds to the hydrostatic term and $p_{t}$ being the pressure at the top of the aquifer. The background flow into the domain can be introduced by providing an increment of pressure ($\mathrm{P_{il}}$ or $\mathrm{P_{ir}}$) at any boundary on $\Omega_{D}$.

These two fields can either be specified by the user prior to model execution or can be initialized by the inbuilt random field generator. After obtaining the random $k-\phi$ field and knowing the value of the viscosity and density (obtained either from phase module or from predefined user specified values), the flow equation can be solved to obtain the pressure and Darcy flux.

To construct the global linear system from equation 7 for both uniform and non-uniform grids, we adopt the 2nd order accurate scheme (Scheme-II of Wang et al. (2011)). We introduce new indexing notation based on the five point stencil, where we denote the center cell $ij$ by $c$ and the neighbouring cells as $n$, $e$, $s$, $w$ corresponding to ($i,j+1$), ($i+1,j$), ($i,j-1$), ($i-1,j$) cell indices. We define weights, $w$’s, used for linear interpolation from cell-centers to cell-interfaces which are given by the volume ratio between the current cell and the neighbouring cells according to:

Following this and assembling the terms on the left and right hand side, we get the final form of our pressure equation:

where,

The $P$’s are the unknown pressure terms and $\rho$ is the density of the mixture which is a function of the local concentration.

Assembling equation 9 results in the linear system of equations in the form:

where, $A$ is the coefficient matrix, $x$ is the vector of unknowns and $b$ is the vector of known coefficients. For optimization purpose, we take special care of solving the linear system of equations. If the fluid viscosity is constant, the pressure coefficient matrix $A$ is invariant. So, we perform a LU decomposition on A and precompute its inverse, $A^{-1}=(LU)^{-1}$, only once at the start of the simulation to be used at successive iterations. On the other hand, if the viscosity is concentration dependent, then the pressure coefficient matrix cannot be precomputed. Hence, in this case, we use a sparse linear solver using the previous iteration pressure values as the initial guess for faster convergence. This approach proves to be really effective even if we are dealing with very high medium heterogeneity.

To solve the non-linear Darcy equation (1) requires us to know the density field locally, i.e. cellwise. To estimate any phase attribute at any point in space, we should know the pressure and temperature value associated with that point. For our case we have an isothermal model with a global predefined temperature and pressure which is insufficient to get a density field. So we configure an initial condition for our problem assuming the local pressure everywhere to be same as the pressure at the top and based on this we can get the density field which can then be used to solve for the pressure field which in turn can be used in successive iterations. Figure 5 shows one such generated initial condition from a predefined global pressure of $30$ MPa and temperature of $50$ ⁰C.

In order to solve the flow and transport equations, it is imperative to estimate the phase attributes of the $\mathrm{CO_{2}}$–brine mixture, namely the solubility, density, viscosity and $\mathrm{CO_{2}}$ diffusivity at the ambient pressure and temperature. This is done with the help of the thermodynamic models that are developed extensively over the past years in the context of $\mathrm{CO_{2}}$ storage at geological conditions. The associated complexity of estimating all the four above mentioned phase parameters demands a detailed description of each individually.

A particle tracker, as its name implies, is a Lagrangian approach that inherently tracks discrete particles representing some physical quantity. For our case that physical quantity being tracked is the $\mathrm{CO_{2}}$–brine mixture concentration. The tracks of those individual particles are stochastic in nature and are governed by the generalized stochastic differential equation (GSDE) as given by the integration scheme (Salamon et al., 2006):

where $dX=X_{p}(t+\Delta t)-X_{p}(t)$ represents the displacement of the particle in an interval $\Delta t$, $A(X_{p},t)=v(X_{p},t)$ is the deterministic drift vector, $B(X_{p},t)$ is the displacement matrix and $\zeta(t)\sim\mathcal{N}(0,1)$ is a white noise parameter. The only caveat is that equation 21 solves the Fokker-Planck equation and not the advection-dispersion equation (ADE) represented in equation 3. The equivalence between the ADE and Fokker-Planck lies in rearranging the terms in ADE such that we end up with an augmented drift vector $A(X_{p},t)=(v+\phi\nabla\cdot D+\frac{1}{\phi}D\cdot\nabla\phi)(X_{p},t)$ for a heterogeneous case (Henri and Diamantopoulos, 2022). For a homogeneous case we drop the tailing term and end up with the two leading terms $v$ and $\nabla\cdot D$. For an isotropic case in 2D, the displacement matrix, $B(X_{p},t)$, is represented as (Salamon et al., 2006):

where $\alpha_{L}$ is the longitudinal dispersion coefficient, $u_{x}$ and $u_{y}$ are the horizontal and vertical velocity components respectively and $|u|$ is the Euclidean norm of the vector field $u$. This aforementioned scheme is also popularly known as random walk particle tracking (RWPT). Although RWPT methods are easier to implement to track individual particle trajectories, additional care must me taken such that the estimated density respects mass conservation and boundary conditions. As the particles we consider are point masses that are dispersed in the continuum, it does not possess a metric to compute field values like concentration. This makes the concentration estimated from particle locations to be a delta function at the Lagrangian node where the particle resides and is given as:

where $c(x,t)$ is the spatially varying concentration field, $x$ is any point within the domain. But to solve the non-linear problem we need to estimate the concentration values at the nodes where the flow velocity is estimated and therefore an accurate mapping from Lagrangian nodes to Eulerian nodes is required. For this reason, the metric to be estimated, i.e. the mass of the particles, are considered to be distributed symmetrically or asymmetrically around each particle giving rise to a particle support volume. This is facilitated by using a projection function in the form of a kernel that acts as an operator to estimate the concentration at any point in the domain and is represented as (Bagtzoglou et al., 1992):

where $W(x-X_{p}(t))$ is the kernel having a support volume whose evaluation is a function of the position of the particle, $X_{p}$ and the point of estimation, $x$, within a predefined interval governed by the kernel bandwidth. This method is also popularly known as kernel smoothing as we find a smoothed continuous distribution of a field from its discrete representation. The kernel contains a bandwidth parameter, $h$, responsible for the quality of density estimation and can be either considered to be a globally uniform value or a locally adaptive one (Sole-Mari and Fernàndez-Garcia, 2018). In order to be a valid approximation of the true concentration field and ensure mass conservation, the kernel function should satisfy the condition that $\int_{\Omega}W(x-X_{p}(t))dx=1$, where $\Omega$ is the domain of interest. For simplicity and computational efficiency we choose the standard simple binning method with the support volume being the dimensions of the grid cells. This makes the estimated concentration, $\hat{c}=\frac{N_{p}}{\phi\Delta V}$ which is basically the number of particle in each grid cell divided by the available cell volume (Wright et al., 2018).

In addition to this, we also need to model the particle influx from the top boundary where we have a constant $\mathrm{CO_{2}}$ solubility computed at each iteration from the solubility model. It is important to note here that as the kernel is also locally conservative, we can enforce the Dirichlet boundary condition simply by reverse engineering. We divide the reservoir into discrete grid cells that acts as mirror reflections of the interior cells. After the binned concentration is computed at each step, based on the $\mathrm{CO_{2}}$ solubility at the boundary, we can define the reservoir concentration such that the boundary conditions are strictly satisfied. After estimating the concentration values, we create particle clouds inside each reservoir control volume which is basically the grid to particle projection, contrary to the particle to grid projection that we described previously. Now once we have those particle clouds created, we apply the simple Fickian diffusion physics to diffuse particle inside domain with a characteristic diffusion coefficient, $D$. The top boundary which the reservoir shares with the computational domain is held symmetric w.r.t diffusion in order to properly impose the diffusive physics, which means that we increase particle count as more particles diffuses from the reservoir to the domain and at the same time delete particles that leave the domain from the top boundary. But this boundary is closed w.r.t dispersion. As we have different boundary conditions for the diffusive and dispersive physics, we decompose the random walk algorithm by first diffusing particles, correcting particle count by applying diffusion boundary conditions and then apply the 2D random walk represented by equation 21 respecting the dispersion boundary conditions.

Figure 6 shows the plume structure as represented by particle cloud and the corresponding density and concentration fields obtained after an elapsed simulation time of around 53 years. Addition to estimating concentration, density, pressure and velocity fields which are stored in binary format, HDF5, to be accessed and used by the user for different post-processing applications, we also store another important attribute, the vorticity field. In 2D, vorticity represents the magnitude of convection which can be used to understand the circulation inside an enclosed domain. For a 2-component vector field $\vec{q}=[u,v]$, vorticity, $\omega$, is given by the relation:

For our case,

From the above expression it can be seen that once the fingers starts propagating, the convection cells starts to develop due to the horizontal density gradients between the fingers and its magnitude increases with increasing density finger propagation. The convective cells are indicated by the velocity plot in Fig. 6 showing the magnitude that advective velocities attain during convection.

The user input to the model is accepted through a JSON file present at the top level of the working directory, that holds the physical attributes required for the simulation. A JSON file represents structured data as key value pairs, with the keys holding the attributes names and the values being their corresponding quantification. Figure 7 gives the detailed insight into the JSON file with the physical parameters along with their units. All the physical parameters are grouped under four categories depending on the parameter type. domain_params groups all the parameters required for defining the computational domain. Grid refinement is achieved by two input parameters — refine_levels, which takes into account the total number of cells from the top where refinement has to be applied and res_levels, which holds the factor by which the initial resolution in the y-direction should be increased for cells affected by refinement. Transition zones due to refinement are handled internally in order to avoid considerable jumps in the cell aspect ratios between the refined and unrefined zones. kfield_params holds the necessary meta-parameters, mean_k, var and corr_length, required to create the heterogeneous random permeability field. This field, kf, is created by invoking the Field class function Field.KField(x, y, mean_k, var, corr_length) where x and y represents the spatial coordinates. For simulations requiring heterogeneous porosity fields, the Field class attribute Field.PHIField(kf) is used. boundary_params contains the dirichlet pressure boundary conditions required for solving the Darcy flux. reservoir_params contains the lateral extent of the Sc-$\mathrm{CO_{2}}$ source, horizontal_extent and molesperparticle, which is a quantity carried by individual particles. The user-defined phase attributes for the model are stored in phase_params, which contains the salinity (NaCl), saturated concentration (CO2SaturatedConcentration), $\mathrm{CO_{2}}$ diffusion coefficient in brine (DiffusionCoefficient) and end-member viscosity (CO2SaturatedBrineViscosity, BrineViscosity) and density (CO2SaturatedBrineDensity, BrineDensity). All user-specified attributes in the JSON file are converted to Python objects and are initialized in the Python file V.py. File TH.py holds the features required for computing phase attributes if the inbuilt phase module is used. These files are declared in the formulate module and can be accessed, modified and used internally whenever required.

The components of the numerical software along with their contents are shown in Fig. 9. Proper execution of the numerical solver requires that the different software components interact efficiently with one another. This is leveraged by the systematic coupling of the different parts of the module as shown in Fig. 8. The Simulate class handles the simulation process and is designed to establish the overall flow of control. To create a Simulate object, the user should pass the arguments — filename, the JSON file which contains the physical parameters used in the simulation, UsePhaseModule, a boolean variable indicating whether to use the inbuilt phase module, datafile, having the name of the binary HDF5 file to store output results and the user supplied permeability and porosity field values, kf and phif, if the inbuilt random field module is not used. Internally it passes the filename to the initialization class as IP.ModelInitialization(filename) creating the attribute repository file V.py which contains domain information along with all other necessary physical parameters. This class also initializes a dictionary, Simulate.attr used for storing the values of concentration, density, viscosity and diffusion coefficient computed at each iteration and making it available to the user for access. If the boolean variable UsePhaseModule corresponds to false, the Simulate class checks for the end-member values of viscosity and density from V.py and based on local concentration, c, declares interpolation functions, Simulate.mu_func(c) and Simulate.rho_func(c), to get local cell-wise values. Conversely, if UsePhaseModule is true, the Simulate.attr dictionary is initialized with attributes computed from the inbuilt phase module. The Simulate class also initializes the random fields either supplied by the user as its arguments, or generated from the inbuilt field class, Field, as described in the previous section. It also links to another class SpeciesInfo which contains information about the solubility of $\mathrm{CO_{2}}$ and mole fraction of different species in the brine. To start the simulation, the user should call the class function ParticleTracker(steps,realizations,intervals,params) in Simulate, where the arguments respectively are number of simulation steps, spatio-temporal evolution of plume dynamics for a single $k-\phi$ field, alternate periods where the user wants to store the data in years and user-specified thermodynamic parameters belonging to Simulate.attr list which is computed form phase equilibrium models other than the one used here. The initial condition is configured by calling the _configure_init_condition() function at the start of every simulation each time the ParticleTracker function is invoked. The initial condition computes the flow field for the Lagrangian particle simulation to work. Solving the flow field is handled by the FlowSolver(field) class which takes in the random permeability field as the only argument and contains the attributes — pressure, vorticity and horizontal and vertical flux components. Simulations pertaining to constant viscosity coefficients and variable viscosity coefficients are handled separately. Functions _AssemblePressureCoefficientMatrix_cv() and _AssembleCoefficientMatrix_rhs(r, PL, PR) are designed to handle cases arising from constant $\mathrm{CO_{2}}$–brine viscosity. The former assembles the left hand side of the pressure coefficient matrix and the later assembles the right hand side depending on the density, r and surface pressure values at the boundary, PL and PR. The left hand side of the pressure coefficient matrix, for a constant permeability and viscosity field, is invariant and assembled only once and inverted at the start of the simulation. This inversion is carried out by a LU-decomposition technique as discussed earlier. For a variable viscosity case, the pressure coefficient matrix can no longer be precomputed and hence has to be iteratively solved at each step. This is carried out by the function _GlobalCoefficientMatrix_vv(). The solve(r, mu) method is invoked by passing in arguments density, r and viscosity, mu, to solve for the pressure field and Darcy flux. At every step the particle tracker defines a particle reservoir, uses the flow field to advect the particle cloud according to the random walk method described before, computes the local concentration by binning and re-computes the flow field for the successive iterations. The RWPT(c) class is responsible for solving the transport equation, where the argument c represents the concentration field. Advection of the particles requires the knowledge of local dispersion tensor, which is computed inside CellwiseDispersion() function present inside FlowSolver(field) class. To disperse particle cloud as per equation 16, the disperse(Dxx, Dxy, Dyx, Dyy, phi, vx, vy, plst) function is invoked which solves the SDE by taking into account the cell-wise disperson tensor (Dxx, Dxy, Dyx, Dyy), local Darcy flux (vx, vy) and porosity (phi) values interpolated at particle locations (plst).

At each iteration or specified intervals, the simulation results are exported in binary format and stored in the memory. As indicated in figure 8, these datasets are divided into two categories – temporally variable, consisting mainly of computed scalar and vector fields like concentration, density, pressure, velocity, vorticity and particle plume configuration that changes at each simulation step and temporally invariable attributes like the domain information consisting of Eulerian nodes and the random $k-\phi$ field that are initialized only once at the start of the simulation.