## Abstract

In this paper, we present a second-order accurate adaptive algorithm for solving multi-phase, incompressible flow in porous media. We assume a multi-phase form of Darcy’s law with relative permeabilities given as a function of the phase saturation. The remaining equations express conservation of mass for the fluid constituents. In this setting, the total velocity, defined to be the sum of the phase velocities, is divergence free. The basic integration method is based on a total-velocity splitting approach in which we solve a second-order elliptic pressure equation to obtain a total velocity. This total velocity is then used to recast component conservation equations as nonlinear hyperbolic equations. Our approach to adaptive refinement uses a nested hierarchy of logically rectangular grids with simultaneous refinement of the grids in both space and time. The integration algorithm on the grid hierarchy is a recursive procedure in which coarse grids are advanced in time, fine grids are advanced multiple steps to reach the same time as the coarse grids and the data at different levels are then synchronized. The single-grid algorithm is described briefly, but the emphasis here is on the time-stepping procedure for the adaptive hierarchy. Numerical examples are presented to demonstrate the algorithm’s accuracy and convergence properties and to illustrate the behaviour of the method.

## 1. Introduction

Multi-component and multi-phase flows in the subsurface are often characterized by localized phenomena such as steep concentration gradients or saturation fronts. Accurately resolving these types of phenomena requires high resolution in regions where the solution is changing rapidly. For this reason, the development of some type of dynamic gridding capability has long been of interest in the porous media community.

Heinemann (1983) and Ewing *et al.* (1989) have considered dynamic local grid refinement approaches. More recent papers by, for example, Sammon (2003) and Christensen *et al.* (2004) discuss development of adaptive techniques in the context of unstructured grids. An alternative approach to local refinement is based on structured-grid adaptive mesh refinement (AMR). This type of approach, based on the strategy introduced for gas dynamics by Berger & Colella (1989), was first applied to porous media flow by Hornung & Trangenstein (1997) and by Propp (1998). Additional developments are discussed in the work by Trangenstein (2002), Trangenstein & Bi (2002) and Hoang & Kleppe (2006). In this approach, regions to be refined are uniformly subdivided in both space and time. Related approaches were developed by Nilsson *et al.* (2005*a*,*b*) who use a spatially anisotropic refinement strategy with no temporal refinement and by Edwards (1996) who uses a cell-by-cell refinement strategy.

The focus of this paper is on developing a structured-grid AMR algorithm for porous media flow. Our approach is similar to the approach introduced by Hornung & Trangenstein (1997) and to the approach discussed by Propp (1998). In these approaches, the solution is represented on a hierarchical sequence of nested grids with successively finer spacing in both time and space. Increasingly finer grids are recursively embedded in coarse grids until the solution is sufficiently resolved. An error estimation procedure based on user-specified criteria evaluates where additional refinement is needed and grid generation procedures dynamically create or remove rectangular fine-grid patches as resolution requirements change. The method presented here uses subcycling in time so that all levels are advanced at the same Courant–Friedrichs–Levy (CFL) number, thus reducing the numerical dissipation of the explicit upwind advection scheme used to advance the solution. The major difference between the approach adopted here and that of Hornung & Trangenstein (1997) is that the current method does not require a global, multi-level pressure solve at each fine-grid time step. Instead, when advancing a given level, we solve the pressure on that level only with boundary conditions obtained from the coarser levels. This approach avoids the computational expense of the global solve but introduces additional complexity into the synchronization step of the algorithm in which inconsistencies between different levels are corrected. The synchronization approach used here is based on the algorithm developed by Almgren *et al.* (1998) for incompressible flows. The methodology is implemented in parallel using the BoxLib framework discussed in Rendleman *et al.* (2000) and Crutchfield (1991). A similar approach was used by Propp (1998); however, his algorithm used a different approach to synchronization and was limited to two dimensions.

Before describing the adaptive algorithm we will briefly review the total velocity splitting approach and discuss our basic fractional step scheme for a single grid. In §3, we describe, in detail, the recursive time-stepping procedure for the adaptive algorithm and other aspects of the adaptive algorithm. Section 4 shows the convergence results and presents computational examples illustrating both the performance and parallel scaling of the method.

## 2. Total velocity-splitting algorithm

Here we consider multi-phase, multi-component incompressible flow in heterogeneous porous media. The multi-component mixture is composed of *N* components (or lumped pseudo-components). We define **n**≡(*n*_{1},…,*n*_{N}) as the vector of component densities per unit pore volume, and let **n**_{α} represent the portion of **n** in phase *α*, where Greek subscripts refer to mobile phases. Thus, is the total component density in the combined fluid system. In a general setting, the separation of components into phases requires a computation of thermodynamic equilibrium for the mixture. For the problems considered here, the phases are incompressible, we assume no volume change on mixing and each component only appears in a single phase. The basic void fraction of the medium is referred to as the porosity, *ϕ*, and the fraction of that void occupied by a particular phase is referred to as the phase saturation, *u*_{α}, which, for the models considered here, can be determined from the composition of the phase, **n**_{α}, and the pure component densities.

The equations that describe the flow represent mass conservation and Darcy’s law. Darcy’s law expresses the volumetric flow rate, *v*_{α}, of each phase in terms of the phase pressure, *p*_{α}, namely,
2.1
where *λ*_{α}≡*K**k*_{r,α}/*μ*_{α} is the phase mobility. Here *K* is the permeability distribution of the medium; *k*_{r,α} is the relative permeability, which is a function of *u*_{α} that expresses modification of the flow rate due to multiphase effects; *μ*_{α} and *ρ*_{α} are the phase viscosity and phase density, respectively, which depend on the phase composition; and ** g** is the gravitational acceleration. The pressure in each phase is related to a reference pressure,

*p*, by a capillary pressure,

*p*

_{c,α}=

*p*

_{α}−

*p*, which is a function of saturation.

For this system, conservation of mass for each component is given by
2.2
where represents diffusive terms that include multi-phase molecular diffusion and dispersion and is the reaction source term. The diffusive term associated with capillary pressure is implicitly included in the definition of the phase velocities. The diffusion terms can be quite complex, depending on the particular problem; see the work of Xu *et al.* (2004) for a detailed discussion of these terms. For purposes of exposition, we will set and for the remainder of the discussion. However, we note that within this framework additional diffusion terms can be added analogously to the capillary pressure term. If reactions are included in the system, they can be included using the operator-split formalism discussed in Day & Bell (2000).

The component conservation equations specify the change in total mass of each component resulting from the transport and diffusion of that component distributed across the phases. The phase behaviour of the system specifies how the components are apportioned into phases and the volume occupied by those phases. In this paper, we consider only simple incompressible systems without mass transfer between phases. We can then define a total velocity
2.3
since *p*_{α}=*p*+*p*_{c,α}. The total velocity is divergence free so that
2.4
This leads to a second-order elliptic equation for pressure
2.5
We can rewrite equation (2.3) to express ∇*p* in terms of *v*_{T}, then using equation (2.1) we can express *v*_{α} in terms of the total velocity,
where is the total mobility. Writing the conservation equations in terms of the total velocity yields the fractional flow form of the component conservation equations given by
2.6
where
and

In the absence of diffusive terms, the total velocity defines a splitting that decomposes the dynamics into an elliptic pressure equation and a system of hyperbolic conservation laws. We note that we are assuming here that the system *ϕ*∂**n**/∂*t*+∇⋅*F*(**n**,*v*_{T})=0 is hyperbolic with *v*_{T} viewed as a function of space and time. This is a condition on the model specification; some models for three-phase flow use relative permeability models that lead to non-hyperbolic behaviour and the resulting problems are ill-posed in the zero-viscosity limit (Bell *et al*. 1986).

### (a) Single-grid algorithm

In this section, we discuss the single-grid algorithm, giving an overview of the time-stepping procedure. The discretization uses a volume of fluid approach in which denotes the average value of **n** over cell (*i*,*j*,*k*) at time *t*^{n}; **n** and *p* are defined on cell centres while *F* and *v*_{T} are defined on cell edges. The temporal discretization fits into the basic framework of implicit pressure–explicit saturation (IMPES)-type methods in which the pressure equation is solved implicitly and the saturation (component) equations are solved explicitly. Here the basic algorithm is modified so that the overall splitting approach is second-order accurate in time and we treat the diffusion terms semi-implicitly so that diffusive terms, *H*∇*P*_{c}, do not limit the time step.

The outline of the algorithm is as follows.

(i)

*Step 1.*Solve the pressure equation (2.5) rewritten in the form for*p*with properties evaluated using**n**^{n}. We then use equation (2.3) to define a total velocity . Here*D*and*G*are second-order accurate discretizations of the divergence and gradient operators, respectively. The divergence operator returns a cell-centred divergence from face-centred values; the gradient operator differences cell-centred values to return normal gradients on faces. For example, in two dimensions, the discretization of at cell (*i*,*j*) would be where Δ*x*and Δ*y*are the mesh spacings in the*x*and*y*directions, respectively.(ii)

*Step 2.*Use equation (2.6) to advance the solution from time*t*^{n}to time*t*^{n+1}using We use an unsplit second-order Godunov scheme to compute the hyperbolic fluxes using the methodology described in Bell*et al.*(1989). (See Almgren*et al.*(1998) for additional details of the multi-dimensional aspects of the advection scheme.) The Godunov discretization is coupled to a Crank–Nicolson discretization of the diffusive terms, so that 2.7 with a suitable linearization of the coefficients of the diffusion term; here we use*H*^{n+1,*}and to denote*H*and*P*_{c}, respectively, which are functions of**n**, evaluated at**n**^{n+1,*}. In this step,*F*^{n+1/2,*}denotes time-centred fluxes computed using the Godunov procedure but with the total velocity evaluated at*t*^{n}.(iii)

*Step 3.*Solve the pressure equation (2.5), with properties evaluated using**n**^{n+1,*}, to compute a new total velocity from equation (2.3). We then define so that we can time centre the dependence of the flux on*v*_{T}.(iv)

*Step 4.*Use equation (2.6) to re-advance the solution from time*t*^{n}to time*t*^{n+1}, this time using to obtain values of**n**^{n+1}, 2.8 again with a suitable linearization of the coefficients of the diffusion term. In this step,*F*^{n+1/2}denotes time-centred fluxes computed using the Godunov procedure but with the total velocity evaluated at*t*^{n+1/2}.

We note here that Step 1 is not always necessary; it can be acceptable to replace with the final *v*_{T} computed in Step 3 of the previous time step without significant loss of accuracy, thus allowing us to avoid the computational cost associated with Step 1. In this case, is computed with **n**^{n,*} instead of **n**^{n}. We shall demonstrate in §4 that for our examples this has no noticeable effect on the solution. For clarity, we shall denote a time step that uses Step 1 as a *full solve* and one without as a *partial solve*.

## 3. Adaptive mesh refinement

In this section, we present the extension of the algorithm described above to an adaptive hierarchy of nested rectangular grids. First we describe the creation of the grid hierarchy and the regridding procedure used to adjust the hierarchy during the computation. Next we discuss the adaptive time step algorithm with subcycling in time, focusing on the synchronization between different levels of refinement.

### (a) Creating and managing the grid hierarchy

The grid hierarchy is composed of different levels of refinement ranging from coarsest (ℓ=0) to finest (). Each level is represented as the union of rectangular grid patches of a given resolution. In this implementation, the refinement ratio is always even, with the same factor of refinement in each coordinate direction, i.e. Δ*x*^{ℓ+1}=Δ*y*^{ℓ+1}=Δ*z*^{ℓ+1}=Δ*x*^{ℓ}/*r*, where *r* is the refinement ratio. (We note here that neither isotropic refinement nor uniform base grids are requirements of the fundamental algorithm.) In the actual implementation, the refinement ratio, either 2 or 4, can be a function of level; however, in the exposition, we will assume that *r* is constant. The grids are properly nested, in the sense that the union of grids at level ℓ+1 is contained in the union of grids at level ℓ for . Furthermore, the containment is strict in the sense that, except at physical boundaries, the level ℓ grids are large enough to guarantee that there is a border at least one level ℓ cell wide surrounding each level ℓ+1 grid. (Grids at all levels are allowed to extend to the physical boundaries so the proper nesting is not strict there.)

The initial creation of the grid hierarchy and the subsequent regridding operations in which the grids are dynamically changed to reflect changing flow conditions use the same procedures as those used by Bell *et al.* (1994) for hyperbolic conservation laws. The construction of the grid hierarchy is based on error estimation criteria specified by the user to indicate where additional resolution is required. The error criteria are currently based on tracking component density gradients for one of the components; however, more sophisticated criteria based on estimating the error can be used (e.g. Berger & Colella 1989). Given grids at level ℓ, we use the error estimation procedure to tag cells where the criteria for further refinement are met. We then tag a buffer region *n*_{buf} cells wide around the originally tagged cells so that the features of interest are safely contained within the newly created fine level. The tagged cells are grouped into rectangular patches using the clustering algorithm given in Berger & Rigoutsos (1991). These rectangular patches are refined to form the grids at the next level. The process is repeated until either the error tolerance criteria are satisfied or a specified maximum level is reached. The proper nesting requirement is imposed at this stage.

At *t*=0, the initial data are used to create grids at level 0 through . (Grids have a user-specified maximum size; therefore, more than one grid may be needed to cover the physical domain.) As the solution advances in time, the regridding algorithm is called every *k*_{ℓ} (also user specified) level ℓ steps to redefine grids at levels ℓ+1 to . The parameter *k*_{ℓ} should reflect the constraint that we do not want the tagged feature to move off level ℓ+1 during *k*_{ℓ} level ℓ+1 time steps, so typically *k*_{ℓ}≤*n*_{buf}. Level 0 grids remain unchanged throughout the calculation. Grids at level ℓ+1 are only modified at the end of level ℓ time steps, but because we subcycle in time, i.e. Δ*t*^{ℓ+1}=Δ*t*^{ℓ}/*r*, level ℓ+2 grids can be created and/or modified in the middle of a level ℓ time step if *k*_{ℓ+1}<*r*.

When new grids are created at level ℓ+1, the data on these new grids are copied from the previous grids at level ℓ+1 if possible or interpolated in space from the underlying level ℓ grids otherwise. After regridding, we always recalculate *v*_{T} from the new data, thus the first fine time step after regridding is always a full solve rather than a partial solve. We note here that while there is a user-specified limit to the number of levels allowed, at any given time in the calculation there may not be that many levels in the hierarchy, i.e. can change dynamically as the calculation proceeds, as long as it does not exceed the user-specified limit.

### (b) Overview of time-stepping procedure

The adaptive time-stepping approach uses temporal subcycling so that each level is advanced independently at its own time step (Δ*t*^{ℓ+1}=Δ*t*^{ℓ}/*r*), requiring no interlevel communication other than the supplying of Dirichlet data from a coarse level to be used as boundary conditions at the next finer level. When coarse and fine data reach the same point in time, the data at the different levels are synchronized. The synchronization approaches are based on the same set of algorithmic ideas as those developed in Almgren *et al.* (1998).

The adaptive time-step algorithm can most easily be thought of as a recursive procedure, in which to advance level ℓ, , the following steps are taken.

(i) Advance level ℓ in time as if it is the only level. Use boundary conditions for component densities and pressure from level ℓ−1 if level ℓ>0, and from the physical domain boundaries.

(ii) If

(a) Advance level (ℓ+1)

*r*times with time step Δ*t*^{ℓ+1}=Δ*t*^{ℓ}/*r*. Use boundary conditions for component densities and pressure from level ℓ and from the physical domain boundaries.(b) Synchronize the data between levels ℓ and ℓ+1 and interpolate corrections to higher levels if .

Before describing the steps of the synchronization in detail, we first discuss, in general terms, how to synchronize the data at different levels so that the solution as computed on each level sequentially can most closely approximate the solution which would be found using composite solves. During the advance of each level, for each operator, we supply Dirichlet boundary data for the fine grids from the next coarser grid. This implies that the values at both levels are consistent, but the computed fluxes at the coarse/fine interfaces are not. The synchronization solves correction equations that account for discrepancies in fluxes between levels. The correction equations reflect the type of operator being corrected. For hyperbolic equations, the correction of flux discrepancies is a simple explicit flux correction as discussed in Berger & Colella (1989). For a self-adjoint elliptic operator, the discrepancy in the fluxes represents a discontinuity in normal derivative at the coarse/fine boundary and the correction equation takes the form of a discrete layer potential problem (Almgren *et al*. 1998).

After the level ℓ+1 data have been advanced to the same point in time as the level ℓ data, there are three sources of discrepancy in the composite solution that need to be corrected in the synchronization step.

(D1) The data at level ℓ that underlie the level ℓ+1 data are not synchronized with the level ℓ+1 data.

(D2) The composite total velocity computed from the pressure equation, defined as the time-averaged (over a level ℓ time step) level ℓ+1 advection velocity on all level ℓ+1 faces including the ℓ|(ℓ+1) interface, and the level ℓ advection velocity on all other level ℓ faces, does not satisfy the composite divergence constraint at the ℓ|(ℓ+1) interface.

(D3) The advective and diffusive fluxes from the level ℓ faces and the level ℓ+1 faces do not agree at the ℓ|(ℓ+1) interface, resulting in a loss of conservation.

The aim of the synchronization steps is to correct the effects of each mismatch. Discrepancy (D1) is easily corrected by averaging **n** at level ℓ+1 onto level ℓ. We denote this correction by (S1). The third discrepancy (D3) is also easily corrected in the case without diffusive terms; we simply correct the solution in the coarse cells immediately adjacent to the fine level by the divergence of an edge-based correction field that is non-zero only on the coarse/fine interface where it contains the difference between the coarse flux and time- and space-averaged fine flux. In combination with (S1), this ensures overall conservation.

The second discrepancy (D2) is discretely manifest as a non-zero difference between the coarse-grid total velocity and the effective time- and space-averaged fine-grid total velocity at the coarse/fine interface. This difference results from not having satisfied the elliptic matching conditions at the coarse/fine interface during the pressure solve. An elliptic solve is necessary to correct for the discrepancy. We perform a level ℓ ‘pressure sync solve’ (S2) with the right-hand side defined as the divergence of the mismatch between the level ℓ and the time-averaged level ℓ+1 total velocity. From this solve, we define a total velocity correction field that is used to modify the explicit dependence of the hyperbolic flux terms on the total velocity. In the case of zero diffusive terms, these ‘re-advection corrections’ are explicitly added to the new-time solution in all cells at level ℓ and are interpolated to the new-time solution at all higher levels.

In the case of non-zero diffusive terms, the modification of the solution by the refluxing (S2) and re-advection (S3) corrections requires solving additional elliptic equations. These will be described in more detail in §3*c*.

### (c) Details of time-stepping procedure

Assume now that we are advancing level ℓ, , one level ℓ time step. We now add the level index, ℓ, to the superscript of each quantity, defining, for example, **n**^{n,ℓ} as the component density on level ℓ at time *t*^{n}. We define Δ*t*^{ℓ} as the time step of level ℓ.

#### (i) Advancing a single level

To advance the data on level ℓ one level ℓ time step, we follow the time-stepping procedure as described for the single-grid algorithm in the previous section. For the level ℓ advancement, we advance all the grids at that level simultaneously. Explicit advection work is decoupled except for the exchange of boundary conditions. When a coarse/fine boundary does not coincide with a physical domain boundary for the level ℓ advection step, level ℓ−1 data are interpolated linearly in time to specify Dirichlet boundary conditions at the coarse/fine boundary. The procedures used for these interpolations are the same as those discussed in Almgren *et al*. (1998).

#### (ii) Computing the coarse–fine discrepancy

Over the course of a level ℓ time step, we must accumulate several quantities at the ℓ|(ℓ+1) interface in order to correctly capture the flux discrepancies at the end of the level ℓ time step. We refer to the face-based data structures that contain these quantities as registers. The total velocity and flux registers accumulate the discrepancies between the level ℓ and level ℓ+1 face-based total velocity *v*_{T} and fluxes *F*, respectively.

These registers are defined only on the ℓ|(ℓ+1) interface and are indexed by level ℓ indices. Note that in *d* dimensions, one level ℓ face contains *r*^{d−1} level ℓ+1 faces; the sums over faces below should be interpreted as summing over all level ℓ+1 faces that are contained in the level ℓ face. The sums over *k* should be understood as summing over the *r* level ℓ+1 time steps contained within a single level ℓ time step.

At the end of the level ℓ time step, the velocity register () holds the difference between the total velocity at level ℓ and the time average over one level ℓ time step of the space average over the area of the level ℓ face of the total velocity at level ℓ+1: We note that the velocity included in the total velocity register is the time-centred velocity, , that is defined in Step 3 as the time average of and because this is the velocity field used in advection.

The flux registers, and contain the differences between the hyperbolic and diffusive fluxes calculated at level ℓ and the time average over the level ℓ time step of the space average over the area of the level ℓ face of the advective fluxes at level ℓ+1: where as computed in Step 4 of the algorithm.

We note here that the signs of the quantities added to the flux registers actually depend on the orientation of the normal facing away from the fine grid. We follow the convention below that the signs are given for the faces at which the fine grid is in the direction of the lower coordinate indices.

#### (iii) Synchronization of data

The first synchronization step (S1) has already been described in §3*b*. Here we give the details of the other two steps of the synchronization.

The discrepancy in the total velocity (D2) is captured in the divergence of defines the right-hand side for the level ℓ elliptic sync solve (S2). We solve
on all grids at level ℓ for the correction *δ**e*^{ℓ}. Recall is defined only at the coarse/fine interface; here is defined to be the discrete divergence operator evaluated only on the level ℓ cells adjacent to the interface but not underlying any level ℓ+1 grids and including contributions only from the coarse/fine interface. Boundary conditions on physical no-flow boundaries are homogeneous Neumann (∂(*δ**e*)^{ℓ}/∂*n*=0); on outflow *δ**e*^{ℓ}=0. If ℓ>0, the boundary conditions for *δ**e*^{ℓ} are given as homogeneous Dirichlet conditions on the level ℓ−1 cells outside the level ℓ grids. We then define the correction velocity field from *δ**e*^{ℓ}

We now use to define new fluxes on all level ℓ faces and define as the difference between the original fluxes and these newly computed fluxes. In the case of non-zero diffusive terms, we must diffuse the re-advection and refluxing corrections before adding them to the new-time solution, so we do not add them directly to the solution. Rather, the divergence of the re-advection flux corrections is added to the advective and diffusive flux mismatches to define the cell-centred right-hand sides for the refluxing solves (S3): 3.1 Then, we solve for the correction to the solution, 3.2 where and .

If ℓ>0, we must now modify the level ℓ−1 velocity registers and flux registers to account for the corrections to the solution owing to the re-advection corrections as well as the diffused corrections. This is analogous to the accumulation of advective and diffusive fluxes during the advance of a single level. To do this, we set These will enter into the correction performed between levels ℓ−1 and ℓ at the end of the level ℓ−1 time step if ℓ>1.

The corrected solution at level ℓ is then given by
If , we interpolate the correction onto the fine grids at all finer levels, *q*, using conservative interpolation
This completes the synchronization steps for scalar quantities.

## 4. Numerical results

In this section, we examine the numerical properties of the method described in the previous section. In §4*a*, we examine the convergence behaviour of the method. While the basic discretization presented here is similar to a traditional IMPES fractional step discretization, we have modified the method such that the total velocity is centred in time. This makes the method formally second-order accurate, in contrast to the standard IMPES approach. In our first example, we use a stable one-dimensional problem based on a two-component single-phase system to demonstrate the method’s second-order rate of convergence and to illustrate the loss of accuracy associated with the standard IMPES-type fractional step scheme. We then demonstrate the second-order convergence rate of the algorithm with a two-dimensional problem and show that the partial solve algorithm maintains the overall accuracy of the method relative to the full solve version. We also examine the effects of centring the total velocity and compare the current approach with the standard IMPES discretization when the flow is unstable. In §4*b*, we compare solutions obtained on a uniform grid with those computed with AMR. For this purpose, we use a more complicated example involving a three-component two-phase system. In §4*c*,*d*, we simulate a large three-dimensional problem using the parallelized AMR algorithm, and present the scaling behaviour. We note that for the numerical results presented here, we neglect capillary pressure.

### (a) Convergence

For our first study, we consider a single-phase two-component system where components 1 and 2 are completely miscible and form a mixture with viscosity given by
4.1
where *c*=*n*_{1}/(*n*_{1}+*n*_{2}) and *M*=10. In our algorithm as presented, both the pressure discretization and the component conservation equations are discretized to second-order accuracy. In the standard IMPES formulation, the total velocity is not centred in time, which introduces a small first-order temporal error. In most realistic problems, this error is typically dominated by spatial truncation errors, so we examine the problem in the context of one dimension where we can use fine spatial resolution. We take as initial data *c*(*x*)=0.5(1+tanh((0.3−*x*)/0.05)) on the interval [0,2] and impose a constant pressure drop. The time step, which is fixed for each resolution, is set so that the CFL is initially approximately 0.35. In addition, we turn off the flux limiting in the advection scheme.

The *L*^{1} norm of the errors for the standard IMPES discretization and the full solve version of the new algorithm are presented in table 1. The new algorithm shows clear second-order convergence, whereas the standard IMPES algorithm is only first-order accurate asymptotically. At the coarsest resolution, the two approaches have nearly the same error and the IMPES discretization shows nearly second-order convergence between the two coarsest discretizations. At these coarser resolutions, the temporal discretization error is dominated by the spatial truncation error, which is second-order accurate in both cases.

Next, we examine the convergence properties of the single-grid algorithm in two dimensions. For this problem, we consider a domain [0,*l*_{x}]×[0,*l*_{y}], where *l*_{x}=16 m and *l*_{y}=4 m. At time *t*=0, the domain contains only component 2. For time *t*>0, component 1 flows into the domain from the left edge of the domain with a composition given by *n*_{1}(0,*y*;*t*)=*τ**ρ*_{1} and *n*_{2}(0,*y*;*t*)=(1−*τ*)*ρ*_{2}, where . The densities of components 1 and 2 are 1000 and 800 kg m^{−3}, respectively. We impose no-flow conditions on the top and bottom edges and a pressure difference of 2 atm between the inlet and the outlet. We define *M*=2. The porosity is uniformly 1 while the permeability function *κ*(*x*,*y*) is given by
We note that *τ*(*t*) and *κ*(*x*,*y*) are chosen so that the solution will be smooth. Gravity points in the negative *y* direction. The calculation is run for 1.2×10^{7} s with Δ*x*=Δ*y*=1/2^{n} m for *n*=3,4,5,6. Uniform grids and a fixed time step are used for the purpose of evaluating the convergence of the solution.

Table 2 shows the discrete , *L*^{1} and *L*^{2} norms of the difference between the solution *n*_{1} obtained on each grid and that obtained on the next finer grid and the resulting convergence rates. The rate between the two columns of error norms is defined as , where *ϵ*_{l} and *ϵ*_{r} are the errors shown in the left and right columns. Table 2 clearly demonstrates the second-order convergence property of the algorithm when measured in the *L*^{1} and *L*^{2} norms. The convergence rate of the error is less than 2; this can be attributed to limiters used in the Godunov scheme. Table 3 shows the same data but for calculations done using the full solve approach. By comparing the two tables, we see that the errors obtained from the different approaches are very similar. We conclude that the omission of Step 1 does not adversely affect the convergence rate of our method.

The next set of calculations compares the solution obtained using our method with that obtained using a first-order method when the solution is no longer smooth. The goal is to show that even though the flow is no longer in the asymptotically second-order regime of the algorithm, the accuracy properties of the formally second-order method still have important consequences for the behaviour of the flow. We construct the first-order method by eliminating Steps 3 and 4 from the method presented in §2 and setting **n**^{n+1}=**n**^{n+1,*}.

However, the error in the solution cannot be directly measured by taking the norm of the difference between two solutions. We instead compare the mixing length given by the length of the fingered zone. Following the work of Manickam & Homsy (1995), we define the mixing length, *L*_{δ}, as
(4.2)
where
and *δ*=1×10^{−5}. We first look at a stable fingering profile that allows us to establish a hypothesis on how mixing length is related to accuracy. We then examine an example where the flow is highly unstable.

For the first example, we again look at the two-component single-phase system described in the previous section. We consider here a uniform *κ*(*x*,*y*)=200 mD and a pressure difference of 0.5 atm between the inlet and the outlet. We see in figure 1 that a stable finger extends from the inlet to the outlet at the bottom of the simulation domain. We then measure the mixing length for three different grid sizes and three different approaches (table 4). The grids we consider are 512×128, 1024×256 and 2048×512. The three approaches are first-order method, second-order method with the partial solve algorithm and second-order method with the full solve algorithm. The solution obtained from the 2048×512 grid using the second-order method with the full solve algorithm is used as the reference solution. The simulation is run to a final time of 8×10^{6} s.

To amplify the difference between the first-order and the second-order methods, we change the permeability distribution function *κ*(*x*,*y*) to
4.3

In addition, *M* is increased to 10, dramatically increasing the inherent instability of the problem. Figure 2 shows that the viscous fingers resulting from the first-order and the second-order methods have some notable differences. Slight variations in the solution can develop into distinct fingering profiles owing to instabilities in the solution. To further illustrate this point, we examine a two-component single-phase system in a domain where *κ*=100(1+10^{−8}*ε*) mD and *ε* is a random number between 0 and 1. Figure 3 shows that two different realizations of *ε* lead to two different fingering profiles. This suggests that any minor perturbation can lead to variation in the details of the fingering. In spite of this sensitivity, we see that the omission of Step 1 does not have significant effects on the solution. In table 5, we list the mixing lengths for grid sizes 512×128 and 1024×256 computed with the three approaches. The minor differences in mixing length for the two second-order methods further demonstrate that the additional work required for the full solve algorithm is unnecessary. Results of the subsequent sections are based on the second-order method with the partial solve algorithm.

### (b) Comparison between adaptive mesh refinement and uniform grid

The next set of calculations compares the composite grid solution with the uniform fine grid solution for a two-phase, three-component system representing a simple polymer flood example. We have an aqueous phase consisting of two components, water (*n*_{w}) and polymer (*n*_{p}), that is injected into a domain filled with an oleic phase composed of a single component (*n*_{o}). We impose the same type of boundary conditions as in §4*a*. For this problem, we have *ρ*_{w}=*ρ*_{p}=1000 kg m^{−3}, *μ*_{aqueous}=0.175+0.35(*n*_{p}/(*n*_{w}+*n*_{p})) cP, *ρ*_{o}=800 kg m^{−3} and *μ*_{oleic}=0.35 cP. The porosity is again uniformly 1. The permeability function *κ*(*x*,*y*) is given by equation (4.3) and the injection composition is given by
*n*_{p}=1000−*n*_{w} kg m^{−3}, *t*>0, and *n*_{o}=0 kg m^{−3}, *t*>0. Again, the flow is not in the asymptotically second-order regime of the algorithm. The goal here is to show that with the adaptive algorithm one can achieve accuracy comparable to that on a uniform fine grid. The regridding criterion, which flags coarse grid cells for refinement, is based on the magnitude of |∇*n*_{w}|.

We note that although the composition for *t*<5×10^{5} s will lead to a stable front, the composition for *t*≥5×10^{5} s will lead to a highly unstable front. A direct comparison between the solution after *t*=5×10^{5} s is then not meaningful for reasons elucidated in the previous section. It nevertheless allows us to examine the behaviour of both propagation modes of the aqueous phase.

For the uniform grid, we use a grid size of 1024×256. For the AMR simulation, we have a base grid of 256×64, with two additional levels of refinement with a factor of two refinement between each level. A particular cell is tagged for refinement if |∇*n*_{w}|>500. For both simulations, the time step is computed with a CFL number of 0.4. Figure 4 shows that stable solution fronts for the uniform and the adaptive simulations are identical. However, while the unstable fronts exhibit similar features, there are notable differences in the fingering structure. Unstable fronts can thus develop differently even though captured at the same effective resolution, as shown by the uniform and adaptive solutions shown in figure 5.

In table 6, we compare the computational time needed to solve the current problem on a uniform grid and an adaptive grid. We note that at 7×10^{5} s, AMR reduces the computational time by 88 per cent. However, as the percentage of the domain that is refined increases, the performance benefit decreases. At 10^{6} s, when the percentage of the domain refined at level 2 increases to 22.5 per cent, the computational saving from using AMR is reduced to 80 per cent. As expected, the efficiency gains from using AMR will be problem dependent.

### (c) A three-dimensional example

We consider the same three-component two-phase system for a three-dimensional problem. The domain is now given by [0,16] m×[0,4] m×[0,2] m. We impose no-flow boundary conditions along the faces with normals in the *y* and *z* directions and a pressure difference of 2 atm in the *x* direction. The properties of the components are the same as in §4*b*, and porosity is again uniformly 1. The permeability function *κ*(*x*,*y*,*z*) here is a log-random distribution with a mean of 150 mD and a variance of 0.1 mD. In addition, *κ* is highly correlated in the *x* and *y* directions but not in the *z* direction, as shown in figure 6.

We performed the simulation on a Cray XT4 supercomputer. At the coarsest level, the grid size is 128×64×16 with up to three additional levels of refinement with refinement ratio 2 between adjacent levels. The resulting grid and component density of component water at time 1.6×10^{6} s are shown in figure 7.

### (d) Parallel performance

By adopting a block-structured form of AMR, the solution at each level in the hierarchy is naturally represented in terms of data defined on a collection of logically rectangular grid patches each containing a large number of points. Thus, the data are represented by a modest collection of relatively large, regular data objects when compared with a point-by-point refinement strategy. This type of approach allows us to amortize the irregular aspects of an adaptive algorithm over large regular operations on the grid patches. This organization of data into large aggregate grid patches also provides a model for parallelization of the AMR methodology.

Our adaptive methodology is embodied in a hybrid C++/Fortran software system. In this framework, memory management and flow control are expressed in the C++ portions of the program and the numerically intensive portions of the computation are handled in Fortran. The software is written using a layered approach, with a foundation library, BoxLib, that is responsible for the basic algorithm domain abstractions at the bottom, and a framework library, AMRLIB, that marshals the components of the AMR algorithm, at the top. Support libraries built on BOXLIB are used as necessary to implement application components such as interpolation of data between levels, the coarse/fine interface synchronization routines and linear solvers used in the pressure and diffusion solves.

The fundamental parallel abstraction is the `MultiFab`, which encapsulates the Fortran-compatible data defined on unions of `Box`s; a `MultiFab` can be used as if it were an array of Fortran-compatible grids. The grids that make up the `MultiFab` are distributed among the processors, with the implementation assigning grids to processors using the distribution given by the load balance scheme described in Crutchfield (1991) and in Rendleman *et al*. (2000). This load balance scheme is based on a dynamic programming approach for solving the knapsack problem: the computational work in the irregularly sized grids of the AMR data structures is equalized among the available processors. After the initial allocation of grids, additional changes to the grid distribution can be performed to reduce communications between processors. For non-reacting flows, the number of cells per grid is often a good work estimate. `MultiFab` operations are performed with an *owner computes* rule with each processor operating independently on its local data. For operations that require data owned by other processors, the `MultiFab` operations are preceded by a data exchange between processors.

Each processor contains *meta-data* that are needed to fully specify the geometry and processor assignments of the `MultiFab`s. At a minimum, this requires the storage of an array of boxes specifying the index space region for each AMR level of refinement. In the parallel implementation, meta-data also include the processor distribution of the Fortran-compatible data. The meta-data can thus be used to dynamically evaluate the necessary communication patterns for sharing data among processors, enabling us to optimize communications patterns within the algorithm.

To measure the parallel performance of the algorithm, we use a weak scaling study based on a replicated problem strategy as discussed in Colella *et al*. (2007). The case considered here is a two-component single-phase system with a layered permeability function. By replicating the problem in the *y* and *z* directions, we are able to scale the problem size without modifying the problem characteristics, particularly with regard to how adaptive criteria and grid generation impact the overall problem. In figure 8, we present scaling data compared with the ideal behaviour for a range of processors from 4 to 256. We see a slight deviation from ideal scaling, which is primarily attributable to increases in time spent in the elliptic solver.

## 5. Conclusions and future work

In this paper, we have presented a structured AMR algorithm for incompressible flows in porous media. The method is based on a total velocity splitting of the equations into an elliptic pressure equation and a system of hyperbolic conservation laws for the component densities. The resulting system is discretized using cell-centred differencing for the pressure equation and a second-order Godunov scheme for the component conservation equations. The adaptive algorithm uses subcycling in time, i.e. the ratio of the time step to the mesh spacing is constant across levels of refinement. Rather than solve a global composite equation for the pressure at every time step at the finest level, we perform single-level solves at each fine step and then use an elliptic correction solve to synchronize the solution across levels. This leads to a considerable improvement in computational efficiency. The method has been implemented in parallel and shows excellent scalability up to 256 processors.

Our goal in future work will be to extend this approach to more realistic problems. As a first step in that direction, we will extend the current approach to include realistic geochemical reactions. Beyond that, our target will be to extend the approach to more realistic fluid models that include compressibility, interphase mass transfer and thermal effects.

## Footnotes

One contribution of 10 to a Theme Issue ‘Mesh generation and mesh adaptation for large-scale Earth-system modelling’.

- © 2009 The Royal Society