Computational Geosciences 4 (2000) 141–164

141

Node-centered finite volume discretizations for the numerical simulation of multiphase flow in heterogeneous porous media ∗ Ralf Huber ∗∗ and Rainer Helmig Institut f¨ur ComputerAnwendungen im Bauingenieurwesen, TU Braunschweig, Pockelsstr. 3, D-38106 Braunschweig, Germany E-mail: {r.huber;r.helmig}@tu-bs.de

Received 23 February 2000

Three node-centered finite volume discretizations for multiphase porous media flow are presented and compared. By combination of these methods two additional discretization methods are generated. The ability of these schemes to describe flows at textural interfaces of different geologic formations is investigated. It was found that models with nonzeroentry pressures for the capillary pressure–saturation relationship in conjunction with the Box discretization may give rise to spurious oscillations for flows around low permeable lenses. Furthermore, the applicability and sensitivity of the discretization methods with regard to the used computational grids is discussed. The schemes are used for the numerical study of two-phase flow in porous media with zones of different material properties. Keywords: finite volume method, heterogeneity, capillary effects AMS subject classification: 76S05, 76T99, 65M99

1.

Introduction

Subsurface contaminations by organic liquids, e.g., solvents, fuels and oils, which are widely used in the industrialized world, represent a danger to many drinking water resources. In the soil they typically form liquid phases which are referred to as nonaqueous phase liquids (NAPLs). By dissolution of particles into the ambient groundwater, NAPLs act as long-term contamination sources in the subsurface and can affect large volumes of potential drinking water. The relatively low aqueous solubilities of various organic chemicals still lie considerably above drinking water standards. NAPLs are distinguished upon their relative density to water into dense NAPLs (DNAPLs) and light NAPLs (LNAPLs). ∗

This research was supported by the Deutsche Forschungsgemeinschaft, SFB 404 Mehrfeldprobleme in der Kontinuumsmechanik. ∗∗ The author was supported by a scholarship from the foundation “Besinnung und Ordnung”.  Baltzer Science Publishers BV

142

R. Huber, R. Helmig / Finite volume discretizations

Numerical modeling can be used in order to gain a better understanding of the migration and fate of these organic contaminants in the subsurface. Two-phase and three-phase models can be used to describe the movement of NAPLs in the saturated and unsaturated ground zones. Numerical models find their application in prognoses of contaminated site conditions and risk assessment studies of industrial facilities which typically entail Monte Carlo simulations [26]. Furthermore, numerical simulators are an important key for the design of new in-situ remediation technologies, e.g., alcohol and surfactant flooding [21,27], steam injection [7], or soil vapor extraction and air sparging [22,29]. Typically, the subsurface is structured into geologic formations which differ in their texture, e.g., the permeability and grain size distribution. These zones are characterized by their porosities, intrinsic permeabilities and capillary–hydraulic properties. An accurate description of capillary and gravity forces is important since they strongly influence the flow behavior of a multiphase fluid system in heterogeneous environments. At transitions from coarse-grained to fine-grained zones the so-called “entry pressure effect” leads to a pooling of the nonwetting phase until a threshold capillary pressure is overcome. On the other hand, at transitions in the reverse direction, the coarse-grained regions may act as “capillary barriers” and may prevent a penetration of the wetting phase into these formations. Experimental and numerical studies of multiphase flow processes in heterogeneous porous media were conducted by several authors. Kueper et al. [18] investigated the migration of tetrachloroethylene (PCE) under saturated conditions in a sand-box which was structured by layers of different sands. Kueper and Frind [19] simulated this experimental study. The main focus in their study lay in the entry pressure effect of lower permeable zones. Since the vertical extents of the different formations were relatively small, this problem seems not very suitable for the investigation of grid orientation effects. Illangasekare et al. [17] conducted several experimental studies with DNAPL moving under saturated and unsaturated conditions through zones of higher and lower permeability. They pointed out that their results could be used for model validation purposes, but they did not, to our knowledge, perform numerical simulations of their experimental studies. Forsyth et al. [12] presented several test problems which entailed the infiltration of water from the top into initially dry soils with zones of different material properties. They investigated the computational costs associated with the use of a passive air phase formulation versus a full two-phase formulation and the use of a variable substitution algorithm. Furthermore, they compared upstream and central weighting of mobility terms. For the solution of their problems they employed the control-volume finite element discretization. Schroth et al. [28] studied the influence of an inclined coarse-grained formation on water and LNAPL flow. The textural interface acted as a capillary barrier to water flow. However, NAPL which had a lower wettability than water penetrated the coarse-grained zone in regions where the water saturations were sufficiently high.

R. Huber, R. Helmig / Finite volume discretizations

143

They also performed simulations of the experimental study by use of the integrated finite difference method in conjunction with a rectangular mesh which was aligned to the textural interface. Their predicted water saturations above the interface were in good agreement with measurements. Schroth et al. emphasized the importance of an accurate description of phase saturations at heterogeneity interfaces on fluid flow. One main issue in modeling multiphase flow in heterogeneous environments is the accurate description of the location and extent of heterogeneity interfaces and their associated interface conditions [24,30,31]. Finite element methods, e.g., the Petrov–Galerkin methods, have proven to be not suitable for heterogeneous porous media problems with nonzero entry pressures [14]. In addition they have the drawback that they do not conserve mass locally. For advection-dominated flow problems the mixed and mixed-hybrid finite element methods are popular. Mixed finite elements have the advantage that the velocity and pressure fields have the same order of accuracy [4] and they yield physically realistic flow fields in highly heterogeneous media [6]. However, these methods are in general used in conjunction with the IMPES approach, which implies a sequential solution of an implicit pressure equation and one or several explicit saturation equations. They are not so favorable for capillary pressure-dominated flow problems [15]. Another class of discretization schemes represent the node-centered finite volume methods. They can be easily applied to fully coupled pressure–saturation formulations which have several advantages over other approaches. For instance, capillary pressure and other nonlinear terms can be handled with a Newton–Raphson algorithm. This gives a stable and robust numerical scheme [16]. In this paper several node-centered finite volume schemes are presented and compared with respect to their ability to predict flows along and around zones of different material properties. This is done by using two experimental studies of twophase porous media flow which entail heterogeneities of relatively simple geometries. 2.

Governing equations

Multiphase flow in porous media is described by the mass conservation of each phase α, i.e., ∂(φ%α Sα ) (1) = −∇ · [%α vα ] + qα , ∂t where Sα is the saturation, %α the density, vα the Darcy velocity, and qα the mass source rate of phase α. φ denotes the porosity. The Darcy velocity of phase α is given by vα = −λα K(∇pα − %α g),

(2)

where λα = krα /µα is the mobility of phase α, i.e., the ratio of the relative permeability kr to the dynamic viscosity µ. K denotes the intrinsic permeability tensor, pα is the α-phase pressure, and g is the gravitational acceleration vector. Equation (2) represents

144

R. Huber, R. Helmig / Finite volume discretizations

a generalized form of the Darcy law. The capillary pressure pc is defined as the difference in fluid phase pressures, i.e., pc = pn − pw , where subscripts n and w denote the nonwetting and wetting phase, respectively. By introduction of the concept of relative permeability the interaction of simultaneous flow and obstruction within the porous medium between the individual fluid phases is taken into account. For the description of the capillary–hydraulic properties the most widely used models are the ones of Brooks and Corey (BC) [2] and Van Genuchten (VG) [32]. The BC capillary pressure–saturation relationship has the following representation:
−1/λ , (3) pc = Pd S w where Pd is the entry pressure, i.e., the capillary pressure necessary for the nonwetting phase to displace the wetting phase from a porous medium. λ is a pore size distribution index, and Sw =

Sw − Swr 1 − Swr

(4)

is the effective water saturation. Sw is the actual water saturation whereas Swr is the irreducible water saturation. The Van Genuchten model entails the following functional dependence of capillary pressure on saturation:
−1/m 1/n 1 Sw −1 , α 1 m=1− , n pc =

(5) (6)

where m, n and α are model parameters. The pc –Sw -curves for the first test problem are shown in figure 1. The main difference of the two concepts lies in the distinct entry pressure of the Brooks–Corey model which is missing in Van Genuchten’s model.

Figure 1. NAPL–water capillary pressure curves for the fine sand and coarse sand.

R. Huber, R. Helmig / Finite volume discretizations

145

The relative permeabilities of the Brooks–Corey model have the following form:
(2+3λ)/λ , (7) krw = S w
2
(2+λ)/λ
krn = 1 − S w 1 − S w , (8) whereas the Van Genuchten parameterizations have the following representation:
1/2 
1/m
m 2 krw = S w , (9) 1 − 1 − Sw
1/2 
1/m 2m . (10) krn = 1 − S w 1 − Sw Capillary pressure and the relative permeabilities are strongly nonlinear functions of the saturation. Substitution of the Darcy velocity from (2) into the continuity equation (1) yields a coupled pressure–saturation formulation:   ∂(φ%α Sα ) (11) = ∇ · %α λα K(∇pα − %α g) + qα , ∀α. ∂t The capillary pressure term in the model equations acts as a diffusion term [13]. Equations (11) which represent coupled strongly nonlinear partial differential equations with changing mathematical properties dependent on the dominating flow mechanisms are the basis of many two-phase and three-phase formulations [13]. In the following section, three finite volume discretizations and two combinations are applied to them. 3.

Finite volume discretization

The computational domain is discretized by a finite element grid, the so-called “primary mesh”, with triangles and quadrilaterals in the two-dimensional case, and tetrahedra and hexahedra in three dimensions. A “secondary mesh” which describes the finite volumes is deduced from this finite element mesh. This is done by connecting the centers of gravity of each element with the associated edge midpoints (see figure 2). The extension to three dimensions is straightforward. Thus, each finite volume is associated with a node. Finite volumes are also referred to as control-volumes, cells or boxes. Porous media properties which vary spatially in heterogeneous environments are specified for control-volumes, i.e., they are assigned to nodes. It should be noted that there also exist approaches which identify porous media properties with finite elements, i.e., control-volumes may include regions of different properties. From the method of weighted residuals by use of mass lumping of the accumulation terms and an implicit Euler time discretization the following discrete equations in the general form of a finite volume scheme are deduced from (11):  X n+1 Vi  n+1 − (φ%α Sα )ni = Fα,ij + qα,i Vi ∀i, α, (12) (φ%α Sα )n+1 i ∆t j∈ηi

where Vi is the volume of the ith control-volume or box, ηi the set of neighbor nodes of i, ∆t the time-step size, and superscripts indicate the time-step level. The discrete

146

R. Huber, R. Helmig / Finite volume discretizations

Figure 2. Spatial discretization of a finite volume/box.

. n+1 flux Fα,ij of phase α between finite volumes i and j during the (n + 1)th time-step is given by

Z n+1 Fα,ij

=

Ω

  X n+1 n+1 n+1 Wi ∇ · %α λα K · ψα,j ∇Nj dΩ.

(13)

j∈ηi

Nj denote the shape functions which are the usual Lagrange polynomial C 0 base functions. Wi are test functions which are specified later. ψα,j represents the total fluid potential of phase α at node j and is defined by ψα,j = pα,j + %α gzj .

(14)

g in (14) is the gravity constant and zj denotes the (vertical) z-coordinate of node j. Based on equations (12), in the following three node-centered finite volume discretizations are developed. The Brooks–Corey model employs in the absence of the nonwetting phase instead of definition (14) a fluid phase potential of ψn = pw + Pd + %n gz,

(15)

where Pd is the entry pressure of the associated porous medium. This definition guarantees that a penetration into a zone is only triggered if the associated entry pressure is overcome.

R. Huber, R. Helmig / Finite volume discretizations

147

3.1. Box method (BOX) The Box method uses in equation (13) piecewise constant test functions Wi = χBi , i.e., characteristic functions of the boxes Bi , which yields the following surface integral: Z X n+1 n+1 %n+1 ψα,k ∇Nk · ~n∂Bi ds, (16) α λα K · ∂Bi

k∈ηi

where ~n∂Bi is the outward unit normal vector with respect to box boundary ∂Bi . The surface integral in equation (16) of the Box discretization comprises the flows of phase α over the box boundary ∂Bi from the neighbor boxes into box Bi . The portion which is associated with the box boundary segments shared by boxes Bi and Bj is Z X n+1 n+1 n+1 Fα,ij = %n+1 ψα,k ∇Nk · ~n∂Bi ds. (17) α λα Kij · ∂Bi ∩∂Bj

k∈ηi

In (17) the harmonic mean Kij of the nodal intrinsic permeabilities Ki and Kj is used since this choice yields a consistent approximation of the flux. BOX uses the following approximation of (17): X XZ
BOX,n+1 n+1 n+1 Fα,ij = (%α λα )ups(i,j) ds Kij · ψα,k ∇Nk · ~n∂Bi e,ij , (18) e∈Th

∂Bi ∩∂Bj ∩e

k∈e

where Th is the set of elements, k ∈ e denotes nodes k of element e, and (·)e,ij designates the value at the midpoint of the box boundary segment ∂Bi ∩ ∂Bj inside of element e. ups(i, j) in (18) denotes the node in the physical upstream direction. Upstream weighting of the nonlinear terms of phase mobilities and densities leads to a stable scheme [3], whereas central weighting gives a higher-order approximation which may give rise to spurious oscillations in the solution. Hence, the Box discretization evaluates the gradients of the fluid potential at the midpoints of the linear box boundary segments. 3.2. Control-volume finite element method (CVFE) Forsyth [9] introduced the so-called control-volume finite element method (CVFE) to the modeling of porous media flows. This discretization scheme is derived now. By application of Green’s theorem to the integral in (13) with Wi = Ni , the following volume integral is obtained: Z X n+1 n+1 ψα,j ∇Nj dΩ. (19) − %n+1 α λα ∇Ni · K · Ω

j∈ηi

148

R. Huber, R. Helmig / Finite volume discretizations

With

X

ψα,j ∇Nj =

j∈ηi

X

(ψα,j − ψα,i )∇Nj ,

j∈ηi

(19) can be rewritten as XZ
n+1 n+1 n+1 − %n+1 α λα ∇Ni · K · ∇Nj dΩ ψα,j − ψα,i . j∈ηi

(20)

Ω

(21)

Now, an influence coefficient technique [5,9] is applied to the integral in (21) which entails the following approximation: Z CVFE %α λα ∇Ni · K · ∇Nj dΩ ≈ (%α λα )ups(i,j) γij , (22) Ω

where the influence coefficients, here %α and λα , are upstream weighted according to  i, if γij (ψα,j − ψα,i ) 6 0, ups(i, j) = (23) j, if γij (ψα,j − ψα,i ) > 0 CVFE of the CVFE discretization is and the transmissivity integral/transmissibility γij given by Z CVFE = − ∇Ni · Kij · ∇Nj dΩ. (24) γij Ω

The integral in (24) is evaluated by use of a one-point quadrature rule. The intrinsic permeability tensor Kij in the transmissivity integrals is the harmonic mean of the two associated nodal permeabilities. Note that the upstream direction defined by (23) is dependent on the sign of the transmissivity integral and does not necessarily coincide with the physical upstream direction. Equation (23) ensures that the saturations remain in the physical range [9]. Letniowski and Forsyth [20] investigated the CVFE method and noticed that obtuse triangles which rectangles with a side √ do not obey the Delauney criterium and CVFE length ratio greater 2 lead to negative transmissibilities γij . These give rise to non-physical discrete fluxes and cause an oscillatory behavior in the Newton iterations which usually leads to a failure of the algorithm. Letniowski and Forsyth also suggested avoiding negative γij when discretizing a given model domain and called the associated constraint the “positive transmissibility condition”. In order to alleviate the problem with negative transmissibilities, downstream influence coefficients according to (23) are used. Hence, CVFE is applicable to a wider range of problems. However, there exist many meshes for which this modification is not sufficient. The discrete flux of the CVFE method between finite volumes i and j can be expressed as
CVFE,n+1 n+1 n+1 CVFE = (%α λα )n+1 − ψα,i ψα,j . (25) Fα,ij ups(i,j) γij

R. Huber, R. Helmig / Finite volume discretizations

149

The discrete fluxes are associated with neighboring node pairs, i.e., pairs of nodes or control-volumes (i, j) for which γij 6= 0. Equation (12) with (25) describes a mass balance for a not closer specified control-volume associated with node i. The discrete flux terms (25) are structured in the following way: discrete flux = [influence coefficients] × [transmissivity integral] × [nodal difference]. The transmissivity integral comprises the geometrical information, i.e., the interfacial area divided by the nodal distance, whereas the nodal difference is related to the driving force, here, the gradient of the fluid potential. The discretization of the Neumann boundary condition along the model boundary deviates from the original control-volume finite element method described by Forsyth [9]. The representation of the boundary fluxes from the original method yields a weak form in analogy to finite element methods. Here, flux boundary conditions are specified for each boundary nodal volume and, hence, the same representation as for the BOX method is employed. 3.3. IFDM method In general, the discrete fluxes of the BOX method cannot be expressed in the form of equation (25). In order to obtain a discretization which uses surface integrals as BOX and at the same time fits into the form of (25), the flux across a box boundary segment is approximated by use of the difference in total fluid potentials of the two associated nodes. This is analogous to integrated finite difference approximations. The resulting discrete fluxes of this scheme, which will be denoted here as the IFDM method, are given by Z
IFDM,n+1 n+1 n+1 n+1 Fα,ij = (%α λα )ups(i,j) Kij · (∇Nj )ij · ~n∂Bi ds ψα,j − ψα,j , (26) ∂Bi ∩∂Bj

where ups(i, j) is the node in the physical upstream direction and (·)ij means the evaluation at the midpoint of edge i–j. These discrete fluxes represent for homogeneous media coarser approximations of the surface fluxes than the ones used by the BOX method. The associated transmissibilities have the following form: Z IFDM γij = Kij · (∇Nj )ij · ~n∂Bi ds. (27) ∂Bi ∩∂Bj

IFDM > 0 ∀(i, j), formula (23) also applies to IFDM. The presented scheme Since γij coincides with the integrated finite difference method [23] in the case of rectangular meshes. For other grids, e.g., triangular meshes, the finite volume interfaces are in general not perpendicular to the connecting line between nodes which gives rise to grid orientation effects.

150

R. Huber, R. Helmig / Finite volume discretizations

3.4. Global algorithm A residual-based Newton–Raphson iteration method in conjunction with the BiCGStab matrix solver is used for the solution of the discretized phase conservation equations in order to resolve the strong nonlinearities in the constitutive relationships of capillary pressure and relative permeability. Two-phase flow is solved with the discrete form of the two-phase flow formulation (12) by using a phase pressure and a phase saturation as primary variables. In the following computations pw and Sn are used. Secondary variables, e.g., the constitutive relations and the remaining saturations, are derived from these. 3.5. Properties of the methods BOX and IFDM always yield local mass conservative schemes. If K is a symmetCVFE is symmetric, and as a consequence the associated discretization ric tensor, then γij is also local mass conservative. In the one-dimensional case CVFE, IFDM and BOX coincide. The finite volumes associated with the CVFE method are, in the multidimensional case, not uniquely determined. For a rectangular mesh the BOX and IFDM methods represent five-point schemes in contrast to the control-volume finite element method CVFE which is associated with a nine-point stencil. Panday et al. [25] and Huber and Helmig [15] have already investigated the grid sensitivity of five- and nine-point schemes by use of the five-spot problem. The five-point solutions exhibited for coarse grids a diamond-shape flood whereas the nine-point scheme predicted circle-shape fronts which were considerably closer to the exact solution. The IFDM discretization has no interaction between diagonal nodes of quadrilateral elements, since their associated boxes have no mutual interfacial surface, and the fluxes between two finite volumes are independent of other nodal values. BOX does not represent a pure five-point scheme because the used gradients of the fluid potentials are interpolated inside of each element and thus all element nodes influence the discrete flux between two finite volumes. This should yield a better approximation of the fluxes in homogeneous media than does IFDM. However, for heterogeneous media an entire decoupling of the discrete fluxes from other nodal values should give physically more realistic results. 3.6. Combination of the methods The three discretization schemes, BOX, CVFE, and IFDM, can be arbitrarily combined. Two combinations are investigated: • The first method, denoted as BOXIFDM, is a combination of the BOX and IFDM schemes. If an element is situated at a textural interface of two or more different porous media, the IFDM scheme is used. Otherwise the BOX scheme is employed for this element. This new discretization method has the advantage that it is less

R. Huber, R. Helmig / Finite volume discretizations

151

sensitive to the grid orientation in homogeneous regions and utilizes a decoupling of element fluxes near heterogeneity interfaces. • The second method, denoted as CVBOXIFDM, combines all three basic discretization schemes. If an element has no negative transmissibilities, CVFE is used. Otherwise, BOX is used in homogeneous regions, whereas IFDM is used in heterogeneous regions. 4.

Test problem I

In the following, an experimental study conducted by Barczewski and Koschitzky [1] is used as a test problem in order to assess the performance of the five finite volume methods – BOX, CVFE, IFDM, BOXIFDM and CVBOXIFDM – and their ability in describing multiphase flows in heterogeneous media. The numerical results are examined with regard to predicted NAPL infiltration, accumulation on top of the fine sand lens and vertical propagation. 4.1. Model assumptions and numerical simulation Trichloroethylene (TCE) is infiltrated into a small container of length 90 cm and height 65 cm which is filled with sand. In the center of this box a fine sand lens is located. The infiltration of TCE occurs over a length of 10 cm at the center of the upper surface. The remaining upper boundary is impermeable. Along the two sides a hydrostatic pressure is prescribed. For symmetry reasons, it is sufficient to model only one half of a cross-section. In figure 3 the model domain is depicted. The fine sand zone is defined by x > 0.1875 m; 0.325 m < z < 0.4625 m. The sand properties are given in table 1. Along the left boundary hydrostatic pressure conditions are prescribed. For the BC model, a permanent TCE (%n = 1462 kg/m3 , µn = 5.7 · 10−4 Pa s) saturation of 0.25 is assumed along the infiltration zone. Four computational grids are employed. Mesh 1 is a rectangular grid consisting of 70 nodes and 6 × 9 elements. A uniform refinement of mesh 1 by subdivision of each rectangular element into 3 × 3 equal rectangular elements yields mesh 2 with Table 1 Material properties of test problem I. Parameter Brooks–Corey λ Brooks–Corey Pd Van Genuchten n Van Genuchten α Residual water saturation Swr Intrinsic permeability K Porosity φ

Unit

Coarse sand

Fine sand

[–] [Pa] [–] [Pa−1 ] [–] [m2 ] [–]

2.7 755 5.755 0.001 0.09 6.64 · 10−11 0.40

2.0 2060 4.37 0.00035 0.12 7.15 · 10−12 0.39

152

R. Huber, R. Helmig / Finite volume discretizations

Figure 3. Configuration of test problem I.

Figure 4. Computational grids 2 and 2lref.

532 nodes, and a subsequent application of this refinement strategy results in mesh 3 which employs 4510 nodes. Mesh 2lref entails grid refinements near the textural interface of the two different sands and uses 1000 nodes. This locally refined mesh provides a better resolution of the interfaces of the fine sand formation. Meshes 2 and 2lref are depicted in figure 4. For all four meshes the heterogeneity interfaces are situated along the center axes of rectangular elements and the prescribed infiltration zone is exactly resolved by the boundary nodal volumes.

R. Huber, R. Helmig / Finite volume discretizations

153

Figure 5. Predicted DNAPL saturations after 100 minutes of infiltration (BC model). From left to right: BOX, CVBOXIFDM, IFDM, BOXIFDM for meshes 2, 2lref and 3 (from top).

All five discretization methods are applied to meshes 1, 2, 2lref and 3 in conjunction with the Brooks–Corey model. 100 minutes of TCE infiltration were simulated. A maximum time-step length of 50 s was allowed. The numerical results of the methods for meshes 2, 2lref and 3 are depicted in figure 5. The plotted lines represent the isocurves associated with Sn = 0.05, 0.1, 0.15, . . .. For meshes 1, 2, and 3, the results of CVFE and CVBOXIFDM coincide. CVFE fails for mesh 2lref. 4.2. Applicability of the methods All meshes include negative transmissibilities which occur inside of the low permeable region. The CVFE discretization could be successfully applied to meshes 1, 2 and 3 because only single-phase flow occurs inside of the fine sand lens. The locally refined mesh 2lref has in addition negative transmissibilities in the elements

154

R. Huber, R. Helmig / Finite volume discretizations

of the refinement stripes. Hence, the Newton algorithm for the CVFE method in conjunction with mesh 2lref failed after 1796.7 s of simulation. This was at the time when the NAPL front reached the first elements with negative transmissibilities. All other methods worked for the four different meshes. When anisotropic permeabilities are assumed, or variable mesh sizes and local grid refinements are employed, in general negative transmissibilities are induced. This usually leads to a failure of the Newton algorithm if CVFE is used. The generation of a grid which obeys the positive transmissibility condition and at the same time accurately describes the soil heterogeneities allows the use of CVFE. However, in many three-dimensional problems this is not possible. 4.3. Pooling on top of the lens During the experiment the infiltrated NAPL migrated from the source area downwards until it reached the textural interface of the low permeable formation. Experimental observations indicate that the NAPL did not penetrate the lens. It accumulated and propagated laterally along the interface. Hence, the maximum NAPL saturations Sn,max must occur at the textural interface below the infiltration zone. In table 2 the values of Sn,max are listed which the five discretization schemes in conjunction with the BC model and the used computational grids after an infiltration of 100 minutes predict. The numerical results indicate that the discretization schemes tend to underpredict Sn,max . The predicted values of the schemes increase with finer grids. This can be explained from a better spatial discretization of the pooling area of NAPL with smaller mesh sizes. The CVFE and CVBOXIFDM results indicate already for coarser meshes good estimates of Sn,max . An accurate prediction of this value is crucial for the description of the flow behavior along such textural interfaces, since an exceedance of the threshold value Sn,d would yield a penetration of NAPL into the fine sand formation. For the Brooks–Corey model, the NAPL saturation Sn,d associated with an exceedance of the entry pressure Pd,fine can be calculated for the given heterogeneity interface to have a value of 0.85 (see figure 1). This value lies considerably above the predictions of the numerical simulations. As Helmig and Huber [14] have already pointed out, a full-upstream weighting of mobility terms in conjunction with a finite difference approximation of the nodal fluid potential values is necessary to capture these physical processes correctly, especially for the Brooks–Corey model [2]. 4.4. Description of flows around corners The fluid flows are deduced in the case of all five discretization schemes by using the total fluid potentials ψw and ψn . The CVFE and IFDM discretizations evaluate the gradients of ψ by using pairs of nodal ψ-values, whereas the Box discretization employs an interpolation of all nodal ψ-values of an element. Hence, the ψ-value of

R. Huber, R. Helmig / Finite volume discretizations

155

Table 2 Characteristics of the numerical results after 100 minutes of infiltration (BC model). Method

Parameter

Unit

Mesh 1

Mesh 2

Mesh 2lref

Mesh 3

BOX

Sn,max D0.1 TCE mass Nonlinear iterations CPU time†

[–] [m] [kg] [–] [s]

0.489 0.320 5.100 272 22.8

0.581 0.443 5.075 488 793

0.605 0.451 5.076 444 2059

0.612 0.457 5.021 518 53673

CVFE

Sn,max D0.1 TCE mass Nonlinear iterations CPU time†

[–] [m] [kg] [–] [s]

0.572 0.289 5.045 261 22.1

0.608 0.372 5.104 503 794

– – – – –

0.620 0.428 5.038 510 52779

IFDM

Sn,max D0.1 TCE mass Nonlinear iterations CPU time†

[–] [m] [kg] [–] [s]

0.492 0.318 5.158 270 23.7

0.582 0.441 5.145 385 598

0.606 0.453 5.145 481 2064

0.612 0.461 5.071 522 56323

BOXIFDM

Sn,max D0.1 TCE mass Nonlinear iterations CPU time†

[–] [m] [kg] [–] [s]

0.491 0.311 5.103 275 25.8

0.581 0.439 5.076 550 935

0.605 0.449 5.076 462 2196

0.612 0.457 5.022 525 54982

CVBOXIFDM

Sn,max D0.1 TCE mass Nonlinear iterations CPU time†

[–] [m] [kg] [–] [s]

0.572 0.289 5.045 261 21.7

0.608 0.372 5.104 503 833

0.616 0.419 5.105 424 2029

0.620 0.428 5.038 510 53810

†

on a Pentium MMX 233 MHz.

the heterogeneity influences the flow around the associated corner in the case of the BOX method. Mesh 1 is a too coarse grid to capture the flows around the corners of the lens adequately. The five schemes yielded NAPL profiles which extended slightly farther than the corner of the lens. However, all five discretization schemes generated monotone saturation profiles. The definition of (15) for the fluid potential which is employed by the Brooks– Corey model leads to large nodal differences in ψn , and hence, especially for fine grids, to high gradients in total fluid potentials along heterogeneity interfaces. As a consequence, BOX generated for meshes 2, 2lref and 3 spurious oscillations in the vicinity of the corners of the lens. This is due to the influence of the fine sand corner nodes on the discrete fluxes around the associated corners of the fine sand lens. CVFE and IFDM employ a decoupling of the discrete fluxes from other nodal values and generated no unphysical over- and undershoots around the corners of the lens.

156

R. Huber, R. Helmig / Finite volume discretizations

Figure 6. Predicted DNAPL saturations after 100 minutes of infiltration (VG model with BOX and meshes 2, 2lref and 3).

The CVFE discretization evaluates fluxes between nodal volumes by using the nodal difference in fluid potentials. However, as pointed out before, the geometries of the nodal finite volumes and thus also of the associated interfaces are of a virtual nature. Hence, heterogeneities, and especially corners, cannot be exactly described by using the CVFE discretization. This is especially true for quadrilateral elements which are situated at corners of formations of different material properties. A flow around a corner typically implies that there is no or only little flow occurring diagonally across the associated quadrilateral element. However, the numerical results show that the diagonal flow across the corner dominates the other element fluxes. The effects of the diagonal flows of CVFE at corners are reduced with finer mesh sizes. The IFDM discretization accurately describes the flow around the corners of the fine-grained zone. There are neither diagonal fluxes nor an interaction of the corner nodal ψ-value with the flows around this corner. Since the BOXIFDM and IFDM methods coincide along the textural interface, they show a similar behavior. The BOX method which entailed spurious jumps in the saturations for meshes 2, 2lref and 3 in conjunction with the BC model is applied to the same meshes, however, now with the VG model. In order to have approximately the same infiltration rate of TCE, a NAPL saturation of 0.158 is assumed along the infiltration boundary for the VG model. The Van Genuchten model allows very small amounts of DNAPL to enter the fine-grained zone. This leads to a moderate increase in the associated fluid potentials which control the further flux into this region. As a result, the fluid potential values along the heterogeneity interfaces do not differ significantly. Hence, the deduced gradients of the BOX method at corners are not strongly perturbed. In figure 6 the predicted DNAPL saturation distributions after 100 minutes are shown. They yield monotone saturation profiles. 4.5. Prediction of DNAPL plume propagation The BOX, IFDM and BOXIFDM methods predict for coarser grids a slimmer plume which propagates much faster downwards than the results of CVFE and

R. Huber, R. Helmig / Finite volume discretizations

157

CVBOXIFDM indicate. The saturation profiles of the latter two methods exhibit a stronger horizontal spreading of the advancing DNAPL plume. As a measure for the quality of the numerical solution, D0.1 , the maximum vertical extent of the Sn = 0.1 isocurve, is selected. The D0.1 -values of the numerical results of the different meshes and methods in conjunction with the BC model are given in table 2. The D0.1 -values converge monotonously for all schemes with decreasing mesh sizes. The D0.1 -values of the IFDM, BOX, and BOXIFDM schemes give already for the coarse meshes good estimates of the convergence solution which must be close to 0.465. The VG model predicts higher relative permeabilities. Furthermore, the VG capillary pressure–saturation curves exhibit large gradients near Sn = 1, whereas the gradients of the corresponding BC-curves are diminishing when Sn = 1 is approached (see figure 1). As a consequence, the VG model induces higher mobilities of the fluid phases and stronger effects of “capillary diffusion”. Hence, the predicted saturation profiles of the BC and VG models differ considerably. Schroth et al. [28], who compared the two models of Brooks–Corey and Van Genuchten but used a flux boundary condition along the infiltration zone, noticed deviations of the predicted plume propagations. For their series of experimental studies they concluded that the BC model gave better estimates of the actual flow processes. The runtimes of the different schemes with the BC model did not vary significantly (see table 2). However, the CPU times for the BOX method are, especially for the finer mesh computations, higher if the Van Genuchten model (mesh 2: 445 Newton steps, 802 s; mesh 2lref: 513 Newton steps, 2457 s; mesh 3: 595 Newton steps, 62 323 s) is employed. This might be caused by the varying entry constraints along the textural interface. The fluid potentials are gradually increased in order to control the flow into the fine-grained zone. This may lead to a slower convergence in the Newton iterations. 4.6. Infiltration of NAPL The amount of infiltrated TCE after 100 minutes due to the prescribed saturation boundary condition which the numerical simulations predicted is also given in table 2. The fine grid solution suggests a total TCE mass of 5.063 kg. The predicted total TCE mass of the IFDM scheme converges from above to this value. The other schemes give similiar results, however, they do not show a monotone behavior for the presented meshes with finer grids. In the case of CVFE this is probably due to the incorrect representation of the flows around corners for the coarse grids. Since the Box method generates for the BC model unphysical extrema in the saturation distribution near the corners of the fine sand lens, its numerical solutions do not converge to the exact solution. Hence, the associated infiltrated NAPL mass is not of interest here.

158

5.

R. Huber, R. Helmig / Finite volume discretizations

Test problem II

Test problem II is targeted to give information on how the results of the five discretization schemes differ with regard to the description of NAPL propagation for complex configurations of heterogeneities. This problem represents a typical example where the size of a given problem forces the use of quite coarse approximations of layers of different material properties. In addition the associated experiment is aimed at generating situations where viscous, gravitational and capillary forces, which control the two-phase flow, were locally of different strengths. The test problem corresponds to the initial phase of a larger-scale essentially two-dimensional experiment of a TCE infiltration and redistribution in an initially saturated structured aquifer having dimensions 6.43 m×2.4 m×0.4 m. Barczewski and Koschitzky [1] described this flume experiment, which was conducted at the VEGAS facility (Institut f¨ur Wasserbau, Univ. Stuttgart). 5.1. Model assumptions and numerical simulation An 85-minute long period of the experiment is considered during which TCE was infiltrated at a constant rate of 2.04 l/min. The model domain represents a vertical cross-section of the aquifer. The entire lower and upper model boundaries except for the infiltration zone are modeled as impermeable. Along the left and right sides of the model domain a constant hydrostatic pressure is prescribed. The difference in water pressure between the two sides amounts to 662.75 Pa. The DNAPL source is modeled as a flux boundary with a constant flow-rate of TCE over the entire length (30 cm) of the infiltration zone. A set of fine sand lenses which are inclined at different angles against the induced hydraulic gradient is situated in the model domain. Three different sands are used. Their properties are listed in table 3. The flume is filled with a coarse sand. The lenses consist of a fine sand, except for the base of the U-shape trough which is made of a medium sand. Lens A has an inclination of 1%, lens B of 5% and lens C of 10%. Lens D is horizontal. The configuration of the fine sand lenses is depicted in figure 7. The lower permeable formations are specified by polygons with the following corner coordinates: lens A: (0.37,0.3), (1.83,0.3), (1.83,0.445), (0.37,0.43); lens B: (1.13,0.67), (2.63,0.67), (2.63,0.875), (1.13,0.8); lens C: (1.87,1.37), (4.67,1.37), (4.67,1.73), (4.2,1.73), (1.87,1.5); lens D: (2.63,1.9), (4.87,1.9), (4.87,2.03), (2.63,2.03); trough Table 3 Material properties of test problem II. Parameter

Unit

Coarse sand

Medium sand

Fine sand

Brooks–Corey λ Entry pressure Pd Irreducible water saturation Swr Intrinsic permeability K Porosity φ

[–] [Pa] [–] [m2 ] [–]

2.0 200 0.05 4.60 · 10−10 0.39

2.3 700 0.15 3.1 · 10−11 0.40

3.5 1800 0.18 9.05 · 10−12 0.43

R. Huber, R. Helmig / Finite volume discretizations

159

Figure 7. Coarse computational grid and problem configuration.

E: (3.4,0.47), (5.6,0.47), (5.63,1.13), (5.43,1.13), (5.43,0.67), (3.6,0.67), (3.6,0.97), (3.4,0.97). The coarse computational grid entails 1156 nodes (see figure 7). The thin fine sand lenses are modeled with this grid only as single layers of finite volumes. The mesh contains no elements with negative transmissibilities. A refined mesh with 9955 nodes is generated similarly to mesh 2lref in test problem I, i.e., each coarse element is subdivided into nine elements. Several negative transmissibilities are introduced by using this refinement strategy. The numerical simulations with the coarse and the refined mesh are conducted with a maximum allowable time-step of 60 s. 5.2. Results The NAPL saturation profiles after a period of 85 minutes predicted by the five schemes in conjunction with the coarse mesh are depicted in figure 8. The lines are associated with the isolines Sn = 0.05, 0.1, 0.2, 0.5. All five discretization methods yield monotone solutions for the coarse mesh. The results of CVFE and CVBOXIFDM are identical for this mesh. All five methods predict that the DNAPL did not penetrate lenses A–D, but was able to overcome the entry pressure of the medium sand which made up the base of trough E. The BOX, IFDM and BOXIFDM results suggest that DNAPL after 85 minutes of infiltration had advanced to the bottom of the flume, whereas the other two methods indicate that the DNAPL was at the stage of passing the base of the trough. The shape of the DNAPL front to the right of lens D is very slim for the BOX, IFDM, and BOXIFDM methods. The result of the IFDM method for the coarse mesh apparently exhibits some grid orientation effects in the region above the trough. The CVFE and CVBOXIFDM predictions show a horizontal widening of the DNAPL fronts, such that the top of the right side wall of the trough is hit. Furthermore, diagonal flows across the corners of the lenses are apparent.

160

R. Huber, R. Helmig / Finite volume discretizations

Figure 8. Coarse mesh. Predicted DNAPL saturations after 85 minutes of infiltration (coarse mesh). From left: top row: BOX, CVFE/CVBOXIFDM; below: IFDM, BOXIFDM.

Figure 9. Refined mesh. Predicted DNAPL saturations after 85 minutes of infiltration (refined mesh). From left: top row: BOX, CVBOXIFDM; below: IFDM, BOXIFDM.

The numerical results of the refined mesh are depicted in figure 9. The finer grid solutions predict a faster propagation of the NAPL fronts and a considerably smaller spreading. The saturation contour plots of the BOX solution shows over- and undershoots near the corners of the lenses. The BOX result predicts almost no propagation to the right along lens C, while the other three methods, which employ here locally the IFDM discretization, yield a considerable flow of NAPL along the ascending interface of the lens. BOX, IFDM, and BOXIFDM predict that the NAPL front has reached the bottom of the flume after 85 minutes, while CVBOXIFDM yields a much smaller vertical propagation of the NAPL front and wider NAPL saturation profiles.

R. Huber, R. Helmig / Finite volume discretizations

161

To the right of lens D, IFDM generates some grid orientation effects, while BOX and BOXIFDM exhibit little and CVFE no effect of the used mesh. All four schemes predict a penetration of the trough. CVBOXIFDM predicts that the NAPL does not flow along the entire width of the trough, while the other methods yield a NAPL pool over the entire width. The predicted NAPL fronts of BOX and IFDM show the strongest propagation below the trough. 6.

Discussion

The numerical solutions of test problems I and II indicate that the use of too coarse grids for the discretization of heterogeneous media may lead to inaccurate predictions. In addition to an appropriate geometric description of zones of different material properties, capillary effects such as the entry pressure of low permeable zones require the use of finer grids along heterogeneity interfaces. Furthermore, the flows around corners cannot be captured in an appropriate way by using very coarse computational grids. On the other hand, fine grids entail a large number of nodes and unknowns for which the problem must be solved. This leads to enormous memory requirements for the used machine and long computation times. So, typically a trade-off is done between the goal of having an accurate prediction with a high spatial resolution and an acceptable computational expense. The experimental observations indicated that conceptual models which assume regions of constant material properties and the associated mathematical problems did not adequately describe the physical processes. This is the case, since apparently existing smaller-scale heterogeneities were ignored. The incorporation of heterogeneities on different length scales suggests not a deterministic but a stochastic approach to the problem. Strategies to handle such effects can be the use of anisotropic permeabilities [8], or the use of geostatistic models with random permeability fields [29]. However, these modeling issues are not discussed in greater detail here. In this study, the numerical solution to the formulated mathematical problem is of interest. Here, the discretization schemes are compared to a fine grid solution which gives the “exact” solution to the posed mathematical problem. A good discretization scheme must be characterized by giving already for coarse grids physically realistic and good estimates of the exact solution, especially in the context of heterogeneity problems. The major drawback of the CVFE method is its failure for a wide range of computational grids. Especially in three dimensions is it difficult to generate meshes which obey the positive transmissibility constraint. The BOX, IFDM, BOXIFDM and CVBOXIFDM methods do not entail any restriction to the used grid. However, IFDM exhibits for non-rectangular meshes grid orientation effects which may become severe. The numerical discretization of the fluxes as it is employed by the BOX discretization is not suitable to problems with varying material properties. If the BC model is used with nonzero entry pressures, the BOX method may yield unphysical saturation profiles.

162

R. Huber, R. Helmig / Finite volume discretizations

The two presented test problems revealed that the CVFE discretization introduces a strong numerical cross-diffusion for flows oriented along the edges of rectangular meshes and may considerably underpredict the propagation of plumes. If there exist many heterogeneities on smaller length scales, they typically yield a spreading of propagating NAPL plumes. This assumption is supported by the observations of the two mentioned experimental studies [1]. The cross-diffusive property of CVFE/CVBOXIFDM mimics this physical macrodispersion better than the IFDM and BOXIFDM schemes. Furthermore, if flows occur diagonally across a rectangular mesh (e.g., five-spot problem), the CVFE discretization yields good qualitative results, whereas five-point schemes like BOX and IFDM for coarse meshes generate diamond-shape profiles [25].

7.

Conclusions

Each of the three presented basic discretization methods has some drawbacks. The BOX discretization may generate spurious oscillations at corners of heterogeneities if the Brooks–Corey model is employed. If the CVFE discretization is used in conjunction with meshes which exhibit negative transmissibilities, typically the numerical algorithm fails. This is unfavorable if complex geometries should be modeled. The IFDM discretization exhibits grid orientation effects in regions where the finite volume interfaces are not perpendicular to the connecting line between two nodes. Two combinations of these discretization schemes were suggested. The BOXIFDM method alleviates severe grid orientation effects in homogeneous regions and ensures that flows around corners of lenses are accurately described and no unphysical oscillations in the saturation distribution are predicted. The combination of all three basic schemes, denoted here as CVBOXIFDM, allows the use of CVFE for all elements which have no negative transmissibilities. For all other elements the BOXIFDM method is applied. Both combined methods are applicable to any mesh without generating unphysical saturation profiles at heterogeneities. The basic discretization scheme of CVBOXIFDM is CVFE which represents for rectangular grids a nine-point scheme, whereas BOXIFDM is exclusively based on five-point schemes.

References [1] B. Barczewski and H.P. Koschitzky, in: Groundwater and Subsurface Remediation, The VEGAS Research Facility: Technical Equipment and Research Projects, eds. H. Kobus, B. Barczewski and H.P. Koschitzky (Springer, Berlin, 1996) pp. 129–157. [2] R.H. Brooks and A.T. Corey, Hydraulic properties of porous media, Hydrol. Pap. 3, Colorado State University (1964). [3] M.A. Celia and P. Binning, A mass conservative numerical solution for two-phase flow in porous media with application to unsaturated flow, Water Res. Res. 28 (1992) 2819–2828.

R. Huber, R. Helmig / Finite volume discretizations

163

[4] G. Chavent and J.E. Roberts, A unified physical presentation of mixed, mixed-hybrid finite elements and standard finite difference approximations for the determination of velocities in waterflow problems, Adv. Water Resources 14 (1991) 329–348. [5] V. Dalen, Simplified finite-element models for reservoir flow problems, SPE J. (1979) 333–343. [6] L.J. Durlofsky, Accuracy of mixed and control volume finite element approximations to Darcy velocity and related quantities, Water Res. Res. 30 (1994) 965–973. [7] R.W. Falta, K. Pruess, I. Javandel and P.A. Witherspoon, Numerical modeling of steam injection for the removal of nonaqueous phase liquids from the subsurface 2. Code validation and application, Water Res. Res. 28 (1992) 451–465. [8] C.R. Faust, J.H. Guswa and J.W. Mercer, Simulation of three-dimensional flow of immiscible fluids within and below the unsaturated zone, Water Res. Res. 25 (1989) 2449–2464. [9] P.A. Forsyth, A control volume finite element approach to NAPL groundwater contamination, SIAM J. Sci. Statist. Comput. 12 (1991) 1029–1057. [10] P.A. Forsyth and B.Y. Shao, Numerical simulation of gas venting for NAPL site remediation, Adv. Water Resources 14 (1991) 354–367. [11] P.A. Forsyth and M.C. Kropinski, Monotonicity considerations for saturated–unsaturated subsurface flow, Technical report CS-94-17, Department of Computer Science, University of Waterloo, Ontario. [12] P.A. Forsyth, Y.S. Wu and K. Pruess, Robust numerical methods for saturated–unsaturated flow with dry initial conditions in heterogeneous media, Adv. Water Resources 18 (1995) 25–38. [13] R. Helmig, Multiphase Flow and Transport Processes in the Subsurface (Springer, Berlin, 1997). [14] R. Helmig and R. Huber, Comparison of Galerkin-type discretization techniques for two-phase flow in heterogeneous porous media, Adv. Water Resources 21 (1998) 697–711. [15] R. Huber and R. Helmig, Multiphase flow in heterogeneous porous media: A classical finite element method versus an IMPES-based mixed FE/FV approach, Internat. J. Numer. Methods Fluids 29 (1999) 899–920. [16] R. Huber, Compositional multiphase flow and transport in heterogeneous porous media, Ph.D. thesis, TU Braunschweig, 1999; also report 1/2000, Institut f¨ur ComputerAnwendungen im Bauingenieurwesen. [17] T.H. Illangasekare, J.L. Ramsey, K.H. Jensen and M.B. Butts, Experimental study of movement and distribution of dense organic contaminants in heterogeneous aquifers, J. Cont. Hydr. 20 (1995) 1–25. [18] B.H. Kueper, W. Abbot and G. Farquhar, Experimental observations of multiphase flow in heterogeneous porous media, J. Cont. Hydr. 5 (1989) 83–95. [19] B.H. Kueper and E.O. Frind, Two-phase flow in heterogeneous porous media 2. Model application, Water Res. Res. 27 (1991) 1059–1070. [20] F.W. Letniowski and P.A. Forsyth, A control volume finite element method for three-dimensional NAPL groundwater contamination, Internat. J. Numer. Methods Fluids 13 (1991) 955–970. [21] S.R.D. Lunn and B.H. Kueper, Removal of pooled dense, nonaqueous phase liquid from saturated porous media using upward gradient alcohol floods, Water Res. Res. 33 (1998) 2207–2219. [22] J.E. McCray and R.W. Falta, Numerical simulation of air sparging for remediation of NAPL contamination, Groundwater 35 (1997) 99–110. [23] T.N. Narasimhan and P.A. Witherspoon, An integrated finite difference method for analyzing fluid flow in porous media, Water Res. Res. 12 (1976) 57–64. [24] M.J. de Neef and J. Molenaar, Analysis of DNAPL infiltration in a medium with a low-permeable lens, J. Comput. Geosciences 1 (1997) 191–214. [25] S. Panday, Y.S. Wu, P.S. Huyakorn and E.P. Springer, A three-dimensional multiphase flow model for assessing NAPL contamination in porous and fractured media, 2. Porous medium simulation examples, J. Cont. Hydr. 16 (1994) 131–156. [26] S. Panday, Y.S. Wu, P.S. Huyakorn, S.C. Wade and Z.A. Saleem, A composite numerical model for assessing subsurface transport of oily wastes and chemical constituents, J. Cont. Hydr. 25 (1997)

164

R. Huber, R. Helmig / Finite volume discretizations

39–62. [27] K.D. Pennell, G.A. Pope and L.M. Abriola, Influence of viscous and buoyancy forces on the mobilization of residual tetrachloroethylene during surfactant flushing, Env. Sci. Tech. 30 (1996) 1328–1335. [28] M.H. Schroth, J.D. Istok, M. Oostrom and M.D. White, Multifluid flow in bedded porous media: laboratory experiments and numerical simulations, Adv. Water Resources 22 (1998) 169–183. [29] B.E. Sleep and J.F. Sykes, Compositional simulation of groundwater contamination by organic compounds 2. Model applications, Water Res. Res. (1993) 1709–1718. [30] C.J. van Duijn, J. Molenaar and M.J. de Neef, The effect of capillary forces on immiscible two-phase flow in heterogeneous porous media, Transport in Porous Media 21 (1995) 71–93. [31] C.J. van Duijn and M.J. de Neef, Similarity solution for capillary redistribution of two phases in a porous medium with a single discontinuity, Adv. Water Resources 21 (1998) 451–461. [32] M.T. van Genuchten, A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J. 44 (1980) 892–898.