A Scalar Conservation Law for Plume Migration in Carbon Sequestration

TL;DR

This paper develops a new mathematical model using a quasi-linear hyperbolic PDE with dual flux curves to better understand and analyze the complex migration and trapping of CO2 in geological formations during sequestration.

Contribution

It introduces a novel dual flux model with a new cross-hatch characteristic construction, extending the classical theory of scalar conservation laws with convex flux.

Findings

Characterizes shock and rarefaction waves in the dual flux model

Reveals differences from classical scalar conservation laws with convex flux

Provides a generalized entropy condition for shocks in the dual flux context

Abstract

A quasi-linear hyperbolic partial differential equation with a discontinuous flux models geologic carbon dioxide migration and storage. Dual flux curves characterize the model, giving rise to flux discontinuities. One convex flux describes the invasion of the plume into pore space, and the other captures the flow as the plume leaves carbon dioxide bubbles behind, which are then trapped in the pore space. We investigate the method of characteristics, the structure of shock and rarefaction waves, and the result of binary wave interactions. The dual flux property introduces unexpected differences between the structure of these solutions and those of a scalar conservation law with a convex flux. During our analysis, we introduce a new construction of cross-hatch characteristics in regions of the space-time plane where the solution is constant, and there are two characteristic speeds. This construction is used to generalize the notion of the Lax entropy condition for admissible shocks, and is crucial to continuing the propagation of a shock wave if its speed becomes characteristic.

Equations

$$u_{t}+\sigma\,f(u)_{x}=Pe^{-1}\sigma\big{(}f(u)u_{x}\big{)}_{x},$$

$$\sigma f(u)=\sigma \frac{u(1-u)}{u(\mathcal{M}-1)+1},$$

$$\sigma=\left\{\begin{array}[]{ll}1-\varepsilon,&\mbox{if}\ \ u_{t}>0\,,\\[2.0pt] 1,&\mbox{if}\ \ u_{t}<0\,.\end{array}\right.$$

$$u_t+\sigma f(u{)}_x=0,$$

$$u_{t}=-\sigma f^{\prime}(u)u_{x}.\$$

$${\gamma }^{\prime }(t)=f^{\prime }(u(\gamma (t),t))\left(1+ϵ\frac{u_x^{+}(t)}{u_x^{-}(t)-u_x^{+}(t)}\right),t>0,$$

$$u(\gamma (t{)}^{-},t)=u(\gamma (t{)}^{+},t)$$

$$u_x^{\pm }(t)=\frac{u_0^{\prime }(\gamma -\sigma f^{\prime }(u(\gamma ,t))}{1+\sigma u_0^{\prime }(\gamma -\sigma f^{\prime }(u(\gamma ,t))f^{\prime \prime }(u_0(\gamma -\sigma f^{\prime }(u(\gamma ,t)))t},$$

$${\gamma }^{\prime }(t)=f^{\prime }(u(\gamma (t),t))\left(1+ϵ\frac{u_x^{-}(t)}{u_x^{-}(t)-u_x^{+}(t)}\right),t>0.$$

$$u(x,t)=\left\{\begin{array}[]{ll}u_{0}(x-f^{\prime}(u)t),&x\leq\gamma(t)\\[6.0pt] u_{0}(x-rf^{\prime}(u)t),&x\geq\gamma(t),\end{array}\right.$$

$$u_t+\sigma (u_t)f(u{)}_x=0,\sigma (u_t)=1-ϵH(u_t),$$

$$u\left(x,t\right)=\left\{\begin{array}[]{ll}u_{L},&\quad x<\Lambda t,\\[2.0pt] u_{R},&\quad x>\Lambda t\end{array}\right.$$

$$\Lambda\,=\,\dfrac{\;\sigma\,\big{[}\,f(u_{R})-f(u_{L})\,\big{]}\;}{u_{R}-u_{L}}$$

$$u\left(x,t\right)=\left\{\begin{array}[]{lc}u_{L}\,,&\quad\dfrac{x}{t}<\sigma\,f^{\prime}(u_{L})\,,\\[6.0pt] \check{u}\left(\dfrac{x}{t}\right),&\quad\sigma\,f^{\prime}(u_{L})\leq\dfrac{x}{t}\leq\sigma\,f^{\prime}(u_{R})\,,\\[6.0pt] u_{R}\,,&\quad\sigma\,f^{\prime}(u_{R})<\dfrac{x}{t}\,,\end{array}\right.$$

$$[u_x]=\frac{1}{tϵf^{\prime \prime }(u^{\ast })}=-\frac{\sqrt{\mathcal{M}}}{2tϵ}$$

$$u(x,0)=\left\{\begin{array}[]{ll}u_{L}\,,&\quad x<0\\[1.0pt] u_{R}\,,&\quad x>0.\end{array}\right.$$

$$u(x,0)=\left\{\begin{array}[]{ll}u_{L}\,,&\quad x<x_{1}\\[3.0pt] u_{M}\,,&\quad x_{1}<x<x_{2}\\[3.0pt] u_{R}\,,&\quad x_{2}<x\end{array}\right.$$

$$\frac{f({\tilde{u}}_M)-f(u_M)}{{\tilde{u}}_M-u_M}=(1-\epsilon )f^{\prime }(u_M),$$

$$u_M<u_R\le {\tilde{u}}_M,u_M<u_L,$$

$$\frac{f(u_R)-f(u_M)}{u_R-u_M}\ge {\lambda }_M.$$

$$(1-\epsilon )\frac{f(u_R)-f(u_M)}{u_R-u_M}>f^{\prime }(u_M).$$

$$\left\{\begin{array}[]{l}\overline{u}_{M}=1\,,\ \ \ \mbox{if}\quad(1-\varepsilon)\,\displaystyle\frac{f(1)-f(u_{M})}{1-u_{M}}>f^{\prime}(u_{M})\\[11.0pt] (1-\varepsilon)\,\displaystyle\frac{f(\overline{u}_{M})-f(u_{M})}{\overline{u}_{M}-u_{M}}=f^{\prime}(u_{M})\,,\ \ \mbox{otherwise.}\end{array}\right.$$

$$u_M<u_R\le {\overline{u}}_M,u_M<u_L.$$

$$f^{\prime }(\tilde{u})=\sigma \frac{f(u_R)-f(\tilde{u})}{u_R-\tilde{u}}.$$

$$\sigma f^{\prime }(\overline{u})=\frac{f(u_R)-f(\overline{u})}{u_R-\overline{u}}$$

$$\left\{\begin{array}[]{l}\tilde{u}_{M}=0\,,\ \ \ \mbox{if}\quad(1-\varepsilon)\,\displaystyle\frac{f(u_{M})-f(0)}{u_{M}}<\Lambda_{LM}\\[11.0pt] (1-\varepsilon)\,\displaystyle\frac{f({u}_{M})-f(\tilde{u}_{M})}{{u}_{M}-\tilde{u}_{M}}=\Lambda_{LM}\,,\ \ \mbox{otherwise.}\end{array}\right.$$

$$0\le u_R\le {\tilde{u}}_M<u_M.$$

$$(1-\epsilon ){\Lambda }_{LM}=\frac{f(u_M)-f({\overline{u}}_M)}{u_M-{\overline{u}}_M}.$$

$$0\le u_R\le {\overline{u}}_M<u_M.$$

$$\displaystyle\left\{\begin{array}[]{ll}u_{t}+\big{(}\sigma\,f(u)\big{)}_{x}=0\,,&\quad x\in\mathbb{R}\,,\;t>0\,,\\[6.0pt] u(x,0)=u_{0}(x)\,,&\quad x\in\mathbb{R}\,,\end{array}\right.$$

Full text

1 Introduction

Some 35.7 billion tonnes of carbon dioxide (CO2) were emitted into the atmosphere in 2014 [18], an increase from the previous year’s global CO2 emissions of 32 gigatonnes [10]. In 2000, the Intergovernmental Panel on Climate Change projected a range of estimated emissions from fossil fuel combustion and industrial processes for the year 2020; current emissions are within that annual planning range of 29 to 44 billions tonnes of CO2 [11]. The capture of CO2 before its exodus into the atmosphere seems to be a promising technological solution to reduce the escalating global impact of CO2 emissions. In such a process, gaseous CO2 is collected at industrial sites and power plants, compressed, and injected into geological formations deep underground. Geotechnical evidence suggests that there is a potential subsurface storage capability of 2,000 billion tonnes of CO2 in porous reservoirs worldwide [11]. A goal of future and ongoing carbon dioxide capture and storage projects, such as the Sleipner project beneath the North Sea [23, 24, 25], is to permanently trap CO2 underground [6, 17]. While a wealth of seismic surveys of the Sleipner project have indicated no signs of leakage [3], the possibility of escape of the injected CO2 from brine-filled aquifers remains a concern.

During injection, the captured gaseous CO2 is compressed and becomes supercritical; hence, upon release into the porous rock, the sequestered CO2 behaves like a liquid. Since it is less dense than the ambient brine, the injected plume rises within the aquifer [7, 8, 9]. Appropriate sites for carbon capture and storage projects have an impermeable cap rock in the geological formation that acts as a barrier to hinder the upward migration of the buoyant plume and keep the CO2 beneath the Earth’s surface. Once the plume rises to the impermeable upper boundary, the CO2 travels along inclines in the cap’s lower surface and spreads through the porous rock as a gravity current. As the plume migrates, it deposits bubbles of CO2 that remain in place. The sequestration is successful if all of the CO2 in the plume is deposited before the plume reaches fractures within the cap rock that would allow leakage of the plume from the aquifer [6, 10, 22, 24].

This mechanism to permanently immobilize CO2 within a porous medium is known as residual trapping. Capillary forces between the two fluids (brine and supercritical CO2) stably trap bubbles of CO2 within pore spaces. Hesse et al. [8] formulated a quasi-linear hyperbolic partial differential equation with a discontinuous flux to model geologic carbon dioxide migration and storage through residual trapping. A striking feature of their model is that, due to the discontinuous flux, the entire CO2 plume is deposited as bubbles in a finite time.

In this paper, we explore the model in more detail, approximating solutions of the Cauchy problem using wave-front tracking. In §2 we describe the model of [8], whose key feature is a switch between two flux functions, occuring when the plume changes from propagating into a region of brine to depositing CO2 droplets. In §3 we describe novel features of the method of characteristics, and the construction of fundamental solutions of the equation, namely shock waves and rarefaction waves. To establish the admissibility of shock waves, we introduce the notion of cross-hatch characteristics to address the ambiguity of characteristic speeds due to the twin flux functions. §4 includes a detailed description of wave interactions, including some properties that do not occur in conventional scalar conservation laws. In §5 we construct piecewise constant approximate solutions of the Cauchy problem using expansion shocks in place of rarefaction waves. We conclude the paper in §6 with some remarks.

2 The Two-Flux Model

In this section, we outline several simplifying assumptions about the aquifer and the nature of the flow, then state the model, a first order partial differential equation with a switch in flux depending on whether, at a given location, the CO2 plume is advancing, or depositing bubbles in its wake.

2.1 Model Assumptions

Subsurface geology often has complicated spatial variability, and three-dimensional models of carbon sequestration require unresolved and difficult issues. To simplify matters, we consider a porous aquifer that is locally uniform in the transverse horizontal direction and analyze the two-dimensional propagation of a cross-section of the flow. Consider a porous aquifer of constant thickness $H$ beneath an impermeable cap rock sloped at constant angle $\theta$. A buoyant plume of supercritical carbon dioxide, CO2, with height $h(x,t)$ at position $x$ and time $t$ is introduced to the brine-filled aquifer for storage, as shown in Fig. 2.1(a). As in the figure, the CO2 plume is represented by a sharp interface, beglecting the dissolution of CO2 into the brine [8, 10]. The viscosity contrast between the two fluids propels the CO2 plume to invade available pore space as it migrates as a gravity current [6, 13, 19].

Isolated ganglia of carbon dioxide will be trapped in a region of the permeable aquifer, with residual surface once the plume recedes, Fig. 2.1(b). Thus, the volume of CO2 within the plume decreases, as the plume migrates and becomes disconnected from the immobilized residual bubbles. It is assumed that pressure within the current is hydrostatic since the advection-dominated migration is mainly horizontal. Within the aquifer, volume is conserved, and the multiphase extension of Darcy’s law is applicable in place of conservation of momentum [7, 10, 12, 24]. Combining these assumptions with a hyperbolic limit yields a non-dimensional first order partial differential equation given in [8].

Let u=\frac{h}{H}\in\big{[}\,0\,,1\,\big{]} be the dimensionless height of the CO2 plume, $t$ the non-dimensional advection-dominated time scale, and $x$ the dimensionless spatial variable, based on the initial width $L$ of a typical plume. The mobility ratio, $\mathcal{M}$, between the supercritical carbon dioxide and the brine depends on permeability and viscosity of each phase; for carbon sequestration, the invading CO2 is more mobile than the ambient brine, so that $\mathcal{M}\geq 1$ [8, 17].

The residual surface of immobile CO2 remaining in the wake of the migrating plume is controlled by a residual trapping parameter, \varepsilon\in\big{[}\,0\,,1\,\big{]}. Both $\varepsilon$ and $\mathcal{M}$ are constant material properties [6, 8, 12].

2.2 Governing Equation

Hesse, Orr, and Tchelepi [8] modeled the evolution of a gravity current with residual trapping as a scalar equation

[TABLE]

in which the flux $\sigma f$ is a fractional flow rate obtained by eliminating pressure from a version of Darcy’s law,

[TABLE]

and $Pe$ is the Peclet number, representing the ratio of advective and diffusive time scales. The parameter \sigma\in\big{[}\,0\,,1\,\big{]} depends on the evolution and is a step function given by

[TABLE]

When $u_{t}>0$, the migrating CO2 is invading new pore spaces, whereas when $u_{t}<0$ the plume is draining, no new trapping locations are sought and the brine invades, isolating bubbles of CO2.

In a sloping aquifer, advection dominates diffusion, and the equation reduces (in the limit $Pe\to\infty$) to the nonlinear conservation law

[TABLE]

The switch between migration and deposition represented by the parameter $\sigma$ gives rise to discontinuities in the flux. As shown in Fig. 2.2, the lower flux curve describes the invasion of the plume into pore space, and the upper flux captures the flow as the plume leaves CO2 bubbles behind, which are then trapped by brine in the pore space. The characteristic speed is therefore increased during deposition.

Flux functions with discontinuities in space have been previously studied, [4, 16, 21]; however, the flux in this model depends on the sign of $u_{t},$ a different kind of discontinuity that introduces new phenomena. For $\varepsilon=0$, there is a single flux function; the aquifer has no available pore space to trap CO2, and the plume migrates according to the classical case in which the plume volume remains fixed and would migrate indefinitely with no deposition. Typically \varepsilon\in\big{(}\,0\,,1\,\big{]} in geologic storage [8, 12, 17, 20], and the entire compactly supported plume may be trapped within available pore space after a finite time and within a finite aquifer volume.

3 Characteristics and Shocks

In this section, we consider equation (2.4) with the switch parameter $\sigma$ given by (2.3), and assumptions on the flux $f$ consistent with the Hesse et al model [8]:

(H) $f:[0,1]\to\mathbb{R}$ is $C^{2},\ f(0)=f(1)=0,\ f^{\prime\prime}(u)<0,\ 0\leq u\leq 1.$

The value $u=u^{*}$ with $f^{\prime}(u^{*})=0$ plays a significant role in the construction of admissible shock solutions of (2.4). For the flux function (2.2), we have u^{*}={1}\,\big{/}\big{(}{1+\sqrt{\mathcal{M}}}\,\big{)}\,.

We explain the role of the discontinuous switch function $\sigma$ (see (2.3)), in the construction of continuous solutions, shocks and rarefactions. We also resolve an ambiguity, related to the constant regions of $u$ in the characteristic plane, by introducing cross-hatch characteristics.

3.1 Method of Characteristics

Suppose $u(x,t)$ is a continuous solution of (2.4) with initial data $u(x,0)=u_{0}(x)$ in $C^{2}.$ For short time, the solution should be obtained by the method of characteristics. However, since there are two possible characteristic speeds, $f^{\prime}(u)$ and $(1-\epsilon)f^{\prime}(u),$ we have to choose between them, at least in open regions of the $x,t$ plane where $u(x,t)$ is non-constant. Where $u_{x}(x,t)\neq 0,$ we see that the choice of characteristic speed depends on the sign of $u_{x}(x,t)$ and the sign of $f^{\prime}(u),$ since

[TABLE]

Suppose $u^{\prime}_{0}(\bar{x})=0.$ If $u^{\prime}_{0}(x)$ is constant in a neighborhood of $\bar{x},$ then the solution is continued to $t>0$ as that constant, and we introduce cross-hatch characteristics, meaning that both charactistic speeds apply where $u(x,t)$ is constant in an open $x,t$ region. If $u_{0}$ has an inflection point at $x,$ then the function is either increasing or decreasing at $x,$ and the characteristic speed is uniquely defined. However, if $u_{0}(x)$ has a maximum or minimum at $x=\bar{x},$ then something interesting happens. Suppose for now that $f^{\prime}(u_{0}(x))>0.$

(1) If $u_{0}(x)$ has a minimum at $x=\bar{x},$ with then near $(x,t)=(\bar{x},0),$ the characteristics originating from $x<\bar{x},t=0$ are slower than those originating from $x>\bar{x},t=0.$ Consequently, the solution satisfies $u(x,t)=u_{0}(\bar{x}),$ for $(x,t)$ between the characteristics $x=(1-\epsilon)f^{\prime}(u_{0}(\bar{x}))t+\bar{x}$ and $x=f^{\prime}(u_{0}(\bar{x}))t+\bar{x}.$

(2) If $u_{0}(x)$ has a maximum at $x=\bar{x},$ then near $(x,t)=(\bar{x},0),$ the characteristics originating from $x<\bar{x},t=0$ are faster than those originating from $x>\bar{x},t=0.$ In this case, the solution has a corner along a curve $x=\gamma(t),$ with $\gamma(0)=\bar{x}.$

Lemma 3.1**.**

Suppose $u(x,t)$ is a piecewise $C^{2}$ solution of equation (2.4) satisfying $u(x,0)=u_{0}(x),$ where $u_{0}\in C^{2}(\mathbb{R}).$ If $u_{0}(x)$ has a maximum at $x=\bar{x},$ and $f^{\prime}(u_{0}(\bar{x}))>0,$ the maximum propagates as a corner $x=\gamma(t)$ in the graph of $u(x,t),\ t>0,$ satisfying

[TABLE]

where $u_{x}^{\pm}(t)=u_{x}(\gamma(t)^{\pm},t);$ $\gamma(0)=\bar{x},\ \gamma^{\prime}(0)=(1-\frac{\epsilon}{2})c,\ c=f^{\prime}(u_{0}(\bar{x})).$

Proof:  To derive an ODE for $\gamma(t),$ we differentiate the continuity condition

[TABLE]

with respect to $t,$ and use the identity (3.1). After some manipulation, we establish (3.2) for $t>0,$ where necessarily $u_{x}^{-}(t)>0>u_{x}^{+}(t).$ Note that away from $x=\gamma(t),$ the solution is determined from the method of characteristics. Thus, $u_{x}^{\pm}(t)$ depend implicitly on $\gamma(t):$

[TABLE]

where $\gamma=\gamma(t)$ and $\sigma=1-\epsilon,1$ for $u_{x}^{\pm}$ respectively.

However, equation (3.2) has a singular limit as $t\to 0,$ since $u_{x}^{\pm}(0)=0.$ Let $c=f^{\prime}(u_{0}(\bar{x})),$ and note that $f^{\prime}(u(\gamma(t),t))=c$ to leading order as $t\to 0.$ Without loss of generality, we assume that $c>0.$ Similarly, since $u_{0}(x)$ is $C^{2}$ and has a maximum at $x=\bar{x},$ if $u_{0}(x_{L})=u_{0}(x_{R}),$ with $x_{L}<\bar{x}<x_{R},$ then to leading order, $x_{L}-\bar{x}=\bar{x}-x_{R}.$ Now consider the solution $u(x,t).$ It is determined at the maximum from two different characteristics, that meet at $x=\gamma(t).$ If the two characteristics emanate from $x_{L}<x_{R},$ then $\gamma=f^{\prime}(u)t+x_{L}=(1-\epsilon)f^{\prime}(u)t+x_{R}.$ Thus, to leading order near $t=0,\ \gamma=ct+2\bar{x}-x_{R}=(1-\epsilon)ct+x_{R}.$ Solving the second equation, we have $x_{R}=\frac{1}{2}ct\epsilon+\bar{x}.$ Hence, $\gamma=\bar{x}+(1-\frac{\epsilon}{2})ct$ to leading order. Thus, as $t\to 0,\ \gamma^{\prime}(t)\to(1-\frac{\epsilon}{2})c.$ That is, the initial speed of the corner, at the local maximum of $u(x,t)$ (with respect to $x$) is the average of the two characteristic speeds $c$ and $(1-\epsilon)c.$

Remarks  1.  If $c=f^{\prime}(u_{0}(\bar{x}))<0,$ a corresponding argument applies, but the propagation is to the left. In this case, we have

[TABLE]

2. The functions $u^{\pm}_{x}(t)$ depend implicitly on $\gamma$ as follows.For $f^{\prime}(u_{0}(\bar{x}))>0,$ we have (for $x$ near $\bar{x}$)

[TABLE]

and note that $u(x,t)$ is continuous, at least over some finite time interval $0\leq t\leq T.$ Differentiating with respect to $x,$ we have $u_{x}^{+}(t)=u_{0}^{\prime}/(1+ru_{0}^{\prime}f^{\prime\prime}(u)t),$ where $u_{0}^{\prime}=u_{0}^{\prime}(\gamma(t)-rf^{\prime}(u)t),$ and $u=u(\gamma(t),t),$ and a similar expression for $u_{x}^{-}(t),$ but dropping $r$ from both expressions.

3.2 Cross-hatch Characteristics

Since the switch parameter $\sigma$ is not defined when $u_{t}=0$, the characteristic speed is not well defined in regions of the characteristic plane where the solution is constant. To resolve this, we include characteristics determined by both flux curves at each point where $u_{t}=0$ ; we refer to them as cross-hatch characteristics since they form a cross-hatch pattern in regions where $u$ is constant (see Fig. 3.1(a)).

The two possible characteristic speeds are $\sigma f^{\prime}(u),\ \sigma=1$ or $\sigma=1-\varepsilon.$ We refer to the larger or greater characteristic speed as the faster speed, and the other characteristic speed as the slower speed. *Thus, the faster speed is $f^{\prime}(u)$ (and hence positive) if and only if $u<u^{*}.$ * To clarify further, when $u>u^{*},$ we have $f^{\prime}(u)<0,$ so that the faster speed is $(1-\epsilon)f^{\prime}(u)$ since it is greater than $f^{\prime}(u)$ in this case.

3.3 Shocks

The definition of weak solution for equation (2.4) does not follow the usual pattern of multiplication by a test function and integration by parts. To see that the usual procedure is problematic, we rewrite the equation as

[TABLE]

where $H$ is the Heaviside step function. This form highlights the difficulty of interpreting the equation in the sense of distributions, as both $\sigma(u_{t})$ and $f(u)_{x}$ may be singular. However, if $u(x,t)$ has only jump discontinuities, then although $u_{t}$ is singular, namely a delta function, the definition of $H(u_{t})$ can be extended by $H(a\delta(x))=1$ if $a>0,$ and $H(a\delta(x))=0$ if $a\leq 0.$ In this way, the notion of solution can be extended to piecewise smooth functions.

To define piecewise smooth solutions with jump discontinuities, it is enough to consider a piecewise constant jump discontinuity

[TABLE]

propagating with speed $\Lambda.$ Since $\sigma$ in (3.3) is selected by the sign of $u_{t}$, we set $\sigma=1$ if $u$ jumps down across the shock as time increases; otherwise, if the jump is up, we set $\sigma=1-\varepsilon$. This fixes the value of $\sigma$, and we can write the Rankine-Hugoniot jump condition,

[TABLE]

Hesse et al. [8] justified the choice of $\sigma$ in a slightly different way by including dissipative terms (in (2.1)) that smooth the shock.

For a scalar conservation law with a single flux function, admissible shocks satisfy the Lax entropy condition, requiring characteristics to enter the shock on both sides [14]. Here, with two fluxes, we specify shock admissibility as follows:

Definition 3.2**.**

The shock wave (3.4) is admissible if and only if the faster characteristics enter the shock from both sides.

We argue that (3.4) is an admissible shock if and only if $u_{L}<u_{R},$ just as it would be for a scalar conservation law with a convex flux. As shown in Fig. 3.1(a), $u_{t}<0$ across an admissible forward shock (i.e., with $\Lambda>0$), so that $\sigma=1.$ Consequently, not only is $\Lambda$ determined from the upper flux curve, but also the faster characteristics enter the shock, see Fig. 3.1(b). For an admissible backward shock, with $\Lambda<0$, we have $\sigma=1-\varepsilon,$ and the shock is admissible if and only if the characteristics found on the lower flux curve impinge on the shock on the right, because they are the less negative characteristics, and enter the shock on the left because either they are the fast characteristics (if $u_{L}>u^{*}$), or both families have positive speed (if $u_{L}<u^{*}$), as shown in Fig. 3.2. Once again, this amounts to the condition $u_{L}<u_{R},$ but there is an important point regarding the slower characteristics, which necessarily enter the shock on the right, but may leave on the left.

Lemma 3.3**.**

*The only characteristics that can leave an admissible shock belong to the slower family, and are on the left of the shock. *

Proof: Consider an admissible shock (3.4). If $u_{R}<u^{*}$, the faster characteristic speed is on the upper flux, so $f^{\prime}(u_{R})<\Lambda$ is required for admissibility. Thus, $(1-\varepsilon)\,f^{\prime}(u_{R})<\Lambda$ also. Hence, both characteristics on the right impinge on the shock. If $u_{R}>u^{*}$, the faster characteristic speed is on the lower flux curve, so admissibility requires $(1-\varepsilon)\,f^{\prime}(u_{R})<\Lambda$. Since $u_{R}>u^{*}$, $f^{\prime}(u_{R})<(1-\varepsilon)\,f^{\prime}(u_{R})$, and the slower characteristic on the right also enters the shock. Hence, both characteristics on the right always impinge on an admissible shock. Since the faster characteristics are required to enter the shock on the left, only the slower characteristics on the left can leave the shock.

It is perhaps instructive to understand when the slower characteristics leave an admissible shock. If $u_{L}<u^{*}$ in a backward admissible shock, both characteristics on the left have positive speed but the shock speed is negative, so both characteristics on the left must enter the shock.

When $u_{L}<u^{*}$ in a forward admissible shock, the faster characteristic entering the shock from the left has speed $f^{\prime}(u_{L})>\Lambda\,,$ since $\sigma=1$ in (3.5) for a forward shock. If $\Lambda<(1-\varepsilon)\,f^{\prime}(u_{L})$, the slower characteristics will also impinge on the forward shock; however, it is possible that $(1-\varepsilon)\,f^{\prime}(u_{L})<\Lambda$ , in which case the slower characteristics on the left emanate from the shock. Similarly, if $u^{*}<u_{L}$, an admissible shock requires $0>(1-\varepsilon)\,f^{\prime}(u_{L})>\Lambda$ since $\sigma=1-\varepsilon$ in (3.5). The more negative characteristic speed $f^{\prime}(u_{L})$ may or may not satisfy $f^{\prime}(u_{L})>\Lambda$ , so the slower characteristics on the left can leave the shock.

In summary, since the faster characteristics must impinge on the shock from both sides, the slower characteristics on the right also enter the shock, but the slower characteristics on the left can leave the shock. Fig. 4.1(b) illustrates the latter behavior of the characteristics.

3.4 Expansion Shocks

Expansion shocks are shock wave solutions of (2.4) that are inadmissible. We characterize them here because we will need them in §5 as approximations to rarefactions in wave-front tracking. For scalar conservation laws, expansion shocks have characteristics leaving in forward time on both sides. Here we define a discontinuous function (3.4) to be an expansion shock if it satisfies the Rankine-Hugoniot jump condition (3.5) and the slower characteristics on each side emanate from the shock. The latter condition is equivalent to $u_{R}<u_{L}.$ Then the faster characteristics on the right also leave the shock, but the faster characteristics on the left may or may not enter the shock, as shown in Fig. 3.3.

3.5 Rarefactions

Centered rarefaction fans are continuous weak solutions of (2.4) obtained via the method of characteristics and have the form

[TABLE]

in which, the function $\check{u}$ is given implicitly by $y=\sigma f^{\prime}(\check{u}(y)).$

The rarefaction in Fig. 3.4(a) has both forward and backward characteristics with speeds that depend on the value of $\sigma\,,$ as explained in the figure caption. In Fig. 3.4(b) we show a rarefaction wave approximated by three expansion shocks; from left to right, the expansion shocks have increasing speeds. In Fig. 3.5(a), we show the construction of the rarefaction wave, resolving the initial step down in $u\,,$ using both flux functions, and in Fig. 3.5(b) we show the corresponding plume profile.

The rarefaction solution (3.6) varies continuously from $u_{L}$ to $u_{R}$ in Fig. 3.4(b). In particular, $u(x/t)$ is continuous across $x=0$ even though $\sigma$ in (3.6) has a discontinuity at this position [8]. Correspondingly, there is a discontinuity in the slope of the plume interface due to the jump in $\sigma\,.$ There is a jump $[u_{x}]$ in the derivative $\partial_{x}u$ at $x=0,$ where $f^{\prime}(u)=0,$ and $\sigma$ switches from $\sigma=1$ to $\sigma=1-\epsilon.$ We calculate it assuming $f^{\prime\prime}(u^{*})<0:$

[TABLE]

where the final equality uses the specific flux function (2.2).

4 Wave Interactions

We consider the Riemann problem, consisting of equation (2.4) with jump initial data

[TABLE]

It follows from §3 that the solution is an admissible shock if  $u_{L}<u_{R}$ and a rarefaction fan if  $u_{L}>u_{R}$. While the structure of these individual waves depends on the details of two flux functions and the switch between them, the outcome is, broadly speaking, the same as for a convex scalar conservation law with a single flux.

In this section, we consider pairs of Riemann problems. Each Riemann problem generates a single wave; we are interested in whether the waves interact, and the result of the interaction. The results have significant differences from the corresponding wave interactions for a scalar equation with a single convex flux.

While a detailed classification is complicated, we focus on the main features of solutions of initial value problems with jump initial data of the form

[TABLE]

in which $u_{L}$ and $u_{R}$ are different from $u_{M}.$ Similar to the classical case, if $u(x,0)$ is decreasing, i.e. $u_{L}>u_{M}>u_{R}\,,$ then the solution consists of two rarefaction waves that do not approach. Consequently, since the speed of an approximating expansion shock is between the speeds of the corresponding rarefaction’s trailing and leading characteristics, two expansion shocks will not approach. We treat the three remaining cases in turn, and, if the data has an initial rarefaction, we examine the interactions involving expansion shock approximations.

4.1 Case A: Shock - Rarefaction: $u_{L}<u_{M}$ and $u_{R}<u_{M}$

In this case, we have a shock with speed $\Lambda$ emanating from $x=x_{1}$ at time $t=0,$ and a rarefaction centered at $x=x_{2}>x_{1},t=0.$ To see that the two waves approach, we check that the shock speed is greater than the speed of the trailing characteristic in the rarefaction. There are two cases to consider. In case (i), $\Lambda>0\,,$ the shock admissibility condition requires $f^{\prime}(u_{M})<\Lambda,$ so that the speed $\sigma f^{\prime}(u_{M})$ of the trailing edge of the rarefaction is less than the shock speed, whether $u_{M}<u^{*},$ for which $\sigma=1-\varepsilon,$ or $u_{M}>u^{*},$ for which $\sigma=1.$ In case (ii), $\Lambda<0,$ so that $u_{t}>0$ and $\sigma=1-\epsilon.$ Thus, $u_{M}>u^{*}$ but now shock admissibility requires $\sigma f^{\prime}(u_{M})<\Lambda,$ and the rarefaction, with trailing edge traveling at speed $f^{\prime}(u_{M})<\sigma f^{\prime}(u_{M})<\Lambda<0,$ approaches the shock.

In Fig. 4.1, we illustrate the solution as the interaction proceeds. In this and other figures, we plot exact solutions using the specific flux (2.2) for illustration. On the left we show the track of the rarefaction through the flux curves as the characteristics fan from negative to positive speed. The rarefaction fan provides the values of $u$ on the right of the shock as the evolution proceeds. The shock speed is represented by the slope of the chords in Fig. 4.1(a). As the speed switches from negative to positive, the chord moves from the lower flux graph to the upper, as $u_{t}$ changes sign. The crossover is represented by the horizontal dashed lines. In this example, the construction proceeds until the rarefaction wave has been completely absorbed by the shock. Since the initial data have $u_{R}>u_{L},$ the long-time behavior is a single shock joining $u_{L}$ to $u_{R}.$ On the other hand, if $u_{L}>u_{R}\,,$ then the long-time behavior would be a rarefaction wave, the remnants of the short-time wave joining $u_{M}$ to $u_{R},$ after the interaction with the shock wave has completed.

This interaction of a shock with a rarefaction, illustrated in In Fig. 4.1, appears to be similar to such interactions for a scalar conservation law with convex flux. However, there is a significant difference. While the shock has negative speed, it is calculated from the flux $(1-\varepsilon)f(u).$ The shock is admissible because the characteristics on the left have positive speed, and the faster characteristics on the right have speed $(1-\varepsilon)f^{\prime}(u),$ which is slower than the shock speed, as shown in Fig. 4.1(b). In fact, for the smaller flux (in the lower graph), the shock satisfies the Lax entropy condition. However, as the shock turns and gains positive speed, we switch to the upper flux curve. The characteristics on the right both have negative speed to start with, and hence impinge on the shock. On the left, both characteristics travel faster than the shock. In fact, as the shock turns, it has zero speed, and the characteristics on the left for both fluxes have positive speed, so this property persists for some further time.

However, as the shock continues to accelerate, there is a time, corresponding to shock location , \raisebox{-.9pt} {1}⃝ in Fig. 4.1(b), when the shock moves with the characteristic speed of $(1-\varepsilon)f^{\prime}(u_{L})$ of the smaller flux, see the inclined dashed lines in Fig. 4.1(a) corresponding to $u=u_{\text{graze}}$. Consequently, if we continue to consider only the single slower family of characteristics (that were significant for the shock when it had negative speed), then the shock would fail to satisfy the Lax entropy condition at this time. By including the characteristics of the larger flux (which has already been invoked to calculate the shock speed) we retain admissibility of the shock. This device is consistent with causality, as the constant value of $u$ is carried by both families of characteristics. This example and other similar instances are the reason for including both families of characteristics (hence, cross-hatch characteristics) in open regions of the $(x,t)$ plane where $u$ is constant.

4.2 Case B: Shock - Shock: $u_{L}<u_{M}<u_{R}$

The second case involves a shock from a left state up to a middle state followed by a shock from the middle state up to a right state. Since the flux function is concave, the shock from $u_{L}$ to $u_{M}$ will have a greater shock speed than the shock from $u_{M}$ to $u_{R}$, so the shocks will approach each other and interact at a finite time to yield a single shock from $u_{L}$ up to $u_{R}$ with strength $u_{R}-u_{L}$. If the speeds of the approaching shocks have the same sign, the resulting shock has the same direction; if not, the resulting shock is forward if $f(u_{L})<f(u_{R})$ , stationary if $f(u_{L})=f(u_{R})$ , or backward if $f(u_{L})>f(u_{R})$ . The total variation is unchanged before and after the discontinuities interact, and the middle state is eliminated in finite time.

4.3 Case C: Rarefaction - Shock: $u_{M}<u_{L}$ and $u_{M}<u_{R}$

This case mirrors Case A, in that the short-time solution is a rarefaction wave to the left of a shock wave. However, whereas in Case A the two waves approach, in Case C their approach depends on further restrictions on the data. The reason for this is that the slower characteristics on the left can leave the shock (Lemma 3.3); they are necessarily parallel to the leading edge of the rarefaction. We distinguish two subcases in which the waves do not approach:

(i) If u_{M}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\leq}u^{*}, define $\tilde{u}_{M}$ by

[TABLE]

shown in Fig. 4.2, and let $\lambda_{M}$ denote this speed. Then $\lambda_{M}>0$ is the speed of the leading edge of the rarefaction, and if $u_{R}=\tilde{u}_{M},$ then it is also the speed of the shock, since the shock has a jump up and positive speed. Then for u_{M}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\leq}u^{*} and

[TABLE]

the shock from $u_{M}$ to $u_{R}$ has positive and larger speed:

[TABLE]

Thus, (4.2) is sufficient to guarantee that the shock and rarefaction do not approach.

(ii) Similarly, if $u_{M}>u^{*},$ then the shock speed and speed of the leading edge of the rarefaction wave are both negative. In this case, the rarefaction is backward and uses the larger flux $f(u)$ whereas the shock uses the lower flux (1-\varepsilon){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f(u)}\,. Consequently, the interaction condition becomes

[TABLE]

Define $\overline{u}_{M}>u_{M}$ by

[TABLE]

Then the two waves do not approach if $u^{*}<u_{M}$ and

[TABLE]

In summary, if neither (4.2) nor (4.3) are satisfied by $u_{R}\,,$ then the rarefaction wave and shock wave interact much as in Case A, see Fig. 4.3(a). Otherwise, the shock travels faster than the rarefaction, and there is no interaction, as in Fig. 4.2(b).

Unlike Cases A and B, not all initial conditions in Case C lead to an eliminated initial middle state in finite time. Some solutions in Case C exhibit unusual behavior, due to the flux discontinuity, that does not arise in scalar equations with a single flux: shock speeds determined by one flux curve can equal corresponding characteristic speeds found on the other flux curve. In Fig. 4.3(b), the plume asymptotically approaches a height of \tilde{u}\in\Big{(}\max(u^{*},u_{M})\,,\,\min(u_{L},u_{R})\Big{)} such that

[TABLE]

Hence, if $u_{L}\geq\tilde{u}$, the backward shock does not reach the rarefaction’s trailing characteristic; the shock speed approaches the characteristic speed corresponding to $\tilde{u}$ labeled e in Fig. 4.3(b). The result approaches a rarefaction from $u_{L}$ down to $\tilde{u}$ and a shock from $\tilde{u}$ up to $u_{R}$; since $u_{M}<\tilde{u}$, the total variation of the solution decreases to $u_{L}+u_{R}-2\,\tilde{u}$.

However, if $u_{L}<\tilde{u}$ as in Fig. 4.3(a), the middle state is eliminated in finite time, resulting in a decrease of total variation to $u_{R}-u_{L}$. It is also possible for a middle state to asymptote to a value \overline{u}\in\Big{(}u_{M}\,,\,\min(u^{*},u_{L},u_{R})\Big{)} such that

[TABLE]

since the speed of a forward shock is determined by the upper flux curve, and the characteristic speed to the right of the center of a rarefaction is found on the lower flux curve. Again, the total variation of the solution decreases. Hence, for Case C, if there is an interaction, the total variation always decreases, and the number of outgoing waves is non-increasing.

4.4 Interactions of Shocks and Expansion Shocks

4.4.1 Case a: Shock - Expansion Shock: $u_{L}<u_{M}$ and $u_{R}<u_{M}$

This sub-case corresponds to case A above, in which a shock necessarily interacts with a rarefaction wave. However, when the rarefaction is replaced by a piecewise constant approximation consisting of expansion shocks, the shock from $u_{L}$ up to $u_{M}$ may not meet the slowest expansion shock on the right. The two waves move apart if the shock has positive speed and the expansion shock has larger speed, or if the shock has negative speed and the expansion shock has either less negative or positive speed. To analyze the situation, we consider the expansion shock from $u_{M}$ to $u_{R}.$ Let $\Lambda_{LM}=\dfrac{f(u_{M})-f(u_{L})}{u_{M}-u_{L}}.$ (i) If $\Lambda_{LM}>0$ , a forward shock with speed $\Lambda_{LM}$ connects $u_{L}$ and $u_{M}$ . If $u_{M}\leq u^{*}$, the expansion shock between $u_{M}$ and $u_{R}$ also has positive speed; however, if $u_{M}>u^{*}$ , the expansion shock could have positive, zero, or negative speed. The shock and expansion shock will move apart only if the expansion shock moves faster than the shock, in which case, the expansion shock has speed $\Lambda_{MR}=(1-\epsilon)\dfrac{f(u_{M})-f(u_{R})}{u_{M}-u_{R}}.$ Let $\tilde{u}_{M}<u_{M}$ be such that

[TABLE]

Hence, the shock and approximating expansion shock(s) do not approach if $u_{L}<u_{M}$ with $\Lambda_{LM}>0$  and $u_{R}$ satisfies

[TABLE]

(ii) For the case when $\Lambda_{LM}\leq 0$, the shock speed $(1-\varepsilon)\Lambda_{LM}$ is reduced due to residual trapping. Define $\overline{u}_{M}<u_{M}$ to be such that

[TABLE]

The shock wave and expansion shock wave do not interact if $u_{L}<u_{M}$ has $\Lambda_{LM}\leq 0$ and $u_{R}$ is such that

[TABLE]

Hence, if we have an expansion shock from $u_{M}$ down to $u_{R}$, where $u_{R}$ does not satisfy (4.4) or (4.5), then the shock and expansion shock collide, producing a single admissible shock from $u_{L}$ to $u_{R}.$

4.4.2 Case c: Expansion Shock - Shock: $u_{M}<u_{L}$ and $u_{M}<u_{R}$

This case is analyzed similarly to Case b, with conditions for the approach or separation of the two waves, analogous to Case C, where a rarefaction wave to the left of a shock may fail to approach the shock because the fastest characteristic in the rarefaction is slower than the shock speed. Correspondingly, when a rarefaction from $u_{L}$ to $u_{M}<u_{L}$ is approximated with one (or more) expansion shock(s), the fastest (right-most) expansion shock connects $u_{E}\in(u_{M}\,,\,u_{L}\,]$ down to $u_{M}$  with speed $\Lambda_{EM}=\sigma_{L}\dfrac{f(u_{E})-f(u_{M})}{u_{E}-u_{M}}.$ If this speed is less than $\Lambda_{RM}=\sigma_{R}\dfrac{f(u_{R})-f(u_{M})}{u_{R}-u_{M}},$ then the two waves fail to interact and all the expansion shocks approximating the rarefaction move away from the shock. Here, $\sigma_{L}=1$ if and only if $\Lambda_{EM}<0,$ and $\sigma_{R}=1$ if and only if $\Lambda_{RM}>0.$

If both waves are moving right, then they interact only if $u_{R}>u_{L}.$ If they do interact, then the result is an admissible shock from $u_{L}$ to $u_{R}$. A similar argument applies to left-moving waves: either they separate, or the result is an admissible shock from $u_{L}$ to $u_{R}$. Consequently, the number of waves either remains at two, with no change in the total variation, or is decreased to one, with a corresponding decrease in total variation.

The overall result of binary interactions between shock waves and expansion shocks is that the total variation and number of waves decreases, but adjacent waves may move apart.

5 Initial Value Problems

Plume migration within a porous aquifer depends on the geometry of the carbon dioxide plume at the end of injection [10], [12]. An analytic solution for a specific idealized CO2 plume is constructed by Hesse, Orr, and Tchelepi [8]. In this section we consider the scalar conservation law (2.4) with a general initial plume of supercritical carbon dioxide,

[TABLE]

in which $u_{0}\in L^{1}(\mathbb{R})\,\cap\,BV(\mathbb{R})\,$ with $\,0\leq u_{0}\leq 1$ .

5.1 Wave-Front Tracking

Dafermos [5] introduced wave-front tracking as a method to construct approximate solutions for scalar, nonlinear partial differential equations. The method has since been greatly generalized to systems of hyperbolic conservation laws [1], [2]. In this section, we describe wave-front tracking, following the approach of LeFloch [15].

In the wave-front tracking algorithm, we first approximate the initial plume shape with a sequence of piecewise constant functions, $u_{0}^{h}(x),\ h>0$, such that

[TABLE]

Each approximation $u_{0}^{h}$ is constructed to have a finite number of discontinuities. The construction of a piecewise-constant solution for short time involves solving the Riemann problems associated with each discontinuity in $u_{0}^{h}.$ Rarefaction waves are replaced by a finite number of expansion shocks of magnitude $h$ or less. When waves meet, we refer to the collision as an interaction. Each interaction results in a Riemann problem in which the initial jump may exceed the threshold $h.$ If the resulting solution is an admissible shock, it is propagated forward without change. If the resulting solution is a rarefaction wave, then (as observed in the previous section) the magnitude is necessarily smaller than $h;$ it is approximated by an expansion shock, traveling with the shock speed of that discontinuity. Continuing in this way, we generate a piecewise constant solution of the conservation law.

In §4, we showed that the number of waves and total variation decreased or remained constant at any interaction. Consequently, since there are finitely many discontinuities initially, there are a finite number of interactions and no accumulation points. Thus, the number of wave interactions and resulting wave-fronts in each $u^{h}$ remains finite for all $t>0$, so the approximations are well defined globally in time. [2].

As observed in the previous section, the total variation is non-increasing, and each approximation $u^{h}(x,t)$ is bounded by $u_{0}^{h}(x)$ . It follows from (5.1) that, at any position and time,  $\inf(u_{0})\leq u^{h}(x,t)\leq\sup(u_{0})$ ; hence,  $\left\|u^{h}(x,t)\right\|_{L^{\infty}}\leq 1$ . We also have  TV\big{(}u^{h}(\cdot,t)\big{)}\leq TV\big{(}u_{0}(\cdot)\big{)}  for all $t>0$ .

Since we have established that there are a finite number of waves, there will be a finite number, $k$, of classical and expansion shocks in $u^{h}$ within \big{[}\,t_{1},t_{2}\,\big{]}, any time interval containing no interaction time. For $m=1,\ldots,k$ , let $y_{m}^{\prime}$ be the speed of propagating shock front $x=y_{m}(t)$ in $u^{h}$ for t\in\big{[}\,t_{1},t_{2}\,\big{]} ; by (2.2), $\left|y_{m}^{\prime}\right|\leq\sup\left|f^{\prime}\right|<\infty$ . The approximate solution to the left/right of wave-front $y_{m}$ is u^{h}\big{(}y_{m}(t)^{\mp},t\big{)}. We following estimate is based on the change in area under the graph of $u^{h}(t)$ due to the motion of individual waves.

[TABLE]

We have shown that conditions for both Helly’s Theorem and the time-dependent version ([2],[15]) are satisfied. Hence, there exists a subsequence of $u^{h}$, which we also label $u^{h}$, and a BV function $u:\mathbb{R}\times\mathbb{R}^{+}\rightarrow[\,0\,,1\,]\,$, such that

[TABLE]

for all $x\in\mathbb{R}$ , $t,\,t_{1},\,t_{2}\in\mathbb{R}^{+}$, and some $\kappa>0$ .

Combining (5.1) with the lower semi-continuity property TV\big{(}u(\cdot,t)\big{)}\leq\displaystyle\mathop{{\lim\inf}}_{h\to 0^{+}}\,TV\big{(}u^{h}(\cdot,t)\big{)} we have TV\big{(}u(\cdot,t)\big{)}\leq TV\big{(}u_{0}(\cdot)\big{)} for all $t\geq 0$ . Similarly, since $u^{h}$ converges to $u$ , it follows that $\,\inf(u_{0})\leq u(x,t)\leq\sup(u_{0})$ . We also have  \left[u^{h}(x,t_{2})-u^{h}(x,t_{1})\right]\rightarrow\big{[}u(x,t_{2})-u(x,t_{1})\big{]}  in $L^{1}_{\text{loc}}$  by (5.1), and it follows from the lower semi-continuity property of norms that  \big{\|}u(x,t_{2})-u(x,t_{1})\big{\|}_{L^{1}}\leq\mathop{{\lim\inf}}_{h\to 0^{+}}\big{\|}u^{h}(x,t_{2})-u^{h}(x,t_{1})\big{\|}_{L^{1}}\,. Finally, from the uniform estimate above, we have \,\big{\|}u(x,t_{2})-u(x,t_{1})\big{\|}_{L^{1}}\leq TV(u_{0})\sup\left|f^{\prime}\right|\big{|}t_{2}-t_{1}\big{|}\, for all $t_{1},t_{2}\geq 0$ .

The wave-front tracking approximations $u^{h}$ are exact solutions of  u^{h}_{t}+\big{(}\sigma\,f(u^{h})\big{)}_{x}=0  since the Rankine-Hugoniot jump condition is satisfied across all classical and expansion shocks. However, we are unable to take the limit as $h\to 0+$ for a pair of reasons: First, we do not have a weak formulation of the Cauchy problem, and second, as $h\to 0,$ the value of $\sigma(u^{h}_{t})$ changes and it is not clear how to formulate the limit $\lim_{h\to 0}\sigma(u^{h}_{t})f(u^{h})_{x}$ in the sense of distributions, which should be $\sigma(u_{t})f(u)_{x}.$ If such problems can be resolved, then establishing that the limit is an entropy solution in the appropriate sense generalized to the model is straightforward.

6 Discussion

The tracking a plume of supercritical carbon dioxide after it has been injected into a deep saline aquifer is modeled by a scalar partial differential equation that has unusual features due to the property of deposition of $CO_{2}$ bubbles as the plume migrates. In the model of Hesse et al [8], this is achieved by reducing the flux by a constant scale as the plume migrates away from a region of space, leaving behind bubbles of sequestered $CO_{2}.$ In this paper, we have explored some interesting properties of the model that fall outside the conventional theory of conservation laws.

The method of characteristics has an interesting twist, due to the presence of two characteristic speeds. Since the switch occurs when either $u_{x}$ or $f^{\prime}(u)$ changes sign, tracking maxima and minima of $u(x,t)$ the solution propagates either as a corner, or as an expanding interval in $x$ over which $u(x,t)$ is constant. Similarly, if a rarefaction wave includes values of $u$ that cross $u=u^{*},$ where $f(u)$ has a maximum, then the rarefaction wave includes a corner, where the slope $u_{x}$ jumps as $u$ crosses $u^{*}.$

In order to define shock waves, we have to generalize the Lax entropy condition in that the admissible behavior of characteristics on either side of the discontinuity has to be interpreted appropriately. A consequence is that each value of $u$ has two possible characteristic speeds, namely $f^{\prime}(u)$ and $(1-\epsilon)f^{\prime}(u).$ The choice depends on the direction of propagation of the wave, so that the choice switches if a shock wave changes direction. To accommodate this behavior, we express admissibility in terms of both families of characteristics.

These phenomena associated with characteristics and shock waves appear when describing the interaction of pairs of waves. We find that shock-to-rarefaction interactions can be complete in finite time, leaving a shock wave, or can persist, resulting in a remaining rarefaction and a shock whose speed approaches characteristic speed.

The asymptotic behavior as $t\to\infty$ shown in Fig. 4.3(b) suggests an unusual rarefaction-shock construction, in which $u=u_{L}$ is connected to $\tilde{u}<u_{L}$ by a rarefaction wave, whose fastest characteristic speed (the speed of the right-most characteristic in the rarefaction fan) is the same as the shock speed of a jump from $\tilde{u}$ to $u_{R}>\tilde{u}.$ This composite wave, in which the shock is characteristic on one side and does not decay, is unusual, because the flux functions are convex, whereas shock-rarefactions are expected to appear only when genuine nonlinearity fails; that is, for non-convex flux functions.

It would be interesting to know how the notion of weak solution can be formulated for equation (2.4). Although it is clear how to treat piecewise smooth solutions, the convergence result from wave front tracking does not guarantee that the limit is a piecewise smooth function, even if the initial data are smooth. In terms of the application, it would be interesting to know whether compactly supported initial data collapses to zero in finite time, signifying the desirable property of complete sequestration in a finite time and over a finite distance. Of particular significance would be an estimate of the maximum time over which this would occur, and the corresponding maximum distance any given plume would migrate before giving up all its $CO_{2}$ to sequestered bubbles.