Efficient modelling of solute transport in heterogeneous media with discrete event simulation
Introduction
Understanding the transport of solutes (e.g., chemicals, tracers, minerals or contaminants) in porous geological formations is essential for a broad range of geo-engineering applications, including enhanced recovery of oil and gas from hydrocarbon reservoirs [1], remediation of groundwater pollution in aquifers [2], and storage of nuclear waste or CO2 in the subsurface [3]. The high degree of heterogeneity in geological systems creates various numerical and computational challenges when solving fluid pressure equations and advection–diffusion equations describing the solute transport [4]. Various numerical methods have been proposed to discretise the governing equations and simulation variables in space and time. Typical examples are finite difference methods (FDM) [5], [6] and finite element methods (FEM) [7], [8], with the latter having more flexibility to represent complex geometries. Combinations of FEM and finite volume method (FVM) have become increasingly popular [4], [9], [10], [11], [12]. In this approach, a node-centred finite volume mesh is superimposed on the finite element mesh. The elliptic fluid pressure equation is solved by FEM and the hyperbolic transport equation by FVM. This combined approach shows improved accuracy and efficiency over the fully coupled FEM, while retaining the geometric flexibility of finite elements [10], [11].
How to efficiently integrate these spatially discretised equations over time is still a fundamental challenge. Traditionally time advance is based on a global time step, which synchronously evolves the system at spatially uniform time intervals. These time-driven simulation (TDS) techniques fall into two categories: explicit and implicit schemes. In explicit schemes, unknown variables are computed at the current time level from quantities already available from previous times. The size of the integration time step is restricted by the most stringent Courant–Friedrichs–Levy stability criterion (also known as the Courant or CFL condition) in the system, which guarantees the capture of the fastest transient change, and ensures numerical stability. This becomes problematic in simulations of solute transport in highly heterogeneous porous media where flow properties, such as permeability, can vary over many orders of magnitude [11], [13], [14], leading to dramatic variations in transport velocities. This problem is further exacerbated in models constructed with variable mesh sizes where the smallest elements and control volumes tend to represent areas of highest permeability (e.g. rock fractures) and therefore coincide with regions of enhanced fluid flow [4]. The resultant infinitesimal global time step leads to prohibitively slow TDS simulations.
Stability of implicit schemes is not restricted by the severe CFL constraint permitting larger time increments. However, this is at the expense of accuracy and an increased computational burden from solving a system of equations at each iteration step. Matthäi, Nick, Pain and Neuweiler [4] found that when simulating solute transport through a 3D stochastic discrete fracture model containing 2000 disc-shaped fractures, a single implicit transport step required approximately 20 times the CPU time needed for a single explicit step, even though an efficient algebraic multigrid solver, SAMG [15], was used. Furthermore, implicit schemes introduce numerical smoothing and tend to be too diffusive [16]. Omelchenko and Karimabadi [17] argue that implicit integrations conducted with large time increments are unable to correctly represent local physical phenomena with characteristic process frequencies higher than the inverse of the time-step size.
To get around global time step restrictions, a number of local time-stepping approaches have been proposed to allow the time step to vary spatially satisfying a local, rather than global CFL condition. In these approaches which are also known as adaptive or multiple time-stepping schemes, the solutions at different elements or cells are updated with different time increments determined from local CFL conditions [18], [19], [20], [21]. Because of the asynchronous update of different elements, it is difficult for these methods to correctly and efficiently calculate numerical fluxes across element interfaces. Some use time interpolation for flux computation, however, these violate local conservation of fluxes, leading to reduced numerical accuracy and stability [17]. Local time-stepping approaches are often used in combination with adaptive mesh refinement, where the spatial grid is locally refined in order to represent complex geometries and better capture more active regions [22], [23], [24], [25], [26], [27]. Spatial elements on each level of refinement are clustered into a logically rectangular patch, leading to a hierarchy of nested patches with different temporal resolutions that locally satisfy the CFL conditions. In order to preserve conservation laws, global synchronization steps are required to impose flux corrections at patch interfaces. For efficiency, a fixed local time step size, computed at the beginning of each global synchronization step, is used throughout the global time step. This is problematic for highly complex systems as the CFL conditions may change rapidly and significantly during a global synchronization step, making the pre-determined time integration sequence inaccurate and unstable.
Discrete event simulation (DES), also called event-driven simulation, provides an alternative to conventional time-stepping schemes. Unlike TDS where the simulation progresses in pre-defined global or local time steps, in DES temporal evolution of the global system is driven by discrete events. This forces transition from one state to another at irregular (or asynchronous) points in time [28], [29]. Such discrete events, representing effective units of information in DES, are simulation objects characterised by a process function responsible for updating the object state, and a time stamp indicating when the process function is scheduled to be executed. A DES simulation begins with scheduling events based on local rates of change in individual elements. All scheduled but not yet executed events are listed into an event queue sorted by their time stamps in increasing order. The top event always corresponds to an event with the earliest time stamp. The DES simulation then progresses by repeatedly removing the top event from the priority queue, executing it by calling its process function and adding a new event with an updated time stamp to the queue. Consequently, DES avoids using pre-determined time steps and integration sequence as in the TDS schemes, but allows individual parts of the simulation domain to evolve based on their own physically determined temporal scales, which may vary in space and over the progress of the simulation. This effectively removes the global CFL restriction and, at the same time, ensures numerical accuracy and stability. DES focuses updates on active parts of the system, while eliminating the computational overhead associated with idling processors.
DES was originally developed for naturally discrete systems, as found in operations research, management science, games and telecommunications [28], [29], [30]. It was extended to modelling continuous systems [31], [32] with reported speedups by a factor of 300 in some cases [33]. The DES scheme has been applied to solve flux-conservative partial differential equations with source terms [17] by introducing a self-adaptive predictor–corrector scheme. The predictor schedules events with a time delay estimated from the local variation rate, and the corrector ensures that pending events are processed in time, if their causality constraints undergo a significant change. This self-adaptive paradigm, in combination with a synchronization mechanism for updating neighbouring events, effectively preserves causality and guarantees numerical stability. Later work by Stone, Geiger and Lord [34] associated an event with the mass transfer at element interfaces, rather than the evolution within an element. In order to improve event processing by preventing unnecessary inter-element synchronisations in parts of the computational domain where the solution properties are fairly homogeneous, Omelchenko and Karimabadi [35] introduced a preemptive-event-processing (PEP) technique, which automatically enforces synchronous execution of events with sufficiently close timestamps. So far the advances in DES simulation of continuous systems have been demonstrated on low-dimensional uniform meshes. In this paper we will investigate the use of DES to complex and realistic problems at larger scales. We will also make use of the advantage that the PEP technique can be used to easily parallelise DES at least for shared memory architectures – a fact which was not realised in the original publication.
The objective of this paper is to develop a parallel DES algorithm for a finite-element node-centred finite volume scheme. The implementation uses the Complex Systems Modelling Platform (CSMP++) [36], an object-oriented application programmer interface designed for simulation of complex geologic processes and their interactions [37]. It provides an efficient simulation framework for solute transport through heterogeneous porous media on high-dimensional, complex and unstructured meshes. The rest of this paper is structured as follows. In section 2, the governing equations for solute transport in porous media are presented, followed by a brief description of the discretization of these equations in the CSMP++ FEM-FVM framework. We then describe in detail the algorithmic and programming aspects of the parallel DES-PEP for the node-centred FV transport scheme. In section 3, the proposed DES algorithm is verified with analytical solutions of the 1-D advection equation with spatially variable coefficients. In section 4, we apply it to tracer advection through a heterogeneous 3D model of fractured rock, and compare its results and performance with conventional TDS simulations. We also conduct a scalability analysis for multiple threads on this model. Section 5 concludes this paper with a summary of the most important results obtained.
Section snippets
Governing equations
Ignoring diffusion which is readily modelled using the finite element approach, advective transport of a non-reactive solute through a porous medium is governed by the general advection equation where c (kg/m3) represents the solute concentration, ϕ refers (m3/m3) to the porosity of the medium, t (s) is time, and q (kg/m3 s) describes external sinks and sources of the solute. The transport velocity v (m/s) is given by Darcy's law in which k (m2) refers to the
Verification of DES scheme for 1D models
Analytical techniques for the solution of the advection equation are generally restricted to simple problems with specific initial and boundary conditions. Here we adopt the analytical solutions to the one dimensional advection equation with spatially variable coefficients [42], where the velocity field is set as a linear function of distance, . Based on different initial and boundary conditions, we have conducted an analysis for two scenarios, a steady point source and a
Performance evaluation on a complex 3D fracture model
The previous simple 1D tests verified the accuracy and correctness of the implemented DES algorithm. In this section we further evaluate its computational performance on a complex field data-based model.
Conclusions
In this paper we present a new parallel DES algorithm for hybrid FEM–FVM simulations performed on unstructured grids. It is implemented in the node-centred FV framework of the CSMP++ modelling platform, for the efficient simulation of solute transport through heterogeneous porous media in a stationary flow velocity field. Compared with conventional synchronous time-driven simulation, the implemented DES scheme has the following advantages: (1) efficiently removing the global CFL restriction via
Acknowledgements
The presented research is funded by the Carbon Capture and Storage Research Development and Demonstration Fund (CCS49356) “Australian subsurface carbon sequestration Simulator” awarded by the Department of Industry, Innovation and Science, Australian Government. Dr. Yuri Omelchenko is thanked for providing additional technical context about his publications on the subject and for sharing his enthusiasm about the potential of DES as a fundamentally new and powerful computational method. The
References (44)
A triangle based mixed finite element – finite volume technique for modeling two phase flow through porous media
J. Comput. Phys.
(1993)- et al.
An asynchronous scheme with local time stepping for multi-scale transport problems: application to gas discharges
J. Comput. Phys.
(2007) - et al.
Self-adaptive time integration of flux-conservative equations with sources
J. Comput. Phys.
(2006) - et al.
Adaptive mesh refinement for hyperbolic partial differential equations
J. Comput. Phys.
(1984) - et al.
Adaptive mesh refinement and multilevel iteration for flow in porous media
J. Comput. Phys.
(1997) - et al.
Multi-level adaptive simulation of transient two-phase flow in heterogeneous porous media
Comput. Fluids
(2010) - et al.
A new asynchronous methodology for modeling of physical systems: breaking the curse of Courant condition
J. Comput. Phys.
(2005) - et al.
Asynchronous discrete event schemes for PDEs
J. Comput. Phys.
(2017) - et al.
A time-accurate explicit multi-scale technique for gas dynamics
J. Comput. Phys.
(2007) - et al.
River flow forecasting through conceptual models, part I — a discussion of principles
J. Hydrol.
(1970)
Enhanced Oil Recovery
Intrinsic remediation in natural-gradient systems
On thermohydrologic conditions near high-level nuclear wastes emplaced in partially saturated fractured tuff, 1: simulation studies with explicit consideration of fracture effects
Water Resour. Res.
Simulation of solute transport through fractured rock: a higher-order accurate finite-element finite-volume method permitting large time steps
Transp. Porous Media
A High-Accuracy Finite-Difference Technique for Treating the Convection–Diffusion Equation
An integrated finite difference method for analyzing fluid flow in porous media
Water Resour. Res.
The Computational Methods in Subsurface Flow
Groundwater Modeling by the Finite Element Method
Multiphase flow in heterogeneous porous media: a classical finite element method versus an implicit pressure–explicit saturation-based mixed finite element–finite volume approach
Int. J. Numer. Methods Fluids
Combining finite element and finite volume methods for efficient multiphase flow simulations in highly heterogeneous and structurally complex geologic media
Geofluids
Hybrid finite element–finite volume discretization of complex geologic structures and a new simulation workflow demonstrated on fractured rocks
Geofluids
Solute transport in heterogeneous porous formations
J. Fluid Mech.
Cited by (14)
Modelling CO<inf>2</inf> plume spreading in highly heterogeneous rocks with anisotropic, rate-dependent saturation functions: A field-data based numeric simulation study of Otway
2022, International Journal of Greenhouse Gas ControlAdaptive total variation stable local timestepping for conservation laws
2022, Journal of Computational PhysicsCitation Excerpt :The scalability of these approaches is limited by the amount of work in the timestepping group with the smallest timestep. Other parallelization approaches have treated the adaptive local timestepping problem as a discrete event simulation [30–35]. While these approaches achieve tremendous speed-ups, these methods typically provide only heuristics and examples as proof of robustness.
Adaptive conservative time integration for transport in fractured porous media
2022, Advances in Water ResourcesCitation Excerpt :Kulka and Jenny (2022) extended the method to unsteady compressible flow. Related approaches include discrete event simulation (DES), which is a totally asynchronous time stepping method (Nutaro et al., 2003; Shao et al., 2019), and multirate methods which integrate each component of the system using a different time step (Constantinescu and Sandu, 2007; Delpopolo Carciopolo et al., 2019, 2020). This work extends the method of Jenny (2020) to tracer transport in fractured porous media and applies it to two-dimensional models where we compare accuracy and computational cost of ACTI with those of global time stepping.
Temporally adaptive conservative scheme for unsteady compressible flow
2022, Journal of Computational PhysicsCitation Excerpt :It is based on re-evaluation of the time step sizes after each global step in order to capture high wave speeds. Nutaro and Zeigler [7] introduced discrete event simulation for solving simple conservation laws, and Shao et al. [8] developed a parallel implementation thereof. Impressive speedup could be demonstrated for flow in fractured porous media.
Simulation of two-phase flow in porous media with sharp material discontinuities
2020, Advances in Water ResourcesCitation Excerpt :Nick and Matthäi (2011a) To alleviate this effect, the idea of the interface treatment, which was introduced in section (2.2), could be combined with the discrete event simulation (DES) technique. DES updates only a small number of FV cells at any one time while excluding slow-flowing cells from computations (Shao et al. (2019)). In this work, we present a capillary interface coupling scheme for immiscible two-phase flow and transport in heterogeneous porous media, taking account discontinuities of material properties.