Efficient modelling of solute transport in heterogeneous media with discrete event simulation

https://doi.org/10.1016/j.jcp.2019.01.026Get rights and content

Highlights

  • Discrete event simulation (DES) is applied to solute transport for the first time.

  • DES is innovatively implemented in a hybrid FE node-centred FV framework.

  • Parallel computations are incorporated for efficient simulation at large scales.

  • A speedup of 670 over a conventional time-stepping method is achieved.

  • Presented DES algorithm is generic and can be applied to other problems.

Abstract

To address the problem of different time scales present in the simulation of solute transport through systems with a complex permeability structure such as fractured porous rocks, we propose a parallel discrete event simulation (DES) algorithm based on local time stepping criteria, specifically developed for the hybrid finite-element node-centred finite volume (FV) framework. A preemptive-event-processing (PEP) approach is applied to synchronise discrete events with sufficiently close time stamps, thereby facilitating the parallelisation for shared memory architectures. The accuracy of the presented DES-PEP scheme is first verified against the analytical solution of a 1D advection equation with spatially variable coefficients. The DES scheme is then applied to simulate tracer advection through a 3D model of highly heterogeneous fractured rock represented by an unstructured adaptively refined mesh with over 1 million elements. DES produces results comparable to those of a conventional time-driven simulation (TDS), but uses less than 1% of the execution time. Analysis of event distributions shows that updates occur almost exclusively in a small number of FV cells marked by order-of-magnitudes faster fluid flow and advection-dominated transfer, while slow-flowing cells remain inactive and excluded from computations. This focusing of the computational effort leads to high simulation efficiency while simultaneously diminishing round-off errors. Scalability tests with a parallel version of DES on shared memory demonstrate further computational speedups mirroring the increased number of threads. With the use of 20 threads, execution time is reduced from 42.5 days (with TDS) to only 1.5 hours, equivalent to a speedup of over 670. This parallel DES algorithm therefore enables efficient multi-core simulation of solute transport in heterogeneous geologically realistic systems.

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ϕct+(vc)=q, 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 lawv=kμ(pρg), 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, u(x,t)=u0x. 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)

  • L.W. Lake

    Enhanced Oil Recovery

    (1989)
  • O.A. Cirpka

    Intrinsic remediation in natural-gradient systems

  • K. Pruess et al.

    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.

    (1990)
  • S.K. Matthäi et al.

    Simulation of solute transport through fractured rock: a higher-order accurate finite-element finite-volume method permitting large time steps

    Transp. Porous Media

    (2010)
  • D.D. Laumbach

    A High-Accuracy Finite-Difference Technique for Treating the Convection–Diffusion Equation

    (1975)
  • T.N. Narasimhan et al.

    An integrated finite difference method for analyzing fluid flow in porous media

    Water Resour. Res.

    (1976)
  • P.S. Huyakorn et al.

    The Computational Methods in Subsurface Flow

    (1983)
  • J. Istok

    Groundwater Modeling by the Finite Element Method

    (1989)
  • R. Huber et al.

    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

    (1999)
  • S. Geiger et al.

    Combining finite element and finite volume methods for efficient multiphase flow simulations in highly heterogeneous and structurally complex geologic media

    Geofluids

    (2004)
  • A. Paluszny et al.

    Hybrid finite element–finite volume discretization of complex geologic structures and a new simulation workflow demonstrated on fractured rocks

    Geofluids

    (2007)
  • G. Dagan

    Solute transport in heterogeneous porous formations

    J. Fluid Mech.

    (1984)
  • Cited by (14)

    • Adaptive total variation stable local timestepping for conservation laws

      2022, Journal of Computational Physics
      Citation 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 Resources
      Citation 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 Physics
      Citation 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 Resources
      Citation 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.

    View all citing articles on Scopus
    View full text