Search

We study in detail the pressure stabilizing effects of the non-iterated fixed-stress splitting in poromechanical problems which are nearly undrained and incompressible. When applied in conjunction with a spatial discretization which does not satisfy the discrete inf-sup condition, namely a mixed piecewise linear - piecewise constant spatial discretization, the explicit fixed-stress scheme can have a pressure stabilizing effect in transient problems. This effect disappears, however, upon time step refinement or the attainment of steady state. The interpretation of the scheme as an Augmented Lagrangian method similar to Uzawa iteration for incompressible flow helps explain these results. Moreover, due to the slowly evolving solution within undrained seal regions, we show that the explicit fixed-stress scheme requires very large time steps to reveal its pressure stabilizing effect in examples of geologic CO$_2$ sequestration. We note that large time steps can result in large errors in drained regions, such as the aquifer or reservoir regions of these examples, and can prevent convergence of nonlinear solvers in the case of multiphase flows, which can make the explicit scheme an unreliable source of pressure stabilization. We conclude by demonstrating that pressure jump stabilization is as effective in the explicit fixed-stress setting as in the fully implicit setting for undrained problems, while maintaining the stability and convergence of the fixed-stress split for drained problems.

The interaction of multiple fluids through a heterogeneous pore space leads to complex pore-scale flow dynamics, such as intermittent pathway flow. The non-local nature of these dynamics, and the size of the 4D datasets acquired to capture them, presents challenges in identifying key fluctuations controlling fluid connectivity. To address these challenges, this work employs Dynamic Mode Decomposition (DMD), a data-driven algorithm that decomposes complex nonlinear systems into dominant spatio-temporal structures without relying on prior system assumptions. We present a workflow that identifies critical spatio-temporal regions exhibiting intermittent flow dynamics. This workflow is validated through three test cases, each exploring the impact of viscosity ratio on flow dynamics while maintaining a constant capillary number. Our findings demonstrate DMD's potential in analyzing extensive experimental datasets and identifying crucial intermittent flow structures, offering a powerful tool for understanding complex fluid behaviors in heterogeneous pore spaces. Using our method, we can quickly identify the timescales and locations of interest in an objective manner, providing a valuable diagnostic tool for analysing large synchrotron datasets.

We consider the numerical behavior of the fixed-stress splitting method for coupled poromechanics as undrained regimes are approached. We explain that pressure stability is related to the splitting error of the scheme, not the fact that the discrete saddle point matrix never appears in the fixed-stress approach. This observation reconciles previous results regarding the pressure stability of the splitting method. Using examples of compositional poromechanics with application to geological CO$_2$ sequestration, we see that solutions obtained using the fixed-stress scheme with a low order finite element-finite volume discretization which is not inherently inf-sup stable can exhibit the same pressure oscillations obtained with the corresponding fully implicit scheme. Moreover, pressure jump stabilization can effectively remove these spurious oscillations in the fixed-stress setting, while also improving the efficiency of the scheme in terms of the number of iterations required at every time step to reach convergence.

The injection of CO$_2$ into underground reservoirs provides a long term solution for anthropogenic emissions. A variable injection method (such as ramping the flow rate up or down) provides flexibility to injection sites, and could increase trapping at the pore-scale. However, the impact of a variable injection method on the connectivity of the gas, and subsequent trapping has not been explored at the pore-scale. Here, we conduct pore-scale imaging in a continuum-scale sample to observe the role of a variable flow rate on residual trapping. We show that the injection method influences how much of the pore space is accessible to the gas, even when total volumes injected, and total flow rates remain constant. Starting at a high flow rate, then decreasing it, leads to a larger amount of the pore space accessed by the gas. Conversely, starting at a low flow rate, and increasing it, leads to a larger role of heterogeneity of the pore space. This can promote trapping efficiency because channelling of the two fluids can occur, but less gas is trapped overall. Overall, a high-to-low injection scenario is optimum for residual trapping in the pore space due to increases in pore space accessibility.

Deep-learning-based surrogate models provide an efficient complement to numerical simulations for subsurface flow problems such as CO$_2$ geological storage. Accurately capturing the impact of faults on CO$_2$ plume migration remains a challenge for many existing deep learning surrogate models based on Convolutional Neural Networks (CNNs) or Neural Operators. We address this challenge with a graph-based neural model leveraging recent developments in the field of Graph Neural Networks (GNNs). Our model combines graph-based convolution Long-Short-Term-Memory (GConvLSTM) with a one-step GNN model, MeshGraphNet (MGN), to operate on complex unstructured meshes and limit temporal error accumulation. We demonstrate that our approach can accurately predict the temporal evolution of gas saturation and pore pressure in a synthetic reservoir with impermeable faults. Our results exhibit a better accuracy and a reduced temporal error accumulation compared to the standard MGN model. We also show the excellent generalizability of our algorithm to mesh configurations, boundary conditions, and heterogeneous permeability fields not included in the training set. This work highlights the potential of GNN-based methods to accurately and rapidly model subsurface flow with complex faults and fractures.

Successful deployment of geological carbon storage (GCS) requires an extensive use of reservoir simulators for screening, ranking and optimization of storage sites. However, the time scales of GCS are such that no sufficient long-term data is available yet to validate the simulators against. As a consequence, there is currently no solid basis for assessing the quality with which the dynamics of large-scale GCS operations can be forecasted. To meet this knowledge gap, we have conducted a major GCS validation benchmark study. To achieve reasonable time scales, a laboratory-size geological storage formation was constructed (the "FluidFlower"), forming the basis for both the experimental and computational work. A validation experiment consisting of repeated GCS operations was conducted in the FluidFlower, providing what we define as the true physical dynamics for this system. Nine different research groups from around the world provided forecasts, both individually and collaboratively, based on a detailed physical and petrophysical characterization of the FluidFlower sands. The major contribution of this paper is a report and discussion of the results of the validation benchmark study, complemented by a description of the benchmarking process and the participating computational models. The forecasts from the participating groups are compared to each other and to the experimental data by means of various indicative qualitative and quantitative measures. By this, we provide a detailed assessment of the capabilities of reservoir simulators and their users to capture both the injection and post-injection dynamics of the GCS operations.

We present a Parametrization of the Physics Informed Neural Network (P-PINN) approach to tackle the problem of uncertainty quantification in reservoir engineering problems. We demonstrate the approach with the immiscible two phase flow displacement (Buckley-Leverett problem) in heterogeneous porous medium. The reservoir properties (porosity, permeability) are treated as random variables. The distribution of these properties can affect dynamic properties such as the fluids saturation, front propagation speed or breakthrough time. We explore and use to our advantage the ability of networks to interpolate complex high dimensional functions. We observe that the additional dimensions resulting from a stochastic treatment of the partial differential equations tend to produce smoother solutions on quantities of interest (distributions parameters) which is shown to improve the performance of PINNS. We show that provided a proper parameterization of the uncertainty space, PINN can produce solutions that match closely both the ensemble realizations and the stochastic moments. We demonstrate applications for both homogeneous and heterogeneous fields of properties. We are able to solve problems that can be challenging for classical methods. This approach gives rise to trained models that are both more robust to variations in the input space and can compete in performance with traditional stochastic sampling methods.

This work focuses on the development of a two-step field-split nonlinear preconditioner to accelerate the convergence of two-phase flow and transport in heterogeneous porous media. We propose a field-split algorithm named Field-Split Multiplicative Schwarz Newton (FSMSN), consisting in two steps: first, we apply a preconditioning step to update pressure and saturations nonlinearly by solving approximately two subproblems in a sequential fashion; then, we apply a global step relying on a Newton update obtained by linearizing the system at the preconditioned state. Using challenging test cases, FSMSN is compared to an existing field-split preconditioner, Multiplicative Schwarz Preconditioned for Inexact Newton (MSPIN), and to standard solution strategies such as the Sequential Fully Implicit (SFI) method or the Fully Implicit Method (FIM). The comparison highlights the impact of the upwinding scheme in the algorithmic performance of the preconditioners and the importance of the dynamic adaptation of the subproblem tolerance in the preconditioning step. Our results demonstrate that the two-step nonlinear preconditioning approach-and in particular, FSMSN-results in a faster outer-loop convergence than with the SFI and FIM methods. The impact of the preconditioners on computational performance-i.e., measured by wall-clock time-will be studied in a subsequent publication.

The "Workshop on Machine learning in heterogeneous porous materials" brought together international scientific communities of applied mathematics, porous media, and material sciences with experts in the areas of heterogeneous materials, machine learning (ML) and applied mathematics to identify how ML can advance materials research. Within the scope of ML and materials research, the goal of the workshop was to discuss the state-of-the-art in each community, promote crosstalk and accelerate multi-disciplinary collaborative research, and identify challenges and opportunities. As the end result, four topic areas were identified: ML in predicting materials properties, and discovery and design of novel materials, ML in porous and fractured media and time-dependent phenomena, Multi-scale modeling in heterogeneous porous materials via ML, and Discovery of materials constitutive laws and new governing equations. This workshop was part of the AmeriMech Symposium series sponsored by the National Academies of Sciences, Engineering and Medicine and the U.S. National Committee on Theoretical and Applied Mechanics., Comment: The workshop link is: https://amerimech.mech.utah.edu

In subsurface multiphase flow simulations, poor nonlinear solver performance is a significant runtime sink. The system of fully implicit mass balance equations is highly nonlinear and often difficult to solve for the nonlinear solver, generally Newton(-Raphson). Strong nonlinearities can cause Newton iterations to converge very slowly. This frequently results in time step cuts, leading to computationally expensive simulations. Much literature has looked into how to improve the nonlinear solver through enhancements or safeguarding updates. In this work, we take a different approach; we aim to improve convergence with a smoother finite volume discretization scheme which is more suitable for the Newton solver. Building on recent work, we propose a novel total velocity hybrid upwinding scheme with weighted average flow mobilities (WA-HU TV) that is unconditionally monotone and extends to compositional multiphase simulations. Analyzing the solution space of a one-cell problem, we demonstrate the improved properties of the scheme and explain how it leverages the advantages of both phase potential upwinding and arithmetic averaging. This results in a flow subproblem that is smooth with respect to changes in the sign of phase fluxes, and is well-behaved when phase velocities are large or when co-current viscous forces dominate. Additionally, we propose a WA-HU scheme with a total mass (WA-HU TM) formulation that includes phase densities in the weighted averaging. The proposed WA-HU TV consistently outperforms existing schemes, yielding benefits from 5\% to over 50\% reduction in nonlinear iterations. The WA-HU TM scheme also shows promising results; in some cases leading to even more efficiency. However, WA-HU TM can occasionally also lead to convergence issues. Overall, based on the current results, we recommend the adoption of the WA-HU TV scheme as it is highly efficient and robust.

In this study, we develop a probabilistic approach to map the parametric uncertainty to the output state uncertainty in first-order hyperbolic conservation laws. We analyze this problem for nonlinear immiscible two-phase transport in heterogeneous porous media in the presence of a stochastic velocity field. The uncertainty in the velocity field can arise from the incomplete description of either porosity field, injection flux, or both. The uncertainty in the total-velocity field leads to the spatiotemporal uncertainty in the saturation field. Given information about the spatial/temporal statistics of the correlated heterogeneity, we leverage method of distributions to derive deterministic equations that govern the evolution of single-point CDF of saturation. Unlike Buckley Leverett equation, the equation for the raw CDF function is linear in space and time. Hereby, we give routes to circumventing the computational cost of Monte Carlo scheme while obtaining the full statistical description of saturation. We conduct a set of numerical experiments and compare statistics of saturation computed with the method of distributions, against those obtained using the statistical moment equations approach and kernel density estimation post-processing of high-resolution Monte Carlo simulations. This comparison demonstrates that the CDF equations remain accurate over a wide range of statistical properties, i.e. standard deviation and correlation length of the underlying random fields, while the corresponding low-order statistical moment equations significantly deviate from Monte Carlo results, unless for very small values of standard deviation and correlation length.

We present our progress on the application of physics informed deep learning to reservoir simulation problems. The model is a neural network that is jointly trained to respect governing physical laws and match boundary conditions. The methodology is hereby used to simulate a 2-phase immiscible transport problem (Buckley-Leverett). The model is able to produce an accurate physical solution both in terms of shock and rarefaction and honors the governing partial differential equation along with initial and boundary conditions. We test various hypothesis (uniform and non-uniform initial conditions) and show that with the proper implementation of physical constraints, a robust solution can be trained within a reasonable amount of time and iterations. We revisit some of the limitations presented in previous work \cite{Fuks2020} and further the applicability of this method in a forward, pure hyperbolic setup. We also share some practical findings on the application of physics informed neural networks (PINN). We review various network architectures presented in the literature and show tips that helped improve their convergence and accuracy. The proposed methodology is a simple and elegant way to instill physical knowledge to machine-learning algorithms. This alleviates the two most significant shortcomings of machine-learning algorithms: the requirement for large datasets and the reliability of extrapolation. The principles presented can be generalized in innumerable ways in the future and should lead to a new class of algorithms to solve both forward and inverse physical problems., Comment: 19 pages, 17 figures. arXiv admin note: text overlap with arXiv:2001.05172 Thanks to Xin Ju for review and comments

We describe a second-order accurate approach to sparsifying the off-diagonal blocks in the hierarchical approximate factorizations of sparse symmetric positive definite matrices. The norm of the error made by the new approach depends quadratically, not linearly, on the error in the low-rank approximation of the given block. The analysis of the resulting two-level preconditioner shows that the preconditioner is second-order accurate as well. We incorporate the new approach into the recent Sparsified Nested Dissection algorithm [SIAM J. Matrix Anal. Appl., 41 (2020), pp. 715-746], and test it on a wide range of problems. The new approach halves the number of Conjugate Gradient iterations needed for convergence, with almost the same factorization complexity, improving the total runtimes of the algorithm. Our approach can be incorporated into other rank-structured methods for solving sparse linear systems., Comment: 26 pages

We propose MeshfreeFlowNet, a novel deep learning-based super-resolution framework to generate continuous (grid-free) spatio-temporal solutions from the low-resolution inputs. While being computationally efficient, MeshfreeFlowNet accurately recovers the fine-scale quantities of interest. MeshfreeFlowNet allows for: (i) the output to be sampled at all spatio-temporal resolutions, (ii) a set of Partial Differential Equation (PDE) constraints to be imposed, and (iii) training on fixed-size inputs on arbitrarily sized spatio-temporal domains owing to its fully convolutional encoder. We empirically study the performance of MeshfreeFlowNet on the task of super-resolution of turbulent flows in the Rayleigh-Benard convection problem. Across a diverse set of evaluation metrics, we show that MeshfreeFlowNet significantly outperforms existing baselines. Furthermore, we provide a large scale implementation of MeshfreeFlowNet and show that it efficiently scales across large clusters, achieving 96.80% scaling efficiency on up to 128 GPUs and a training time of less than 4 minutes., Comment: Supplementary Video: https://youtu.be/mjqwPch9gDo. Accepted to SC20

In this work, the Minkowski functionals are used as a framework to study how morphology (i.e. the shape of a structure) and topology (i.e. how different structures are connected) influence wall adsorption and capillary condensation under tight confinement. Numerical simulations based on classical density functional theory (DFT) are run for a wide variety of geometries using both hard-sphere and Lennard-Jones fluids. These DFT computations are compared to results obtained using the Minkowski functionals. It is found that the Minkowski functionals can provide a good description of the behavior of Lennard-Jones fluids down to small system sizes. In addition, through decomposition of the free energy, the Minkowski functionals provide a good framework to better understand what are the dominant contributions to the physics of a system. Lastly, while studying the phase envelope shift as a function of the Minkowski functionals it is found that topology has a different effect depending on whether the phase transition under consideration is a first- or a second-order transition.

Accurate numerical simulation of coupled fracture/fault deformation and fluid flow is crucial to the performance and safety assessment of many subsurface systems. In this work, we consider the discretization and enforcement of contact conditions at such surfaces. The bulk rock deformation is simulated using low-order continuous finite elements, while frictional contact conditions are imposed by means of a Lagrange multiplier method. We employ a cell-centered finite-volume scheme to solve the fracture fluid mass balance equation. From a modeling perspective, a convenient choice is to use a single grid for both mechanical and flow processes, with piecewise-constant interpolation of Lagrange multipliers, i.e., contact tractions and fluid pressure. Unfortunately, this combination of displacement and multiplier variables is not uniformly inf-sup stable, and therefore requires a stabilization technique. Starting from a macroelement analysis, we develop two algebraic stabilization approaches and compare them in terms of robustness and convergence rate. The proposed approaches are validated against challenging analytical two- and three-dimensional benchmarks to demonstrate accuracy and robustness. These benchmarks include both pure contact mechanics problems and well as problems with tightly-coupled fracture flow.

We present a new hybrid physics-based machine-learning approach to reservoir modeling. The methodology relies on a series of deep adversarial neural network architecture with physics-based regularization. The network is used to simulate the dynamic behavior of physical quantities (i.e. saturation) subject to a set of governing laws (e.g. mass conservation) and corresponding boundary and initial conditions. A residual equation is formed from the governing partial-differential equation and used as part of the training. Derivatives of the estimated physical quantities are computed using automatic differentiation algorithms. This allows the model to avoid overfitting, by reducing the variance and permits extrapolation beyond the range of the training data including uncertainty implicitely derived from the distribution output of the generative adversarial networks. The approach is used to simulate a 2 phase immiscible transport problem (Buckley Leverett). From a very limited dataset, the model learns the parameters of the governing equation and is able to provide an accurate physical solution, both in terms of shock and rarefaction. We demonstrate how this method can be applied in the context of a forward simulation for continuous problems. The use of these models for the inverse problem is also presented, where the model simultaneously learns the physical laws and determines key uncertainty subsurface parameters. The proposed methodology is a simple and elegant way to instill physical knowledge to machine-learning algorithms. This alleviates the two most significant shortcomings of machine-learning algorithms: the requirement for large datasets and the reliability of extrapolation. The principles presented in this paper can be generalized in innumerable ways in the future and should lead to a new class of algorithms to solve both forward and inverse physical problems., Comment: 21 pages, 15 figures

In numerical simulations of multiphase flow and transport in fractured porous media, the estimation of the hydrocarbon recovery requires accurately predicting the capillary-driven imbibition rate of the wetting phase initially present in the fracture into the low-permeability matrix. In the fully implicit finite-volume scheme, this entails a robust methodology that captures the capillary flux at the interface between the matrix and the fracture even when very coarse control volumes are used to discretize the matrix. Here, we investigate the application of discrete interface conditions at the matrix-fracture interface to improve the accuracy of the flux computation without relying on extreme grid refinement. In particular, we study the interaction of the upwinding scheme with the discrete interface conditions. Considering first capillary-dominated spontaneous imbibition and then forced imbibition with viscous, buoyancy, and capillary forces, we illustrate the importance of the interface conditions to accurately capture the matrix-fracture flux and correctly represent the flow dynamics in the problem.

Many applications involving porous media--notably reservoir engineering and geologic applications--involve tight coupling between multiphase fluid flow, transport, and poromechanical deformation. While numerical models for these processes have become commonplace in research and industry, the poor scalability of existing solution algorithms has limited the size and resolution of models that may be practically solved. In this work, we propose a two-stage Newton-Krylov solution algorithm to address this shortfall. The proposed solver exhibits rapid convergence, good parallel scalability, and is robust in the presence of highly heterogeneous material properties. The key to success of the solver is a block-preconditioning strategy that breaks the fully-coupled system of mass and momentum balance equations into simpler sub-problems that may be readily addressed using targeted algebraic methods. Numerical results are presented to illustrate the performance of the solver on challenging benchmark problems.

The sequential fully implicit (SFI) method was introduced along with the development of the multiscale finite volume (MSFV) framework, and has received considerable attention in recent years. Each time step for SFI consists of an outer loop to solve the coupled system, in which there is one inner Newton loop to implicitly solve the pressure equation and another loop to implicitly solve the transport equations. Limited research has been conducted that deals with the outer coupling level to investigate the convergence performance. In this paper we extend the basic SFI method with several nonlinear acceleration techniques for improving the outer-loop convergence. Specifically, we consider numerical relaxation, quasi-Newton (QN) and Anderson acceleration (AA) methods. The acceleration techniques are adapted and studied for the first time within the context of SFI for coupled flow and transport in porous media. We reveal that the iterative form of SFI is equivalent to a nonlinear block Gauss-Seidel (BGS) process. The effectiveness of the acceleration techniques is demonstrated using several challenging examples. The results show that the basic SFI method is quite inefficient, suffering from slow convergence or even convergence failure. In order to better understand the behaviors of SFI, we carry out detailed analysis on the coupling mechanisms between the sub-problems. Compared with the basic SFI method, superior convergence performance is achieved by the acceleration techniques, which can resolve the convergence difficulties associated with various types of coupling effects. We show across a wide range of flow conditions that the acceleration techniques can stabilize the iterative process, and largely reduce the outer iteration count.

Conservation laws in the form of elliptic and parabolic partial differential equations (PDEs) are fundamental to the modeling of many problems such as heat transfer and flow in porous media. Many of such PDEs are stochastic due to the presence of uncertainty in the conductivity field. Based on the relation between stochastic diffusion processes and PDEs, Monte Carlo (MC) methods are available to solve these PDEs. These methods are especially relevant for cases where we are interested in the solution in a small subset of the domain. The existing MC methods based on the stochastic formulation require restrictively small time steps for high variance conductivity fields. Moreover, in many applications the conductivity is piecewise constant and the existing methods are not readily applicable in these cases. Here we provide an algorithm to solve one-dimensional elliptic problems that bypasses these two limitations. The methodology is demonstrated using problems governed by deterministic and stochastic PDEs. It is shown that the method provides an efficient alternative to compute the statistical moments of the solution to a stochastic PDE at any point in the domain. A variance reduction scheme is proposed for applying the method for efficient mean calculations.

In reservoir simulation, solution of the coupled systems of nonlinear algebraic equations that are associated with fully-implicit (backward Euler) discretization is challenging. Having a robust and efficient nonlinear solver is necessary in order for reservoir simulation to serve as the primary tool for managing the recovery processes of large-scale reservoirs. Here, we develop a continuation method based on the use of a dissipation operator. We focus on nonlinear two-phase flow and transport in heterogeneous formations in the presence of viscous, gravitational, and capillary forces. The homotopy is constructed by adding numerical dissipation to the coupled discrete conservation equations. A continuation parameter is introduced to control the amount of dissipation. Numerical evidence of multi-dimensional models and detailed analysis of single-cell problems are used to explain how the dissipation operator improves the nonlinear convergence of the coupled system of equations. An adaptive strategy to determine the dissipation coefficient is proposed. The dissipation level is computed locally for each cell interface. We demonstrate the efficiency of the dissipation-based continuation (DBC) nonlinear solver using several examples, including 1D scalar transport and 2D heterogeneous problems with fully-coupled flow and transport. The DBC solver has better convergence properties compared with the standard damped-Newton solvers used in reservoir simulation.

In the context of stochastic two-phase flow in porous media, we introduce a novel and efficient method to estimate the probability distribution of the wetting saturation field under uncertain rock properties in highly heterogeneous porous systems, where streamline patterns are dominated by permeability heterogeneity, and for slow displacement processes (viscosity ratio close to unity). Our method, referred to as the frozen streamline distribution method (FROST), is based on a physical understanding of the stochastic problem. Indeed, we identify key random fields that guide the wetting saturation variability, namely fluid particle times of flight and injection times. By comparing saturation statistics against full-physics Monte Carlo simulations, we illustrate how this simple, yet accurate FROST method performs under the preliminary approximation of frozen streamlines. Further, we inspect the performance of an accelerated FROST variant that relies on a simplification about injection time statistics. Finally, we introduce how quantiles of saturation can be efficiently computed within the FROST framework, hence leading to robust uncertainty assessment., Comment: 34 pages

We developed a sharp interface level-set approach for two-phase immiscible flow with moving contact lines. The Cox-Voinov model is used to describe the moving contact line. A piecewise linear interface method is used to construct the signed distance function and to implement the contact angle boundary condition. Spurious currents are studied in the case of static and moving fluid interfaces, which show convergence behavior. Pressure and the surface tension force are balanced up to machine precision for static parabolic interfaces, while the velocity error decreases steadily with grid refinement when the interface is advected in a uniform flow field. The moving contact line problem is studied and validated through comparison with theory and experiments for an advancing interface and capillary rise in a tube.

Monte Carlo (MC) simulations of transport in random porous networks indicate that for high variances of the log-normal permeability distribution, the transport of a passive tracer is non-Fickian. Here we model this non-Fickian dispersion in random porous networks using discrete temporal Markov models. We show that such temporal models capture the spreading behavior accurately. This is true despite the fact that the slow velocities are strongly correlated in time, and some studies have suggested that the persistence of low velocities would render the temporal Markovian model inapplicable. Compared to previously proposed temporal stochastic differential equations with case specific drift and diffusion terms, the models presented here require fewer modeling assumptions. Moreover, we show that discrete temporal Markov models can be used to represent dispersion in unstructured networks, which are widely used to model porous media. A new method is proposed to extend the state space of temporal Markov models to improve the model predictions in the presence of extremely low velocities in particle trajectories and extend the applicability of the model to higher temporal resolutions. Finally, it is shown that by combining multiple transitions, temporal models are more efficient for computing particle evolution compared to correlated CTRW with spatial increments that are equal to the lengths of the links in the network.

The generalized polynomial chaos method is applied to the Buckley-Leverett equation. We consider a spatially homogeneous domain modeled as a random field. The problem is projected onto stochastic basis functions which yields an extended system of partial differential equations. Analysis and numerical methods leading to reduced computational cost are presented for the extended system of equations. The accurate representation of the evolution of a discontinuous stochastic solution over time requires a large number of stochastic basis functions. Adaptivity of the stochastic basis to reduce computational cost is challenging in the stochastic Galerkin setting since the change of basis affects the system matrix itself. To achieve adaptivity without adding overhead by rewriting the entire system of equations for every grid cell, we devise a basis reduction method that distinguishes between locally significant and insignificant modes without changing the actual system matrices. Results are presented for problems in one and two spatial dimensions, with varying number of stochastic dimensions. We show how to obtain stochastic velocity fields from realistic permeability fields and demonstrate the performance of the stochastic Galerkin method with local basis reduction. The system of conservation laws is discretized with a finite volume method and we demonstrate numerical convergence to the reference solution obtained through Monte Carlo sampling.