ABSTRACT. In this talk, we address some of the challenges in mathematicalmodelling and numerical simulationof complex nonlinear coupled systems. The focus is on variationalinequalities, mixed-dimensional formulations and non-localdifferential operators in space and time. We discuss several aspectssuch as the proper analytical framework, the different role offractional operators andefficient numericalalgorithms. For each aspect, we give an application relevant example andprovide numerical simulation results illustrating the flexibility ofthe proposed techniques. Among others,  we considercomplex ice-water phase transitions in soil freezing models and thegrowth of micro-vascular networks based on the coupling of bloodflow and oxygen transport.

ABSTRACT. The characterization of thermal transport in fractured rocks is crucial for understanding numerous systems of environmental, geological and industrial importance. Fracture wall roughness exhibits long range spatial correlations and induces a heterogeneous aperture field, thus promoting the formation of preferential flow channels within fracture planes. Here, we develop a modelling approach to identify geometrical parameters of individual fractures that control the heat exchange at the fluid/rock interface and the advection of heat along the fracture.

ABSTRACT. Here, we present a numerical analysis of available Equivalent Continuous Porous Media (ECPM) upscaling methods for use in groundwater flow modelling and subsequent radionuclide transport simulations. The exercise is based on rock domains of the Forsmark groundwater ow model. Results are interpreted in terms of statistical distributions of groundwater travel time and transport resistance as well as radionuclide breakthrough curves. The sensitivity of the ECPM solution to grid refinement and to the selected upscaling method is discussed and empirical guidelines for use in safety assessment (SA) applications are provided.

ABSTRACT. Simulating the propagation and (bio)degradation of contaminants in heterogeneous fractured porous media is considered by developing a new modeling strategy that relies on the advantages of multi-scale multi-continuum random walk approaches. The method is adapted to deal with non-linear bio-chemical reactions occurring in heterogeneous immobile domains which interact with an explicit representation of the fracture network. Avoiding a costly discretization of the matrix, this method is well suited for running parametric studies and Monte Carlo simulations.

ABSTRACT. We present a model-based investigation of the transport of tracers and pesticides through fractured clayey tills. To characterize flow and transport properties, we conducted well-controlled flow-through experiments in large undisturbed columns (LUC, diameter 0.5 m, height 0.5 m, [1]) of fractured clayey tills from two field sites in eastern Denmark. Hydraulic tests yielded matrix hydraulic conductivity and porosity values and allowed the determination of fracture hydraulic apertures. Solutes with different sorption and diffusion characteristics were injected in the LUCs at a constant and a variable flow rate: bromide as a conservative solute and tebuconazole, a pesticide that sorbs to clayey tills. Finally, preferential flow paths were made visible by injection of the dye tracer Brilliant Blue followed by segmentation of the columns.

A detailed 3D discrete-fracture model comprising the observed preferential flow through fractures and macropores, the interaction with the clayey-till matrix (advection, diffusion, sorption) and the physical setup of the LUC was used to interpret the experiments and to shed light on the transport behavior in the two columns. Different transport mechanisms prevailed; the first experiment with clayey till from about 2 m bgs was dominated by fast advective transport through larger fractures and macropores with apertures of about ~100 µm, whereas the transport in the second column from 5 m bgs was happening much slower through smaller fractures (~10 µm) and by advection and diffusion through the matrix being more prevalent. The model could reproduce the breakthrough curves from the flow-through experiments. The sorption behavior of tebuconazole to the clayey tills could be captured by introducing non-equilibrium sorption description kinetics in the governing transport equation. Furthermore, the flow and transport parameters determined by the experiments and the calibrated 3D models served as basis for a modeling investigation, where the properties of the LUCs were applied to larger vertical cross sections and the leaching behavior of conservative and sorbing solutes/pesticides was analyzed under the influence of depth-varying fracture properties. This allowed us to study possible bottlenecks for vertical solute transport at the considered field sites.

References [1] P.R. Jørgensen, K. Mosthaf and M. Rolle (2019). A Large Undisturbed Column Method to Study Flow and Transport in Macropores and Fractured Media. Groundwater, 57, 951–961, (2019).

ABSTRACT. Available field data has shown interference between parallel, hydraulically-fractured horizontal wells in gas reservoirs, which indicates potential deviations with the prevailing single-well simulations. Three-dimensional, multi-phase, multi-well, multi-fracture simulation of gas production from unconventional reservoirs are performed to investigate the interference mechanism as functions of fracture characteristics. Realistic well configuration, reservoir condition and production data are used in the modeling scenarios. Results will provide insights on the significance of multi-well interference in practical productions.

ABSTRACT. This work presents a framework to analyze geometrical, topological and hydraulic properties of 3D Discrete Fracture Networks (DFN). A set of efficient algorithms have been developed to perform geometrical and topological analysis upon 3D networks of planar fractures with various shapes (mainly circular and ellipsoidal fractures). The present set of algorithms is capable of (i) calculating all possible intersections in the 3D networks, as well as the resulting trace lengths, the fracture areas inside a parallelepipedic domain, and other geometrical attributes of the DFN; (ii) extracting the percolating clusters and eliminating dead end clusters; (iii) constructing the corresponding graph of the 3D network of planar fractures; and (iv) solving the 3D flow on the corresponding DFN graph (case of a fractured rock with impervious rock matrix). All the calculations implemented in the algorithms have been strictly validated by direct numerical simulations, which increase the confidence in our algorithm methodology. The innovations of the present approach concern mainly (a) the adaptation of the multiple labelling techniques for the search of clusters to the case of 3D DFN’s of planar fractures, and (b) the use of efficient algorithms to eliminate dead end clusters in DFN’s. The graph approach can be used to simulate the flow in the DFN by solving the Laplacian of the corresponding graph of the DFN (algebraic approach with rectangular matrices). This method enables, based on the detailed flow calculation, to obtain the upscaled permeability of the DFN in a computationally efficient way, which avoids the direct numerical simulation of flow on the DFN. Although the graph approach does not give exact equivalent permeability of the DFN, it enables a reliable estimation of the permeability while gaining several orders of magnitudes in terms of CPU time. Different options for the graph representation (nodes and links) are considered, with the nodes centered either on the planar fractures, or on the intersecting traces. Appropriate transmissivity measures between connected nodes are analyzed. Conclusions are drawn about the physical relevance of each approach, and their accuracy in estimating the equivalent permeability of the DFN.

ABSTRACT. In the research to be presented, we express concentration or flux of solutes as a distribution over their velocity. We then derive an integrodifferential equation that governs the plume evolution and multi-component reactions solute distribution over velocity at given times and locations for a particle ensemble, based on a presumed velocity correlation structure and an ergodic cross-sectional velocity distribution. This way, the spatial evolution of breakthrough curves away from the source and the effective reaction rates are predicted based on cross-sectional velocity distribution and the connectivity, which is expressed by the velocity transition probability density. The transition probability can be specified via a copula function that can help construct a joint distribution with a given correlation and given marginal velocities or or in the special case of Gaussian copula, via an Ornstein-Uhlenbeck process. Using this approach, we analyze the breakthrough curves depending on the velocity distribution and correlation properties.

ABSTRACT. We develop efficient algorithms based on micro-macro models to simulate mineral dissolution and precipitation processes, analyzing potentially degenerating bulk properties of the medium such as porosity, diffusivity and permeability. Our model comprises transport equations at the scale of the porous medium (macroscale) including convection, diffusion and reaction. They feature averaged time- and space-dependent coefficients explicitly computed by means of auxiliary cell problems (microscale). We validated our approach against the dissolution of an array of dolomite grains in the micro-macro context.

ABSTRACT. Simulations with well-calibrated models offer a unique way to predict the multifaceted behavior of environmental surface or subsurface systems. Due to the lack of available data and high computational costs of the numerical simulation, this class of problems is still very challenging for uncertainty quantification. However, prediction uncertainty must be quantified through stochastic simulations and parameter inference. We offer varies strategies (conventional, sparse and adaptive) based on arbitrary polynomial chaos to quantify uncertainty in environment systems incorporating the observation data.

ABSTRACT. Localizing groundwater divides helps to define catchment boundaries and groundwater protection zones. Usually, groundwater divides can be derived from hydraulic head observations; however, the installation of observation wells is usually financially limited so that the localization is highly uncertain. We present a methodology to determine optimal monitoring strategies that minimize the uncertainty in locating groundwater divides with a limited number of observation wells, with application to a real catchment.

ABSTRACT. This work provides a prediction-driven uncertainty quantification framework for local decisions in groundwater management. We first design a multi-scale model parameterization through Local Principal Component Analysis, a new embedding method that generated consistency in the uncertainty updates across scales. Then we develop a comprehensive understanding of how model variables in different scales impact the local prediction from global sensitivity analysis over all scales. This understanding guides the design of a prediction-driven uncertainty reduction method using Bayesian methods. Our proposed framework is applied to quantify the uncertainty of local well catchment predictions, a real case study in Denmark.

ABSTRACT. There is an increased interest in managed aquifer recharge projects in California after the implementation of SGMA. Understanding and modeling subsurface heterogeneity is essential for determining areas of efficient recharge. The goal of this project is to determine how to statistically account for the stochastic nature of geostatistical modeling. Using 100 generated geologic realizations of our study site, statistical analysis of preliminary recharge results will be used to determine confidence in results of these preliminary recharge calculations.

ABSTRACT. In this work we present a sparse approximate Bayesian inference method for model inversion of partial differential equation (PDE) models with heterogeneous parameters. In our approach we construct a probabilistic representation of model parameters in terms of "pilot" values of said parameters evaluated at a finite set of "pilot" points. Inference is performed via variational inference algorithms for empirical Bayes. The proposed method provides an accurate and cost-effective alternatives to Markov Chain Monte Carlo simulation for model inversion.

ABSTRACT. We consider the problem of inferring the basal sliding coefficient field for an uncertain Stokes ice sheet forward model from surface velocity measurements. The uncertainty in the forward model stems from unknown (or uncertain) auxiliary parameters (e.g., rheology parameters). This inverse problem is posed within the Bayesian framework, which provides a systematic means of quantifying uncertainty in the solution. To account for the associated model uncertainty (error), we employ the Bayesian Approximation Error (BAE) approach to approximately premarginalize simultaneously over both the noise in measurements and uncertainty in the forward model. We also carry out approximative posterior uncertainty quantification based on a linearization of the parameter-to-observable map centered at the maximum a posteriori (MAP) basal sliding coefficient estimate, i.e., by taking the Laplace approximation. The MAP estimate is found by minimizing the negative log posterior using an inexact Newton conjugate gradient method. The gradient and Hessian actions to vectors are efficiently computed using adjoints. Sampling from the approximate covariance is made tractable by invoking a low-rank approximation of the data misfit component of the Hessian. We study the performance of the BAE approach in the context of three numerical examples in two and three dimensions. For each example the basal sliding coefficient field is the parameter of primary interest, which we seek to infer, and the rheology parameters (e.g., the flow rate factor, or the Glen’s flow law exponent coefficient field) represent so-called nuisance (secondary uncertain) parameters. Our results indicate that accounting for model uncertainty stemming from the presence of nuisance parameters is crucial. Namely our findings suggest that using nominal values for these parameters, as is often done in practice, without taking into account the resulting modeling error, can lead to overconfident and heavily biased results. We also show that the BAE approach can be used to account for the additional model uncertainty at no additional cost at the online stage.

ABSTRACT. Simulation of underground fluid flow is challenging because it involves solving a multi-physics problem in which multiphase flow and transport are tightly coupled. In this work, we develop an algebraic framework based on multigrid reduction (MGR) to construct preconditioners for the discrete linear systems resulting from a fully implicit (monolithic) approach for problems in porous and fractured media. We demonstrate the applicability of the framework when used as a “black-box” solver, and show that the preconditioners are efficient and scalable.

ABSTRACT. Crosshole GPR full-waveform inversion (FWI) has shown a high potential to characterize the near surface at a decimeter scale which is crucial for flow and transport. GPR FWI provide high-resolution tomograms of dielectric permittivity and electrical conductivity, which can be linked lithological properties. This study tests the potential of time-lapse GPR FWI to monitor tracers of different geophysical properties (salt, heat, ethanol). Synthetic and preliminary field results show that both properties can resolve major transport processes.

ABSTRACT. A common discretization strategy for coupled poromechanics is a continuous finite element scheme for the momentum balance and a finite volume scheme for the mass balance equations. When applied within a fully-implicit solution strategy this discretization is not intrinsically stable. We propose a new stabilization technique in which the balance of mass equations are supplemented with stabilizing flux terms. Results demonstrate that the stabilization treats spurious pressure oscillations and prevents the degradation in the convergence rate of iterative linear solvers.

ABSTRACT. Exchange across the interfaces separating atmospheric free-flows and multiphase porous medium flows are sensitive to the near interface flow dynamics in both flow regimes, as well as to characteristics of the interface that separates them. Experimental and numerical evaluations of wind speeds, interfacial forms, and heat fluxes each affect the exchange between these flow domains. Further, data-driven model expansions including sub-scale processes at the REV scale, and model comparisons including physical processes for large scale evaluation are investigated.

ABSTRACT. We present our results regarding the optimization of the fixed-stress splitting scheme applied to the quasi-static linear Biot equations. Both numerical schemes for finding the optimal stabilization term and theoretical convergence properties will be discussed.

ABSTRACT. I was one of the co-investigators in The Center for Geological Storage of CO2, an Energy Frontier Research Center at the University of Illinois from 2014 to 2019. With prior expertise in pore-scale direct numerical simulation of single phase flow, and having studied the Digital Rocks Physics literature, I was confident that my research team could use micro-CT scans of rock samples from a field injection site along with lattice Boltzmann models of two-phase flow to readily compute core-scale relative permeability and capillary pressure relations. At the time, the University of Illinois was home to one of the most powerful supercomputers hosted at an academic institution. We had a talented group of graduate students and post-docs, plus another team doing cutting edge microfluidics experiments for validation. We had a real field site with active CO2 injection, monitoring and sampling. What could go wrong?

We learned that there were challenges at every step: technical limits and artifacts in the xray scans, the highly heterogeneous Mt. Simon sandstone, daunting computational requirements and numerical instabilities in conventional lattice Boltzmann models, unforeseen complications due to the small viscosity of supercritical CO2, and others. In this talk I will give an overview of these challenges, and discuss development, validation and applications of a lattice Boltzmann code. I will also report some results showing that core-scale relative permeability curves computed by computationally simple porenetwork models are often comparable to those computed by lattice Boltzmann models. While still a neophyte with only six years of experience in multiphase pore scale simulation, I will offer some perspectives on future opportunities.

ABSTRACT. Many subsurface flow applications involve components where physical laws are well understood and other components where the physical laws are either poorly understood or not applicable. Numerical modeling excels at the former whereas interpolating data with machine learning (ML) excels at the latter, but neither approach can tackle these components simultaneously. Existing ML ap- proaches (often called physics-informed ML, or PIML) to handling these types of components simultaneously are minor tweaks to standard ML methods (e.g., PIML might use physics data to train or a loss function that encourages the ML to obey an equation without any accuracy guarantees). Tweaking black-box ML models is fundamentally limited because "big data does not interpret itself"–meaningful, interpretable structure in models is a necessity to improve predictability, enable human understanding, and maximize the impact of small data. We show how Differentiable Programming (DP) enables us to to meld trustworthy numerical modeling with trainable ML to produce fast models that can thrive on small data. This will be illustrated with an open-source, differentiable subsurface flow simulation, DPFEHM, that is available at https://github.com/OrchardLANL/DPFEHM.jl

ABSTRACT. Permeability levels in aquifers are often computed via inverse analysis from hydraulic head measurements. This can require discrete optimization, a task ill-suited to classical computers. In this talk we show how a widely accessible cloud-based quantum computer, called a quantum annealer, can efficiently perform discrete optimization for the hydrological inverse analysis. We give examples for 1D and 2D aquifers, compare with classical software, and discuss the technical challenges of working with quantum computing hardware and software.

ABSTRACT. In this work, the d3f software is used for numerically solving the PDEs describing the Elder problem. We adapted this software to the conditions of a Spark cluster to implement the mass parallel runs of d3f software and efficient further analysis of results using modern Big Data technologies. Having this powerful tool for numerical analysis of PDEs, we achieved the following scientific goals. 1. Investigate the steady-state solutions of the Elder problem with regards to different parameters. 2. Quantify the uncertainty of the investigation. 3. Develop predictive models for forecasting system’s steady states basing on early-time observations.

ABSTRACT. Oil and gas recovery rates from unconventional reservoirs are very low (<10% of resources are extracted). The current industry practices for reservoir development and management are generally ad-hoc and are mainly based on field experience. In addition, the physics processes related to hydrocarbons storage and recovery from unconventional reservoirs are not well understood. Approaches based machine learning are currently actively developed and applied to better understand the subsurface processes and enhance oil production. Here we present machine learning analyses based on synthetic and real-world datasets representing oil and gas production from unconventional reservoirs. The analyses are based on recently developed unsupervised and physics-informed machine learning method. Unsupervised methods utilize novel matrix and tensor factorization techniques. In the more general case of tensors, the factorization of a given tensor (high-dimensional array) is typically performed by minimization of the Frobenius norm representing the discrepancies between the original tensor and its approximation. Using Tucker decomposition, the original tensor is reconstructed into a mixing “core” tensor and a series of “feature” factors (vectors or matrices; see the figure). The size of “core” tensor defines the number of extracted features (signals) in each of the tensor dimensions. The factor matrices represent the extracted features (signals) in each of the tensor dimensions. Sparsity constraints on the elements of “core” tensor reduce the number of features and their mixing. Nonnegativity enforces parts-based representation of the original data which also allows the results to be easily interrelated. An algorithm called NTFk for tensor deconstruction which uses custom k-means clustering has been recently developed (https://github.com/orgs/TensorDecompositions). The machine learning methods based on NTFk are capable to reveal the temporal and spatial hidden (latent) features associated with the physical processes embedded in the analyzed datasets. The NTFk machine learning analyses of synthetic and real-world datasets related to oil and gas production demonstrated the applicability of the developed methodology to extract features characterizing differences in the obtained oil and gas extraction rates at different production wells. Physics-informed machine learning methods directly embed physics constraints, laws and simulators in the deep neural networks. Traditional deep neural networks have no information about the analyzed problem before the training is initiated. Physics-informed neural networks contain preconceived knowledge about the processes governing the analyzed problem. We have applied both unsupervised and physics-informed machine learning method to train, validate and verify machine learning models predicting oil and gas production from unconventional reservoirs in Texas.

ABSTRACT. An upscaling-based Bayesian geostatistical approach is developed for solving large-dimensional inverse problems. Instead of estimating the unknown parameter field, the approach estimates the truncated principal component projections. Most computations are implemented on an upscaled coarse grid, significantly saving computational costs.

ABSTRACT. Solving efficiently the large-size linearized systems of equations, that stem from the numerical modelling of multiphase flow by using the Mixed Hybrid Finite Element Method, is still a major computational challenge. In this work, a specifically designed preconditioning framework, founded on an original family of preconditioners, has been developed with the aim of improving the performance of Krylov subspace solvers, by reducing the total solution time and limiting the memory requirement.