INTEGRATED RESPONSE P L O T DESIGNS FOR INDICATORS OF DESERTIFICATION *

R. O. KUEHL*, R. P. BRECKENRIDGE** and M. PANDA* *The University of Arizona, Tucson, AZ 85721 and Lockheed Environmental Systems and Technologies, Las Vegas, Nevada 89119, USA; **Idaho National Engineering Laboratory, P O. Box 1625, Idaho Falls, Idaho 83415-2213, USA

Abstract. The improvement of land management practices on lands susceptible to desertification requires information on the status and condition of the existing resources as well as any change occurring in the resource condition over time. The Environmental Monitoring and Assessment Program (EMAP) of the U.S. Environmental Protection Agency has developed a statistical survey design for monitoring the condition of ecological resources on large spatial scales. EMAP-Rangelands used a uniformity sampling study in 1993 to evaluate response plot designs for three categories of indicators (soils, vegetation, and spectral reflectance) to be used for monitoring ecological condition of a site. The response plot design study was developed to integrate on-site measurements for the three indicator categories. The study was conducted on the Colorado Plateau in southern Utah in three rangeland resource classes (grassland, desertscrub, and conifer woodland) of differing productivity levels in an attempt to develop a common plot design for all three resource classes. Basic measurement units were developed to facilitate integration of data collection. Preliminary spatial analysis of the sampling study found considerable differences in variation patterns among the study sites and measurement categories for the indicator classes used by EMAP-Rangelands. Evidence of substantial trends in the indicator measurements on monitoring sites relative to regional trends leads to the conclusion that nonstationary spatial models for biological processes on a monitoring site may be needed to fulfill the requirements for developing plot designs and indicator criteria.

1. Introduction

Occurring in arid and semiarid areas, the process of desertification is the result of human activities or climatic change (Hellden, 1991; Primack, 1993). One of the main difficulties in assessing and trying to control desertification is trying to determine the threatened areas before they reach a threshold beyond which huge amounts of resources are needed to correct the problem area. This problem is compounded by the different schools of thought within rangeland science. Should these areas be managed following a more traditional range succession-retrogression model prompted by Clementsian ecology or following a state and transition model that looks at changes resulting from excessive erosion conditions or from major changes in weather such as drought (Westoby et al., 1989). To address these concerns and collect data for assessing issues related to desertification, assessment * The U.S. Environmental Protection Agency, through the Office of Research and Development, funded the research described here. This paper has been subjected to the Agency's peer and administrative review and has been approved as an EPA publication. The U.S. Government has the right to retain a nonexclusive, royalty-free license in and to any copyright covering this article. Environmental Monitoring and Assessment 37:189-209, 1995. (~) 1995 Kluwer Academic Publishers. Printed in the Netherlands.

190

R.O. KUEHL, R. R BRECKENRIDGE AND M. PANDA

and survey methods are needed that accurately characterize the resource and can give managers information in a timely manner. The assessments need to be made at a scale large enough to capture the variability of the site and represent the condition of the resource within similar ecological classes of the region. This paper addresses the variability of arid resources in the context of regional monitoring. The plot design study reported here is important to desertification because without understanding the patterns and site-specific processes that govern change at a site, scientists cannot make associations between human- or climate-induced change and select management options appropriately. Surveys are a good way to estimate the status of natural resources. Traditionally surveys have been applied to resources where commodities were removed to estimate the amount of merchantable wood volume by tree species and diameter class (Schreuder and Czaplewski, 1993), to evaluate crop production, or to gain data on the condition of water or soil resources. In the late 1980s the Environmental Protection Agency (EPA) and collaborating federal agencies, universities, research institutes, and private organizations established the Environmental Monitoring and Assessment Program (EMAP) to survey the status and determine trends in the condition of the ecological resources in the United States. EMAP was divided into resource groups to focus the science on different areas. Indicators of resource condition were tailored to the ecological unit and to resources being surveyed (e.g. fish harvest in lakes, vegetation structure and composition on rangelands, or soils quality for forests). The primary objective of EMAP is to provide information, with known confidence, on status and trends in the condition of the Nation's ecological resources as well as estimates, with known confidence, of their extent. EMAP developed a statistical sampling design for the survey of natural resources as a basis for statistical inferences about the condition and extent of ecological resources. The design furnishes a probability sample of ecological resources that forms the foundation for statistical estimation with known confidence levels and the framework to detect changes and trends in ecological indicators. The monitoring plan and survey design are crafted with regional or national assessment as goals (e.g. making statements about changes in desertification or regional shifts in vegetation communities). Thus monitoring and assessment are on a much different scale than that required to address issues associated with localized perturbations of the environment. The EMAP-Rangelands resource group (EMAP-Rangelands) used a uniformity sampling study to evaluate response plot designs for a group of preselected indicators, spectral .properties, soil properties, and vegetation composition and abundance. These indicator categories of arid ecosystem condition were previously tested during the summer of 1992 in the southeastern Utah portion of the Colorado Plateau (Kepner et al., 1994). These indicator categories were selected through workshops and peer reviews (Breckenridge et al., 1993, in press; Mouat et al., 1992; Kepner and Fox, 1991) and appear to meet all the criteria for indicator development as suggested by Hunsaker and Carpenter (1990), such as being

INTEGRATED RESPONSE PLOT DESIGNS FOR INDICATORS OF DESERTIFICATION

191

applicable and interpretable on a regional scale. Thus these indicator categories are of high priority for EMAP-Rangelands evaluation. Although they appear to demonstrate the highest potential or capability for diagnosing ecosystem change (i.e. the ability to be merged with other data sets to make integrated assessments of ecosystem condition at the regional level), they must first be further field tested for confirmation of regional diagnostic ability before their incorporation into long-term implementation.

2. Study Objectives and Rationale The objective of the 1993 EMAP-Rangelands pilot study was to determine the optimum sampling support area and plot size for measuring the indicator attributes for vegetation, soils, and spectral properties in EMAP-Rangeland's extensive resource classes of the Colorado Plateau biogeographical region, where the sampling support area is that area at the monitoring site that provides an adequate sample representation of EMAP-Rangelands subpopulations. The extensive EMAP-Rangelands resource classes included in the 1993 study were the desertscrub, grassland, and woodland resource classes in the Colorado Plateau biogeographical region (Kepner et al., 1993). The development of indicators of the condition of rangeland ecosystems depends on the quality of measurements needed to quantify the indicators. Effective ecological interpretation of information gathered on monitored sites requires measurements that adequately describe the biological communities at the sites. To make these measurements consistent with the biology of the EMAP-Rangelands resource class requires a sample plot configuration, within a sampling support area, that is sufficient to capture the characteristics (i.e. degree of variability) of the site's biological communities. A distinction must be made between the sample plot design and the sample survey design that will be used for the EMAP-Rangelands monitoring survey. The monitoring sites for EMAP-Rangelands are point samples in the context of the EMAP survey design (Overton et al., 1990; White et al., 1992). Point samples for measurements such as temperature do not need a support area around the point to acquire the measurement. But point samples for characteristics such as vegetation cover and soil properties require observation in some area surrounding the point to describe the characteristic identified by the point sample. The plot design at a site may be considered a "response" design for indicator measurement at a monitoring site, whereas the collection of monitoring sites constitutes the survey design. The measurements on the sample plot at a site will be used as one observation of the indicator in the survey design. The survey sample is complete only when that indicator observation is included with the observations collected at all other sites combined by resource class within the landscape or

192

R.O. KUEHL, R. E BRECKENRIDGE AND M. PANDA

region. The observations from the survey design are those used for regional and landscape interpretations. Interpretation of ecological condition at EMAP-Rangelands monitoring sites requires an adequate representation of the vegetation and soil communities at a site. The ecological concept of minimal area (Dietvorst et al., 1982) is usually defined as the smallest area on which the species composition of a plant community is adequately represented. Two renderings of the minimum area concept have been described, the biological and methodological. Biological minimum area refers to the smallest area on which the species composition of a plant community is adequately represented or, equivalently, to the size of stand required to be well developed. Methodological minimum area within a stand refers to the size of sample plots needed to complete at least an adequate description of that stand (Barkman, 1989). The methodological minimum area has been further divided into qualitative and quantitative parts. The qualitative area is that plot size above which the number of species does not increase at all or increases only insignificantly within the same stand. The quantitative area is distinguished as that above which the quantitative shares of all species do not change significantly (Barkman, 1989). For example, Barkman (1989) observed that the gain in information is less each time the sample plot area is doubled, particularly with respect to the cost of sampling. Thus EMAPRangelands concluded that a pragmatic definition of methodological minimum area should be that plot size whose further enlargement produces an insufficient gain in information for the added cost. Researchers have used several criteria to determine the minimum plot size and shape needed to adequately describe an area. Species area curves, similarity analysis, frequency area curves, species representation, and pattern representation have had various applications. The Braun-Blanquet cover abundance scale (Bonham, 1989) is commonly used as a measurement to evaluate the minimum area. Historically, monotonic relationships have been observed between the size of a sampled area and the criteria used to determine the methodological minimum area such as species area curves and similarity measurements (Dietvorst et al., 1982; Barbour et al., 1980). Monotonic relationships between size of sampled areas and variances of measured variables have also been observed in a variety of settings, including biomass of grasses and forbs (Wiegert, 1962), tree volume (Tardif, 1965), basal area (Bormann, 1953), agricultural experimental yield trials (Smith, 1938), surveys for plant disease incidence (Proctor, 1985), and agricultural acreage and yield sample surveys (Cochran, 1977). One standard method to assess efficient plot sizes and configuration considers an empirical model developed by Smith (1938), which relates variance to size of plot. The method has been used for decades to evaluate plot sizes for agricultural field experiments and for natural systems (Brummer et al., 1994). A parameter estimated for the Smith model from the plot data is used in conjunction with a cost function to obtain an efficient size of sample plot. Extensions and modifications of

INTEGRATED RESPONSE PLOT DESIGNS FOR INDICATORS OF DESERTIFICATION

193

the Smith variance model include terms to account for variation between plants, (Pearce, 1955; Freeman, 1963) and plot shape, (Reddy and Chetty, 1982). The use of spatial correlation analysis in place of or in conjunction with the Smith model by Modjeska and Rawlings (1983) included a generalized method of plot optimization. Brewer and Mead (1986) extensively reviewed models that have been developed by themselves and others for analyses of field crop experiments to evaluate the efficiencies of plot designs. Other methods to determine optimum plot size include the tangent to the species area curve used by Hyder et al. (1963) for plant frequency data in sagebrushbunchgrass communities and sample size requirements for specific degrees of precision for confidence interval estimation with cost constraints (Bormann, 1953; O'Regan and Arvanitis, 1966). The measurement of variables for indicators on EMAP-Rangelands monitoring sites is intended to reflect the status of biotic and abiotic communities within the sampled resource community. Thus the sampling support area must be large enough to adequately characterize indicator variables for the plant and soil communities under consideration. The literature previously cited on plot size has been concerned primarily with vegetation communities. In addition, plot size in natural systems most often refers to quadrat for which the entire quadrat area is censused for vegetation community characteristics, such as biomass and species abundance or composition. The agricultural studies have concentrated on plot size studies for biomass, grain, or fruit yield in comparative experimental trials. The uniformity sampling trial is a commonly used method to provide data for several types of studies related to sampling area. Uniformity trials for agricultural crops use harvest data from a field planted to a single cultivar divided into small contiguous areas known as basic units. These data are used to find optimum plot sizes from among aggregations of the basic units. Contiguous rows of quadrants are used in ecological studies to evaluate the occurrence of contagious patterns. Increasing quadrat sizes are built up by blocking adjacent quadrants in pairs, fours, eights, and so forth, and the variances of different block sizes are related to block size in a graph to detect different scales as peaks in the graph (Greig-Smith, 1952). Contiguous quadrants are measured over a designated area in ecological studies to determine biological minimum areas (Barkman, 1989). The contiguous quadrants are combined into larger quadrants of different sizes, and various measures are calculated from the nested sets of quadrants to ascertain biological and methodological minimum areas. The EMAP-Rangelands survey requirements for sampling plots differed from previous studies in several facets. The sampling plots have to serve as monitoring sites for a suite of indicator measurements with the potential for repeated visits to the site over decades. A complete census of the plots is neither feasible nor desirable under these circumstances. Plots must be set up with minimal disturbance

194

R.O. KUEHL, R. E BRECKENRIDGE AND M. PANDA

~a

~z

[]

[]

[][][]

~w

~m

Fig. 1. MacroplotDesign for soils, vegetation,and spectral propertiesindicators.

to the measurement areas. The number of plots or replicate samples at the site for monitoring purposes will be the number required for specified confidence estimates in the survey. Thus questions of replication blocking with different treatments on the plots within blocks as in experimental trials are not relevant to the survey. Finally, for indicator attributes in this study EMAP-Rangelands is concerned not only with vegetation but also with soils and spectral measurements, which provide a much stronger preponderance of evidence for assessing desertification. Thus EMAP-Rangelands must construct a sample plot design that integrates the sampling support area and optimum plot size required for the measurements made for the three indicator categories.

INTEGRATED RESPONSE PLOT DESIGNS FOR INDICATORS OF DESERTIFICATION

3.

195

Approach

The field sampling areas required for this study were coordinated and integrated by three indicator groups. Large macroplots estimated to be larger than the indicator sampling support area were established in relatively homogeneous resource classes of either desertscrub, woodland, or grassland. The macroplots were the basis for a uniformity sampling trial at each site. Each macroplot (Figure 1) was subdivided into a rectangular array of basic sampling units. The size of igasic units was such that a sample measurement could be made for each of the selected measurements in the three indicator categories in each of the biomes. This study evaluated five combinations of resource class and productivity levels. Productivity was defined as the potential mean annual vascular vegetation biomass as described in existing ecological (range) site descriptions (USDA, 1976; USDI, 1990). Within a biogeographic province, vegetation community physiognomy and cover characteristics may differ significantly between naturally low production sites and high production sites. Physiognomy and cover characteristics, however, are expected to exhibit greater similarity between communities with similar production potential when subjected to similar environmental and anthropogenic stresses. Arid environments in the United States typically produce natural communities that range from a mean annual productivity of less than 178 kg/ha to more than 890 kg/ha because of the natural variance in the interrelationships between soil, climate, and vegetation community. Ecological site correlation procedures developed by the Bureau of Land Management (Leonard et al., 1992) suggest natural breaks at various productivity levels. Grouping of sites less than 445 kg/ha (low), 445 to 890 kg/ha (medium), and greater than 890 kg/ha (medium-high) are used for comparison. Also described are other criteria, such as dominant species and community composition, that relate to resource class characteristics at the larger scale considered by EMAP. Study sites for 1993 were selected to represent readily observable differences in resource class and productivity combinations, but not every possible combination was represented because of time and cost constraints. The desertscrub resource class is the most extensive vegetation class in the Great Basin biogeographic region. Low, medium, and (medium) high-producing sites were selected to represent this resource class. A medium-producing grassland and high-producing pinyon-juniper site were also selected. Low-producing grassland and pinyon-juniper sites are rare in the Great Basin biogeographic province. The size of the macroplot for the 1993 plot study (Figure 1) was selected to be more than large enough to provide an adequate sample thought to characterize the local scale of variation at the site. Literature and professional judgment were used to select the macroplots for soils, vegetation, and spectral indicators, but they all had to be of a size that could be supported with existing funding and staff. The final macroplot size (180 m x 180 m) was adjusted to maximize integration of data between the indicator groups and to ensure that all data could be collected within

196

R.O. KUEHL, R. R BRECKENRIDGE AND M. PANDA

1 week. The overall size of the soil macroplot was selected to be 180 m by 180 m to allow adequate soil samples to characterize variability. Previous research (Campbell, 1978; Warrick et al., 1986; Wilding and Drees, 1983) suggests that a spacing of 30 m would result in spatially independent samples. Thus the 180 m by 180 m macroplot contained 36 basic units (30 m by 30 m). A total of 36 other 10 m by 10 m basic units were nested within the 180 m by 180 m macroplot to have spatially dependent basic units in order to measure spatial correlations. Soil properties have different orders of variation and spatial correlation. To measure the variation and spatial correlation, observations and samples were taken at different densities. The 10 m spacing was chosen to estimate the spatial correlation of surface properties that vary more at a small scale and have higher coefficients of variation. These surface properties were organic carbon, total nitrogen, and hydraulic conductivity (Mausbach et al., 1980; Wilding and Drees, 1983). Properties that are variable at a larger scale and have coefficients of variation from 15-30% were sampled at 30 m intervals. These properties were the texture and structure of the A horizon. Properties that are less variable and have coefficients of variation of less than 15%, such as the soil color and thickness of the A horizon and soil classification, were sampled at 60 m intervals. Spectral reflectance properties of biotic and abiotic objects tend to be wavelength dependent, and determining the relationships and the support area is critical for characterizing or discriminating the object. The plot size selected for the spectral part of this study was based on the requirement to collect both ground-based and satellite spectral measurements to evaluate their relationship with estimated vegetation and soil features on the site. The spectral properties plot requirements were developed from findings of other researchers (Gholz, 1982; Waring et al., 1978) about the relationship between ecosystem structure (e.g. biomass, LAI) and functional properties (e.g. net primary productivity). Satellite sensors (Landsat Thematic Mapper) typically collect data that is integrated or mixed into small segments (pixels) representing the spectral properties of an area about 30 m by 30 m on the ground. The ground spectral data is collected using a portable spectrometer at a smaller scale, typically 30 cm when the instrument is held 1 m above the ground. Ground spectral measurements were spatially coincident with vegetation measurements (collected in 5 m by 5 m quadrats) and combined so that the relationship between vegetation composition and abundance and surface spectral attributes could be determined at the basic unit of 10 m by 10 m. Ground spectral measurements were collected for the entire macroplot of 120 m by 120 m so that several satellite pixels could be combined to develop a spectral signature for the entire macroplot. Vegetation composition, structure, and abundance are measured because of their importance as indicators of changes in arid and semiarid ecosystems. Several sampling options were considered, and a decision was made to use a modified Daubenmire cover class method within a quadrat along a line transect. Similar methods

INTEGRATED RESPONSE PLOT DESIGNS FOR INDICATORS OF DESERTIFICATION

197

are used in other inventory, monitoring, and evaluation programs conducted by the USDA Forest Service in the intermountain West (USDI, 1985; O'Brien and Van Hooser, 1983). For plant sampling, rectangular plots have been found to yield better results than other shapes. A rectangle with sides in a 1:2 ratio works well and was selected for grass and shrub quadrats. Tree quadrats were increased to 5 m by 5 m to account for the larger spacing between the pinyon-juniper plants and modified to a square to be integrated into overall macroplot design. Vegetation samples were coordinated with spectral and soils sampling so that relationships between the data could be evaluated. Vegetation line transects of 120 m were used, and plots were located every 5 m (two 1 m x 2 m quadrats were within each basic unit) to provide adequate coverage to characterize the variability of the more widely spaced organisms in arid ecosystems. Centered within the soil macroplot, the macroplot for vegetation and spectral properties was 120 m by 120 m. The macroplot size for each indicator group was selected to provide enough replication to supply adequate degrees of freedom to statistically evaluate the spatial variability of the resources across the site.

4. Results and Discussion

The observations for plot design analysis are spatially related whether they come from ecological studies of natural systems, agricultural experiments, or the current study on plot design for EMAP-Rangelands monitoring. The success of any method of analysis on spatially related data depends considerably on some assumptions about statistical properties of spatial data. More specifically, assumptions are required about the spatial process under study. These properties are discussed in detail by Ripley (1981), Haining (1990), and Cressie (1991). For the uniformity sampling study the data Y(sl), Y(S2) . . . . . Y(sn) collected at spatial locations si = {XliXZi} for i = 1,2,... ,n are considered observations from some biological process at the scale of the macroplot somewhat larger than 1 ha at any of the study locations. One of the desired features of the data for plot design analysis is that of stationarity in the spatial setting. Informally, stationarity refers to variational properties of the data that do not change throughout the sampled region of interest. More formally, a second order stationary process, Y(s), will have a constant mean, E[Y(s)] = # for all locations s in the region of interest, and a constant covariance function Coy[Z(8),

Y(t)] = e{z(s)

- .]{Y(t)

-

for the distance s to t anywhere in the region of interest for a constant direction. If the covariance depends only on the distance and not the direction, the process Y(s) is considered isotropic, a common assumption in many spatial statistics studies.

198

R.O. KUEHL, R. E BRECKENRIDGE AND M. PANDA

TABLE I Spatial correlation matrix for total grass-forb cover on a desertscrub site k 1

-5

-4

-3

-2

-1

0

1

0 -0.18 0.19 0.00 -0.08 0.18 1.00 0.18 1 -0.11 -0.06 0.22 0.08 -0.03 0.03 0.52 2 0 . 1 6 -0.25 -0.22 0.43 -0.02 -0.02 -0.01 3 -0.19 0.24 0.03 -0.09 0.00 0.40 0.11 4 0 . 0 4 0.09 0.01 0.14 0 . 0 0 -0.06 0.77 5 0.41 -0.14 -0.12 0.61 0.12 -0.02 -0.03

2

3

4

5

-0.08 -0.02 -0.06 -0.11 -0.01 0.02

0.00 -0.13 -0.04 -0.06 -0.16 -0.15

0.19 -0.08 -0.15 -0.01 -0.09 -0.08

-0.18 0.06 -0.12 -0.11 0.01 -0.09

TABLE II Standardized mean square errors for total grass-forb cover on a per-unit basis for plots a x b units in size on the grassland site b a

1

2

3

4

1

1.00

1.23

2 3 4 5 6

1.36 1.68 2.04 2.39 2.76

1.84 2.46 3.08 3.72 4.35

1.39 2.15 2.96 3.76 4.61 5.43

1.57 2.50 3.48 4.48 5.53 6.56

5 1.71 2.78 3.92 5.08 ' 6.29 7.47

6 1.77 2.89 4.11 5.34 6.62 7.87

A number of methods may be used to explore the stationarity or nonstationarity, as the case mail be, of the process, Y(s), with the spatially observed data. For the uniformity sampling trials, the exploration included estimates of spatial correlations, variograms, trend across the sites, and variances for a variety of plot configurations. But nonstationarity of the covariance and of the mean are difficult to jointly evaluate, and completely satisfactory methods are not fully developed (Handcock and Wallis, 1994). The initial analyses conducted on the 1993 plot data included the estimate of spatial correlations along with variances for sample plot configurations at the five sites for all of the indicator measurements. The computational algorithm for the covariances, correlations, and plot variances follow those of Modjeska and Rawlings (1983). Spatial correlations, r(l,k), are illustrated with total grass-forb cover measures on a desertscrub site (Table I) where r(l,k) is the correlation between observations

INTEGRATED RESPONSE PLOT DESIGNS FOR INDICATORS OF DESERTIFICATION

199

separated by l rows and k columns in the 12 × 12 array of data. Because of symmetry of correlations, r(l,k) = r(-l,-k) and r(l,-k) = r(-l,k). Under many of the common stationary models one would expect to see some patterns in the correlation matrix. But no regular patterns were apparent with these correlations, and in general we found no regular patterns for most of the indicator measurements and sites. Standardized mean squares, STMSE(ab), for sampling plots of size a rows by b columns computed for total grass-forb cover on the grassland macroplot (Table II) differed considerably among plot shapes. The mean square MSE(ab) for any configuration was divided by the mean square for the unit size plot, a = I and b = I, to achieve the standardized value. The relative efficiency of two plot shapes of the same size can be evaluated by their ratios. Consider plots of six basic units, for example, STMSE(61) = 2.76, STMSE(32) = 2.46, STMSE(23) = 2.15, and STMSE(16) = 1.77. The ratio of standardized mean squares for a 6 x 1 plot to 1 x 6 plot is 2.76/1.77 = 1.56, indicating the 1 row by 6 column plots are 56% more efficient than the 6 row by 1 column plots. Likewise, the efficiencies of 3 x 2 and 2 x 3 plots relative to the 6 x 1 are 1.12 and 1.29 respectively. These efficiencies show that plots were more efficient as they became elongated in the west-to-east direction. Upon viewing the first row and first column of Table II, one can see that any long and narrow plot in the north-to-south direction (i.e. 2 x 1,3 x 1. . . . . 6 x 1) is less efficient than a long and narrow plot in the west-east direction, (i.e. I x 2,1 x 3 . . . . . I x 6). Thus the plots elongated in the west-to-east direction captured considerably more variation on the site than those oriented in the northto-south direction. Similar differences in efficiencies can be seen for other plots of same size but different configurations. Again, similar results were observed for all measurements at all sites. Because nonstationary processes can produce trends in the spatial extent of the data, trends across rows and columns were evaluated for the data on the 12-row by 12-column sampling grid. Patterns of row effects and column effects in two-way arrays have commonly been analyzed with traditional analysis of variance methods. Underlying the two-way analysis of variance is the additive decomposition data = constant + row + column + residual and trends in the two-way array are revealed by the magnitudes and signs of the row and column effects. The row and column effects for the current study were evaluated with a median polish (Tukey, 1977) on the two-way array. The median polish resists the influence of outliers and reduces their effect on the evaluation of row and column contributions to the variation in a two-way table. The row effects from the median polish for total vascular plant cover in the north-to-south direction on a desertscrub site (Figure 2) show definite trends, as does the plot of column effects in the west-to-east direction on the site (Figure 3). Other variables on the sites also exhibited row and column patterns.

200

R . O . K U E H L , R. P. B R E C K E N R I D G E A N D M. P A N D A

O

ID

o

o

It3 !

I

I

I

2

4

6

I

I

I

8

10

12

l'OW

Fig. 2. Row effectsfor total vascularplant cover.

The stability o f the mean and variance on the macroplots can be evaluated with a moving window similar to the moving average technique used for time series. The mean and relative variation, percent coefficient of variation (%CV), for the normalized difference vegetation index (NDVI) was computed in 2 × 2 moving windows of basic units across the macroplot. The moving average for NDVI on a desertscrub plot (Figure 4) appears to be fairly stable, but the %CV of NDVI measurements (Figure 5) are quite divergent across the macroplot, suggesting the possibility of a nonstationary variance structure. For all variables at all sites these nonstationarity diagnostics for the most part exhibited patterns that discourage the use of the simple stationary spatial models most often used for plot design studies. In addition, the patterns or behavior of the diagnostics showed no consistency among sites or among variables. If the response plot design were of concern only for a specific site, the nonstationarity is less troublesome. The orientation of transects or shapes of sampling plots can be determined conveniently from the standardized mean squares and costs as suggested by Modjeska and Rawlings (1983). Of course, each site in that case would have its unique sampling requirements.

INTEGRATED RESPONSE PLOT DESIGNS FOR INDICATORS OF DESERTIFICATION

cq

-

o

e.j

!

I

I

I

I

I

I

2

4

6

8

l0

Fig. 3. Column effects for total vascular plant cover.

J

Fig. 4. Moving avarage for NDVI in a 2 × 2 window.

12

201

202

R. O. KUEHL, R. R BRECKENRIDGE AND M. PANDA

~,

~o~t~to sox~t~

Fig. 5. Moving coefficient of variation for NDVI in a 2 x 2 window.

TABLE III Averages of site means and coefficients of variation for basic units from uniformity sampling Variable

Mean

% Coefficient variation

Vegetation Cover (n = 144) Total vascular plant Total shrub Total grass-forb Dominant shrub Dominant grass-forb Total tree (1 site)

23.3 15.5 3.7 13.7 1.7 12.6

61 60 142 69 282 77

Spectral Reflectance (n = 144) NDVI SAVI Albedo Soils (n = 36) Organic carbon Nitrogen C/N ratio Water content

0.16 0.11 0.19

22 21 21

0.54 0.04 10.75 5.42

40 43 61 18

INTEGRATED RESPONSE PLOT DESIGNS FOR INDICATORS OF DESERTIFICATION

203

But the objective for EMAP-Rangelands requires developing a sampling strategy that applies generally to all monitoring sites. Thus the requirements for stationarity to evaluate the covariance structure across sites and variables become an important consideration for the plot design development. Historically high sampling variability has been experienced on rangeland sample sites, resulting in large sample-size requirements for precision in the estimates of means and other parameters. The average site means and coefficients of variation from the basic units at the five sites for some of the measurements (Table III) differ considerably in variation among the different measurements for indicator classes in indicating that sample size requirements will vary considerably. For some indicator measurements the needed sample information appears to increase with the detail required in the indicator. The variation for species-level information, such as dominant species cover, tends to be greater than that for total vegetation cover, thus requiring larger sample sizes for the same level of information. Likewise, carbonto-nitrogen ratios in soil samples tended to be more variable than the individual carbon or nitrogen samples. If biological processes are operating at a scale equal to or smaller than the sampled area in such a way to produce trends with local minima and maxima on the site, then the usual sample statistics are inadequate to characterize the site. The sample variances fail to separate the variation in trends across the site from the natural sampling variation about the trend. And the site means may be inadequate to describe the site characteristics. Thus these simple stationary statistical models are insufficient to characterize the range site. Under the supposition that biological and geological or soil processes are operating at the site scale, the data can be modeled as Y(s) = #(s) + Z(s) where #(s) is the deterministic mean structure or trend over the spatial extent of the sampled plot and Z(s) is the stationary correlated error process. The challenge with this broader view model for the range site is to achieve a sensible decomposition of the variation into its trend and random components. The component #(s) models the large scale variation on the site as a consequence of biological processes in operation at the scale of the sampled site. The Z(s) component models the random sampling errors about the mean structure, #(s). This representation is analogous to that used in time series analysis to model trend and error processes. Since the model recognizes spatial variability at the large scale through #(s) and at the small scale through Z(s), the objective for an analysis of the model is to estimate the trend on the site, #(s), and subsequently, from the estimates of the residuals, estimate the correlations and variances associated with the stationary process, Z(s). Interpolated surface plots of the 144 observations for total vascular plant cover (Figure 6), total grass-forb cover (Figure 7), and total shrub cover (Figure 8) on a desertscrub site reveal substantially irregular surfaces. The nonstationary mean model presumes that the roughness of the surface plot for total vascular plant

204

R.O. KUEHL, R. R BRECKENRIDGE AND M. PANDA

O

%

Fig. 6. Total vascular plant cover.

O

o

--;

6~3to so,~t~

~o Fig. 7. Total grass-forb cover.

and shrub cover and the obvious trend in the grass-forb cover results from more than simple random sampling variation. Specifically this model presumes that the surfaces can be explained as the aggregate of an unknown deterministic mean structure plus a correlated error Structure. Providing an informative description of the range site through the large-scale trend #(s) and the attendant correlated error

INTEGRATED RESPONSE PLOT DESIGNS FOR INDICATORS OF DESERTIFICATION

205

O

20 eo

z9

60o ~ogtkx xxo~t~

Fig. 8. Total shrub cover.

structure for Z(s) is a requirement placed on model estimation. The desired result is a parsimonious partition of the variation between the parameterization for #(s) and a suitably stationary error structure Z(s). Parameterization of the mean function #(s) and the subsequent estimates for the stationary error process are subjects of much current research in spatial statistics. Ripley (1981) and Cres sie (1991) extensively discuss some of the approaches. The challenge posed by this investigation requires a decision as to the need of nonstationary models for these local-scale monitoring plots when inferences are ultimately based on the results of these point samples over a much broader regional landscape. One consideration in making that decision is the range of variation experienced at the two scales. For example, the range in basic unit total vascular plant cover values was typically 2% to 70% at any one of the sites. The averages of the five sites ranged from 17% to 29%. Because the local scale of variation considerably exceeds the larger scale amount, the nonstationary model at the local site may be a reasonable choice. Another consideration in deciding to use a nonstationary model is whether the potential exists for local processes to occur at these monitoring sites. For example, the total vascular plant cover (Figure '6) is an aggregation of the total grass-forb cover (Figure 7) and total shrub cover (Figure 8). If the substantial differences in spatial patterns for vegetation cover of different types observed over the 1 ha sites are a function of interactions among and within the different plant communities or between the plant and soil communities over the site, then local processes are in effect, and the nonstationary model is appropriate. An overview of criteria and indicators of rangeland health by the National Research Council (1994) Committee

206

R.O. KUEHL, R. P. BRECKENRIDGE AND M. PANDA

on Rangeland Health recognized many of the biotic and abiotic processes on range sites as indicators of rangeland condition. Additionally, the existence of processes on the site may substantially affect the criteria used as indicators of ecological condition for the monitoring site. For example, are simple sample means of the measurements adequate to judge the ecological condition of the site? Perhaps not if the mean varies across the site and the signal is difficult to distinguish from the noise. Other metrics derived from an estimate of the nonstationary mean function #(s) may be more reasonable indicators of ecological condition. The challenge then is to determine those metrics. The scale of the processes under study will dictate the scale of measurement distances and the observed spatial pattern on the site. The challenge posed here is to determine the scale of activity for those processes identified for the monitoring process and then match the scale of measurement to the scale of the process for the indicator under consideration in order to have an effective monitoring indicator. Finally, the measurement tools or methodologies may figure prominently in the modeling process. For example, the use of quadrats, line transect intercepts, or point frames for collecting vegetation data may greatly influence the spatial patterns of the measurements because the cost and measurement errors differ among the methods. Thus the choice of methodology may considerably affect decisions on the most useful metrics for monitoring.

5. Conclusions Evaluating desertification requires a suite of indicators and methodologies that of necessity must be integrated into an overall monitoring program. The EMAPRangelands pilot study addressed several facets of this integration process at the field-sampling level, including the logistical requirements for multiple indicator classes, integrating measurement methodologies, plot traverse, and patterns of responses for the measurements to evaluate plot design strategies. This paper discussed the exploration of spatial relationships in the measurements and the consequences of these relationships for plot design strategies in the monitoring program. Preliminary analyses clearly found the need to expand current thinking on strategies for plot design and indicator criteria. To properly assess the condition of an arid area, the plot design needs to account for the larger spatial variability associated with natural resources on rangelands. This plot design study showed that some of the measurements exhibit a pattern of lower variability in a north-tosouth orientation than east to west. To properly capture this variability, areal cover designs appear to have advantages over long and narrow transects for properly characterizing a site. The complexities of indicator measurement patterns as they vary among monitoring sites beg the question of more complex site models that incorporate non-

INTEGRATED RESPONSE PLOT DESIGNS FOR INDICATORSOF DESERTIFICATION

207

stationary mean structure over the monitoring sites. These changes from current practices in modeling site-specific characteristics in turn promote new thinking about the type of metrics needed for indicator criteria to judge ecological condition. The scale of measurement necessarily must match the scale of the process being used as an indicator, and measurement tools or methodologies may figure prominently in the response plot design development. The selection of indicators for evaluating desertification on rangelands must reflect the biotic and abiotic processes that are producing observed patterns on a sample site. An oversized plot design study is a useful approach for evaluating the scale of processes and their patterns and should be employed in different ecological resource classes (i.e. grassland versus desertscrub) to ensure the resource is accurately assessed before making management decisions.

References Barbour, M.G., J.H. Burk and W.D. Pitts. 1980. Terrestrialplantecology, Benjamirt/Cumming, Menlo Park, California. Barkman, J.J. 1989. A critical evaluation of minimum area concepts. Vegetatio 85: 89-104. Bonham, C.D. 1989. Measurements for terrestrial vegetation, Wiley, New York. Bormann, EH. 1953. The statistical efficiency of sample plot size and shape in forest ecology. Ecology 34: 474--487. Breckenridge, R.P., W.G. Kepner and D.A. Mouat. 1993. Monitoring for assessing sustainability of arid ecosystems. Paper presented at the 25th International Symposium, Remote Sensing and Global Environmental Change. Graz, Austria, 4-8 April, 1993. Breckenridge, R.P., W.G. Kepner and D.A. Mouat. In press. A process for selecting indicators for monitoring conditions of rangeland health. Journal of Environmental Monitoring and Assessment. Brewer, A.C. and R. Mead, 1986. Continuous second order models of spatial variation with application to the efficiency of field crop experiments. Journal of the Royal Statistical Society A 149: 314348. Brummer, J., J. Nichols, R. Engel and K. Eskridge. 1994. Efficiency of different quadrat sizes and shapes for sampling standing crop. Journal of Range Management 47: 84-89. Campbell, J.B. 1978. Spatial variation of sand content and pH within single contiguous delineations of two soil-mapping units. Soil Science Society of America Journal 42: 460--464. Cochran, W.G. 1977. Sampling techniques. 3rd Ed. Wiley, New York. Cressie, N. 1991. Statistics for spatial data, Wiley, New York. Dietvorst, P., E. van der Maarel and H. van der Putten. 1982. A new approach to the minimal area of a plant community. Vegetatio 50:77-91. Freeman, G.H. 1963. The combined effect of environmental and plant variation. Biometrics 19: 273-277. Gholz, H. 1982. Environmental limits on above ground net primary production, leaf area, and biomass in vegetation zones of the Pacific Northwest. Ecology 63:469--481. Greig-Smith, P. 1952. The use of random and contiguous quadrats in the study of structure of plant communities. Annals of Botany 16:293-316. Haining, R.P., 1990. Spatial data analysis in the social and environmental sciences, Cambridge University Press, New York. Handcock, M.S and J.R. Wallis. 1994. An approach to statistical spatial-temporal modeling of meteorological fields. Journal of the American Statistical Association 89: 368-378. Hellden, U. 1991. Desertification - Time for an assessment? Ambio 20(8): 372-383. Hyder, D.N., C,E. Conrad, P.T. Tueller, L.D. Calvin, C.E. Poulton and EA. Sneva. 1963. Frequency sampling in sagebrush-bunchgrass vegetation. Ecology 44: 740-746.

208

R.O. KUEHL,R. P. BRECKENRIDGEAND M. PANDA

Hunsaker, C.T. and D.A. Carpenter (editors). 1990. Ecological indicators for the Environmental Monitoring and Assessment Program. EPA 600/3-90/060. U.S. Environmental Protection Agency, Office of Research and Development, Research Triangle Park, North Carolina. Kepner, W.G. and C.A. Fox (editors). 1991. Environmental Monitoring and Assessment Program. Arid ecosystems strategic monitoring plan. EPA 600/4-91/018, U.S. Environmental Protection Agency, Washington. Kepner, W. G., R.O. Kuehl, R.P. Breckenridge, J.R. Baker, D. O'Leary, J.M. Lancaster, S.G. Leonard, D.A. Mouat, T.G. Reinsch, M.A. Weltz, D.W. Sutton, R.L. Tidwell, A.C. Neale. G.E. Byers and R.L. Slagle. 1993. Environmental Monitoring and Assessment Program Arid Ecosystems 1993 implementation plan Colorado Plateau plot design pilot study. EPA/620/R-93/016. U.S. Environmental Protection Agency, EMSL-LV, Las Vegas, Nevada. Kepner, W.G., R.O. Kuehl, R.P. Breckenridge, J.M. Lancaster, S.G. Leonard, D.A. Mouat, T.G. Reinsch, K.B. Jones, A.C. Neale, T.B. Minor and N. Tallent-Halsell. 1994. Environmental Monitoring and Assessment Program Arid Ecosystems 1992 pilot report, EPA/620/R-94/015. U.S. Environmental Protection Agency, EMSL-LV, Las Vegas, Nevada. Leonard, S.G., G. Staidl, J. Fogg, K. Gebhardt, W. Hagenbuck and D. Prichard. 1992. Riparian area management: Procedures for ecological site inventory - with special reference to riparianwetland sites. U.S. Department of the Interior, Bureau of Land Management, Service Center, Denver. Mausbach, M.J., B.R. Brasher, R.D. Yeck and W.D. Nettleton. 1980. Variability of measured properties in morphologically matched pedons. Soil Science Society of America Journal 44: 358-363. Modjeska, J.S. and J.O. Rawlings. 1983. Spatial correlation analysis of uniformity data. Biometrics 39: 373-384. Mouat, D.A., C.A. Fox and M.R. Rose. 1992, Ecological indicator strategies for monitoring arid ecosystems, in D.H. McKenzie, D.E. Hyatt and V.J. McDonald (Eds.), Ecological indicators (vol. I and II), Elsevier Applied Science, New York. National Research Council. 1994. Rangeland health: New methods to classify, inventory, and monitor rangelands. Committee on Rangeland Classification, Board on Agriculture. National Academy Press, Washington. O'Brien, T. and D.D. Van Hooser. 1983. Understory vegetation inventory: An efficient procedure. Research paper INT-323. U.S. Department of Agriculture Forest Service, Intermountain Forest and Range Experiment Station, Ogden, Utah. O'Regan, W.G. and Arvanitis, L.G. 1966. Cost-effectiveness in forest sampling. Forest Science 12: 406--414. Overton, W.S., D. White and D.L. Stevens, Jr. 1990. Design report for EMAP, Environmental Monitoring and Assessment Program. EPA/600/3-91/053. U.S. Environmental Protection Agency, ERL, Corvallis, Oregon. Pearce, S.C. 1955. Some considerations in deciding plot size in field trials with trees and bushes. Journal of the Indian Society of Agricultural Statistics. 7: 23-26. Primack, R.P. 1993. Essentials of conservation biology, Sinauer Associates, Sunderland, Massachusetts, pp. 129-130. Proctor, C.H. 1985. Fitting H.E Smith's empirical law to cluster variances for use in designing multistage sample surveys. Journal of the American Statistical Association 80: 294-300. Reddy, M.N. and C,K.R. Chetty. 1982. Effect of plot shape on variability in Smith's variance law. Experimental Agriculture 18: 333-338. Ripley, B.D. 1981. Spatial statistics, Wiley, New York. Schreuder, H.T. and R.L., Czaplewski. 1993. Long-term strategy for the statistical design of a forest health monitoring system. Environmental Monitoring and Assessment 17:81-94. Smith, H.E 1938. An empirical law describing heterogeneity in the yields of agricultural crops. Journal of Agricultural Science 28: 1-23. Tardif, G. 1965. Some considerations concerning the establishment of optimum plot size in forest survey, Forestry Chronicle 41: 93-102. Tukey, J.W. 1977. Exploratory data analysis, Addison-Wesley, Reading, Massachusetts. U.S. Department of Agriculture (USDA). 1976. National range handbook. USDA Soil Conservation Service, Washington.

INTEGRATED RESPONSE PLOT DESIGNS FOR INDICATORSOF DESERTIFICATION

209

U.S. Department of the Interior (USDI). 1985. Rangeland monitoring trend studies. Technical Reference 4400-4. USDI Bureau of Land Management, Service Center, Denver. U.S. Department of the Interior (USDI). 1990. National range handbook. BLM Manual Handbook H-4410-1. USDI Bureau of Land Management, Washington. Waring, R.H., W.H. Emmingham, H.L. Gholz and C.C. Grief. 1978. Variation in maximum leaf area of coniferous forests in Oregon and its ecological significance. Forest Science 24: 131-139. Warrick, A.W., D.E. Myers and D.R. Nielsen. 1986. Geostatistical methods applied to soil science, in A. Klute (Ed.), Methods of soil analysis. Part 1. Physical and mineralogical methods, 2rid ed., Agronomy 9: 53-82. American Society of Agronomy, Madison, Wisconsin. Westoby, M., B. Walker and I. Noy-Meir. 1989. Opportunistic management for rangelands not at equilibrium. Journal of Range Management 42: 265-274. White, D., A.J. Kimberling and W.S. Overton. 1992. Cartographic and geometric components of a global sampling design for environmental monitoring. Cartography and Geographical Information Systems 19: 5-22. Wiegert, R.G., 1962. The selection of an optimum quadrat size for sampling the standing crop of grasses and forbs. Ecology 43: 125-129. Wilding, L.E and L.R. Drees. 1983. Spatial variability and pedology, in L.E Wilding, N.E. Smeck and G.E Hall (Eds.), Pedogenesis and soil taxonomy: Concepts and interactions, Elsevier Applied Science, New York, pp. 83-116.