Three simulation models specifically designed for use on mined and reclaimed lands are discussed. The first model attempts to assess the impact of mining and reclamation on site biomass productivity. The second model describes potential erosion at a mined site, draws most likely rill and gully pathways, and delineates contributing zones and depositional areas. The third and final model considers ordinary and bacterially catalyzed pyrite oxidation and leaching of acid and acid products from a reclaimed site. This site specific model can be programmed to evaluate effects of placement, fragment size, spoil bulk density and porosity as well as long- and short-term consequences of revegetation and use of limestone.
Much as the primary goal of any mining operation must remain the profitable recovery of a deposit, the need to restore land to previous or better use following mining cannot be overemphasized. Traditionally, exploration and exploitation of mineral deposits has been the province of mining engineers and geologists. Only recently have soil and plant scientists been called upon to make a contribution in this field. The nation has realized that it can no longer afford to exploit its resources without meaningful restoration. Such restoration constitutes a great challenge to ingenuity, because there may be an opportunity to create productive agricultural complexes where none existed before. Equally if not more important is the preservation of existing agricultural land by restoring it to its original productivity. Novel approaches to reclamation have demonstrated that a mined area can be handsomely reclaimed, revegetated, and subsequently used for a variety of purposes, ranging from farming and ranching to recreation and wildlife habitat. Reclaimed mined land can produce food and fiber for years to come despite a temporary mining disturbance. Although initially the new "soil" will need to be supplemented with fertilizer and organic matter, eventually a full-fledged soil profile will form, perhaps not unlike the original profile it replaced.
In the first section we discuss different approaches to topsoil handling based on comparative biomass productivity criteria subject to time and the existing climatic constraints (Rogowski and Weinrich, 1983). Although the reconstituted topsoil is not immediately like the material it replaces, it can be made into an adequate medium for plant growth through proper handling. With time, it should also develop its own particular layering (horizons) and characteristics.
Raindrops that fall in a storm can be both useful and damaging. While water is a necessary input to plant growth, it can also act as a primary source of energy for soil erosion especially in newly reclaimed humid areas. Search for the appropriate guidelines in identifying and evaluating the nature and extent of erosion and deposition on mined and reclaimed watersheds has recently been intensified. A properly posed erosion and deposition model can be used as a design tool to develop plans for control of erosion and siting of sediment ponds. In the second section, we discuss briefly a mathematical soil erosion and deposition model based on the fundamental mechanics of erosion (Khanbilvardi et al., 1983), and show how this model is applied to a large watershed, where data, except for those readily available from Soil Conservation Service or University Extension, are largely lacking.
Acid drainage from reclaimed coal-stripmines is a severe problem in the humid Eastern United States. Although the amount of pollutants leaving a mined site can be reduced by proper reclamation, the development and evaluation of best management techniques is difficult because of the great expense and time required to implement and assess the value of a specific practice. The final section describes a mathematical model of acid production and leaching from stripmined lands (Jaynes et al., 1983a,b,c,d) recently developed at our Center. This model is used to assess the potential of acid mine drainage at a given site and to test the probable effects of different management practices.
I. ASSESSMENT OF BIOMASS PRODUCTIVITY POTENTIAL
A taxonomic soil description may contain much pertinent information about a soil, yet it says very little in a quantitative way of how productive a given soil is prior to mining and how productive it will be following mining. Consequently, it gives us little or no information on how the soil should be handled for optimum results.
Depending on whether the area mined is in the East, Midwest, West, or even South, it would be desirable to show objectively how mining affects soil productivity. It would also be desirable to propose two or more alternate ways of topsoil handling which attempt to minimize the impact of mining operations. To do so we need to resort to modeling of soil productivity before and after a simulated mining operation.
Some currently available biomass productivity models approach the whole issue of productivity strictly from the perspective of climate. In Table 1, for example, biomass productivity (B) is calculated using mean annual values of temperature (Eq. 1) and precipitation (Eq. 2) in the so-called Miami Model (Lieth, 1975, P. 246),
The site productivity is taken as the lesser of temperature (B ) or precipitation (Bp) dependent values (these are the starred values in Table T Other models approach the productivity potential through soil effects alone. Kiniry et al. (1983) assumed that the productivity index was a function of certain physical and chemical soil properties which are or may become limiting. The primary soil properties of interest were bulk density, aeration porosity, available water, pH, and electrical conductivity. The model was weighted with depth by root distribution.
The Kiniry et al. (1983) approach has been modified here, and combined with Lieth's (1975) Miami Model for assessment of site productivity potential before and after mining. It will allow the user to choose which soil horizons he wishes to save and what soil properties he needs to maintain, improve, or amend in the reconstituted spoil. Coupled with the soil moisture budgeting procedure, the model can hopefully make a realistic appraisal of the site's reclamation potential. One must of course realize that values used in the computations of moisture budgets and productivity factors such as bulk density, pH, water availability--to mention just a few--as well as mean annual values of temperature and precipitation needed for the Miami Model of biomass productivity, are in fact distributions in space and time, spanning a range of values with different probabilities of occurrence. We should also bear in mind that for a given year or a given point within a site results may differ considerably from those predicted, although on the average, the model will probably be correct. Thus, this model should be used as a relative rather than an absolute predictor of productivity. In this sense it is particularly suitable for use in mining operations where the operator generally wants to know how pre- and post-mining conditions relate to one another.
Productivity Index
Relative productivity of a site (P) in g/m2 /yr can be written,
The undisturbed field soil profile usually consists of several distinct horizons in a solum of different depths. Similarly, plant root distributions vary. Some plants are shallow rooted, others have roots going down to considerable depths. However, in reconstituted topsoil, root distribution, at least initially, should be quite similar and reasonably constant for most profiles. To simplify the procedure we are considering here a 100-cm. deep soil profile, assuming that such a profile approximates an average realistic depth of reclaimed topsoiled profile which is available for plant growth. For comparison purposes the 100-cm depth will be used throughout, realizing that some topsoiled profiles may be shallower and others deeper than the value discussed here. Taylor and Terrell (1982) compiled information on rooting depth of over 200 plant species. Their data show that, on the average, 1.5 ± 1 meter is the depth reached by many roots, or their branches to which the absorption of water and nutrients takes place.
The structure of our model given in Eq. (4) is fairly simple. Productivity factors based on soil data are computed individually for each horizon by considering the minima, or at times the maxima required for plant growth. These values are then multiplied together to show where the problems may arise. The products are subsequently multiplied by an appropriate weighting factor for each horizon, approximating plant root distribution, summed, and multiplied by a site biomass productivity potential (Table 1) giving an overall site productivity index.
To compute the roots distribution weighting factor we used the method proposed by Horn (1971) and adapted by Kiniry et al. (1983). Horn (1971) computed root distribution from depletion measurements of water holding capacity in a 3.5 m soil profile under maple trees. Kiniry et al. (1983) assumed that the same approach can be used in soil profiles other than 3.5 m deep, provided that the rooting depth of growing plants was known. In here we postulate a 1 m rooting depth profile and assume that the reconstituted profile will have no physical or chemical barriers to root growth and no shallow water table. This appears to be a reasonable assumption for engineered profiles where we have considerable latitude in selection of desirable horizons and in handling the soil in such a manner as to minimize adverse effects of compaction.
The principal assumption of this model is that the biomass production is a function of root growth which in turn depends primarily on available water, bulk density, pH, and where applicable, on conductivity, sodium adsorption ratio, and nutrient status of the soil.
Equation (5) gives a fractional depletion of profile moisture (W) at a depth r,
where
R = rooting depth (cm)
r = depth of layer in the profile
When integrated with depth and divided by total profile depletion the equation approximates the root distribution weighting function (W).
Available water capacity (AWC) in mm should realistically describe the actual amount of water available for plant use during the growing season. This figure will vary depending on the climate, type of soil, and water use efficiency of a particular plant species. Realizing this we must choose an appropriate value of a critical content of soil water that will likely limit the plant growth under particular site conditions. Working in Missouri, Kiniry et al. (1983) chose a value of 0.20 cm/cm as the limiting available water content and represented PAWC (available water capacity productivity factor) as,
PAWC = AWC/0.20
with AWC values being always equal to or less than 0.20. Thus, for AWC '>' 0.20 cm/cm, PAWC = 1.0. For a 100-cm profile this would amount to 200 mm of available water. Table 2 shows available water capacity values (AWC) based on texture
Figure 1. Estimated distribution of roots with depth in a dimensionless profile.
Table 2. Potential available water capacity estimated from soil texture
(Neal, 1979; Peterson et al., 1968). It has also been customary to express AWC as the difference between water contents of 33 and 1500 kPa. Such values may be satisfactory for use in the model, provided they reflect actual soil conditions during a growing season from planting to harvest. In the humid East, where the annual percolation is usually in excess of 180 mm, this generally is the case; however, as we move further west climatic consideration may override textural properties and the precipitation may simply not be enough to fill the available storage capacity. We therefore suggest that AWC be chosen as the lesser of the values in Table 2 (or a value computed from a desorption curve at 33 and 1500 kPa), and the actual measured or estimated field values. Profile moisture contents at planting may be taken as starting values. The difference between field values at planting and the moisture content at 1500 kPa will give an approximate amount of water available to plants during a growing season for a reconstituted 1 m profile.
Both bulk density and aeration porosity are modified during topsoil handling. Assuming that reserves of moisture and nutrients are adequate, soil productivity depends largely on how well plant roots can access these reserves. It has long been established that compacted layers impede root growth by preventing root elongation, limiting respiration, and at times contributing to water logging. Soil bulk density (BD) is commonly used as an index of compaction. Its effects, however, should be evaluated relative to soil texture, moisture and moisture content at the time the soil is handled. Soil moisture content is particularly relevant to consider on reclaimed profiles where dense layers can often be produced if the topsoil is handled at or near so-called optimum moisture content. It is a well-known fact in soil mechanics that compaction at the optimum moisture content will result in maximum attainable values of bulk density for a given soil layer (Terzaghi and Peck, 1948). To account for these limitations, we have used in the model a bulk density productivity factor (D) suggested by Kiniry et al. (1983), Figure 2,
For certain specific conditions either lower or higher values of critical bulk density may be indicated, but in general BD suitable for plant growth will range from 1.30 to 1.80 g/cm3. A good discussion of the effects of bulk density is given by Pearson (1965) and again by Bowen (1981).
In the model presented here, most physical and hydraulic effects of bulk density are assumed to be incorporated into the factor D. other effects, such as adequate aeration of root growth medium so necessary for proper respiration and functioning of the roots, are grouped under productivity factor aeration porosity. The two are related through the equation,
Figure 2. Productivity factor bulk density.
where P a is the aeration porosity of a fully recharged profile, 6 is the moisture content at field capacity, and 2.65 is the particle density of a mined soil. Field capacity will vary with soil texture, structure, and the amount of organic matter present in the profile. In general, it should be somewhere between moisture contents measured at 10 and 33 kPa. Critical aeration porosity (Pcrit = 0.10) when root growth may become restricted ranges from 0.05 to 0.15 pore space by volume (Cannel and Jackson, 1971; Pearson, 1965). Realistically, aeration porosity effect should also include (but does not) a built-in dependence on time, a geometry factor to describe degree of continuity between air-filled pores and concentration level of CO2. In some of the mine spoils where root respiration may compete for pore oxygen with oxidation of pyrite or iron in the profile, or when heavy additions of organic material (such as sludge) place an additional demand on pore oxygen, it may well be advisable, particularly for deeper layers, to set Pcrit higher than the recommended value.
Productivity factor aeration porosity (S) was estimated in the model as an integral with depth of aeration porosity reciprocal divided into an integral with depth of Pcrit reciprocal
Using Eq. (8) it is of course possible to derive values of porosity at planting from field measured values of bulk density and moisture content. Workable estimates of Pa can be obtained by computing anticipated moisture storage in the spring for a given profile and by substituting its volumetric equivalent into Eq. (8). To be practical, however, an operator should ask himself if at any time during the growing season, particularly when soil is wet, or following a heavy incorporation of organic matter, the soil oxygen concentration is likely to drop near critical level for even a short period of time (Cannel and Jackson, 1981). In the event of water logging, for instance, factor S for a particular profile may become quite critical.
At times, however, we may wish to minimize organic matter conversion or nutrient leaching in the topsoil that is stockpiled before being spread. Under these conditions we may want to manipulate bulk density or aeration porosity so as to make values of D in Eq. (7) and S in Eq. (9) as small as possible.
Soil reaction (pH) values appear particularly well suited to characterize productivity response of reconstituted minesoil profiles in the East, but the approach proposed here is rather simplistic. Critical pH (Spurway, 1941) varies among soils and plant species, and with time. The response to pH on acid soils may result from H toxicity, Al toxicity, Mn toxicity, Ca deficiency, or Mo deficiency. Thus, soils at the same value of pH could have limited yields for different reasons and the limiting factors would operate at different intensities in time and space (Adams, 1981; Pearson, 1965). Consequently, the proposed model, which follows Neal's (1979) original formulation, should be used with caution and can be adjusted if sufficient information about a particular site or plant species is readily available. In the meantime, the model will provide sufficient guidance for the potential user in differentiating between the layers that may cause him problems and those that will not.
The pH productivity factor (pHF) is shown in Figure 3 and can be written,
where pH denotes a measured value of pH in the 1:1 aqueous solution.
An argument similar to that offered above for pH applies as well to the effects of salinity on plant growth (Hoffman, 1981). Salinity effects may vary spatially with soil type, texture and moisture status while different plant species will exhibit different tolerances. Salinity associated problems are almost certain to be present if reclamation is carried out in arid regions when original soils contain sufficient soluble salts derived either from marine deposits or from weathering. The problems can also arise in semiarid regions when rainfall is equal to evapotranspiration; as a result of upward artesian flow from aquifers; from overirrigation and from resulting saline seeps in adjacent areas; or if high water tables are present. Consequently, particular attention needs to be paid to the water regime and potential flow pathways in reconstituted profiles. Hoffman (1981) rates plants (his Table 9.3.1, p. 315) according to their salt tolerance, and suggests an appropriate form of the productivity factor, ECF
where A is salinity threshold value, B is the yield reduction per unit of salinity increase and EC is the electrical conductivity of soil saturation extract in ds/m. Selected values abstracted from Hoffman (1981) table are shown in our Table 3. Equation (11), which separates moderately sensitive (MS) and moderately tolerant (MT) plant species with A = 3.0 and B = 0.0769, can be written,
Figure 3. Productivity factor pH.
Table 3. Salt tolerance of some agricultural crops (Hoffman, 1981)
Cursory inspection of Table 3 will tell us that ECF can vary considerably depending on the crop used. Nevertheless, the model will alert the user that the salinity problem may exist if the profile is reclaimed in a similar form to original soil.
Soils or soil horizons that have excess sodium on their exchange complex are known as sodic soils or sodic horizons. Such soils when leached with low electrolyte content waters may show a marked decrease in permeability (Reeve and Bower, 1960; Frenkel et al., 1978). Many studies (Rhoades, 1982) have dealt with reclamation of sodic soils and with evaluation of the irrigation water quality (Rhoades, 1972; Oster and Rhoades, 1976). Guidelines, regarding the suitability of irrigation waters for agriculture, based on the type of predominant clay mineral, have been proposed (Ayers and Westcot, 1976) and questioned by more recent findings (Shainberg et al., 1981; Frenkel et al., 1978, and Suarez et al., 1983). Currently the general consensus appears to be that the permeability of sodic soils depends on the electrolyte level a soil maintains--substantial decreases having been observed with low electrolyte contents--(Shainberg et al., 1981; Rhoades, 1982), and may decrease with increasing pH (Suarez et al., 1983). Such variations in permeability would have a significant effect on infiltration and subsequently on the amount of water available to plants for a given soil.
Figure 4 from Rhoades (1983) summarizes the present situation for some of the more sensitive arid land soils and is used here to derive a productivity correction factor for high sodium soils. The SAR values in Figure 4 are the adjusted SARa values (Oster and Rhoades, 1976; Bower et al., 1968) in the topsoil and the electrical conductivity is that of infiltrating water, here assumed to be in equilibrium with the soil solution (ECe). Accordingly, the curve in Figure 4 is broken into two straight line segments at SARa ~ 10 and ECe = 1,
On the average, Frenkel et al. (1978) have observed an 83 percent reduction in permeability of montmorillonitic, kaolinitic and vermiculitic soils when
Figure 4. Threshold values of adjusted sodium adsorption ratio of topsoil and electrical conductivity of infiltrating water (assumed to be in equilibrium with soil solution) for maintenance of soil permeability (from Rhoades, 1982).
leached with distilled water as compared to 1N NaCl-CaCl solution. In here we have assumed a productivity reduction factor (SARF) proportional to the ECe/EC ratio, where EC is the electrical conductivity of the soil,
SARF = 1, for any of its values greater or equal to 1 and SARF = 0 for any of its values less or equal to 0. At this stage no pH dependence was incorporated, and the program subroutine is only activated when a problem is thought to exist.
The adjusted SARa values are calculated following the procedure outlined in Ayers and Westcot (1976),
pHc --a theoretical pH of water in contact with lime and in equilibrium with soil CO2 -- is computed from Eqs. (20), (21) and (22) fitted to curves in Figure 5. The input required is the concentration in meq/1 of Ca, Mg, Na, CO3 and HCO3 in the soil solution extract.
The computations outlined here will help to identify the problem if one exists. In the event of attempted irrigation procedures given by Rhoades (1982) should be followed.
The user should by now have realized the great complexity of a soil system we are working with. We have attempted to deal with this complexity in two dimensions, by considering the distributions of the productivity index identified with a particular soil series. However, it is a well known fact in pedology (Jenny, 1980) that local relief, aspect, and drainage, are among the most significant modifiers of the soil profile, particularly on a small scale. Thus, south facing slopes may have a different vegetation, moisture or temperature regime from the north facing ones. Soils developed on the hilltops are likely to be shallow compared to those developed near the base. Properties such as cation exchange capacity or clay content may vary significantly with elevation and drainage (Figure 6) and available moisture may range from excessive to impaired depending on position relative to the slope and degree of profile anisotropy (Zaslavsky and Rogowski, 1969).
Considerations such as the ones outlined above should be incorporated into the productivity assessment primarily through a representative sampling of the area to be reclaimed and through use of parameter values that adequately reflect the area heterogeneity. This is of paramount importance when attempting to reclaim land in Appalachia and elsewhere where mountainous, or rapidly changing conditions prevail.
Table 4 illustrates a typical computer program output sheet for a soil, in this case a Mollisol, from Bowman County, North Dakota. The lower values in Table 4 contains input parameters such as bulk density, available water content, pH and electrical conductivity as well as the mean annual temperature, precipitation, and clay content in each horizon. A place for user's soil number, and a space for a three-digit soil mapping unit name although not utilized here are also available. In computing available soil water holding capacity, values higher than 0.20 cm/cm were set equal to 0.20. The actual water holding capacity for this profile was taken as the customary difference between water content at 33 and 1500 kPa. When a profile is or can be fully recharged in spring, the above procedure may be correct. Soils can retain water at tensions lower than 33 kPa (1/3 bar), and some plants can utilize water at tensions higher than 1500 kPa (15 bar). Thus, a good practical estimate of soil water holding capacity is the difference between water held at 0.9 saturation and 1500 kPa (15 bar). In areas with insufficient rainfall, such as Bowman County, North Dakota, example in Table 4, a more nearly correct estimate is to subtract, as is done in Table 5, water content at 1500 kPa from an estimated profile water content in the spring.
Figure 5. Plots of algorithms used for calculating the components (Eq. 19) of adjusted sodium adsorption ration, pK'2-pK'c (A), p(Ca+Mg) (B), and p(Alk) (C). 2
Figure 6. Distribution of clay and cation exchange capacity in a semiarid toposequence (from Nettleton et al., 1968), numbers 1, 2, 3 and 4 denote successive soil series which have developed at different elevations.
Table 5. Computations of available water capacity
Productivity factor values given in the upper part of Table 4 constitute the output which is pertinent to choosing an appropriate procedure for topsoil handling. The magnitude of the five productivity factors is listed by depth and horizon for this soil in columns 6 through 10. The factors range from 0 (critical value exceeded) to 1 (no soil limitations to root growth). Their product in column 11 shows the quality of each layer. Thus in Table 4, the bottom layer (>97 cm), a 30-to-53-cm layer, and a 71-to-97-cm layer show low product values. In the bottom layer critical density is exceeded (1.84 g/cm3) while the 30-to-53-cm layer exhibits high density throughout. These layers, because of relatively higher clay contents, may compact too much during reclamation and should be handled with care. Even a more serious condition are the high electrical conductivity values which reduce the EC productivity factor in all layers below 30 cm.
The analysis suggests that at this site the 0 to 30 cm layer is best for plant growth, therefore it should be segregated and handled with care.
Column 12 in Table 4 gives the estimated profile root distribution which when multiplied by the values in column 11 (product) and summed gives the cumulative profile productivity index. The total profile productivity value of 0.6960 when multiplied by the Biomass Productivity (Table 1) gives the profile Biomass Productivity Factor of 467 g/m2/yr.
Because of the limited recharge at this site only part (84/216) of the available water capacity will be filled. Consequently, available water should be reduced by that amount (Table 5 column marked 3). The profile productivity will then be reduced in the ratio of column (3) to column (2) as given by column (4). When this correction is applied to the column 11 (Table 4) and subsequent steps are carried out as before, the Profile Biomass Productivity drops to 236 g/m2/yr.
Table 6 summarizes the output of the Comparative Biomass Productivity Model for some typical profiles listed in Table 1. The input information for these soils was chosen from states and counties near major mining areas. The results illustrate the likely range of values, as well as principal site limitations. During the reclamation process some of these limitations can be overcome to produce a better or at least the same kind of a site as was there originally.
II. EROSION AND DEPOSITION ON LARGE WATERSHEDS
In the Khanbilvardi et al (1983a) model the erosion process is partitioned into rill and interrill components. Soil particles detached by rainfall impact on the interrill areas are assumed to be transported by sheet flow into the rill microchannels. The amount and extent of rill erosion depends on a balance between rill flow detachment capacity and rill flow transport capacity.
The large 600 ha watershed considered here as an example is first divided into a grid of 100 x 100 m subareas each represented by a node in the center. These 1 ha subareas are assumed to act as homogeneous subwatersheds, and the parameters at each node are assumed constant for the whole subarea. While this assumption certainly holds for small areas of few square meters, it may also hold for larger areas. Rogowski et al. (1983) have also shown that most USLE parameters are correlated for points within 400 m of each other. The values of infiltration appear more variable (Rogowski, 1980), but can be assumed correlated for points within 100 m of each other.
Table 6. Biomass productivity and limitations for selected locations near or in mining areas
To estimate the amount of runoff available for a storm, the model first computes cumulative infiltration at each node from Philip's (1957) infiltration equation. Subsequently, the model computes rainfall excess at each node and checks for possible rill sources. Rill patterns which start at the rill sources and end at the area boundaries, are generated by computer, based on node elevation and the steepness and direction of slope. Laminar sheet flow regime is assumed on interrill areas. If an area is identified as a contributing area, this flow can contribute soil to rill microchannels. Soil loss on a contributing area is computed in the model by the Universal Soil Loss Equation (USLE) (Wischmeier and Smith, 1978). At each node m the soil transported by the rill flow up to that point plus the soil eroded from interrill contributing area are compared to rill flow transport capacity to determine the net amount of soil deposited, or carried away from node m. The computer then examines the new adjusted rill flow transport capacity and computes the amount of soil to be detached by scour by the time the rill flow reaches the following node (m + 1). The routing procedure beginning at each of the rill sources and ending at a deposition basin or at one of several rill outlets computes the net amount of erosion and/or deposition at each node. The model can estimate the pattern of sediment movement as well as the amount of soil lost in an area under study.
The general theory of infiltration summarized by Philip (1969) defines cumulative infiltration (I) as,
where S is the sorptivity (cm/sec 1/2 ), t is time (sec), and A is a parameter depending on soil water content (cm/sec). These characteristic properties of the soil material can be determined experimentally in the laboratory or in the field (Rogowski, 1980). However, when characterizing a large watershed with many different soil series estimates of S and A derived from a suitable mathematical model appear especially advantageous.
In the erosion-deposition model presented here sorptivity S was computed from Parlange's (1972) approximation written as (Engman and Rogowski, 1974),
where:
The procedure requires a soil moisture characteristic curve, hydraulic conductivity and the value of initial water content. When reliable experimental curves of moisture characteristic are not available, reasonable estimates can be obtained (Rogowski, 1971) from moisture content and pressure at 1500 kPa, at air-entry (Bouwer, 1966), and at saturation. The moisture content (215 ) values at 1500 kPa (R15) are generally well known for many soils; when tensiometer pressure at saturation (RO) is set equal to zero, saturation moisture content (2o) can be taken as equal to total pore space and the values of moisture content (2e and pressure at air-entry Re) can be estimated. For practical purposes 2e = 0.90 2o~ and in general it is likely to fluctuate between 0.80 20 and 20 (Rogowski, 1971), while the magnitude of Re can be found by assuming the existence of a linearly declining moisture content between 33 kPa and saturation (Figure 7).
The hydraulic conductivity (K) written as (Rogowski, 1972) is,
Figure 7. Schematic representation of soil moisture characteristic.
Appropriate substitutions in Eq. (24) provided we know the initial water content 2i will now give Philip's sorptivity (S) values. Philip (1969) also stated that A-values in Eq. (23) vary between 0.38 Ks and 0.66 Ks, where KS = saturated conductivity. Thus, A-values can be estimated from Ks values by assuming that A equals the conductivity Ke (Ke = 1/2 Ks) (Rogowski, 1972).
The above methodology was incorporated into the original model, enabling us to determine S and A-values in Philip's infiltration equation from readily available values of water content at 33 and 1500 kPa, soil bulk density (used as estimator of total porosity), and soil permeability (assumed equal to Ks).
Largely because of the subarea size (1 ha) and the type of rain simulated, the distinction between the laminar and turbulent flow regime used previously as criteria for selection of rill sources no longer applied. Intuitively, the selection of an appropriate rill source depends not only on rainfall properties, but is also influenced by infiltration. Both parameters were used in the adjusted model to compute excess rainfall at a point and are expressed as a runoff/ infiltration ratio. Two other important factors are soil erodibility (KE) and the presence of positive slope. Thus, some type of a general relationship between runoff/infiltration ratio and soil erodibility suggests itself (Figure 8). This relationship could also be used to represent the conditions necessary for initiating a rill. Specifically, if a point representing the runoff/infiltration ratio and soil erodibility (KE) for a subarea falls above the curve in Figure 8 the rill will start to develop on that subarea, if not, there will be no rill initiation.
To illustrate the model applicability we have selected as a site a 600 ha area in Central Pennsylvania near a village of Pine Glen which is currently being stripmined. In the model the watershed is represented by a 20 x 30 node regular square grid of 100 x 100 m subareas (1 ha each).
Assuming that the whole area has just been mined and reclaimed the C and P factors were set equal to 1.0. Values of soil erodibility (KE) were entered into the computer, based on the soil name. The storm EI value was calculated and soil names and slopes were entered from soils map while elevations and slope lengths (upslope from the grid node) were estimated at each grid point from a 7.5' topographic map for each subarea. Potential soil loss for each of the 600 subareas was computed with the Universal Soil Loss Equation (USLE) (Wischmeier and Smith, 1978).
To compute erosion or deposition an estimate runoff and infiltration was needed and to do so the knowledge of initial moisture content (2i) was necessary. The initial moisture content was estimated (Engman and Rogowski, 1974) from values measured in previous years on a similar area situated close to the experimental site. Those measurements were compiled, ranked in ascending order and numbered.
Figure 8. Suggested relationship between runoff/infiltration ratio and soil erodibility.
n = 1,2,3……………., N
The probability of occurrence corresponding to each value was computed as,
Figure 9 shows the nearly identical probability distributions of moisture content for the three months in 1979 and in 1980. A 50 percent probability value in July was selected as the initial water content ei. for use in the model.
The model is executed for a 75 mm-6 hr storm event. After determining the cumulative infiltration on each subarea, the program checks in Figure 8 if the combination of runoff/infiltration ratio and soil erodibility is sufficient for a node to become a rill source. The model then generates a pattern of rills from source nodes to nodes at the boundary (numbers 1 to 4 in Figure 10) or to incipient ponding areas within a watershed (letters a to g in Figure 10), while simultaneously delineating the sediment contributing interrill areas. In our example, because of the high intensity and long duration of rainfall, the rills and contributing areas cover almost the entire watershed.
Using the values of potential soil loss at each node, the model computes for all rills the total amount of soil available for transport at each rill node. This amount includes the soil transported to the rill by sheet flow from the interrill areas and the soil detached in the rill itself by scour. The program then considers the transport capacity of the rill flow with reference to the load carried, with reference to the amount of soil eroded from adjacent contributing area, and with reference to the potential rill scour before the flow arrives at the next node. By comparing the transport capacity of rill flow with the total amount of soil available for transport the program determines net erosion or deposition at each node. Figure 11 gives the distribution mosaic of net erosion and deposition over the area. This shows where the zones of largest erosion or deposition are likely to be. Such information suggests optimal areas to install conservation structures or to locate holding ponds.
Although this model was developed to predict soil loss from watersheds, it can equally well be used to model the transport of sediment sorbed pollutants from the same sites.
III. ACID MINE DRAINAGE ON RECLAIMED AREAS
The chemistry of pyrite oxidation is extremely complex and only partially understood. At least two mechanisms for pyrite oxidation are possible. One possibility is that oxygen can react directly with pyrite to form sulfate and acid,
Figure 9. Distribution of relative moisture contents on minesoil 10 km south of the experimental site in 1979 and 1980.
Figure 10. Distribution patterns of rills, of contributing areas (grey), of sinks (a ... g) and outlets (1 ... 4), at.the 3 x 2 km experimental site.
Figure 11. Distribution mosaic of soil loss (mm) at the 3 x 2 km experimental
stripmined site.
Alternatively, ferric iron can replace oxygen as the direct oxidant,
It has been suggested that in stripmine spoil, the only important source of ferric iron for (28) is the in situ oxidation of ferrous iron (Lau et al., 1970; Singer and Stumm, 1968),
While the oxidation of ferrous iron is thermodynamically favorable, the kinetics are extremely slow at the normal pH's of stripmine waters (pH < 4); (Singer and Stumm, 1970). However, certain chemoautotrophic bacteria (Thiobacillus ferrooxidans are known to use the energy released by Eq. (29) as their energy source and can significantly increase the oxidation rate (Lundgren, 1975; Beck, 1960). Thus, Eq. (29) and consequently Eq. (28) are thought to be bacterially catalyzed since ferric oxidation of pyrite is significant only when bacteria are active.
As a final step in the pyrite oxidation process, the ferric iron produced in Eq. 29) may precipitate as some form of ferric hydroxide, which we'll represent as,
Precipitation often takes place after the iron has been leached from the stripmine site, with the iron hydroxide precipitating in surface streams and ponds. Summing Eq. (28) and Eq. (29) results in Eq. (27), thus regardless of mechanism the result is the same; two moles of acid are produced for every mole of pyrite and three and a half moles of oxygen consumed.
Figure 12 represents a conceptual model used in developing the acid mine drainage model. The reclaimed-stripmine environment can be described as consisting primarily of coarse (>2 mm) fragments (Pedersen et al., 1980; Ciolkosz et al., 1977). Pyrite oxidation, represented by Eqs. (28) and (29), takes place at the surfaces of the pyrite grains contained within the coarse fragments (left side of Figure 12). Oxidation reactants and products diffuse through the fragment between the pyrite grains and the spoil solution surrounding the fragments. Oxidation of ferrous iron is primarily by bacteria located in the
Figure 12. Schematic diagram of the acid mine drainage model. Reactions are described in text. Ferric iron complexes are represented by Fe3+.
spoil solution. Pyrite oxidation, taking place within the coarse fragments, and iron oxidation by autotrophic bacteria are linked to each other by their interactions with the spoil solution. Oxygen diffusion occurs through the air-filled pores (upper left corner), with the oxygen in equilibrium with the spoil solution. Side reactions, that may influence acid mine drainage, also occur within the spoil solution. Complexation of ferric iron is assumed to occur within the spoil solution and helps modify the total ferric iron concentration. Precipitation of ferric hydroxide species removes ferric iron from the solution and produces H+ Reactions between H+ in solution and the spoil matrix or gangue may also take place; removing H+ from solution, raising the pH, and producing a variety of reaction products, HC (lower left corner). Finally, leaching of the soluble species by percolating water removes these species from the reclaimed-spoil profile.
Expressions for each of these processes are developed. The processes can be divided into those that are relatively slow and thus may serve as rate controlling steps and those that are relatively fast and can be treated as equilibrium reactions (Ohio State University Research Foundation, 1971). In this model, ferrous iron oxidation, pyrite oxidation, and oxygen diffusion are assumed to be kinetically controlled, whereas iron complexation and precipitation or dissolution, along with acid-spoil reactions, and solution leaching are considered to be equilibrium processes.
To compute the oxygen balance within the spoil program solves the following nonlinear, nonhomogeneous, second-order partial differential equation in two independent variables,
D02 is the effective diffusion coefficient of O2 in an O2, CO2, N2 atmosphere and is found from the binary diffusion coefficients and the molar flux ratios, r1 = NC/NO and r2 = NN/NO (Jaynes and Rogowski, 1983).
This equation is solved by rewriting it as an implicit finite difference equation, using the central difference for the space derivative and the backward difference for the time derivative. The nonlinearity is removed by estimating the 02 and C02 mole fractions and fluxes at the t + )t time and then calculating the diffusion coefficient for these values. The harmonic average was used to calculate the transmissibility term,
was calculated by using the geometric mean of the source term at t time and the estimated value at t + )t. The geometric mean was used to average over the time step of the source/sink reactions and yet prevent over consumption of O2 in a layer.
The Q term was calculated for the separate processes that contribute to it, such as autotrophic and heterotrophic bacterial respiration, pyrite oxidation and chemical oxidation of ferrous iron. The solution scheme involved estimating the mole fractions of O2 and CO2, the fluxes, and the concentration of the chemical constituents in solution. The Q term was then calculated based on these estimates and the diffusion equation was solved. The values from this solution were then used to solve the various source/sink reactions. Results for these reactions were adjusted for changes in iron complex concentrations and acid neutralization by the rock matrix and carbonates. The entire series of equations were then solved again using the new values for mole fractions and fluxes of O2 and CO2 the rate of leaching of the dissolved constituents, and the concentrations of Fe2+, Fe3+, H+ and SO42-. This procedure was repeated until convergence was reached, which completed one sweep of the time step. After convergence, the estimates of K B (first-order reaction rate coefficient for bacterial oxidation of ferrous iron sec-1) for all the layers were updated on the basis of the values from the previous sweep and the entire series was again solved cyclically until convergence. Up to 25 sweeps were performed for each time step or until convergence between sweeps of all the values was reached. This procedure was stable but involved considerable computer time. However, the procedure was necessary due to the extreme nonlinearity of the K B function and the complex interrelationships of all the variables needed to calculate Q. Large time steps are possible (>1/2 year) provided the gradients of the leaching constituents (Fe3+, Fe2+, H+, S042-, and acidity, HC) are not too steep.
After convergence at a new time, the program outputs the requested results and then proceeds to the next time step or stops if the end of the simulation has been reached and dumps to an interact file if computing time has expired.
The oxidation rate of pyrite is controlled in the model by the chemical reaction rates of the oxidation process and the diffusion rate of the reactants and products within the coarse spoil fragments. These relationships are expressed by Eqs. (32), (33), and (34).
Chemoautotrophic bacteria can accelerate the pyrite oxidation rate by increasing the ferric iron concentration. Their activity is considered to depend on the quality of their energy substrate and on the deviation of the spoil environment from ideal growing conditions. As environmental conditions deviate further from these ideal conditions, the energy requirements of the bacteria increase, causing their role in pyrite oxidation to decrease. These relationships expressed in the model are illustrated in Figure 13.
The computer program POLS was used in each simulation. A 10 meter deep, reclaimed profile, was divided into 20 equal horizontal layers, each layer was assumed to contain 75 percent (kg/kg-1), coarse fragments. The saturated fragments were 4-cm thick, had a pyrite content of 0.25 percent (kg/kg-l), with a 20 percent porosity. The bulk density of each layer was set at 1800 kg/m3, with the total pore space 70 percent water-filled.
At the beginning of each simulation the profile air-space contained 0.21 mole fraction O2 and 0.0003 CO2. The water in the profile contained 0.5 mmol/L Fe2+, 50 mmol/L sulfate, no ferric iron, and was at a pH of 5.0. No soluble amorphous ferric hydroxide was present in the spoil. Water infiltrating the surface of the profile was given a pH of 5.0 and contained no dissolved iron and 50 mmol/L sulfate. In simulations where iron-oxidizing bacteria were active, a 32 day "inoculation" period was assumed.
Figure 14, used here as one example of model application, shows that the oxidation of pyrite increased significantly with the increase in the effective diffusion coefficient. Run 5, without active bacteria, showed a significant increase over Run 3, with 40 percent of the pyrite consumed after 10,000 days. When the bacteria were allowed to be active (Run 6), a very large increase in pyrite oxidation occurred with 70 percent of the pyrite consumed after 10,000 days. This represents a 133 percent increase in pyrite oxidation over Run 3 and a 75 percent increase over Run 5.
In general, simulation results showed that in systems without iron-oxidizing bacteria, oxygen is the only important oxidizer of pyrite. Chemical oxidation of ferrous to ferric iron is not sufficiently fast to affect the overall oxidation rate and is not important in these systems. For the combination of factors used and expected to be found on reclaimed sites, the oxygen diffusion rate is the primary factor in controlling the rate of pyrite oxidation. The pyrite oxidation rate is decreased by a decreased air-filled porosity and increased diffusion path tortuosity, by increased size of the spoil fragments, and by burying the pyritic material at deeper depths. In zones where oxygen is not limiting, autotrophic activity can greatly increase the pyrite oxidation rate. Whether or not bacteria are important in these zones depends primarily on the solution pH. Only if the solution pH can be maintained in the range between reduced bacterial efficiency and reduced ferric iron solubility (2.0 < pH < 3.0) will the oxidation rate of pyrite be affected by bacterially produced ferric iron. The interaction between H+ , produced by pyrite oxidation, and the spoil appears to be crucial in establishing the pH of the spoil solution.
Acknowledgements: Contribution from the U.S. Department of Agriculture, Agricultural Research Service, in cooperation with the Civil Engineering Department, the Pennsylvania Agricultural Experiment Station, The Pennsylvania State University, University Park, PA 16802. Supported by the EPA-ARS Inter- agency Agreement Funds: EPA-IAG-D5-E763.
REFERENCES
Adams, F. 1981. Alleviating chemical toxicities p.269. In G. F. Arkin and H. M. Taylor (ed.) Modifying the root environment to reducing stress. ASAE Monograph No. 4, Am. Soc. Ag. Eng., 2950 Niles Rd., St. Joseph, MO.
Ayers, R. S. and D. W. Westcot. 1976. Water quality for agriculture. Irrigation and Drainage Paper No. 29, Food and Agriculture Organization of the United Nations, Rome, Italy.
Beck, J. V. 1960. A ferrous-ion-oxidizing bacteria. I. Isolation and some general physiological characteristics. J. Bacteriol. 79(4):502-509.
Belly, R. T. and T. D. Brock. 1974. Ecology of iron-oxidizing bacteria in pyrite materials associated with coal. J. Bacteriol. 117(2):726-732.
Bouwer, H. 1966. Rapid field measurement of air-entry value and hydraulic conductivity of soil as significant parameters in flow system analysis. Water Resour. Res. 2(4):729-738.
Bowen, H. D. 1981. Alleviating mechanical impedance. p. 2111. In G. F. Arkin and H. M. Taylor (ed.) Modifying the root environment to reducing stress. ASAE Monograph No. 4, Am. Soc. Ag. Eng., 2950 Niles Rd., St. Joseph, MO.
Bower, C. A., G. Ogata, and J. M. Tucker. 1968. Sodium hazard of irrigation waters as influenced by leaching fraction and by precipitation or solution of calcium carbonate. Soil Sci. 106:29-34.
Cannell, R. Q. and M. B. Jackson. 1981. Alleviating aeration stress. p. 191. In G. F. Arkin and H. M. Taylor (ed.) Modifying the root environment to reducing stress. ASAE Monograph No. 4, Am. Soc. Ag. Eng., 2950 Niles Rd., St. Joseph, MO.
Ciolkosz, E. J., R. L. Cunningham, G. W. Peterson, and R. Pennock, Jr. 1977. Characteristic interpretations and uses of Pennsylvania minesoils. Pa. Agric. Exp. Stn. Prog. Rep.
Engman, E. T. and A. S. Rogowski. 1974. A partial area model for strom flow synthesis. Water Resour. Res. 19(3):464-472.
Frenkel, H., J. 0. Goertzen, and J. D. Rhoades. 1978. Effects of clay type and content, exchangeable sodium percentage, and electrolyte concentration on clay dispersion and soil hydraulic conductivity. Soil Sci. Soc. Am. J. 42(l):32-39.
Hoffman, G. J. 1981. Alleviating salinity stress. p. 305-341. In G. F. Arkin and H. M. Taylor (ed.) Modifying the root environment to reducing stress. ASAE Monograph No. 4, Am. Soc. Ag. Eng., 2950 Niles Rd., St. Joseph, MO.
Horn, F. W. 1971. The prediction of amounts and depth distribution of water in a well-drained soil. M.S. Thesis, Univ. of Missouri, Columbia, MO.
Jaynes, D. B. and A. S. Rogowski. 1983a. Applicability of Fick's Law to gas diffusion. Soil Sci. Soc. Am. J. (In press).
Jaynes, D. B., A. S. Rogowski, H. B. Pionke, and E. L. Jacoby, Jr. 1983b. Spoil atmosphere and temperature in a reclaimed coal-stripmine. Soil Sci. (In press).
Jaynes, D. B., A. S. Rogowski, and H. B. Pionke. 1983c. Acid mine drainage from reclaimed coal-stripmines. I. Model description. Water Resour. Res. (In press).
Jaynes, D. B., H. B. Pionke, and A. S. Rogowski. 1983d. Acid mine drainage from reclaimed coal-stripmines: II. Simulation results of model. Water Resour. Res. (In press).
Jenny, H. 1980. The soil resource: ecological studies 37. Springer-Verlag, New York.
Khanbilvardi, R. M., A. S. Rogowski, and A. C. Miller. 1983a. Modeling upland erosion. Water Resour. Bull. Vol. 19(l).
Khanbilvardi, R. M., A. S. Rogowski, and A. C. Miller. 1983b. Predicting erosion and deposition on a stripmined and reclaimed area. Water Resour. Bull. (In press).
Kiniry, L. N., C. L. Scrivner, and M. E. Keener. 1983. A soil productivity index based upon predicted water depletion and root growth. Missouri Agric. Exp. Stn. Bull. 1051 (In press).
Lieth, H. 1975. Modeling the primary productivity of the world. In. H. Lieth and R. H. Whitiaker (ed.) Primary productivity of the biosphere. Springer-Verlag, New York.
Lundgren, D. G. 1975. Microbiological problems in strip mine areas: Relationship to the metabolism of Thiobacillus ferrooxidans. Ohio J. Sci. 75(6):280-287.
Malouf, E. E. and J. D. Prater. 1961. Role of bacteria in the alteration of sulfide minerals. J. Metals 13(5):353-356.
Neill, L. L. (L. N. Kiniry). 1979. An evaluation of soil productivity based on root growth and water depletion. Unpublished M.S. Thesis, Univ. of Missouri, Columbia, MO.
Nettleton, W. D., K. W. Flach, and G. Borst. 1968. A toposequence of soils in Tonalite grus in the southern California Peninsula Range. USDA, SCS Report 21, Washington DC.
Ohio State University Research Foundation. 1971. Acid mine drainage formation an abatement. abatement. Water Pollut. Control Res. Series Program 14010 FPR, 81 pp., Environ. Prot. Agency, Washington, DC.
Oster, J. D. and J. D. Rhoades. 1976. Various indices for evaluating the effective salinity and sodicity of irrigation waters. Proc. Int. Salinity Conference, Texas Tech. Univ., Lubbock, Texas, p. 1-14.
Parlange, J. Y. 1972. Analytical theory of water-movement in soils. Joint IAHR and ISSS Symposium on Fundamentals of Transport Phenomena in Power Media, Vol. 1, p. 222-236, Huddleston and Barney Lts., Woodstock, Ont.
Pearson, R. W. 1965. Soil environment and root development. p. 9511. In N. H. Pierre (ed.) Plant environment and efficient water use. Am. Soc. of Agronomy, 677 South Segoe Rd., Madison, WI.
Pedersen, T. A., A. S. Rogowski, and R. Pennock, Jr. 1980. Physical characteristics of some minesoils. Soil Sci. Soc. Am. J. 44(2):321-328.
Peterson, G. W., R. L. Cunningham, and R. P. Matelski. 1968. Moisture characteristics of Pennsylvania soils: Moisture retention as related to texture. Soil Sci. Soc. Am. Proc. 32:271-275.
Philip, J. R. 1957. The theory of infiltration: 4. Sorptivity and algebraic infiltration equation. Soil Sci. 84(3):257-264.
Philip, J. R. 1969. Theory of infiltration. p. 215-296. -In V. T. Chow (ed.) Advances in Hydroscience, Academic Press, New York.
Reeve, R. C. and C. A. Bower. 1960. Use of high salt waters as a flocculent and source of divalent cations for reclaiming sodic soils. Soil Sci. 90:139-144.
Rhoades, J. D. 1972. Quality of water for irrigation. Soil Sci. 113(4):277-284.
Rhoades, J. D. 1982. Reclamation and management of salt-deffected soils after drainage. Soil and Water Management Seminar, Nov. 29-Dec. 2, Lethbridge, Alberta, Canada.
Rogowski, A. S. 1971. Watershed physics: Model of the soil moisture characteristic. Water Resour. Res. 7(6):1575-1582.
Rogowski, A. S. 1972. Estimation of soil moisture characteristic and hydraulic conductivity: Comparison of models. Soil Sci. 114(6):423-429.
Rogowski, A. S. 1980. Hydrologic parameter distribution on a mine spoil. In Proceedings, Symposium on Watershed Management '80, July 21-23, Bo' Idaho, p. 764-780.
Rogowski, A. S., R. M. Khanbilvardi, and R. J. DeAngelis. 1983. Estimating erosion on a plot field and watershed. International Conference on Soil Erosion and Conservation, Honolulu, Hawaii, Jan. 10-22.
Rogowski, A. S. and B. E. Weinrich. 1983. A biomass productivity approach to topsoil handling. In Bruce Kennedy (ed.) Surface Mining, Edition II. Soc. of Mining Engineers, AIME, Littleton, Colorado (In press).
Shainberg, I., J. D. Rhoades, D. L. Suarez, and R. J. Prather. 1981. Effect of mineral weathering on clay dispersion and hydraulic conductivity of sodic soils. Soil Sci. Soc. Am. J. 45(2):287-291.
Silverman, M. P. and D. G. Lundgren. 1959. Studies of the chemoautotrophic iron bacterium Ferrobacillus ferrooxidans. II. Manometric studies. J. Bacteriol. 78(3):326-331.
Singer, P. C. and W. Stumm. 1968. Kinetics of the oxidation of ferrous iron. In Proceedings, Second Symposium on Coal Mine Drainage Research, Pittsburgh, Pennsylvania, May 14-15, pp. 12-34, Coal industry Advisory Comm. to Ohio River Valley Sanitation Comm.
Soil Survey Staff. 1975. Soil taxonomy. Agriculture Handbook No. 436, U.S. Govt. Printing Office, Washington, DC 20402.
Spurway, C. H. 1941. Soil reaction (pH) preferences of plants. Michigan Agric. Exp. Sta. Spec. Bull. 306.
Suarez, D. L., J. D. Rhoades, R. Lavado, and C. M. Grieve. 1983. Effect of pH on saturated hydraulic conductivity and soil dispersion.
Taylor, H. M. and E. E. Terrell. 1982. Rooting pattern and plant productivity. In Miloslav Rechligl, Jr. (ed.) Handbook of Agricultural Productivity, Vol 17. CRC Press, Inc., Boca Raton, FL.
Terzaghi, V. and R. B. Peck. 1948. Soil mechanics in engineering practice. p. 5611. John Wiley and Sons, New York.
Wischmeier, W. H. and D. D. Smith. 1978. Predicting rainfall erosion losses - A guide to conservation planning. Science and Education Administration, U.S. Department of Agriculture in cooperation with Purdue Agricultural Experiment Station, Agriculture Handbook No. 537, p. 1-58.
Zaslavsky, D. and A. S. Rogowski. 1969. Hydrologic and morphologic implications of anisotropy and infiltration in soil profile development. Soil Sci. Soc. Am. Proc. 33(4):594-599.