# Numerical Validation of Analytical Estimates for Propagation of Thermal Waves Generated by Gas-Solid Combustion.

1. Introduction

There is a renewed interest in using combustion to recover medium- or high-viscosity oil. In situ combustion (ISC) is the oldest thermal recovery technique and it has been considered a potential method for the recovery of heavy oil [1, 2]. In this method, air is injected in the porous medium. Heavy and immobile components of the crude oil are used as fuel, producing in-place heat and consequently reducing the viscosity of the oil. Most of the oil is driven toward the producers. An operational advantage of this technique is the abundance of injection gas regardless of location.

Modeling the combustion in porous media presents a number of research challenges, mainly because of the physical complexity of the process. The crude oil is composed of many substances which directly affect the combustion process. Furthermore, the combustion reaction and other thermal reactions, like pyrolysis, change the chemical structure of the components along the process. Due to these reasons, a good model for this technique generally consists of many equations, which hinders its analytical and numerical study. That is why a number of works simplify ISC process to one of its main components, gas-solid combustion; see [3-7] and references therein.

Analytical solutions and estimates are simple and cheap tools which can be used for validation of existing numerical solutions and to obtain information about some aspects of oil recovery before investing in computational codes. They also help to better understand the physical process. Usually the combustion front is regarded as a traveling wave. For one-dimensional models there are two main methods used to study such wave: (1) the first method explores the strong nonlinearity of the Arrhenius factor in the reaction rate, which allows neglecting the reaction rate as soon as the temperature decreases [8]. Other works use this method to obtain estimates on the combustion wave speed, temperature, and oxygen concentration; see [3-5, 9] and references therein. (2) The second method considers that the reaction is active for all temperatures but heat losses are negligible; see [7] and references therein. The last method also allows obtaining physically consistent process parameters such as combustion wave speed and combustion front temperature. The previously mentioned works consider one-dimensional models that do not take the pressure equation into account. None of these is directly applicable to two-dimensional models, nor to models containing the pressure equation. The main goal of this work is to verify whether the analytical estimates proposed in [7] can be applied to more general two-dimensional nonhomogeneous models.

Other works address two-dimensional combustion. For example in [10], two-dimensional simulations were performed and it was observed that the combustion front assumes a curved profile when the domain has strongly spatially correlated characteristics, whereas it is almost linear if heterogeneity is random. In [11] the compositional flow in porous media was considered to model forward ISC. The viscous fingering effect was observed for heavy oil oxidation. There is a number of other works that study the filtration combustion in forward (e.g., [12]) or reverse (e.g., [13, 14]) regimes. Although many papers deal with multidimensional combustion models, there is a lack of analytical studies on this topic. Similar to [10], fingering instabilities observed in the present work appear only for correlated nonhomogeneous medium, when fingers increase in length over time.

In this work we follow [3,7] and consider the combustion of immobile fuel. In Section 2 we present a two-dimensional gas-solid model containing an equation for pressure, balance laws for energy, oxygen content, fuel molar mass, and total gas mass, where the gas density is described by the ideal gas law and the average gas speed by Darcy's law. In Section 3 we show that the model proposed in this work is indeed a generalization of the one-dimensional model proposed in [3]. In this section we also review and generalize some analytical results obtained in [7]. In order to simulate the gas-solid process we use the finite element method described in Section 4. Section 5 presents numerical results, which are compared with analytical estimates. Further numerical examples are performed in order to validate the numerical model. Finally, Section 6 summarizes final conclusions.

2. Model

Air is injected into a porous medium initially filled with a nonreactive gas and immobile fuel, which can be solid or liquid at low saturations. We consider a combustion reaction that consumes oxygen and fuel, generating heat. We assume that only a small part of the void space is occupied by the fuel, so its consumption does not affect the porosity of the matrix. We consider the local thermal equilibrium resulting in equal temperatures for all phases (fuel, gas, and rock). Heat losses are neglected, which is a reasonable hypothesis for ISC in field conditions.

Some other assumptions are considered: heat transfer due to radiation and energy source due to pressure variation are neglected; the work related to surface and body forces is also neglected; the ideal gas law is the equation of state for the gas phase; gas heat capacity is negligible when compared to the heat capacity of the rock. To simplify the model, we neglect changes in the effective molecular weight of the gas phase due to the reaction by approximating this quantity by a constant. We do not consider the effect of gravity in the flow and the dispersion terms in the diffusion-dispersion tensor.

Thermodynamic and transport properties such as thermal conductivity of the medium, heat capacity of the rock, and specific heat capacity of the gas are considered constants. Generally, the porosity is related with permeability through polynomial relations, it also changes within the combustion reaction. For simplicity, we consider isotropic porous medium with constant porosity and simplified combustion reaction: C + [O.sub.2] [right arrow] C[O.sub.2].

Let t [member of] [R.sup.+] and x [member of] [R.sup.2] be the time and space coordinates. The primary dependent variables are the temperature, T(x, t) (Kelvin), the oxygen concentration in terms of molar fraction, Y(x,t) (moles of oxygen per moles of gas), the molar concentration of fuel, [[rho].sub.f] (x, t) (moles of fuel per cubic meter), and molar density of gas, [rho] (x, t) (moles of gas per cubic meter).

The system of equations extends the one-dimensional model from [6] by considering variable pressure, two-dimensional domain, and Darcy's law, based on [15]. Parameter names and typical values are given in Table 1. We present the following equations under the assumptions described above. The energy balance law is

[C.sub.m][[partial derivative].sub.t] T + [c.sub.g][nabla] x ([rho]Tu) = [k.sub.T][DELTA]T + QW, (1)

where [C.sub.m] is the volumetric heat capacity of the rock, [c.sub.g] is the specific heat capacity of the gas, u is the average gas speed, [k.sub.T] is the thermal conductivity, Q is the heat released due to combustion of each mole of fuel, and W is the combustion rate. The molar mass balance of oxygen in the gas phase is

[phi][[partial derivative].sub.t] (Y[rho]) + [nabla] x (Y[rho]u) = [phi][d.sub.m][nabla] x ([rho][nabla]Y) - W, (2)

where [phi] is the porosity of the medium and [d.sub.m] is the molecular diffusion coefficient. The molar mass balance of fuel is

[[partal derivative].sub.t][[rho].sub.f] = -W. (3)

The molar mass balance of total gas is

[phi][[partial derivative].sub.t] [rho] + [nabla] x ([rho]u) = 0. (4)

The reaction rate is described by the Arrhenius law combined with the linear Law of Mass Action [16].

W = [k.sub.p]Y[[rho].sub.f] p/[p.sub.atm] exp (-E/RT), (5)

where [p.sub.atm] is the reference pressure, E is the activation energy, [k.sub.p] is the preexponential factor, and R is the ideal gas constant. The gas phase equation of state is the ideal gas law:

P = [rho]RT. (6)

The average gas speed is given by Darcy's law:

u = - ([kappa]/[mu]) [nabla]p, (7)

where [mu] is the dynamic gas viscosity and [kappa] is the absolute permeability of the porous medium.

2.1. Dimensionless Equations. In order to simplify the mathematical analysis, we introduce dimensionless dependent and independent variables (denoted by tildes)

[mathematical expression not reproducible] (8)

and some reference quantities (denoted by stars)

[mathematical expression not reproducible] (9)

where [T.sub.res], [[rho].sup.res.sub.f], and [p.sub.res] are the initial reservoir conditions for temperature, molar density of fuel, and gas pressure, respectively; [Y.sub.inj] is the oxygen fraction in the injected air; and [kappa] and [bar.[mu]] are the average permeability of the porous medium and average dynamic gas viscosity, respectively.

Using (8), (9) and omitting tildes, (1)-(7) are written in dimensionless form as

[[partial derivative].sub.t][THETA] + [V.sub.T][nabla] x ([rho]([THETA] + [[THETA].sub.d])u) = 1/[Pe.sub.T] [DELTA][THETA] + [PHI], (10)

[[partial derivative].sub.t] (Y[rho]) + [sigma][nabla] x (Y[rho]u) = 1/Pe [nabla] x ([rho][nabla]Y) - [sigma][PHI], (11)

[[partial derivative].sub.t][[rho].sub.f] = -[PHI], (12)

[[partial derivative].sub.t][rho] + [sima][nabla] x ([rho]u) = 0, (13)

p = [p.sub.d] ([rho] (1 + [THETA]/[[THETA].sub.d]) - 1), (14)

u = - [lambda][nabla]p, (15)

[PHI] = [t.sup.*]/[[rho].sup.*.sub.f] W = [[rho].sub.f] Y (1 + p/[p.sub.d]) exp (E/[[THETA].sub.d] - E/[THETA] + [[THETA].sub.d]), (16)

where the dimensionless constants are

[mathematical expression not reproducible] (17)

where [[THETA].sub.d] and [p.sub.d] are the scaled reservoir temperature and pressure, [Pe.sub.t] and Pe are the Peclet numbers for thermal and mass diffusion, [V.sub.t] is the dimensionless thermal wave speed, a is the fuel to oxygen concentration rate, and E is the scaled activation energy. In this work we consider that permeability [kappa] depends on x and viscosity [micro] depends on [THETA]; thus the scaled gas mobility function (see [[17], Sec. 4.9]) is defined as [lambda]([THETA], x) = ([bar.[mu]]/[mu] ([THETA])) ([kappa](x)/[bar.[kappa]]). In this work [mu] is considered constant, the general case was left for the future work.

3. Analytical Estimates for a Simplified Model

In this section we simplify system (10)-(16) in order to obtain analytical estimates. Following [6, 7] we restrict ourselves to the one-dimensional case and consider small pressure variation when compared to reservoir reference pressure. Mathematically, this is equivalent to considering the constant prevailing pressure [bar.p] and removing Darcy's law from the system. To simplify notation we denote p = 1 + [bar.p]/[p.sub.d]. The new model considers gas density as a function of temperature given by the ideal gas law. Under these hypotheses, the dimensionless system can be rewritten as

[[partial derivative].sub.t][THETA] + [V.sub.T][[partial derivative].sub.x] ([rho] ([THETA] + [[THETA].sub.d]) u) = [Pe.sup.-1.sub.T] [[partial derivative].sub.xx] [THETA] + [PHI], (18)

[[partial derivative].sub.t] (Y[rho]) + [sigma][[partial derivative].sub.x] (Y[rho]u) [Pe.sup.-1] [[partial derivative].sub.x] ([rho][[partial derivative].sub.x]Y) - [sigma][PHI], (19)

[[partial derivative].sub.t][[rho].sub.f] = -[PHI], (20)

[[partial derivative].sub.t][rho] + [sigma][[partial derivative].sub.x] ([rho]u) = 0, (21)

[rho] = P[[THETA].sub.d]/([THETA] + [[THETA].sub.d]), (22)

[PHI] = P[[rho].sub.f]Y exp (E/[[THETA].sub.d] - E/[THETA] + [[THETA].sub.d]). (23)

When the prevailing pressure is considered to be equal to the atmospheric pressure, we obtain P = 1. In this case the model is equivalent to ones studied in [3, 6, 7]. Differently from [7] we do not restrict ourselves to the limit case with large values of a. We consider

Initial conditions: t = 0, x [greater than or equal to] 0: [THETA] = Y = 0

u = [[rho].sub.f], (24)

Injection conditions: t > 0, x = 0: [THETA] = [[rho].sub.f] = 0,

Y = 1, u = [u.sub.inj]. (25)

The solution of system (18)-(22) with general boundary conditions can be described as a sequence of waves: fuel concentration wave, thermal wave, combustion wave, and gas composition wave separated by constant states; see [6]. Considering conditions (24)-(25) after a sufficiently large time only thermal and combustion waves are presented. We consider the case in which the thermal wave is slower than the combustion wave, which is reasonable in ISC.

We consider the combustion wave as a traveling wave solution of system (18)-(22) moving with constant speed V > 0. We indicate the combustion wave neighbor constant states as burned (indicated by letter b) and unburned (indicated by letter u) and define the traveling variable [xi] = x - Vt. The left (burned) limit states for the combustion wave are

[xi][right arrow] -[infinity]: [[THETA].sup.b] > 0. [Y.sup.b] = 1, [[rho].sup.b.sub.f] = 0, [u.sup.b] = f ([[rho].sup.b], [u.sub.inj]), (26)

where f is a known function that can be obtained from the analysis of the thermal wave; see Appendix A. Variable [[rho].sup.b] is obtained by substituting the unknown [[THETA].sup.b] into (22). The right (unburned) states corresponds to

[xi] [right arrow] +[infinity]: [[THETA].sup.u] = [Y.sup.u] = 0,

[[rho].sup.u.sub.f] = 1, [u.sup.u] > 0. (27)

Only [u.sup.u] is unknown. The variable values at these states can be determined by analysis of the initial and boundary conditions; see Appendix A.

Under condition (24) the sequence of waves moves in a positive direction of the x-axis; then the conditions to the left of the thermal wave are given by the injection conditions. Since the thermal wave is slower than the combustion wave, the conditions to the left of the combustion wave are equal to the conditions to the right of the thermal wave.

The equations describing the traveling wave solution are obtained from (18) to (22) substituting the derivatives [partial derivative]x and [[partial derivative].sub.t] by d/d[xi], and -Vd/d[xi]:

-V[THETA]' + [V.sub.T]P[[THETA].sub.d]u' = [Pe.sup.-1.sub.T] [THETA]" + [PHI], (28)

-V ([Y.sub.[rho]])' + [sigma] (Y[rho]u)' = [Pe.sup.-1] ([rho]Y')' - [sigma][PHI], (29)

-V[[rho]'.sub.f] = -[PHI], (30)

-V[rho]' + [sigma] ([rho]u)' = 0, (31)

where the prime (') means derivative in [xi]. We substitute (30) into (28) and (29) obtaining

[mathematical expression not reproducible]. (32)

Integrating these equations in (-[infinity], +[infinity]), substituting the limit states (26)-(27), and considering that derivatives vanish in the boundary [THETA]' = Y' = 0, we obtain the estimates for the combustion wave temperature [[THETA].sup.b], combustion wave speed V, and gas speed [u.sup.u]:

[mathematical expression not reproducible]. (33)

Analytical formulas in (33) are generalizations of equivalent ones present in [7], where was considered atmospheric pressure [bar.p] = 0 and [sigma] [much greater than] 1. In Section 5, the values of parameters [[THETA].sup.b], V, and [u.sup.u], from (33), are compared with values obtained from numerical simulations using the general model described by (10)-(16). These parameters are also used in Appendix B to validate the numerical model from next section.

4. Numerical Model

In this section we present the numerical model used to solve system (10)-(16) in an open bounded domain [OMEGA] [member of] [R.sup.2] with boundary [GAMMA]. The spatial discretization uses FEM; time discretization uses Crank-Nicolson method and the resulting implicit scheme is solved using Newton's method.

In order to allow the implementation of the numerical method we use (13) to obtain a nonconservative form of (11):

[[partial derivative].sub.t]Y + ([sigma]u - 1/Pe [nabla][rho]/[rho]) x [nabla]Y = 1/Pe [DELTA]Y - [sigma] [PHI]/[rho]. (34)

In what follows we shorten the notation by grouping the variables as U = ([THETA], Y, [[rho].sub.f], [rho]) and consider [OMEGA] boundary as union of Dirichlet type boundary [[GAMMA].sub.1] and zero-flow Neumann type boundary [[GAMMA].sub.2] resulting in r = [bar.[[GAMMA].sub.1] [union] [[GAMMA].sub.2]]. We always use boundary conditions compatible with the initial condition:

U[|.sub.t=0] = [U.sub.0] = ([U.sub.01], [U.sub.02], [U.sub.03], [U.sub.04]) = ([[THETA].sub.0], [Y.sub.0], [[rho].sub.f0], [[rho].sub.0]). (35)

Implemented boundary conditions for t > 0 are

U (x, t) = g (x), [for all] x [member of] [[GAMMA].sub.1], (36)

[nabla]U (x, t) x n (x) = 0, [for all]x [member of] [[GAMMA].sub.2]. (37)

In this work we apply the standard FEM on system (10)-(16). This simple approach meets our goal, although the use of specific methodologies could improve the precision and convergence of the numerical procedure. In Appendix B.2 we show that the numerical instabilities appearing in the simulations do not affect the structural stability of the combustion front and do not compromise present results.

4.1. Weak Formulation and Discretization. Let [mathematical expression not reproducible] be the space of test functions. We use <u, v> to denote the inner product between u, v [member of] [L.sup.2] ([OMEGA]). In the weak formulation, solving the system described in (10)-(16) is equivalent to find U such that for all v [member of] V and t > 0 it satisfies

[mathematical expression not reproducible]. (38)

together with boundary conditions (36) and initial conditions

<[U.sub.i][|.sub.t=0], [v.sub.i]> = <[U.sub.0i], [v.sub.i]>, i = 1, 2, 3, 4. (39)

Let [V.sup.h] [member of] V be a finite element space with basis [beta] = [{[[phi].sub.i]}.sup.m.sub.i=1]. After applying Standard Galerkin and Crank-Nicolson methods to system (38)-(39), the matrix formulation of the weak problem reads: for each k = 0, 1, ..., n, determine [[alpha].sup.k] [member of] [R.sup.m] such that

[mathematical expression not reproducible], (40)

where M [member of] [R.sup.mxm] and F([alpha]), [F.sub.0] [member of] [R.sup.m] are described in Appendix C and

[U.sup.k.sub.h] = [m.summation over (i=1)] [[alpha].sup.k.sub.i] [[phi].sub.i] + g. (41)

For each time step k > 0 we use Newton's method to obtain zeros of the function [f.sup.k] considering the previous time step [[alpha].sup.k-1]. We use the adaptive time step to speed up the simulations, where the size of next step was based in the convergence of Newton's method with tolerance equal to [10.sup.-6] for the relative error using the Euclidean norm.

5. Numerical Experiments

Using parameter values given in Table 1 and (9)-(17) we obtain the value of the characteristic quantities and dimensionless parameters: [[rho].sup.*] = 32.94 mol/[m.sup.3], [[delta].sup.*.sub.T] = 74.40 m, [t.sup.*] = 1.470 x [10.sup.6] s, u = 4.306 x [10.sup.-4] m/s, [x.sup.*] = 11.77 m, [[THETA].sub.d] = 4.97, [p.sub.d] = 1, [Pe.sub.T] = 216.56, Pe = 47.10, [V.sub.T] = 0.0243, [sigma] = 179.27, and E = 93.77.

System (10)-(16) was solved using FEM in a rectangular domain discretized using square elements with 200 divisions in horizontal x axis and 50 divisions in the vertical y axis. The algorithm uses an adaptive time step with the initial time step value equal to [DELTA]t = 4 x [10.sup.-7]. The right and left sides are modeled using the Dirichlet boundary conditions [g.sub.i] = [U.sub.0i], for each i = 1, 2, 3, 4. The upper and lower boundaries are modeled using Neumann zero-flow conditions. The initial conditions are

[mathematical expression not reproducible], (42)

and the initial value of p is obtained by substituting (42) into (14). Air is injected from the left side corresponding to zero on the horizontal axis.

Generally, the permeability function [kappa] depends on space coordinates. In this paper we study two possibilities for [kappa]. The first case corresponds to a homogeneous porous medium with constant permeability [kappa] = [bar.[kappa]] given in Table 1. The second kind of permeability field is more realistic from the physical point of view and represents typical reservoir conditions [18]. A number of works use similar permeability fields in numerical experiments [11, 19, 20]. This case corresponds to a nonhomogeneous porous medium with log-normal permeability field, similar to the one used in [11]. Permeability values vary from 0.40 x [10.sup.-12] [[m.sup.2]] to 9.97 x [10.sup.-12] [[m.sup.2]] and the mean value is equal to k, given in Table 1. This field is given by

k (x) = [e.sup.a+bZ(x)], (43)

where a = -27.693 and b = 0.35297 are the mean and standard deviation of the permeability's natural logarithm, respectively and Z is a standard normal distribution with constant values in each element; see Figure 1. Figure 2 shows the histogram and probability density function for this distribution.

5.1. Parameter Values from Numerical Simulations. In order to obtain the combustion wave speed for two-dimensional simulations we assume that the combustion wave moves only in the x direction. This is exactly the case for homogeneous porous media. For the nonhomogeneous case, the x component of gas velocity [u.sub.x] is more than 11 times greater than [absolute value of ([u.sub.y])] as can be seen in Figure 9.

At time t we track the position of the combustion wave for each value of y by the value of x for which [p.sub.f] = 0.2. This approach is accurate since there are no diffusion nor advection terms in the fuel balance equation (12). The combustion wave speed [V.sup.s] (t) is the average of the velocities of each y and [sigma] ([V.sup.s])(f) denoting the standard deviation of these velocities.

Analogously to the combustion wave speed, at time t for each value of y, we consider the combustion wave temperature [[THETA].sup.b.sub.S] [(y, t) as the maximum value of [THETA] varying v. In this case, the combustion wave temperature, [[THETA].sup.b.sub.s] (t), is the average of the values [[THETA].sup.b.sub.s] (y, t) for all y and the standard deviation of such values is denoted by [sigma] ([[THETA].sup.b.sub.s]) (t). Finally, we consider that the gas speed downstream the combustion wave, [u.sup.u.sub.S](t), can be obtained by the average of u over the right side of the spatial domain.

5.2. Parameter Values for Analytical Formulas. In order to use (33) we need the value of the prevailing pressure, [bar.p], and the injection gas speed, [u.sub.inj], which are not constant for system (10)-(15). Here we study three different choices for the prevailing pressure:

(1) Minimum pressure in the reservoir, which coincides with the native reservoir pressure.

(2) Maximum pressure in the reservoir, which coincides with the injection pressure.

(3) The average pressure measured at the top of the combustion front, which is a numerical estimate for the pressure in the burned state of combustion wave.

We estimate the parameter [u.sub.inj] in two different ways:

(A) Injection gas speed is equal to the average injection gas speed ([u.sup.s.sub.inj]):

[u.sub.inj] (t) = [u.sup.s.sub.inj] (t). (44)

(B) Injection gas flux is equal to the average injection gas flux:

[u.sub.inj] (t) = [p.sup.s.sub.inj] (t) + [p.sub.d] / [bar.p] (t) + [p.sub.d] [u.sup.s.sub.inj] (t). (45)

where [p.sup.s.sub.inj] (t) is the average pressure of injection in the simulation. This equation was obtained using (14) and (22).

Both [u.sup.s.sub.inj] and [p.sup.s.sub.inj] (t) are obtained as the [u.sup.u.sub.s] parameter described in Section 5.1. In the simulations we use constant injection pressure, so estimates (A2) and (B2) are equal.

5.3. Simulation Results: Homogeneous Porous Medium. In the homogeneous case, the solution is constant in the y direction. In order to improve visibility we plot the numerical result corresponding to the cross-section at y = 0. The simulation results for dimensionless variables [THETA], Y, [[rho].sub.f], [rho], p, and [u.sub.x] corresponding to 10 and 50 days are plotted in Figures 3 and 4, respectively. The analytically estimated parameter [[THETA].sup.b], given in (33), is also plotted showing good agreement with the value of [THETA] obtained numerically. The combustion and thermal waves are present for both plots. After 50 days one can observe small numerical oscillations arising in variables Y and [[rho].sub.f] as shown in Figure 4. In Appendix B.2 we show evidence that these oscillations do not affect the structural stability of the combustion front.

In Figures 5, 6, and 7 we compare the values of parameters [[THETA].sup.b.sub.s], [V.sup.S], and [u.sup.u.sub.S] obtained from numerical simulations described in Section 5.1 with analytical estimates from (33). Parameters [u.sub.inj] and [bar.p] were chosen as described in Section 5.2. Combustion wave temperature in the simulated model is the same as in the simplified one. This means that parameter [[THETA].sup.b] can be used to predict combustion wave temperature with good accuracy. Notice that this estimate does not depend on either [u.sub.inj] or [bar.p].

In Figure 6, we see that the estimates corresponding to case (B) are more accurate than the ones corresponding to case (A). All estimates (B1), (B2), and (B3) give good approach to the combustion wave velocity. Among all considered cases, only case (B1) presents good estimates for [u.sup.u], as shown in Figure 7. We conclude that the approximation (B1) provides a better description of the simulations presented in this work.

5.4. Simulation Results: Nonhomogeneous Porous Medium. The simulation results for the nonhomogeneous case at 10 days and 50 days are shown in Figure 8 for variables [THETA], Y, [[rho].sub.f], [rho], and p and Figure 9 for variables u = [absolute value of (u)], [u.sub.x], and [absolute value of ([u.sub.y])]. One can observe the fast propagating combustion wave, corresponding to sudden variation in variables [THETA], Y, and [p.sub.f] and the slow thermal wave corresponding to variation in variable [THETA]. Interface instabilities (fingering) can be seen in the combustion front of the solution for [THETA], Y, and [[rho].sub.f]. We stress that since we simulate single phase flow, these instabilities are not related to the process called viscous fingering, which occurs when a less viscous fluid meets a more viscous one. The fingering instabilities seem to increase in length from 10 to 50 days, as can be seen in Figure 8. This effect is expected to be even greater in more realistic models since gas speed increases from upstream to downstream of the combustion wave. This relation between gas speed and temperature is due to the dynamic gas viscosity, which increases as temperature decreases; see Sutherland's formula in [21].

We can see that horizontal velocity [u.sub.x] is higher than vertical velocity [u.sub.y] indicating that the flow occurs mainly in the horizontal direction. The results present numerical instabilities near the top of the combustion wave. This kind of error is due to the heterogeneity of the medium and diminishes if we use a more refined mesh. Although it cannot be clearly seen in Figure 8, numerical instabilities similar to those observed for the homogeneous porous medium appear near the right side of the domain. The numerical results from Appendix B.2 show evidence that these instabilities do not compromise the solution.

As done in the previous section, in Figures 10, 11, and 12 we compare parameter values obtained from the numerical simulations with the corresponding analytical estimates obtained from (33), where [u.sub.inj] and [bar.p] are chosen as described in Section 5.2.

In Figure 10 we plot the analytical estimate of [[THETA].sup.b], numerically obtained temperature [[THETA].sup.b.sub.S], and its standard deviation [mu] ([[THETA].sup.b.sub.S]). The relation [absolute value of ([[THETA].sup.b] - [mu] ([[THETA].sup.b.sub.S]))] < [sigma] ([[THETA].sup.b.sub.S]) is achieved and after some time the average of combustion wave temperature remains near [[THETA].sup.b].

In Figure 11 we see that combustion wave speed estimates corresponding to case (B) are more accurate than ones corresponding to case (A). All estimates (B1), (B2), and (B3) approach the combustion wave velocity. Furthermore, relation [absolute value of (V(B) - [V.sup.S])] < [sigma]([V.sup.S]) is satisfied for all cases (B) and (A2). Figure 12 shows that only case (B1) presents good accuracy estimating [u.sup.u.sub.S]. As in homogeneous case the approximation (B1) describes better the simulation results.

We present the relative error for main parameters and cases in Table 2. The errors were calculated at each time step for t > 10 days. Each set of errors of the average and standard deviation were evaluated. This data evidences the conclusion that (B1) is a good parameter choice.

6. Conclusions

Previous work analyzed traveling wave solutions and obtained analytical formulas describing combustion wave temperature, velocity, and gas velocity for the one-dimensional gas-solid combustion model. In order to apply these estimates from [7] some assumptions on the reservoir prevailing pressure and gas injection speed needed to be considered. The best fit was obtained using average gas injection flux and minimum reservoir pressure measured on the production side. For numerical validation we presented a two-dimensional generalization of the gas-solid model studied in [3, 6, 7]. The numerical implementation used the finite element method for space discretization and Crank-Nicolson scheme for time discretization.

For simulations using the nonhomogeneous porous medium, fingering instability phenomena were observed. In order to verify that these phenomena are not generated by numerical errors simulations using perturbed homogeneous porous medium were performed.

Appendix

A. Noncombustion Waves

From the simulations presented in Section 5 one can clearly see that noncombustion waves appear in the solution. Studying these waves is important to obtain the analytical estimates presented in Section 3. Here we follow [6], where such an analysis was performed for system (18)-(22) in the hyperbolic framework, i.e., in the absence of the reaction term ([PHI] = 0). Here we derive the relation between parameters [[THETA].sup.b] and [u.sup.b] used in Section 3.

For a large space scale, diffusion effect can be disregarded; see [6]. Using this assumption and disregarding the reaction term we arrive at a simplified version of system (18)-(22) (here [theta] = ([THETA] + [[THETA].sub.d])/(P[[THETA].sub.d])):

[[partial derivative].sub.t][theta] + [V.sub.T][[partial derivative].sub.x] ([rho][theta]u) = 0, (A.1)

[[partial derivative].sub.t] (Y[rho]) + [sigma][[partial derivative].sub.x] (Y[rho]U) = 0, (A.2)

[[partial derivative].sub.t][[rho].sub.f] = 0, (A.3)

[[partial derivative].sub.t][rho] + [alpha][[partial derivative].sub.x] ([rho]u) = 0, (A.4)

[rho][theta] = 1. (A.5)

Similar to what is done in Section 4, we consider solution U = ([theta], Y, [[rho].sub.f], u). Thus, the characteristic velocities of system (A.1)-(A.4) and corresponding characteristic vectors are

[mathematical expression not reproducible]. (A.6)

It is easy to see that the characteristic speeds are constant along the integral curves defined by characteristic vector fields. Thus the corresponding waves are contact discontinuities; see [22]. In ascending order of speed, these waves are the fuel composition wave (only [[rho].sub.f] varies), the thermal wave (only [theta] and u vary), and the gas composition wave (only Y varies); see Figure 13.

For each fixed injection state [U.sub.inj] = ([[??].sub.inj], [Y.sub.inj], [[rho].sup.inj.sub.f], [u.sub.inj]) the integral curve y corresponding to the thermal wave is a line because

[mathematical expression not reproducible]. (A.7)

and, by the definition of integral curves [22], [V.sup.[THETA]] is tangent to [gamma]. State [U.sup.b] is connected to [U.sub.inj] through thermal wave if [U.sub.inj] [member of] [gamma], i.e., [U.sup.b] = [U.sub.inj] + t[V.sup.[THETA]]. Using (A.5) and (A.7),

[u.sup.b] = 1/[[rho].sup.b] + [V.sub.T]/[sigma] / 1/[[rho].sub.inj] + [V.sub.T]/[sigma] [u.sub.inj]. (A.8)

In Section 3, this relation corresponds to f ([[rho].sup.b], [u.sub.inj]) for [[rho].sub.inj] = P, where [[rho].sup.b] is defined as function of [[THETA].sup.b] by (22).

B. Validation of the Numerical Implementation

B.1. Experiment with Constant Pressure. For this test, the program was modified to disregard Darcy's law. The balance equation for variable u has the accumulation term that does not depend on u, making direct simulations difficult. We follow [7], where the same difficulty was observed. We combine (18), (21), and (22) yielding

[mathematical expression not reproducible]. (B.1)

For each time step, the value for variable u is obtained from (B.1) using other variable values from the previous time step. Thus, only (18)-(20) are used by the finite element method described in Section 4.

For this simulation we used the same parameter values as in [7]: [[THETA].sub.d] = 4.97, [Pe.sub.T] = 6180, Pe = 1344, [V.sub.T] = 0.0243, [sigma] = 179, E = 93.8, r = 1, and [bar.p] = 0. The simulated domain has dimensionless length L = 10 with elements of size [increment of x] = 0.01 and adaptive time steps with initial value [DELTA]t = [10.sup.-4]. The simulations results are plotted in Figure 14, where we used [x.sup.*] = 62.9 m and [t.sup.*] = 1.470 x [10.sup.6] s. In [7] the same model was solved numerically using finite difference schemes and analytically using the singular perturbation technique. These results are visually indistinguishable from ones presented in [7].

For the parameter values from the simulation we consider [[THETA].sup.b] as the biggest value of [THETA]; the value of V was estimated using wave speed corresponding to the value [[rho].sub.f] = 0.2 and [u.sup.u] was taken as the value of u at the right boundary. Parameter values given in (33) almost coincide with those obtained in numerical simulations as shown in Table 3.

B.2. Structural Stability and Sensitivity Analysis. Since the value of parameter a in (11) is big for gas injection, the influence of the advection term is increased, which impairs the use of standard FEM with Dirichlet condition in outflow boundary. This is the source of numerical instabilities observed in the simulations presented in previous sections. Next, we show evidence that such instabilities do not influence the structural stability of the combustion front.

We perform sensitivity analysis in all relevant parameters over the homogeneous simulation described in Section 5 in order to validate the numerical implementation. For large times the solution behaves similar to the one described in Section 5.3. The combustion front profile stays approximately constant along the y-direction as expected.

Perturbation on the Initial Conditions in Variables [p.sub.f], [theta], and Y. We considered perturbation in initial conditions in variables [mathematical expression not reproducible] and Y as plotted in Figure 15. Simulations with perturbations in initial conditions of similar shape as plotted in Figure 15 and different size were realized always resulting in stable combustion wave.

The simulation results at 10 and 50 days for perturbed [mathematical expression not reproducible] are shown in Figure 16 for variables [THETA], Y, [[rho].sub.f], [u.sub.x], and [absolute value of ([u.sub.y])]. In each point of the domain, [absolute value of ([u.sub.y])] is less than a tenth of [u.sub.x]. The perturbation of the initial condition creates small perturbation in the combustion front propagation, which becomes smoother in time and approaches the unperturbed solution. The results for perturbed initial data in Y are similar.

Numerical solution corresponding to different initial reservoir temperatures was made for T = 355 K and T = 385 K. When compared to reservoir temperature considered in the paper (T = 370 K), we observe that for lower temperature (T = 355) the reaction is slower. As a result the oxygen penetrates deeper in the reservoir, combustion front speed decreases and the problem shows more hyperbolic behaviour resulting in small numerical oscillations in variables [[rho].sub.f] and Y. For higher temperature (T = 385 K), the reaction rate is higher and the combustion front speed increases.

Perturbation on the Boundary Condition in Variable p. Consider the function [p.sub.0] (x, y) = (4 - x) (1 + cos (4[pi]y)/100)/2 as the initial condition for p. Since we use Dirichlet boundary conditions compatible with the initial conditions, the left boundary condition for p is [p.sub.0] (0, y). It represents a relative perturbation of 1% in the boundary condition of the homogeneous case studied in Section 5.3.

The simulation results at 10 and 50 days are shown in Figure 17 for variables [THETA], Y, [[rho].sub.f], [u.sub.x], and [absolute value of ([u.sub.y])]. In this case, [u.sub.y] is approximately zero in the entire spatial domain, except in a layer near the left boundary. Since the combustion wave moves in the positive direction of the x-axis, the effect of great velocities uy is reduced for large times. The combustion front becomes smoother in time and approaches the not perturbed solution. Simulations with similar perturbations in p were performed resulting in analogous results.

Perturbation of the Permeability Field. In this simulation we consider the heterogeneous permeability field with k close to the homogeneous [bar.[kappa]]. We use the log-normal permeability field with values varying from 0.97 x [10.sup.-12] [[m.sup.2]] to 1.10 x [10.sup.-12] [[m.sup.2]] and the mean value equals to [bar.[kappa]], given in Table 3. In this case, (43) holds with a = -27.631, b = 0.014241, and Z as in Figure 1.

The simulation results at 10 and 50 days are shown in Figure 18 for variables [THETA], Y, [[rho].sub.f], [u.sub.x], and [absolute value of ([u.sub.y])]. The combustion front profile of the obtained solution is almost constant along y direction and the position of the wave is very similar to the one observed in Figures 3 and 4.

C. Discretization Terms

This appendix shows expressions for M, F([alpha]), and [F.sub.0] from Section 4. To obtain decoupled equations we choose integers 0 = [m.sub.0] < [m.sub.1] < [m.sub.2] < [m.sub.3] < [m.sub.4] = m and use a basis [beta] whose functions satisfy

[[phi].sub.i] = [[phi].sub.i][e.sub.k], if [m.sub.k-1] < i [less than or equal to] [m.sub.k], k [member of] 1, 2, 3, 4. (C.1)

Substituting (C.1) in (41), we define

[mathematical expression not reproducible], (C.2)

where the dependence of (x, t) was omitted. Using (14)-(16) we define

[mathematical expression not reproducible]. (C.3)

We note that [U.sub.h] satisfies (38) for all [[phi].sub.i] [member of] [beta]; in other words

[mathematical expression not reproducible], (C.4)

that give us the vector F([alpha]). Observe that if [m.sub.0] < i [less than or equal to] [m.sub.1] we have

[mathematical expression not reproducible] (C.5)

and similar equations for other ranges of i. Using matrices [mathematical expression not reproducible] we define matrix M = diag [{[M.sub.i]}.sup.4.sub.i=1] [member of] [R.sup.mxm].

Since [U.sub.h] satisfies (39) for all [[phi].sub.i] [member of] [beta], then for [m.sub.0] < i [less than or equal to] [m.sub.1] we have

[[m.sub.1].summation over (j=[m.sub.0])] [M.sub.ij] [[alpha].sub.j] (0) = <[[THETA].sub.h][|.sub.t=0], [[phi].sub.i]> = <[[THETA].sub.0], [[phi].sub.i]> = [F.sub.0i] (C.6)

and similar equations for other ranges of i, which defines the vector [F.sub.0].

http://dx.doi.org/ 10.1155/2017/1806052

Competing Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors wish to thank Professor Johannes Bruining for his help in modeling of the problem. The first author was supported in part by CNPq. The second author was supported in part by FAPEMIG under Grant APQ 01377/15.

References

[1] T. C. Boberg, Thermal Methods of Oil Recovery, John Wiley & Sons, New York, NY, USA, 1988.

[2] M. Prats, Thermal Recovery, SPE, New York, NY, USA, 1982.

[3] I. Y. Akkutlu and Y. C. Yortsos, "The dynamics of in-situ combustion fronts in porous media," Combustion and Flame, vol. 134, no. 3, pp. 229-247, 2003.

[4] A. A. Bayliss, B. J. Matkowsky, and D. A. Schult, "Traveling waves in natural counterflow filtration combustion and their stability," SIAM Journal on Applied Mathematics, vol. 58, no. 3, pp. 806-852, 1998.

[5] A. P Aldushin, I. E. Rumanov, and B. J. Matkowsky, "Maximal energy accumulation in a superadiabatic filtration combustion wave," Combustion and Flame, vol. 118, no. 1-2, pp. 76-90, 1999.

[6] A. J. de Souza, I. Y. Akkutlu, and D. Marchesin, "Wave sequences for solid fuel adiabatic in-situ combustion in porous media," Computational & Applied Mathematics, vol. 25, no. 1, pp. 27-54, 2006.

[7] G. Chapiro, A. A. Mailybaev, A. J. de Souza, D. Marchesin, and J. Bruining, "Asymptotic approximation of long-time solution for low-temperature filtration combustion," Computational Geosciences, vol. 16, no. 3, pp. 799-808, 2012.

[8] Y. B. Zeldovich, G. I. Barenblatt, V. B. Librovich, and G. M. Makhviladze, The Mathematical Theory of Combustion and Explosions, Springer, New York, NY, USA, 1985.

[9] A. A. Mailybaev, J. Bruining, and D. Marchesin, "Recovery of light oil by medium temperature oxidation," Transport in Porous Media, vol. 97, no. 3, pp. 317-343, 2013.

[10] C. Lu and Y. C. Yortsos, "Dynamics of forward filtration combustion at the pore-network level," AIChE Journal, vol. 51, no. 4, pp. 1279-1296, 2005.

[11] S. Lovett, F. Monmont, and N. Nikiforakis, "An experimentally-based in-situ combustion model with adaptive meshing," Combustion and Flame, vol. 162, no. 4, pp. 960-977, 2015.

[12] Z. Zhu, M. Bazargan, A. Lapene, M. G. Gerritsen, L. M. Castanier, and A. R. Kovscek, "Upscaling for field-scale in-situ combustion simulation," in Proceedings of the Western North American Region Meeting, SPE, Anchorage, Ala, USA, May 2011.

[13] O. Zik and E. Moses, "Fingering instability in combustion: an extended view," Physical Review E, vol. 60, no. 1, pp. 518-531, 1999.

[14] E. R. Ijioma, A. Muntean, and T. Ogawa, "Effect of material anisotropy on the fingering instability in reverse smoldering combustion," International Journal of Heat and Mass Transfer, vol. 81, pp. 924-938, 2015.

[15] Z. Chen, G. Huan, and Y. Ma, Computational Methods for Multiphase Flows in Porous Media, SIAM, Philadelphia, Pa, USA, 2006.

[16] L. D. Landau and E. M. Lifshitz, Statistical Physics, Elsevier, Oxford, UK, 1980.

[17] L. P Dake, Fundamentals of Reservoir Engineering, Elsevier, Oxford, UK, 1983.

[18] K. M. Hess, S. H. Wolf, and M. A. Celia, "Large-scale natural gradient tracer test in sand and gravel, cape cod, massachusetts: 3. Hydraulic conductivity variability and calculated macrodispersivities," Water Resources Research, vol. 28, no. 8, pp. 2011-2027, 1992.

[19] E. Abreu, J. Douglas Jr., F. Furtado, D. Marchesin, and F. Pereira, "Three-phase immiscible displacement in heterogeneous petroleum reservoirs," Mathematics and Computers in Simulation, vol. 73, no. 1-4, pp. 2-20, 2006.

[20] M. R. Correa and M. R. Borges, "A semi-discrete central scheme for scalar hyperbolic conservation laws with heterogeneous storage coefficient and its application to porous media ow," International Journal for Numerical Methods in Engineering, vol. 73, no. 3, pp. 205-224, 2013.

[21] Crane, "Flow of uids: through valves, fittings and pipe," Tech. Rep., Crane Co., 1982.

[22] R. J. LeVeque, Numerical Methods for Conservation Laws, Birkhauser, Basel, Switzerland, 1990.

Weslley da Silva Pereira (1) and Grigori Chapiro (2)

(1) National Laboratory for Scientific Computing (LNCC), Av. Getulio Vargas, 333, 25651-070 Petropolis, RJ, Brazil

(2) Department of Mathematics, Federal University of Juiz de Fora (UFJF), Rua Jose Lourenco Kelmer, s/n, Campus Universitario, 36036-900 Juiz de Fora, MG, Brazil

Correspondence should be addressed to Grigori Chapiro; grigori@ice.ufjf.br

Received 2 July 2016; Revised 16 September 2016; Accepted 27 October 2016; Published 12 January 2017

Caption: FIGURE 1: Standard normal distribution Z correlated by formula (43) to the log-normal permeability field in a 200 x 50 grid.

Caption: FIGURE 2: Histogram of the log-normal permeability field [kappa] and respective probability density function [f.sub.[kappa]]. The "Freq." function shows the frequency of elements that has permeability value occurring in each range of permeability.

Caption: FIGURE 3: Simulation results for the homogeneous permeability field at 10 days. Dimensionless variables values for [THETA], Y, and [[rho].sub.f] are on (a) and for [rho], p, and [u.sub.x] are on (b). The figure corresponds to the cross-section y = 0. The analytically estimated parameter [[THETA].sup.b] is plotted on (a). The space variable x corresponds to the horizontal axis.

Caption: FIGURE 4: Simulation results for the homogeneous permeability field at 50 days. Dimensionless variable values for [THETA], Y, and [[rho].sub.f] are on (a) and for [rho], p, and [u.sub.x] are on (b). The figure corresponds to the cross-section y = 0. The analytically estimated parameter [[THETA].sup.b] is plotted on (a). The space variable x corresponds to the horizontal axis.

Caption: FIGURE 5: Combustion wave temperature for the homogeneous porous medium: analytically estimated [[THETA].sup.b] and numerically obtained [[THETA].sup.b.sub.S]; see Section 5.1.

Caption: FIGURE 6: Analytically estimated combustion wave speed V (see Section 5.2) and numerically obtained combustion wave speed [V.sup.S] corresponding to the homogeneous porous medium.

Caption: FIGURE 7: Analytically estimated gas velocity [u.sup.u] (see Section 5.2) and numerically obtained gas velocity [u.sup.u.sub.S] corresponding to the homogeneous porous medium.

Caption: FIGURE 8: Simulation results in a log-normal permeability field for the dimensionless variables from top to bottom: [THETA], Y, [[rho].sub.f], [rho], and p. The solution is shown at 10 days (a) and 50 days (b). Space variables x and y are plotted in meters.

Caption: FIGURE 9: Simulation results in a log-normal permeability field for the dimensionless variables from top to bottom: u = [absolute value of (u)], [u.sub.x], and [absolute value of ([u.sub.y])]. The solution is shown at 10 days (a) and 50 days (b). Space variables x and y are plotted in meters.

Caption: FIGURE 10: Analytically estimated combustion wave temperature [[THETA].sup.b], numerically obtained temperature [[THETA].sup.b.sub.S], and its standard deviation [sigma] ([[THETA].sup.b.sub.S]) in the nonhomogeneous porous medium.

Caption: FIGURE 11: Analytically estimated combustion wave speed V (see Section 5.2) and numerically obtained combustion wave speed [V.sup.S] corresponding to the nonhomogeneous porous medium.

Caption: FIGURE 12: Analytically estimated gas velocity [u.sup.u] (see Section 5.2) and numerically obtained gas velocity [u.sup.u.sub.S] corresponding to the nonhomogeneous porous medium.

Caption: FIGURE 13: Regions separated by immobile fuel wave, thermal wave, combustion wave, and gas composition wave in the Riemann solution. Values of U = ([??], Y, [[rho].sub.f], u) in each region.

Caption: FIGURE 14: Numerical simulation of system (18)-(23) at times t = 2.6 x [10.sup.6] s (a) and t = 5.8 x [10.sup.6] s (b). The value of the parameter [[THETA].sup.b] from (33) is also plotted. Horizontal axis, x [m]; vertical axis, dimensionless quantities [THETA], Y, [[rho].sub.f] and u.

Caption: FIGURE 15: Perturbed initial conditions [mathematical expression not reproducible].

Caption: FIGURE 16: Simulation results for the dimensionless variables [THETA], Y, [p.sub.f], [u.sub.x], and [absolute value of ([u.sub.y])]. The solution is shown at 10 days (a) and 50 days (b). Variables x and y are plotted in meters.

Caption: FIGURE 17: Simulation results for the dimensionless variables [THETA], Y, [[rho].sub.f], [u.sub.x], and [absolute value of ([u.sub.y])]. The solution is shown at 10 days (a) and 50 days (b). Space variables x and y are plotted in meters.

Caption: FIGURE 18: Simulation results for the dimensionless variables [THETA], Y, [[rho].sub.f], [u.sub.x], and [absolute value of ([u.sub.y])]. The solution is shown at 10 days (a) and 50 days (b). Space variables x and y are plotted in meters.
```
TABLE 1: Dimensional parameters for ISC and their
typical values.

Symbol                               Description

R                                Ideal gas constant
[C.sub.m]                     Heat capacity of the rock
[c.sub.g]                 Specific heat capacity of the gas
[k.sub.T]                Thermal conductivity of the medium
Q                          Heat released due to combustion
[phi]                          Porosity of the medium
[d.sub.m]                  Molecular diffusion coefficient
[k.sub.p]                       Preexponential factor
E                                 Activation energy
[p.sub.atm]                     Atmospheric pressure
[T.sub.res]             Initial temperature of the reservoir
[Y.sub.inj]              Molar fraction of oxygen in the air
[[rho].sup.res.sub.f]    Initial molar concentration of fuel
[p.sub.res]                Reference pressure in reservoir
[bar.[kappa]]            Average permeability of the medium
[bar.[mu]]                      Average gas viscosity

Symbol                       Value              Unit

R                            8.314           J/(mol x K)
[C.sub.m]               2 x [10.sup.6]    J/([m.sup.3] x K)
[c.sub.g]                    27.42           J/(mol x K)
[k.sub.T]                    0.87           J/(m x s x K)
Q                       4 x [10.sup.5]          J/mol
[phi]                         0.3                 *
[d.sub.m]               2 x [10.sup.-6]      [m.sup.2]/s
[k.sub.p]                     500                1/s
E                            58000              J/mol
[p.sub.atm]                    1                 atm
[T.sub.res]                   370                 K
[Y.sub.inj]                  0.21                 *
[[rho].sup.res.sub.f]         372           mol/[m.sup.3]
[p.sub.res]                    1                 atm
[bar.[kappa]]            [10.sup.-12]         [m.sup.2]
[bar.[mu]]              2 x [10.sup.-5]        Pa x s

TABLE 2: The relative error (its mean value and standard deviation
(std)) for all estimates for homogeneous and nonhomogeneous porous
medium.

Relative error         [[THETA].sup.b]   V (B1)   V (B2)
of estimate

Homogeneous case:
Mean                      0.05%        4.04%    3.08%
std                       0.07%        0.73%    0.72%
Nonhomogeneous case:
Mean                      4.99%        3.91%    3.03%
std                       1.11%        1.85%    1.72%

Relative error         V (B3)   [u.sup.u] (B1)
of estimate

Homogeneous case:
Mean                 3.35%        2.00%
std                  0.76%        0.45%
Nonhomogeneous case:
Mean                 3.21%        1.70%
std                  1.78%        0.42%

TABLE 3: Comparing analytical estimates for parameter values of
[[THETA].sub.b], V, and [u.sup.u] with numerically obtained values
for times [t.sub.t] = 2.6 x [10.sup.6] s and
[t.sub.2] = 5.8 x [10.sup.6] s.

Parameter                      Numerical at [t.sub.1]

Analytical   Value    Relative error

[[THETA].sup.b]     1.0249     1.0230       0.19%
V                   0.9954     0.9959       0.05%
[u.sup.u]           1.0009     1.0000       0.09%

Parameter         Numerical at [t.sub.2]

Value    Relative error

[[THETA].sup.b]   1.0248       0.01%
V                 0.9950       0.04%
[u.sup.u]         1.0000       0.09%
```