VZJ Download to Citation Manager
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
 QUICK SEARCH:   [advanced]


     


Published online 1 August 2007
Published in Vadose Zone J 6:471-482 (2007)
DOI: 10.2136/vzj2006.0171
© 2007 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Web of Science (3) Citing Articles via Web of Science (3)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Fagerlund, F.
Right arrow Articles by Niemi, A.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Fagerlund, F.
Right arrow Articles by Niemi, A.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Fagerlund, F.
Right arrow Articles by Niemi, A.
Related Collections
Right arrow Multiphase Fluid Flow
Right arrow Experiment Design

ORIGINAL RESEARCH

Nonaqueous-Phase Liquid Infiltration and Immobilization in Heterogeneous Media: 1. Experimental Methods and Two-Layered Reference Case

F. Fagerlunda,*, T.H. Illangasekareb and A. Niemia

a Dep. of Earth Sciences, Uppsala Univ., Villavägen 16, 75236 Uppsala, Sweden
b Center for Experimental Study of Subsurface Environmental Processes (CESEP), Environmental Science and Engineering, Colorado School of Mines, Golden, CO 80401-1887

* Corresponding author (ffagerlu{at}mines.edu).

All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.


Received 24 November 2006.



    ABSTRACT
 TOP
 ABSTRACT
 INTRODUCTION
 Theory and Methods
 Results and Discussion
 Conclusions
 REFERENCES
 
Accurate data to understand the migration and entrapment of nonaqueous-phase liquids (NAPLs) in heterogeneous formations are presently lacking. A series of well-controlled laboratory experiments were conducted to investigate the infiltration and subsequent immobilization of dense NAPLs in saturated heterogeneous media. The focus of this first study was the development of a special experimental methodology for measuring the dynamic evolution of a NAPL plume in space and time. To demonstrate the method, a reference case of a two-layered formation consisting of two homogeneous sands separated by a dipping interface is presented. The dipping formation in the reference case allows the study of NAPL behavior at texture interfaces under the influence of both capillary and gravitational forces. The NAPL-saturation measurement methodology, based on a multiple-energy x-ray attenuation technique, correctly captured the known injected NAPL volume as well as the general spreading and entrapment behavior in space and time. Time-continuous measurements of NAPL saturations allow the study of the history dependence of entrapped saturations. The Land model predicted the observed trend in the entrapment behavior well. The entrapment architecture was parameterized using spatial moments and moments of mass distribution at different saturations. The general features of the NAPL architecture were successfully characterized by a simultaneous interpretation of these moments, while the domination of discontinuous or continuous NAPL was captured by the ganglia/pool ratio.

Abbreviations: DNAPL, dense nonaqueous-phase liquid • GPR, ganglia/pool ratio • NAPL, nonaqueous-phase liquid


    INTRODUCTION
 TOP
 ABSTRACT
 INTRODUCTION
 Theory and Methods
 Results and Discussion
 Conclusions
 REFERENCES
 
Dense nonaqueous-phase liquids (DNAPLs), such as chlorinated hydrocarbons, often cause serious and persistent environmental problems when present in the subsurface. These DNAPLs can penetrate the water table and after immobilization may act as long-term sources of soil and groundwater contamination. Dissolution of contaminants from the DNAPL source controls the delivery of mass to downgradient contaminant plumes and consequently determines the severity and longevity of contamination problems.

Several studies investigating the process of DNAPL dissolution have concluded that the DNAPL entrapment architecture in the source zone strongly influences the rate of mass transfer from the NAPL to the aqueous phase (e.g., Parker and Park, 2004; Phelan et al., 2004; Soga et al., 2004; Bradford et al., 2003; Saenton et al., 2002; Dekker and Abriola, 2000b; Oostrom et al., 1999; Powers et al., 1994, 1998; Mayer and Miller, 1996). Two important factors that contribute to rate-limited mass transfer include the NAPL–water contact area at the pore scale and the water flow velocity in the vicinity of the sources. These two factors depend on the subsurface morphology or "architecture" of the NAPL and its distribution in the nonuniform groundwater flow field. The groundwater flow field in the source zone is affected by the reduction in aqueous-phase permeability due to NAPL entrapment.

Formation of source zone architecture after immobilization of flowing DNAPL has been studied experimentally in the field (Kueper et al., 1993; Poulsen and Kueper, 1992) and in the laboratory (Oostrom et al., 1999; Hofstee et al., 1998; Illangasekare et al., 1995; Compos, 1998; Held and Illangasekare, 1995; Smith and Zhang, 2001) as well as by numerical modeling (Christ et al., 2005; Lemke et al., 2004; Dekker and Abriola, 2000a; Bradford et al., 1998; Rathfelder and Abriola, 1998; Kueper and Frind, 1991). Existing field data (Kueper et al., 1993; Poulsen and Kueper, 1992) do not capture the dynamic behavior of the NAPL, however, because the mapping of immobile NAPL saturations was done by excavation. Most previous laboratory studies (Kueper et al., 1989; Illangasekare et al., 1995; Hofstee et al., 1998; Oostrom et al., 1999), in turn, have investigated relatively simple designs of heterogeneity consisting of one or a series of lenses in a homogeneous surrounding medium. Furthermore, in these studies, the transient flow of the NAPL was only monitored visually and recorded photographically, while saturations were measured only at the end of the experiments using dual-energy gamma-attenuation techniques.

To address some of these limitations in the test configurations and data, a set of well-controlled laboratory experiments were designed and conducted. In the first part of the study presented here, an experimental methodology for continuous monitoring and measurement of fluid saturations in space and time was developed. The methodology was applied to a two-layered formation consisting of two relatively homogeneous zones separated by a dipping interface. A NAPL was injected in the coarser zone and migrated first toward the interface and then along it under the influence of both gravity (buoyancy) and capillary forces. The methodology allows generation of data to carefully study transient flow and final entrapment of NAPLs in source zones. In the second part of the study (Fagerlund et al., 2007), the effects of heterogeneity in each of the two zones of the dipping formation were investigated.

In numerical modeling using hysteretic constitutive relations, a model by Land (1968) has commonly been used to predict entrapped NAPL saturations (e.g., Parker and Lenhard, 1987; Lenhard, 1992; Van Geel and Roy, 2002; Gerhard and Kueper, 2003; Lenhard et al., 2004). While Land (1968) based his entrapment model on entrapped gas saturation data, few experimental studies (e.g., Van Geel and Roy, 2002) have specifically investigated the saturation-history dependence of entrapped NAPL saturations in NAPL–water systems.

The detailed data set presented here is unique because NAPL saturations (not only the front) were monitored continuously in space and time using a multiple-energy x-ray attenuation system. In a study of three-phase flow, Kechavarsi et al. (2005) achieved time-continuous monitoring of fluid saturations by using digital image analysis. With this technique, saturations are measured at the walls of the test tank as visible from the outside. In contrast, the x-ray attenuation method used in this study more accurately measures average saturations across the thickness of the soil in the tank. Furthermore, the methodology presented here is based on measurements of NAPL saturation with time at a large number of points in space. Therefore, in addition to capturing the dynamic behavior of the system, this method allows quantitative study of the effects of NAPL saturation history on the final entrapment.


    Theory and Methods
 TOP
 ABSTRACT
 INTRODUCTION
 Theory and Methods
 Results and Discussion
 Conclusions
 REFERENCES
 
Equivalent Upside-Down Setup
The objective of the experiment was to gain understanding of DNAPL infiltration and immobilization in a water-saturated system. Experimentally, however, there are several advantages in using nontoxic light nonaqueous-phase liquids (LNAPLs) instead of the toxic DNAPLs that are usually of interest (e.g., trichlorethylene, perchloroethylene, or trichloroethane). Nonhazardous DNAPLs that behave well in laboratory experiments are not easily found, whereas LNAPLs of low toxicity are readily available. The use of appropriate LNAPLs minimizes the generation of hazardous laboratory wastes and eliminates health hazards and risks. An approach using an equivalent upside-down setup was therefore evaluated. Here LNAPL injected at the bottom of the tank migrates upward driven by buoyancy forces, thus mimicking the downward migration of a DNAPL. As shown by Fagerlund et al. (2006) an LNAPL and a DNAPL with equal density differences with water (equal |{rho}w{rho}n|, where {rho}w and {rho}n are the densities of water and NAPL, respectively) behave the same in equivalent upside-down systems provided that the other fluid properties (e.g., viscosity, interfacial tension, etc.) are the same. Therefore, conceptually, DNAPL migration can be analyzed by studying LNAPL migration in an upside-down system.

Experimental Flume and Nonaqueous-Phase Liquid Infiltration
The experimental flume, as shown schematically in Fig. 1 , has inner dimensions of 70.5 by 53 by 4.7 cm (71 by 53 by 4.7 cm in the heterogeneous experiments presented in Fagerlund et al. [2007]). The tank walls are made of Plexiglas reinforced with a metal frame and lined with glass. Two distinct zones of sand with a dipping interface were packed in the tank. The dip of the interface between the two zones was set at an angle of 3.25°. The tank was packed wet, i.e., water was added before sand. To produce a smooth and straight interface between the two sand zones, the tank was tilted to the desired angle of the interface during packing. The depth of water was maintained at a few centimeters above the sand and after pouring in sand using a tube, the sand was continuously remixed and compacted using a metal rod to minimize layering.


Figure 1
View larger version (26K):
[in this window]
[in a new window]

 
FIG. 1. Experimental flume, showing the point of light nonaqueous-phase liquid (LNAPL) injection.

 
Despite the reinforcement of the flume walls, slight bulging occurred when the tank was packed with sand and water. The effect was small enough not to have any large influence on the dimensionality of fluid flow in the tank. To correctly calculate NAPL saturations and porosities, however, the exact width of the tank during water loading was measured using the x-ray attenuation system. In the area of main interest with respect to NAPL flow and entrapment in the flume, the width of the tank varied from approximately 4.6 to 4.8 cm, with an average of 4.7 cm.

A known volume of LNAPL was injected at a constant rate in the lower coarser layer and was allowed to migrate up toward the dipping interface. The test LNAPL used was a mixture of Soltrol 220 spiked with iodoheptane (10% w/w) and colored red with Sudan IV (0.1% w/w). The mixture had a density of 836 kg m–3. Measurements of viscosity and interfacial tension for similar mixtures have been reported in the literature. For a mixture of Soltrol 220 and 16% iodoheptane by weight (1:9 v/v), colored using Sudan III, Oostrom and Lenhard (1998) reported a viscosity of 4.7 Pa s and interfacial tension with water {sigma}nw = 0.036 N m–1. Two constant-head reservoirs at each end of the tank were used to control the groundwater flow across the length of the tank.

For the reference case presented here, the head difference across the tank that maintained a small discharge was 1.5 cm, the injected volume was 500 cm3, and the injection rate was 3.23x 10–5 kg s–1. The duration of the injection was 2.63 h, after which the redistribution and immobilization of the NAPL was observed for 22 d using both x-ray techniques and photographs, and then an additional 2 wk using only photographs. The lower coarser zone was homogeneous 0.6-mm (no. 30) sand and the upper finer zone was homogeneous 0.3-mm (no. 50) sand. Intrinsic permeability, k (L2), as well as Brooks–Corey air–water displacement pressure, Pd (M L–1 T–2), and pore-size distribution index, {lambda}, of the sands have been measured by Sakaki and Illangasekare (personal communication, 2006) and are given in Table 1. The values were measured in columns packed with sand in a fashion similar to the flume.


View this table:
[in this window]
[in a new window]

 
TABLE 1. Properties of the sands used in the study; Brooks–Corey pore-size distribution index ({lambda}) and air–water displacement pressure (Pd) refer to primary drainage.{dagger}

 
In all experiments, the origin of the xz coordinate system is defined at the center of a lead collimator in the lower left corner of the left constant-head well. The reason is that the 2- by 2-mm hole in the collimator allows exact verification of the position of the x-ray beam with respect to the flume. For this coordinate system, with x as the horizontal and z as the vertical axis, the injection coordinates for the reference case are (x,z) = (51.5, 5.2). The NAPL is injected along a line source across the depth (y dimension) of the flume, creating a point source in two dimensions.

X-Ray Measurement Methods
The NAPL saturation in the tank was estimated by measuring the equivalent material path lengths when a collimated beam of x-ray traverses the soil containing the NAPL. The photons in the x-ray beam are attenuated by different materials (sand, water, and NAPL) that they encounter along their path. The x-ray attenuation system consists of a tube that generates a beam of photons across a range of energies within the x-ray band and a detector that receives and counts all photons that are not attenuated within the tank containing the soil. The system is schematically shown in Fig. 2 .


Figure 2
View larger version (13K):
[in this window]
[in a new window]

 
FIG. 2. X-ray attenuation measurement system used in the study.

 
At the detector, photons of different energies are first counted as the x-ray beam travels through the flume at a reference state. Subsequent photon counts taken at the exact same locations in space will reveal how photons are attenuated within different energy bands. Thereby, changes in material path lengths can be calculated provided that the attenuation properties of the materials that change are known. The lengths lß of the materials (L) (ß = w for water, n for NAPL, s for solid silica sand) are calculated using a maximum-likelihood-estimator methodology (developed and described in more detail by Hill et al., 2002), based on Lambert–Beer's law:

Formula 1[1]
Here Y{varepsilon} is the number of photons of energy {varepsilon} (M L2 T–2) observed at the detector for the reference state of the system, X{varepsilon} is the number of photons of energy {varepsilon} observed with changes in material lengths {Delta}lß (L) that the photon beam passes through, {rho}ß (M L–3) is the density of material ß, and µß{varepsilon} (L2 M–1) is the attenuation coefficient of material ß at energy {varepsilon}. The attenuation coefficients for the three materials—water, NAPL mixture, and silica sand—as measured for the system, are shown in Fig. 3 . The material length lß is negative if a material is removed and positive if it is added, compared with the reference state. Further description of the x-ray system and the methodology for material path length calculation is given by Ferré et al. (2005), Hill et al. (2002), and Hill (unpublished data, 2001).


Figure 3
View larger version (17K):
[in this window]
[in a new window]

 
FIG. 3. Attenuation coefficients as a function of photon energy for the nonaqueous-phase liquid mixture (red), silica sand (black), and water (blue).

 
The x-ray tube and detector are mounted so that they can be moved in synchronization on each side of the flume using high-precision servo motors. The system is automated to take point measurements at predetermined locations. One measurement, including movement of the x-ray tube and detector, takes approximately 35 s. The diameter of the collimated photon beam is about 2 mm.

The material length measurements are most precise when few materials (preferably two or less) change with respect to the reference system and when these materials have well-separated attenuation coefficients. Porosity, {Phi}, is best measured in two steps. First, the inner width of the water-filled flume, lt(x,z) (L), is determined by using an empty flume as a reference system and then filling the flume with water. As will be shown below, lt varies slightly throughout the flume due to bulging of the flume walls and therefore needs to be measured as a function of the spatial coordinates. In the second step, the water-filled flume is used as reference and the length of solid sand, ls(x,z) (L) (not including voids between grains), is measured after packing the tank, thus replacing water with sand. Porosity is then calculated as

Formula 2[2]
The NAPL path lengths, ln(x,z,t) (L), are measured using the packed flume containing sand and water as a reference, thus measuring the length of water replaced by NAPL at a given time, t. To achieve accuracy in this measurement, the attenuation coefficients of the NAPL and water need to be significantly different. Because pure Soltrol 220 has attenuation properties similar to those of water, iodoheptane (10% w/w) was added to the Soltrol. The value of ln(x,z,t) together with the length of void space available for fluids lv(x,z), given by lt(x, z) – ls(x,z), is used to calculate NAPL saturation Sn(x,z,t):

Formula 3[3]
During the NAPL injection and redistribution, NAPL path lengths were measured continuously. As the movement of the NAPL slowed down, the number of measurement points was progressively increased in several steps. The reason for this measurement scheme was that there is a trade-off between the number of spatial points and the total time it takes to measure each point once. At early times when NAPL saturations change fast as the NAPL spreads and moves into uncontaminated areas, measurements were made at relatively few points to monitor changes in saturations accurately at short time intervals. At later times, when changes in Sn were smaller but the size of the NAPL plume had increased, more measurement points were included to provide better spatial coverage.

For the reference case presented here, 48 spatial points (~25 min per scan cycle) were measured at 3 h after injection, 90 points at 7 h, 112 points at 8 h, and 126 points at 26 h. The increasing number of measurement points is shown in Fig. 4 . To capture the behavior at the interface between the two sand zones, measurements were taken at 1 cm above the interface as well as at 0.2, 1, and 3 cm below it. The relatively low number of spatial points resulted in excellent monitoring of saturations with time also at later times. In experiments with heterogeneous sands, a denser measurement spacing is needed to resolve the small-scale spatial heterogeneities, as shown in Fagerlund et al. (2006, 2007). By returning to take measurements at the same points, data on NAPL saturation as a function of time were recorded. At all points of measurement, Sn can be plotted against time, and for any given time, a value of Sn can be determined by interpolation. Thereby a spatial image of NAPL distribution in the flume can be produced for any time.


Figure 4
View larger version (16K):
[in this window]
[in a new window]

 
FIG. 4. X-ray measurement points: green squares, early time; red filled circles, additional points at intermediate times; blue empty circles, additional points at late times. The injection point is indicated with an x and points shown below in Fig. 10 are circled.

 

Figure 10
View larger version (17K):
[in this window]
[in a new window]

 
FIG. 10. Nonaqueous-phase liquid (NAPL) saturation with time at three points: (x,z) = (54,7.7), circles; (54,12.7), squares; (29,27.4), crosses. Figure 4 shows the locations of these points in the flume.

 
Model of Ganglion Entrapment
Land (1968) proposed that entrapped nonwetting-phase saturation is a function of the maximum nonwetting saturation. This saturation occurs at the reversal point between drainage (drying) and imbibition (wetting). The higher the maximum nonwetting saturation, the broader is the range of pore sizes into which the nonwetting phase enters and subsequently can be entrapped. In a NAPL–water system, where the NAPL is the nonwetting phase, the Land (1968) model takes the following form:

Formula 4[4]
Here Formula 4nt is the effective entrapped NAPL saturation, Formula 4nmax is the effective maximum NAPL saturation previously reached, and

Formula 5[5]
where Formula 5ntmax is the maximum possible effective entrapped NAPL saturation, corresponding to the case where the medium has been fully saturated by NAPL, that is, Formula 5nmax. For a given homogeneous porous medium (e.g., a type of sand), Formula 5ntmax can be seen as a curve-fitting parameter for the Land model, describing the entrapment properties of the medium. Effective saturations, here denoted with an overbar, are defined as saturations scaled to exclude residual or irreducible water saturation, which is thereby regarded as part of the solid medium:

Formula 6[6]
where Swr is the residual water saturation. Note that, following Lenhard et al. (2004), we refer to the immobile wetting phase as residual and the immobilized nonwetting phase as entrapped to make a distinction between the different immobilization mechanisms. The NAPL, being the nonwetting phase in a NAPL–water system, is immobilized as water imbibes entrapping blobs and ganglia of NAPL, which become discontinuous and occluded by water.

A simpler model of entrapment is achieved by letting the entrapped NAPL saturation Snt = Sn for Sn≤Sntmax* and Sn>Sntmax*, where is Sntmax* the maximum, constant entrapped saturation. The model is independent of saturation history and all NAPL is assumed immobile until Sn exceeds the constant Sntmax*. This model is referred to as the constant entrapment model and is implicitly assumed when entrapped or residual NAPL is included as a constant in the NAPL relative permeability (krn) function. Such krn functions have been presented in various textbooks, for example, Helmig (1997).

The time series of NAPL saturations recorded at all points of measurement provide detailed information about the saturation history. Both the maximum saturation and final immobile saturation have been registered for a wide range of peak saturations. For the type of medium investigated (uniform sands), the data can therefore be used to test the suitability of models of NAPL entrapment based on saturation history, such as the Land (1968) model.

Spatial Moments
Spatially dense measurements of NAPL saturations allow accurate calculation of spatial moments, which can be used to describe the average spatial behavior of the NAPL plume with time. In an analogy with moments of solute concentration (Freyburg, 1986), the spatial moments of NAPL volume distribution (around the origin of the coordinate system) Mijk(t) (L3+i+j+k) are defined as

Formula 7[7]
where x0, x1, y0, y1, z0, and z1 are the limits of integration, which here are defined as the (inner) limits of the experimental domain, inside the constant-head wells. The order of the moment is given by i + j + k. By multiplying with {rho}n, which here is a constant, moments of mass distribution are obtained. The moments normalized by the total volume or mass are, however, identical.

The total volume of NAPL in the tank, Vn(t) (L3), is given by the spatial integral of NAPL length, ln, across the tank:

Formula 8[8]
Noting that

Formula 9[9]
the zeroth-order moment M000(t) also equals Vn(t).

Using average properties across the y dimension in a two-dimensional xz system, the first-order moments around the origin, normalized by the zeroth moment, define the location of the center of mass (L) (xc,zc):

Formula 10[10]
The two-dimensional second moments taken around the center of mass define a 2 x 2 covariance matrix with diagonal elements {sigma}xx and {sigma}zz (variances) (L2) and covariance elements {sigma}xz and {sigma}zx (L2) where

Formula 11[11]
The coefficient of skewness, Cs, describes asymmetry of the mass distribution around the center of mass. It is defined as the third moment around the center of mass normalized by the variance to the power of 3/2 (i.e., the standard deviation to the power of 3). In the x and z direction, Cs can be calculated as follows (Råde and Westergren, 1995):

Formula 12[12]
The kurtosis, Ck, describes the peakedness of the distribution. It is defined as the fourth moment around the center of mass normalized by the square of the second moment, and subtracting 3. In the x direction (and similarly for z) it can be calculated as (Råde and Westergren, 1995)

Formula 13[13]


    Results and Discussion
 TOP
 ABSTRACT
 INTRODUCTION
 Theory and Methods
 Results and Discussion
 Conclusions
 REFERENCES
 
Images of Spatial Evolution of Nonaqueous-Phase Liquid Saturation
To calculate NAPL saturations from Eq. [3], the spatial distribution of porosity in the tank was measured before each NAPL infiltration experiment. For the reference case, the average porosity ({Phi}) of the lower coarser layer was 0.452 with a standard deviation of 0.009 and for the upper finer layer, the average {Phi} was 0.467 with standard deviation 0.012.

Figure 5 shows the evolution of spatial NAPL distribution with time. Contours of NAPL saturation, measured with x-ray attenuation methods, have been plotted on top of digital images. A volume of 500 cm3 of NAPL was injected at a constant rate in 2.6 h. At the end of injection, the NAPL had almost reached the interface between the sand zones. At 4.3 h, a pool of NAPL was building up at the interface and was starting to spread along the interface. The dip of the interface was mild enough to permit some spreading also downward to the right. A small effect of the hydraulic gradient pushing the NAPL body to the left can be observed. At 9.6 and 23.9 h, more mass had moved up to the interface, but there was not enough pressure in the NAPL pool to produce penetration into the upper finer zone. Instead, the NAPL migrated along the interface and, at approximately 22 h after the start of injection, it started to move out of the experiment domain into the constant-head well on the left side. After 505.6 h (~3 wk), there was only a thin pool of NAPL left at the interface. Elsewhere in the tank, the NAPL was immobile. After 5 wk (not shown), the pools had disappeared and only entrapped NAPL remained in the tank.


Figure 5
View larger version (85K):
[in this window]
[in a new window]

 
FIG. 5. Evolution of nonaqueous-phase liquid (NAPL) distribution, with x-ray attenuation measurements of NAPL saturations (contours) superimposed on digital images.

 
In the digital images in Fig. 5 , some preferential flow channels and irregularities in the shape of NAPL body can be observed. Although the tank was packed to achieve a homogeneous packing of each zone, some heterogeneity in the form of slight layering existed in places. This produced some preferential flow in the direction parallel to the interface. It can be concluded that even very subtle heterogeneities in a nearly homogeneous medium can have an effect on the NAPL migration. In this case, the effects on the average migration behavior and the spatial moments were deemed minor. Additionally, where layering exists, the NAPL may be entrapped at higher saturations due to capillary-barrier effects. That this effect also was very small can be seen below in the entrapment analysis.

In this experiment, NAPL-saturation measurements were taken on a relatively coarse mesh compared with the following experiments with heterogeneous sands (Fagerlund et al., 2006, 2007), where a finer spacing was used. This means that fewer points in space were measured, but that these points were measured more often in time (Fig. 5a is based on 48 spatial points, Fig. 5b and 5c on 112 points, and Fig. 5d on 126 points). The result is that the behavior with time was well captured while at the same time some of the irregular details were missed in the contour plots. The x-ray attenuation measurements can also "see" NAPL that is not present directly adjacent to the glass wall, however, and can therefore include some NAPL that may not be seen in the digital images. The x-ray measurements captured the average behavior of the NAPL very well.

Spatial Moments of the Nonaqueous-Phase Liquid Distribution
Spatial moments of the NAPL distribution are based on spatial contour plots of the type shown in Fig. 5 and are used to quantify the spatial distribution of NAPL at various times. The zeroth moment, or the total NAPL volume inside the experimental domain, Vn(t), is shown in Fig. 6 as triangles. After approximately 5 h, enough x-ray measurements had been made to record the whole NAPL body. As can be seen, the measured Vn was quite constant during the redistribution until, after ~22 h, the NAPL started to exit the tank. For this period, the measured total volume can be compared with the known injected volume of 500 cm3. Determination of Vn is based on linear interpolation between the data points. In a few places where the outermost measured data points show positive NAPL saturations, extra points of zero saturation were added based on the digital images to correctly define the edge of the NAPL body. Testing showed that the estimation of Vn was not very sensitive to the method of interpolation. The linear method worked well for describing the abrupt edge of the NAPL pool at the interface and was therefore chosen in the analysis.


Figure 6
View larger version (10K):
[in this window]
[in a new window]

 
FIG. 6. Total nonaqueous-phase liquid volume (Vn) inside the experimental domain, based on original measurements (triangles) and after calibration (circles).

 
Because the estimated Vn was quite constant even though the NAPL was redistributing, the number and location of measurement points seem adequate for capturing Vn. The slight overestimation of Vn seen in Fig. 6 is probably a result of the fine-tuning of the x-ray measurement system (live time during photon count, etc.), leading to a slight difference in actual and expected attenuation used in the material-path-length calculations. Because the total injected volume is known, however, the measurements of ln and Sn can be calibrated by using the volume estimate. The application of a (linear) correction factor of 0.97 to all measurements of ln and Sn resulted in the correct volume estimate (marked with circles in Fig. 6), which was stable at 500 cm3 until the NAPL started to exit through the constant-head well. At 515 to 525 h, when the NAPL was immobile everywhere except in a thin pool at the interface (see Fig. 5d), roughly half of the injected NAPL remained in the experiment domain.

Figure 7 shows the first moments around the origin, defining the center of mass (xc,zc), and the second moments around the center of mass, defining the spatial variances ({sigma}xx, {sigma}zz). As more and more mass moved up and pooled at the interface, zc shifted upward. Because the NAPL after a while started to exit the measurement domain, however, eventually a point was reached when immobile NAPL started to dominate and, as the pool at the interface was depleted, zc shifted downward. The point xc moved to the left because of the NAPL movement along the interface. It can be noted that the trend in xc was never reversed because the pool did not exhibit a slug-like movement, but rather became thinner and thinner as its mass was depleted. The variance in the x direction first increased because NAPL spread at the interface, then during depletion of NAPL from the domain, it continued to increase because saturations mostly decreased near xc, increasing the relative importance of NAPL in the periphery. The variance {sigma}zz was much smaller than {sigma}xx and had an early peak when the NAPL was spread out evenly between the injection point and the interface, but then decreased as mass was concentrated along the interface. In the end, {sigma}zz became larger because, as for {sigma}xx, mass near zc was depleted. The small discontinuity at 26 h, seen most clearly in the plot of xc, occurred with an increase in the number of measurement points, when a row of additional measurements near the interface was added. Although this has no significant effect on the recorded total NAPL volume (see Fig. 6), the slightly increased detail in the information about the spatial NAPL distribution near the interface caused a small shift in the estimated location of xc.


Figure 7
View larger version (20K):
[in this window]
[in a new window]

 
FIG. 7. Center of mass of nonaqueous-phase liquid plume (xc,zc), and variances ({sigma}xx, {sigma}zz); x-directional moments (circles) with axis scale shown on the left side, z-directional moments (triangles) with axis scale on the right side.

 
The covariance, shown in Fig. 8 , was negative because the main spreading of the NAPL occurred in the second and fourth quadrants with respect to the center of mass (i.e., the NAPL spread toward horizontal [x] and vertical [z] coordinates of opposite sign), which is a result of the interface dipping up to the left. Because the NAPL spread more horizontally, the trend in magnitude of {sigma}xz was similar to that of {sigma}xx. Just like the variances, the covariance continued to increase (in magnitude) at later times because the relative importance of the mass located in the periphery increased as mass in the pool zone, close to (xc,zc), was depleted.


Figure 8
View larger version (12K):
[in this window]
[in a new window]

 
FIG. 8. Covariance ({sigma}xz) of the spatial distribution of nonaqueous-phase liquid in the horizontal (x) and vertical (z) directions.

 
Figure 9 shows the coefficients of skewness (Cs) and kurtosis (Ck) of the NAPL distribution in the x and z directions. Both Csx and Csz were negative because, in both cases, the plume was more spread out to the negative side (left and below) of the center of mass. In the z direction, the process where mass moved away from the lower regions and the plume at the interface became more pronounced can be seen both in the skewness and kurtosis. The value of Csz became more negative as zc moved upward and Ckz increased as the peak shape of the pool became sharper. When the pool was depleted, these trends were reversed as the distribution became flatter and zc moved back down. In the x direction, the skew became more negative as more mass moved toward the left in the pool. The kurtosis had its maximum when high saturations (<0.5, see Fig. 5) existed in a narrow region. The value of Ckx exhibits a small jump at 26 h with the increase in points of measurement near the interface also mentioned above, indicating that the sharpness of the saturation peak was slightly overestimated at earlier times (i.e., before 26 h).


Figure 9
View larger version (22K):
[in this window]
[in a new window]

 
FIG. 9. Coefficients of skewness and kurtosis of the nonaqueous-phase liquid distribution in the x (circles) and z (triangles) directions.

 
For practical purposes, the vertical high-order spatial moments (Csz, Ckz) seem to contain more useful information than the horizontal ones (Csx, Ckx). For the equivalent upside-down case when a DNAPL moves downward, a strong positive vertical skew (high and positive Csz) can be expected when the NAPL has moved down, leaving a trace of entrapped NAPL behind and then pooled on a horizontally oriented layer. The center of mass is near or at the pool but a substantial amount of NAPL also remains trapped above the pool. The distribution either consists of one pool zone (high saturation) with entrapped ganglia (low saturation) above, or multiple pools and ganglia with one main pool zone at the bottom of the distribution. If, on the other hand, the distribution consists only of homogeneously entrapped ganglia or of one single pool zone, the skew will be close to zero. Positive vertical kurtosis (Ckz) can be expected when a pool of high saturation exists, whereas if the distribution is flat with no dominating pool, Ckz will be negative.

Nonaqueous-Phase Liquid Saturation as a Function Time
Figure 10 shows the development of NAPL saturation (Sn) with time at three representative points, the locations of which are indicated also in Fig. 4. The three points are at different vertical locations and are reached by the NAPL front at different times. At all three points, Sn declined essentially linearly, given that the peak saturation had been reached and injection of NAPL had stopped. At the lowest point (x,z) = (54,7.7) closest to the injection (circles), the peak in Sn was followed by a fast decline. In general, at points where vertical migration toward the interface occurred, the rate of decline in Sn decreased with decreasing distance to the interface. Also generally, Sn stabilized at roughly the same constant value when the NAPL became immobile. The NAPL was first immobilized in the lower parts of the flume. Within 2-cm distance from the interface, where the point (29,27.4) is located, a very thin mobile pool of NAPL still existed during the second measurement period between 505 and 530 h. This is the reason why the final saturation was considerably higher at point (29,27.4) than the other two points shown in Fig. 10.

The differences in the rates of decline of Sn are not surprising and should result from the effects of pooling and change in direction of movement (to a direction with a smaller component of buoyancy) at the interface, which are stronger closer to the interface. It is, however, interesting that the rates of decline were essentially linear everywhere (i.e., at all measured locations) in the flume where the NAPL was mobile, at least once the effect of the interface was felt by the NAPL body. In general, it seems that when NAPL is redistributing homogeneously under gravity forces, dSn/dt can be approximated as piecewise constant after the peak saturation has been reached. The appearance of effects of downstream damming may lower the rate of decline, as seems to be the case for point (54,12.7) (marked with squares in Fig. 10), where the decline was slightly faster between 3 and 7 h than between 7 and 48 h.

Entrapment
Figure 11 shows a comparison of maximum NAPL saturations, Snmax, and the saturations measured at the end of the experiment between 505 and 530 h. Points located within 2 cm of the interface, where a thin pool of NAPL still remained at this time, are marked with triangles. The circles show points located farther than 2 cm from the interface, where the NAPL had already become immobile. It is assumed that, at these points, the NAPL was entrapped as discontinuous blobs and ganglia. Hence, for this zone, the relationship between entrapped NAPL saturations Snt and maximum saturations Snmax can be studied. The Land (1968) model (Eq. [4]; solid line) and a model of constant entrapment (dotted line) have been fitted to this data. The Land model has a better least-square fit with a RMSE of 0.0205 compared with 0.0296 for the constant-entrapment model. The constant model overestimates entrapped saturations for low values of Snmax, but if a high maximum saturation has been reached, it seems to perform well. Hence, depending on the application, it may perform satisfactorily. The Land model was fitted using Swr = 0.06, and the fitted value for Sntmax was 0.17 (Formula 13ntmax = 0.18).


Figure 11
View larger version (18K):
[in this window]
[in a new window]

 
FIG. 11. Late-time nonaqueous-phase liquid (NAPL) saturation as a function of maximum NAPL saturation previously reached. Late time is defined as 505 to 530 h after injection. Points in the pool zone within 2 cm below the interface are indicated by triangles; points below this zone are circles. Models fitted to the data were Land (1968) (solid line) and constant entrapment (dotted line).

 
At the triangles, where saturations are >0.2, continuous NAPL existed at this time. As long as the NAPL in the pool slowly keeps moving away along the interface, saturation will decrease at these points. It can therefore be assumed that at later times more of the triangles will be located along the same trend as the circles and eventually the only remaining trend will be the one seen for the circles. Trends like the one in the triangles can be expected when a capillary barrier controls entrapment and NAPL is immobilized in a continuous state (i.e., as a pool). At the end of the measurement period (3 wk), the pool had still not disappeared; however, visual observations indicated that the pool disappeared after two additional weeks.

Nonaqueous-Phase Liquid Occurrence at Different Saturations
While spatial moments discussed above describe how the NAPL, on average, was distributed in space, they do not contain direct information about the relative occurrence of different NAPL saturations. The relative magnitude of high vs. low NAPL saturations is, however, important for the dissolution characteristics of the NAPL, as NAPL saturation is related both to the interfacial area between the NAPL and water and the aqueous-phase relative permeability. The occurrence of different NAPL saturations can be described by means of frequency histograms, as shown in Fig. 12 , where the volumetric fractions of NAPL contained within 1% intervals of NAPL saturation have been plotted for two different times of the experiment. Because {rho}n is constant, mass fractions equal volume fractions. In the early stage, 20 h after the start of injection, most of the NAPL was contained within a relatively high saturation range, indicating that most of the NAPL was continuous and potentially mobile. The average saturation was high and the distribution has a distinct negative skew. This figure is typical for all times during the first measurement period up to 50 h after the start of injection. At 520 h after injection, and in general during the second measurement period between 505 and 530 h, the shape of the distribution had changed considerably. The bulk of the NAPL left in the tank was now contained in a relatively narrow band around Sn = 0.15. The average saturation was low and the NAPL was largely discontinuous and immobile. The small remaining pool gives the distribution a positive skew that may become negative again if the pool completely disappears.


Figure 12
View larger version (35K):
[in this window]
[in a new window]

 
FIG. 12. Nonaqueous-phase liquid (NAPL) distribution as a function of saturation at 20 and 520 h.

 
The general characteristics of the NAPL distribution across the saturation range (Fig. 12) can also be described by its moments. The mean and variance as well as coefficients of skew and kurtosis as functions of time are shown in Fig. 13 . As the NAPL spread out, the mean saturation decreased, and as the occurrences of higher saturations disappeared, the spread around the mean (the variance) also decreased. The higher order moments, i.e., the coefficients of skew and kurtosis, need to be interpreted simultaneously with the mean. Negative skew and relatively high mean (first 50 h) indicates that the bulk of the NAPL was not entrapped as discontinuous ganglia but was potentially mobile. Positive skew and low mean (505–530 h) indicates that most of the NAPL occurrence was entrapped as ganglia but some NAPL still resided in pools. Should the last occurrences of higher saturations disappear (i.e., the NAPL contained at saturations higher than ~0.2 in Fig. 12 redistributed to lower saturations), the skew will again become negative. High kurtosis occurred both in the beginning and in the end of the measurement period, when a large part of the NAPL was contained within a narrow saturation range. Knowing the mean, the approximate location of this range is known. Therefore from Fig. 13 one can tell that at the beginning, the NAPL distribution was pool dominated, with pools of mean Sn around 0.4 to 0.45, while at the end of measurement period, the distribution was dominated by NAPL at saturations around 0.15. This is confirmed by the ganglia/pool ratio discussed below. At intermediate times when no distinct peak existed and the NAPL was spread out relatively evenly across the entire saturation range, the kurtosis is negative.


Figure 13
View larger version (26K):
[in this window]
[in a new window]

 
FIG. 13. Moments of the nonaqueous-phase liquid (NAPL) saturation (Sn) distribution.

 
Ganglia/Pool Ratio of Nonaqueous-Phase Liquid Occurrences
The characteristics of the NAPL saturation distribution can also be described by the ganglia/pool ratio (GPR; Lemke et al., 2004), defined as the mass of NAPL occurring as ganglia divided by the mass occurring as pools. The most easily implemented way to distinguish between ganglia and pools is to use a single cutoff value, above which the NAPL occurrence is regarded as a pool. Lemke et al. (2004) used the maximum entrapped saturation, Sntmax for this. If a stricter definition based on the potential mobility (state of NAPL continuity) is used, the cutoff value between ganglia and pools varies locally, depending on the properties of the porous medium and the maximum NAPL saturation previously reached. If the medium is homogeneous and the relationship between Snt and Sntmax (see e.g., Fig. 11) is relatively flat, a single, constant cutoff value may be sufficient for a reasonable ganglia and pool characterization. In the two-layered reference-case experiment presented here, the NAPL did not penetrate into the upper finer zone and, therefore, was distributed within one type of sand. For this sand (0.6-mm sand) and NAPL, Sntmax according to the Land (1968) model is fitted to 0.17. As can be seen in Fig. 11, however, immobile NAPL (ganglia) exists also at higher saturations up to Sn = 0.22, and a dividing line between the immobile NAPL (marked with circles in Fig. 11) and the mobile part of the pool zone (triangles in the upper part of Fig. 11) may be better placed at Sn ~ 0.19 to 0.20.

Figure 14 shows the GPR for different cutoff values. As can be seen, the different cutoff values result in a similar general appearance of GPR with time. Especially in the second measurement period (505–530 h), however, GPR is sensitive to the cutoff value and almost tripled between cutoff values 0.15 to 0.19. The reason for this sensitivity is that a large amount of the NAPL mass was contained within a small saturation range from Sn = 0.11 to 0.19, which can be seen in Fig. 12, for 520 h. Even though the total NAPL mass was stable and diminished very slowly (Fig. 6) during this period, the interpolated edge of the mobile part of the pool at the interface changed somewhat from one point in time to the next, thereby causing jumps in the GPR. Still, the trend that the pool zone was diminishing can be seen in the GPR within the range of tested cutoff values shown in Fig. 14. In conclusion, even if the use of a constant cutoff value results in some uncertainty and possible mischaracterization, the GPR provides information about whether pools or ganglia are dominating and the general trend with time.


Figure 14
View larger version (20K):
[in this window]
[in a new window]

 
FIG. 14. Ganglia/pool ratios as a function of time for different cutoff values: 0.15 (green squares), 0.17 (red circles), and 0.19 (yellow triangles).

 
Nonaqueous-Phase Liquid Architecture Characterization Summary
The main findings of the above characterization of NAPL architecture can be summarized as follows:

The spatial moments describe the average characteristics of the distribution of NAPL in space. The first moments define the center of the mass and the second moments provide measures of the spread. The third and fourth moments contain information about how high-saturation (pool) zones and low-saturation (ganglia) zones are distributed in space. For real contamination scenarios, the vertical higher order moments may contain more readily interpretable information than those in the horizontal direction, concerning the positions of pools and ganglia and thereby the state of the plume. Applied to analysis of controlled laboratory experiments, spatial moments are useful for the validation of numerical models striving to reproduce the average behavior observed in the experiments.

The level of NAPL saturation is closely related to the state of continuity of the NAPL, and thereby its mobility. Furthermore, saturation is related both to the interfacial area between the NAPL and the aqueous phase and the aqueous-phase permeability, which are parameters that influence NAPL dissolution. The distribution of NAPL across different saturations (Fig. 12) therefore provides information concerning both the NAPL mobility and the dissolution properties of the source zone.

Moments of the distribution across the saturation range as well as the GPR constitute average measures of how NAPL is distributed across different saturations. A joint interpretation of the higher order moments along with the mean provides a general method for characterizing the saturation distribution. For example, negative skew in combination with high mean indicates high mobility for the bulk of the mass (Fig. 12a), while a positive skew in combination with a low mean (Fig. 12b) indicates a high level of immobility. The GPR is a measure of the state of continuity of the NAPL. To calculate it, the NAPL mass needs to be characterized as either discontinuous (ganglia) or continuous (pool). Despite the fact that the NAPL resided in one single sand type, GPR turns out to be sensitive to the cutoff value selected, especially at the end of the experiment when most of the NAPL was entrapped as ganglia at saturations close to the cutoff value. Nevertheless, GPR clearly shows whether ganglia or pools dominate the distribution for the range of tested cutoff values.

The final entrapped NAPL saturation, when the NAPL has become discontinuous and immobile, depends on the maximum saturation previously reached (Fig. 11). The saturation at which the NAPL becomes immobile is important for characterization purposes and affects, e.g., the GPR. The Land (1968) model describes the observed trend of the SntSntmax relationship well; however, the spread around the trend shows that there are local variations in the NAPL immobilization even for almost homogeneously packed sand. As discussed above, slight heterogeneity exists even though care was taken to pack the sand homogeneously. These heterogeneities may produce the observed variations in entrapped saturations; however, judging also from the following experiments presented in Fagerlund et al. (2007), there will be such variations, making the SntSntmax relationship less distinct regardless of how homogeneously the sand has been packed. We hypothesize that heterogeneity on what could be called the representative elementary volume (REV) scale, commonly used to analyze subsurface transport phenomena and similar to our scale of x-ray measurements, may not be sufficient to explain these local variations in entrapment. Entrapment by snap-off is governed by processes on the pore scale that it may not be possible to upscale so that a fully deterministic relationship between SntSntmax can be found on the REV scale. Therefore, a practical possibility to improve the description of ganglion entrapment could be to add a random term to account for such local, random variations to the Land (1968) model.


    Conclusions
 TOP
 ABSTRACT
 INTRODUCTION
 Theory and Methods
 Results and Discussion
 Conclusions
 REFERENCES
 
This study was the first part of an experimental study on the migration and immobilization of NAPL in saturated heterogeneous media. Here, the methods of measurement and analysis were developed and applied to a two-layered reference case. In the second part (Fagerlund et al., 2006, 2007), the methodology was applied to more complex systems with stochastic heterogeneity, which are compared with the reference case presented here.

We presented a methodology to conduct laboratory experiments where dynamic NAPL behavior is continuously monitored in space and time. Precise measurements of fluid saturations as well as porosity were obtained using a multiple-energy x-ray attenuation technique. The measurements were shown to correctly predict the total injected NAPL volume as well as the dynamic spreading behavior of the NAPL. The advantage of the methodology in comparison with photographic recording of NAPL migration is that the x-ray measurement constitutes an average of the saturation across the width of the tank, whereas the photographic techniques record behavior at the edges only; that is, the x-ray "sees" the fluids inside the sand. Light-transmission techniques constitute an alternative; however, because the visible light is attenuated across much shorter distances than x-rays, a much thicker flume can be used with x-ray techniques, thereby allowing much less influence of edge effects at the sand–flume–wall interface on the average flow and entrapment behavior in the flume. Compared with gamma-ray attenuation, the x-ray technique is in general both more accurate and significantly faster. The approach presented here takes advantage of the speed of the x-ray measurement system; thereby a methodology to accurately capture dynamic NAPL behavior in both space and time is presented. Such methodology is not possible with slower measurement equipment such as gamma-ray attenuation devices.

The methodology was applied to a base-case study of NAPL migration and entrapment under gravitational and capillary forces in a two-layered dipping formation, where the two layers were homogeneous. The time-continuous measurements allowed the study of the dynamic NAPL-plume movement and the evolution of its saturation distribution with time, thereby allowing the study also of the dependence of entrapped saturations on the saturation history. The entrapment model by Land (1968) was shown to predict the trend in the Snt Sntmax relationship well.

Different methods for characterizing the NAPL architecture were tested. Spatial moments are useful for characterizing the spatial distribution of the NAPL plume at different times. The NAPL distribution across different saturations is related to NAPL dissolution parameters and the moments of the distribution across saturations are useful for this purpose. Also, the GPR (Lemke et al., 2004) can describe whether discontinuous (ganglia) or continuous (pool) NAPL dominates, at least for the case when the NAPL is contained in a homogeneous medium. The results showed that the actual value of the GPR is sensitive to the cutoff value chosen, but the GPR nevertheless gives a reliable indication of whether ganglia or pools dominate.

All the data and observations were reported in a detailed manner to allow other investigators to use the data, e.g., as the basis for conceptual or numerical model validations, thereby improving the reliability of model predictions in, e.g., field scenarios.


    ACKNOWLEDGMENTS
 
We thank the Swedish Rescue Services Agency (Räddningsverket), the Swedish Road Administration (Vägverket), and the U.S. National Science Foundation (award no. DMS-0222286) for providing funding for this research. Furthermore, F. Fagerlund thanks Dr. T. Sakaki and Dr. M. Komatsu for measurements of sand permeabilities and capillary pressure–fluid saturation curves, Dr. K. Glover, Dr. D. Rodriguez, and J. Gago for training and assistance with the x-ray attenuation measurement system, and Dr. M. Mathew for assistance with flume packing and modeling.


    REFERENCES
 TOP
 ABSTRACT
 INTRODUCTION
 Theory and Methods
 Results and Discussion
 Conclusions
 REFERENCES
 




This article has been cited by other articles:


Home page
Vadose Zone JHome page
F. Fagerlund, T.H. Illangasekare, and A. Niemi
Nonaqueous-Phase Liquid Infiltration and Immobilization in Heterogeneous Media: 2. Application to Stochastically Heterogeneous Formations
Vadose Zone J., August 1, 2007; 6(3): 483 - 495.
[Abstract] [Full Text] [PDF]


This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Web of Science (3) Citing Articles via Web of Science (3)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Fagerlund, F.
Right arrow Articles by Niemi, A.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Fagerlund, F.
Right arrow Articles by Niemi, A.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Fagerlund, F.
Right arrow Articles by Niemi, A.
Related Collections
Right arrow Multiphase Fluid Flow
Right arrow Experiment Design


HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
The SCI Journals Agronomy Journal Crop Science
Journal of Natural Resources
and Life Sciences Education
Soil Science Society of America Journal
Journal of Plant Registrations Journal of
Environmental Quality
The Plant Genome