Fully Conservative Characteristic Methods for Flow and Transport: - - PowerPoint PPT Presentation

Fully Conservative Characteristic Methods for Flow and Transport:

Part I, Linear Transport Part II, Theoretical Considerations Part III, Nonlinear Two-Phase Flow Todd Arbogast Department of Mathematics

and

Center for Subsurface Modeling, Institute for Computational Engineering and Sciences (ICES) The University of Texas at Austin

Chieh-Sen (Jason) Huang, National Sun Yat-sen University (Taiwan) Thomas F. Russell, U.S. National Science Foundation Wenhao Wang, The University of Texas at Austin This work was supported by

• U.S. National Science Foundation
• KAUST through the Academic Excellence Alliance

Center for Subsurface Modeling Institute for Computational Engineering and Sciences The University of Texas at Austin, USA

Fully Conservative Characteristic Methods for Flow and Transport: Part I, Linear Transport

Todd Arbogast Department of Mathematics

and

Center for Subsurface Modeling, Institute for Computational Engineering and Sciences (ICES) The University of Texas at Austin Chieh-Sen (Jason) Huang, National Sun Yat-sen University (Taiwan)

This work was supported by

• U.S. National Science Foundation
• KAUST through the Academic Excellence Alliance

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 1 The University of Texas at Austin, USA

Outline

• 1. Transport Problems and Local Conservation Principles
• 2. Characteristic Methods for Linear Transport
• 3. Local Volume Conservation in Characteristic Methods
• 4. Some Numerical Results
• 5. Summary and Conclusions

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 2 The University of Texas at Austin, USA

Transport Problems and Local Conservation Principles

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 3 The University of Texas at Austin, USA

Conservative Fluid Flow Suppose ξ is a conserved quantity ξ (mass/volume)

v is the fluid velocity (length/time)

ξv is the flux of ξ (mass/area/time) Q is an external source or sink of fluid (mass/volume/time)

❩❩❩❩ ❩❩❩❩ ✓ ✓ ✓ ✼

A

v

Within a region of space R, the total amount of ξ changes in time by d dt

• R ξ dx
• = −
• ∂R ξv · ν da(x)
• +
• R Q dx
• Change in R

Flow across ∂R Sources/sinks = ⇒ conservation locally on R

• R ξt dx

= −

• R ∇ · (ξv) dx
• +
• R Q dx

Divergence Theorem

❩❩❩❩ ❩❩❩❩ ❇ ❇ ❇ ❇ ❇ ▼ ✓ ✓ ✓ ✼ ❩ ❩ ❩ ❩

∂R

v·ν v

This is true for each region R, so in fact ξt + ∇ · (ξv) = Q

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 4 The University of Texas at Austin, USA

A Transport Problem–1 One incompressible fluid (tracer) flowing miscibly in another incompressible fluid, within an incompressible medium. Velocity of the bulk fluid. Conservation of bulk fluid mass (ξ = φρ) gives ξt + ∇ · (ξv) = Q = ⇒ ∇ · u = q

u is the (unknown) bulk fluid velocity (v = u/φ)

φ is the porosity (constant in time) ρ is the (constant) density q is the source/sink (wells, Q = ρq) Simple Tracer Transport. Conservation of tracer mass (ξ = φc) gives φct + ∇ · (cu) = cIq+ + cq− ≡ qc(c) c is the (unknown) tracer concentration cI is the given concentration of injected fluid q+/q− is q when positive/negative

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 5 The University of Texas at Austin, USA

A Transport Problem–2 However, transport is not the only process occurring! Mass flux.

v = cu − D∇c

(Transported plus Diffusive Flux) D is the diffusion/dispersion coefficient Chemical reactions. q = qc(c) + R(c) (Wells plus Reactions) R is the reaction term Tracer Transport. Conservation of tracer mass gives φct + ∇ · (cu − D∇c) = qc(c) + R(c)

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 6 The University of Texas at Austin, USA

Operator Splitting of Transport Equation—1 φct + ∇ · (cu − D∇c) = qc(c) + R(c) Discretization in time: ∆t > 0 and tn = n∆t. We want to solve the transport and reactive part of the equation explicitly and the diffusive part implicitly. Thus, we want φcn+1 − cn ∆t + ∇ · (cnu) − ∇ · (D∇cn+1) = qc(cn) + R(cn) This is equivalent to the three steps (Reaction) φcr − cn ∆t = R(cn)

• φct = R(c)

(Transport) φcrt − cr ∆t + ∇ · (cnu) = qc(cn)

• φct + ∇ · (cu) = qc(c)

(Diffusion) φcn+1 − crt ∆t − ∇ · (D∇cn+1) = 0

• φct − ∇ · (D∇c) = 0

with some intermediate ˜ c and ˆ c.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 7 The University of Texas at Austin, USA

Operator Splitting of Transport Equation—2 Nonlinear Ordinary Differential Equation part (Reaction) φct = R(c) Linear Hyperbolic part (Transport) φct + ∇ · (cu) = qc(c) Linear Parabolic part (Diffusion/dispersion) φct − ∇ · (D∇c) = 0 We discuss approximations of the transport step only.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 8 The University of Texas at Austin, USA

Locally Conservative Methods A locally conservative numerical method is one for which the approximate solution satisfies the conservation principle, but only over certain discrete regions. Bulk volume conservation:

E ∇ · u dx =

• E q dx

We solve for the flow u using a locally conservative mixed finite element method, conserving it over the grid elements E. Tracer mass conservation: A compatible method would require

• E
• φct + ∇ · (cu)
• dx =
• E qc dx

Characteristic methods do not use the grid elements. They use more general space-time regions E. Question: Is there a numerical incompatibility in the conservation of the two fluids?

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 9 The University of Texas at Austin, USA

Compatibility Tracer mass conservation:

• E
• φct + ∇ · (cu)
• dx dt =
• E qc dx dt

Ambient fluid mass conservation:

• E
• φ(1 − c)t + ∇ · ((1 − c)u)
• dx dt ?

=

• E(q − qc) dx dt

⇐ ⇒

• E ∇ · u dx dt ?

=

• E q dx dt

For mixed finite element methods, ∇ · u = Pq pointwise (P is projection onto the elements E). Thus there is no incompatibility (up to treatment of the wells). Remark: This may not be true for nonconservative Galerkin methods. Remark: The reactive and diffusive steps of the operator splitting must also be solved by locally conservative methods (e.g., mixed finite elements or cell-centered finite differences)!

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 10 The University of Texas at Austin, USA

Characteristic Methods for Linear Transport

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 11 The University of Texas at Austin, USA

Characteristic Tracing of Points The characteristic trace-forward of the point x is denoted ˆ

x = ˆ x(x; t).

It satisfies the ordinary differential equation dˆ

x

dt = u(ˆ

x, t)

φ(ˆ

x)

, tn < t ≤ tn+1 ˆ

x(tn) = x

In the absence of sources/sinks and diffusion, fluid particles simply travel along the characteristics of the equation.

✲

tn tn+1

✻

ˆ

x x

Space Time The concentration is constant along this space-time path, since dc(ˆ

x, t)

dt = ∂c ∂t + ∇c · dˆ

x

dt = ct + ∇c · u φ = 1 φ

• φct + u · ∇c
• = 0

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 12 The University of Texas at Austin, USA

Characteristic Trace-back of Points The characteristic trace-back of the point x is denoted ˇ x = ˇ x(x; t). It satisfies the (time backward) ordinary differential equation dˇ x dt = u(ˇ x, t) φ(ˇ x) , tn ≤ t < tn+1 ˇ x(tn+1) = x In the absence of sources/sinks and diffusion, fluid particles simply travel along the characteristics of the equation.

✲

tn tn+1

✻

x

ˇ x Space Time Again, the concentration is constant along this space-time path.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 13 The University of Texas at Austin, USA

Modified Method of Characteristics (MMOC) (Douglas and Russell, 1982) Key idea: Use a finite difference approximation of the characteristic derivative dc dt ≡ ct(x, t) + u(x, t) φ · ∇c(x, t) ≈ c(x, t + ∆t) − c(ˇ x, t) ∆t This results in the approximation φc(x, t + ∆t) − c(ˇ x, t) ∆t = (cI − c)q+ at each grid point

✲

tn tn+1

✻

c(x) c(ˇ x) Problems: Because the method is based on points, it violates local mass conservation constraints for both the bulk fluid and the tracer.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 14 The University of Texas at Austin, USA

Characteristic Trace-back of Regions To obtain mass conservation, ... Key idea: Trace regions rather than points! The particles in a grid element E trace back to a region ˇ E ˇ E = {ˇ x ∈ Ω : ˇ x = ˇ x(x; tn) for some x ∈ E}. In space and time, we actually trace a region E = E(E) given by E = {(ˇ x, t) ∈ Ω × [tn, tn+1] : ˇ x = ˇ x(x; t) for some x ∈ E}.

✲

tn tn+1

✻

E ˇ E E

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 15 The University of Texas at Austin, USA

Local Mass Conservation of the Tracer φct + ∇ · (cu) = qc(c) Integrate in space-time over E and use the divergence theorem

• E ∇x,t ·
• cu

φc

• dx dt =
• ∂E c
• u

φ

• · νx,t dσ

=

• E φcn+1 dx −
• ˇ

E φcn dx+

• S c
• u

φ

• · νx,t dσ

The last term is integration on the space-time sides S of E, but

• u

φ

• is orthogonal to νx,t there!

The local mass constraint:

• E φcn+1 dx =
• ˇ

E φcn dx

+

• E qc dx dt −
• SB

cu · ν dσ

✲

tn tn+1

✻

E ˇ E E SB

❥

νx,t

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 16 The University of Texas at Austin, USA

Local Mass Conservation of the Bulk Fluid A similar local mass constraint holds for the bulk fluid (c ≡ 1) φt + ∇ · u = ∇ · u = q Since we are dealing with incompressible fluids, we call this the local volume constraint. The local volume constraint:

• E φ dx =
• ˇ

E φ dx +

• E q dx dt −
• SB

u · ν dσ

Let the pore volume of region R be |R|φ =

• R φ dx

Then we have |E|φ = | ˇ E|φ +

• E q dx dt −
• SB

u · ν dσ

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 17 The University of Texas at Austin, USA

Characteristics Mixed Method (CMM) (A., Chilakapati, and Wheeler, 1992; A. and Wheeler, 1995) Use lowest order Raviart-Thomas mixed finite elements. The scalar function is approximated as a constant on each element in space: c(x, tn) ≈ cn

h(x) = cn E ∈ R

when x ∈ E Then cn+1

E

|E|φ =

• ˇ

E φcn h dx +

• E qcn

h dx dt −

• SB

cn

hu · ν dσ

Remark: Practical implementation requires that ˇ E be approximated by ˜ E ≈ ˇ E, a polygon. This is equivalent to modifying the velocity field, so tracer mass is still conserved locally by the above equation.

② ② ② ② ② ② ② ②

ˇ E ≈ ˜ E | ˇ E|φ = | ˜ E|φ

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 18 The University of Texas at Austin, USA

Volume Error Problem: The local volume constraint may be violated for this perturbed velocity! This is because the volume constraint does not enter into the definition of the method. This leads to incorrect densities of the tracer, which leads, e.g., to bad approximation of reaction dynamics.

2.00 0.10 0.01

• 0.01
• 0.10

Relative volume errors around an injection well.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 19 The University of Texas at Austin, USA

Local Volume Conservation in Characteristic Methods

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 20 The University of Texas at Austin, USA

Volume Conservation Key idea: To obtain volume conservation, define ˜ E ≈ ˇ E to be a simple shape that satisfies the volume constraint |E|φ − | ˜ E|φ =

• ˜

E q dx dt −

• ˜

SB

u · ν dσ

Strategy: Suppose that E is a rectangle. Perturb the vertices and midpoints of ˇ E

• nly a little so that we get a polygon ˜

E with 8 sides such that the above constraint is satisfied.

② ② ② ② ② ② ② ②

ˇ E ≈ ˜ E | ˇ E|φ = | ˜ E|φ

② ② ② ② ② ② ② ②

Volume Corrected Characteristics Mixed Method (VCCMM) cn+1

E

|E|φ =

• ˜

E φcn h dx +

• ˜

E qcn

h dx dt −

• ˜

SB

cn

hu · ν dσ

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 21 The University of Texas at Austin, USA

An Example of Unphysical Flow Problem: It is easy to introduce systematic bias into the flow field and thereby produce unphysical flows. We must adjust points very carefully!

10 20 30 40 50 60 10 20 30 40 50 60 Volume not conserved 10 years 10 20 30 40 50 60 10 20 30 40 50 60 30 years 10 20 30 40 50 60 10 20 30 40 50 60 Volume conserved 10 20 30 40 50 60 10 20 30 40 50 60

Biased trace-back adjustment has introduced unphysical flow corresponding to a large, incorrect velocity channel.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 22 The University of Texas at Austin, USA

Forward Trace of Injection Wells Most of the error is near injection wells. Characteristic tracing back in time traces into the well, which is difficult to approximate. Key idea 1: Trace the well forward (out of the well). Adjust region to satisfy the volume constraint (cf. Healy and Russell, 2000). The characteristic trace-forward of the point x is denoted ˆ

x = ˆ x(x; t),

and it satisfies the ordinary differential equation dˆ

x

dt = u(ˆ

x, t)

φ(ˆ

x)

, tn < t ≤ tn+1 ˆ

x(tn) = x

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 23 The University of Texas at Austin, USA

Conservation Near Injection Wells W E ˜ W ˜ E Well volume constraint: Adjust ˜ W so | ˜ W|φ = |W|φ +

• ˜

Ef

W

q dx dt Element volume constraint: Adjust ˜ E so |E|φ = | ˜ E|φ + |E ∩ ˜ W|φ | ˜ W|φ

• ˜

Ef

W

q dx dt Transport: cn+1

E

|E|φ =

• ˜

E φcn h dx + |E ∩ ˜

W|φ | ˜ W|φ

• ˜

Ef

W

qcn

h dx dt Center for Subsurface Modeling Institute for Computational Engineering and Sciences 24 The University of Texas at Austin, USA

Inflow boundaries Like injection wells, inflow boundaries trace back “out of the domain.” Idea: Either trace inflow boundaries forward, or “fold” the time axis down to the xy-plane to create a ghost region. t y x ∂Ω Ω −u · ν φ tn tn+1 t y x ˜ E Volume constraint: Replace φ by u · ν in the ghost region: |E|φ − | ˜ E ∩ Ω|φ +

• ˜

SB

u · ν dσ = |E|φ − | ˜

E|φ =

• ˜

E q dx dt,

Mass constraint: Replace φcn

h by cn I u · ν in the ghost region:

cn+1

E

|E|φ −

• ˜

E φcn h dx =

• ˜

E qcn

h dx dt. Center for Subsurface Modeling Institute for Computational Engineering and Sciences 25 The University of Texas at Austin, USA

Trace-back Point Adjustment Key idea 2: Adjust points in the direction of the flow; that is, along the characteristics in time (cf. Douglas, Huang, and Pereira, 1999). To define ˜

x, for τn ≈ tn, we solve (backwards)

d˜

x

dt = u(˜

x, t)

φ(˜

x)

, τn < t ≤ tn+1 ˜

x(tn+1) = x

We convert space error into time error:

✲

tn tn+1

✻

x

ˇ x

✲

τn ˜

x

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 26 The University of Texas at Austin, USA

Trace-back (or Forward) Point Adjustment—1 Proceed away from injection wells and inflow boundaries by “layers.” For each layer, obtain volume conservation in two steps.

• 1. Volume conservation of the layer. Adjust the exterior contour of the

entire layer along the characteristics until the volume of the layer is correct (within a small tolerance). That is, in the absence of other sources, inflow boundaries, and sinks,

• E in the layer

|E|φ =

• E in the layer

| ˜ E|φ

✇ Adjusted point (fixed) ❣ Points adjusted simultaneously in the

direction of the characteristic (we use a type of “bisection” algorithm) × Points adjusted individually transverse to the flow in Step 2

❅

• Flow

❛❛ ❛✥✥ ✥ ✥✥ ✥❛❛ ❛✥✥ ✥ ❵❵ ❵✥✥ ✥✥✥ ✥ ❛❛ ❛✥✥ ✥ ❵❵ ❵✥✥ ✥ ✦✦ ✦❵❵ ❵✥✥ ✥ ❵❵ ❵✥✥ ✥✥✥ ✥✦✦ ✦ ☎ ☎☎ ☎ ☎☎ ❉ ❉ ❉ ❉ ❉ ❉ ☎ ☎☎ ✂ ✂ ✂ ✂ ❉ ❉ ❉ ✄ ✄ ❉ ❉ ❉ ☎ ☎☎ ✆ ✆ ✆ ✆ ❉ ❉ ❉ ❈ ❈ ☎ ☎☎ ☞ ☞ ☞ ❊ ❊ ❊ ❊ ☎ ☎☎ ❊ ❊ ❊ ❊ ✇ ✇ ✇ ✇ ✇ ✇ ✇ ✇

×

❣

×

❣ ❣ ✇ ✇ ✇ × ❣ ✇ ✇ ✇ ❣ ✲ ✛ ✇ ❣ ✇ ✇ ✇ × ❣

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 27 The University of Texas at Austin, USA

Trace-back (or Forward) Point Adjustment—2

• 2. Element volume conservation. Within the layer, sequentially adjust

the interior midpoint of each element transverse to the flow until the volume of the element is correct (within a small tolerance).

③ Adjusted point (fixed)

× Points individually adjusted transverse to the flow > Flow ×xi,j+1/2

xi,j+1 xi,j

< >

③ ③ ③ ③ ③

×

③ ③ ③ ③ ③ ③

×

❈ ❈ ❈ ❈ ❈ ❈ ❈ ❈ ❈ ✄ ✄ ✄ ✄ ✄✄ ✂ ✂ ✂ ✂ ✂ ✂ ✂ ❈ ❈ ❈ ❈ ❈ ❈ ✄ ✄ ✄ ✄ ✄ ✄ ✄ ✄ ✄ ☎ ☎ ☎ ☎ ☎ ☎ ☎ ✏✏✏✏✏✏ ✏PPPPPP P ✦✦✦✦✦✦✦❛❛❛❛❛❛❛ ✏✏✏✏✏✏ ✏❛❛❛❛❛❛❛

Remark: This is an extremely fast direct algorithm.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 28 The University of Texas at Austin, USA

Some Numerical Results

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 29 The University of Texas at Austin, USA

Numerical Results Darcy’s Law: Complete the equation with

u = −k

µ∇p p is the fluid pressure k is the permeability µ is the fluid viscosity Solve using locally conservative mixed finite elements (or cell-centered finite differences). Coefficient of variation: Measure the variability of k by Cv = 1 Mk

• 1

|Ω|

• Ω
• k(x) − Mk

2 dx 1/2

where the mean is Mk = 1 |Ω|

• Ω k(x) dx

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 30 The University of Texas at Austin, USA

A Nuclear Contamination Problem—1 The permeability is log-normal and fractal. Mk = 2 × 10−10 cm2 (about 20 md) Cv = 0.522 (varies over five orders of magnitude).

1E-08 1E-09 1E-10 1E-11 1E-12 1E-13

Inflow Outflow

✲ ✲ ✲ ✲ ✲ ✲ t Injection well

Domain is 256 m×256 m

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 31 The University of Texas at Austin, USA

A Nuclear Contamination Problem—2 If we use a small time step of ∆t = 1.5 years, we can trace back into the injection well.

2.00 0.10 0.01

• 0.01
• 0.10

CMM volume errors VCCMM (volume errors 10−9)

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 32 The University of Texas at Austin, USA

A Nuclear Contamination Problem—3 Concentration at 30 years on a 64 × 64 grid, with ∆t = 1.5 yr ≈ 2.6 CFL

32 64 96 128 160 64 96 128 160 192 1.1E-05 1.0E-05 9.0E-06 8.0E-06 7.0E-06 6.0E-06 5.0E-06 4.0E-06

CMM

32 64 96 128 160 64 96 128 160 192 1.1E-05 1.0E-05 9.0E-06 8.0E-06 7.0E-06 6.0E-06 5.0E-06 4.0E-06

VCCMM CMM overshoots the maximum concentration of 1E-5 by 34% up to 1.34E-5.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 33 The University of Texas at Austin, USA

The Courant-Friedrichs-Lewy (CFL) Condition Explicit methods in 1-dimension have a time step restriction, known as the Courant-Friedrichs-Lewy (CFL) time-step, given by ∆t ≤ ∆tCFL,1-D = max

x∈Ω

hφ(x) |u(x)| where h is the grid spacing. In 2-dimensions, we should limit ∆t to half this value, ∆t ≤ ∆tCFL,2-D = max

x∈Ω

hφ(x) 2|u(x)| Godunov’s method is a popular method, that is unstable if the CFL condition is violated. In principle, characteristic methods are not subject to this constraint, and large time steps can be used.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 34 The University of Texas at Austin, USA

A Nuclear Contamination Problem—4 Concentration at 30 years on a 64 × 64 grid

32 64 96 128 160 64 96 128 160 192 1.0E-05 9.0E-06 8.0E-06 7.0E-06 6.0E-06 5.0E-06 4.0E-06

Godunov ∆t = 0.586 yr = 1 CFL

32 64 96 128 160 192 160 128 96 64 1.0E-05 9.0E-06 8.0E-06 7.0E-06 6.0E-06 5.0E-06 4.0E-06

VCCMM-TF ∆t = 3 yr = 5.1 CFL

• We use trace-forwarding near the well.
• No overshoot for either method.
• Less numerical diffusion for VCCMM-TF.
• 51 Godunov steps vs. 10 for VCCMM-TF.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 35 The University of Texas at Austin, USA

A Nuclear Contamination Problem—5 Concentration at 30 years on a 128 × 128 grid

32 64 96 128 160 64 96 128 160 192 1.0E-05 9.0E-06 8.0E-06 7.0E-06 6.0E-06 5.0E-06 4.0E-06

Godunov ∆t = 0.146 yr = 1 CFL

32 64 96 128 160 192 160 128 96 64 1.0E-05 9.0E-06 8.0E-06 7.0E-06 6.0E-06 5.0E-06 4.0E-06

VCCMM-TF ∆t = 1 yr = 6.8 CFL

• We use trace-forwarding near the well.
• No overshoot for either method.
• Less numerical diffusion for VCCMM-TF.
• 205 Godunov steps vs. 30 for VCCMM-TF.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 36 The University of Texas at Austin, USA

A Nuclear Contamination Problem—6 Concentration at 30 years on a 256 × 256 grid

32 64 96 128 160 64 96 128 160 192 1.0E-05 9.0E-06 8.0E-06 7.0E-06 6.0E-06 5.0E-06 4.0E-06

Godunov ∆t = .0366 yr = 1 CFL

32 64 96 128 160 192 160 128 96 64 1.0E-05 9.0E-06 8.0E-06 7.0E-06 6.0E-06 5.0E-06 4.0E-06

VCCMM-TF ∆t = 0.5 yr = 13.7 CFL

• We use trace-forwarding near the well.
• No overshoot for either method.
• Less numerical diffusion for VCCMM-TF.
• 820 Godunov steps vs. 60 for VCCMM-TF.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 37 The University of Texas at Austin, USA

A Quarter Five-Spot Problem—1 Geostatistically generated permeability. Mk = 100 md Cv = 2.58 (varies over four orders of magnitude).

5E-12 1E-12 5E-13 1E-13 5E-14 1E-14 5E-15 1E-15

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 38 The University of Texas at Austin, USA

A Quarter Five-Spot Problem—2 Concentration at 3.36 years using ∆t = 0.012 yr = 10.68 CFL

1.00 0.85 0.70 0.55 0.40 0.25 0.10

CMM

1.00 0.85 0.70 0.55 0.40 0.25 0.10

VCCMM-TF

• CMM shows both overshoot and undershoot.
• Very large initial volume imbalances throughout the domain.
• If ∆t = 0.0136 yr = 12.10 CFL initially creates degenerate trace-back

regions, which cannot be used.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 39 The University of Texas at Austin, USA

A Linear Flood Problem—1 A test with an inflow boundary.

• Permeability field of the quarter five-spot problem
• Linear pressure drop across the domain in the x-direction.

Concentration at 25 years using ∆t = 0.18 year = 3 CFL.

1.000 0.800 0.600 0.400 0.200 0.000

CMM

1.000 0.800 0.600 0.400 0.200 0.000

VCCMM

• CMM has severe overshoots up to c = 1.65.
• Initial volume errors exceed 10% in the interior of the domain.
• If ∆t = 9 CFL, initial relative volume errors are around 25%.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 40 The University of Texas at Austin, USA

A Linear Flood Problem—2 Concentration at 25 years using trace-forwarding of the inflow boundary.

1.000 0.800 0.600 0.400 0.200 0.000

VCCMM-TF ∆t = 0.18 yr = 3 CFL

1.000 0.800 0.600 0.400 0.200 0.000

VCCMM-TF ∆t = 0.36 yr = 6 CFL

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 41 The University of Texas at Austin, USA

A Fluvial Domain Problem—1

• Domain 600 × 600 feet2 solved on a 60 × 60 grid.
• Permeability of 3 values, Mk = 4.056 Darcy and Cv = 1.15. φ = 0.2.
• Wells in opposite corners, injecting 1 pore volume every 3 years.
• ∆t = 0.015 years = 12.84 CFL.

0.1 1.0 10.0

The permeability, in Darcies

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 42 The University of Texas at Austin, USA

A Fluvial Domain Problem—2 VCCMM-TF concentration

1.000 0.800 0.600 0.400 0.200 0.000

Time 1.05 years (step 70)

1.000 0.800 0.600 0.400 0.200 0.000

Time 1.65 years (step 110) Remark: CMM alone produces negative concentrations on a 30 × 30 grid with ∆t = 0.01 year, indicating that the trace-back regions self intersect.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 43 The University of Texas at Austin, USA

A Fluvial Domain Problem—3 CMM with trace-forwarding of wells only

1.000 0.800 0.600 0.400 0.200 0.000

Time 1.05 years (step 70)

1.000 0.800 0.600 0.400 0.200 0.000

Time 1.65 years (step 110)

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 44 The University of Texas at Austin, USA

Comparison with an Analytic Solution Comparison of VCCMM-TF concentration with an analytic solution for radial flow from a well in a horizontally infinite, uniform porous medium.

50 100 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 50 100 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Analytic VCCMM-TF

50 100 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 50 100 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

Analytic VCCMM-TF Longitudinal dispersion 20 cm Longitudinal dispersion 100 cm The analytic solution is derived by a code from P.A. Hsieh (1986). VCCMM-TF uses ∆t = 8 CFL.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 45 The University of Texas at Austin, USA

Cellular Flow—1 We consider the divergence-free (solenoidal) velocity

u =

• −0.1 + 50 sin(2πx) cos(2πy)

−50 cos(2πx) sin(2πy)

• ,

The velocity field, with a region

• ccupied initially by a higher

tracer concentration value. Early tracer concentration

• n a 256 × 256 grid, ∆t = 3e–4.

Steps 10, 30, 60, and 80.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 46 The University of Texas at Austin, USA

Cellular Flow—2 Tracer on 256 × 256 grid for CMM (left) and VCCMM (right) at t = 0.3 (3000 steps, top) and t = 0.4 (4000 steps, bottom). Tracer on 512 × 512 grid for CMM (left) and VCCMM (right) at t = 0.3 (6000 steps, top) and t = 0.4 (8000 steps, bottom). Remark: ∆t is 2.56 times the CFL time step Remark: For CMM, tracer fills the center of the flow cells due to a larger numerical diffusion from the volume errors.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 47 The University of Texas at Austin, USA

Cellular Flow—3 VCCMM tracer 256 × 256 grid, ∆t = 3e–4 (left) 512 × 512 grid, ∆t = 1.5e–4 (right) t = 0.3 (top) and 0.4 (bottom). Remark: ∆t is 7.69 times the CFL time step

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 48 The University of Texas at Austin, USA

Cellular Flow—4 Uncorrected volume errors using 256 × 256 grid and ∆t = 1e–4. The variation is from ±1.72e–3%.

✲ ✲ ✲ ✲ ✲ ✲ ✲ ✲ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✻ ✲ ✲ ✲ ✲ ✲

1 2 7 8 12 16 13 3 4 9 10 17 14 15 5 6 11 18 19 20 21 22 23 24 Remark: A simple adjustment is all that is necessary, since there are no

• wells. We use a staircase-like diagonal element-by-element adjustment

strategy from the bottom-left to the upper-right corner. Note that errors cancel along the paths.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 49 The University of Texas at Austin, USA

Summary and Conclusions

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 50 The University of Texas at Austin, USA

Summary and Conclusions

• 1. The goal. A critical aspect of approximating hyperbolic transport

problems is to conserve the mass of the tracer locally. However, it is just as critical to the conserve locally the mass of the bulk fluid.

• 2. The problem. Current characteristic methods can conserve tracer

mass locally. However, they do not conserve bulk fluid mass locally.

• 3. The solution. Adjust trace-back regions so volume is conserved:
• Avoid systematic bias by adjusting along the characteristics;
• Use trace-forwarding around injection wells and inflow boundaries.

The only time step restriction is that the trace-back regions do not degenerate/self-intersect (or trace more points of the boundary ∂E).

• 4. The results. Numerical examples show:
• Large time steps can be taken in practice (such as 14 times the

2-D CFL limited time step);

• Little numerical diffusion;
• No overshoots;
• Most of the problems are near wells.

Center for Subsurface Modeling Institute for Computational Engineering and Sciences 51 The University of Texas at Austin, USA