Abstract

Ecological functioning of the intensive, homogeneous agroecosystems in the Chippewa River Watershed (CRW), MN, USA, can be improved by reducing soil erosion, runoff, and nutrient leaching. These ecosystem services can be achieved through increased perennials in crop rotations to diversify land use and sustain carbon sequestration. We calibrated, validated, and used APSIM software to simulate the effect of 100 yrs each of historical and future climate change scenario (IPCC-A2) on biophysical processes in representative soil types of the predominant farming systems in CRW. The interrelationships between crop rotations, soil types, climate variables, and ecosystem services indicated that not all objectives of sustainable agro-ecosystem are compatible, and tradeoffs among them are necessary. Site-specific and diversified crop rotations that comply with the environmental constraints of climate and soils could lead to more efficient implementation of strategies to improve ecosystem services in the watershed if current management practices of high external inputs and tillage persisted.

1. Introduction

The extent and global scale of the environmental costs associated with intensive agriculture are already documented [1], and their effects on ecosystem services have been extensively researched [2, 3]. Ecological functioning of the intensive agroecosystems in the Upper Midwest, in general [4], and in MN’s Chippewa River Watershed (CRW), in particular, is being compromised by soil erosion, runoff, and nutrient leaching [5, 6]. The land use in CRW is predominately centered on commodity production systems that deliver corn, soybean, and livestock products for domestic consumption and for an expanding export market [6, 7]. Similar to other conventional land-use systems in the Upper Midwest of the Corn Belt (i.e., corn-soybean rotation, conventional tillage, N fertilizer), losses of ~30% of the average 125 kg applied nitrogen (N) per ha to subsurface drainage in CRW are not uncommon [7, 8]. The rise in intensive, homogeneous croplands, especially corn production systems, is responsible for a substantial proportion of the increase in NO3-N values observed in surface and underground waters of the Midwest since the 1970s. River-born nutrients, especially NO3-N from tile-drained fields in the CRW, contribute to environmental pollution within and beyond the watershed [9]. Increased corn production as a feedstock for ethanol would result in a concomitant increase in N loads of ~5 to 18% [10]. This raises concerns about the environmental, agronomic, and economic sustainability of the corn-soybean crop rotation in this and similar watersheds in the Upper Midwest.

Recent research [11] indicated that managed ecosystems often fail to respond smoothly to external changes. Therefore, more efforts are being devoted to studying ecological regime shifts, thresholds, and resilience of ecosystems. This is becoming ever more important for agroecosystems because many rural communities depend on their provisioning services. These services are primarily determined by the structural diversity (i.e., number of crops within a field) and temporal diversity (i.e., complex crop rotations), and their interaction with management and climatic factors, and contribute to diversification of agroecosystems. However, the antagonistic relationship between ecosystems’ capacity to simultaneously sustain the availability of regulating services and agricultural production of an agro-ecosystem should be considered when designing new management practices and complex crop rotations [12].

The CRW supports agroecosystems that flourish in a temperate climate regime that could see significant changes in the next few decades due to climate [13] and land-use changes [13, 14]. This watershed is a source of many ecosystem services including provisioning services such as the production of food, feed, and biofuel; regulating services such as climate regulation, flood regulation, and water quality control; cultural services such as recreational, aesthetic, and other benefits; supporting services such as nutrient cycling [3]. This watershed is among the most altered landscapes in the Upper Midwest. The majority (60–94%) of the land area within the watershed is used for row crops and livestock production, and pressure to increase row crop production is intensifying with elevated demand for biofuels [10]. Fertilizer application, erosion attributed to mechanical tillage and exposed soil, and rapid movement of nutrients and sediment to streams via tile drainage are primary causes for high nitrate and turbidity in this [15] and other watersheds in the Midwest [14].

Historical land-use changes from prairie native grasses to row crops have had a significant effect on the hydrology of the CRW [4, 5]. It is anticipated that future population growth and changing agricultural demands will inevitably force land-use changes and, therefore, land cover in the watershed [15]. Limited land cover during the short growing season (~110 days yr−1) by annual crops in the CRW plays a significant role in influencing the water and energy balance at the land surface via its effect on light interception, transpiration, and evaporation from crop canopies. In addition to the physical properties of soil types in the watershed that can have a strong influence on the nature of the hydrologic response, changes in vegetation cover can alter the surface energy balance and evapotranspiration. Land-use change from native prairie in the CRW to row crops resulted in about 10–15% decrease in evapotranspiration and increases of 20–30% in total runoff [16].

Results of earlier simulation studies [15] suggested that average annual percent change in precipitation, runoff, soil loss, and crop yield, in response to projected climate change in six locations in MN, would increase by 16.1, 34.6, 22.5, and 58.6%, respectively. Similar projections in biophysical processes, except crop yield, were obtained [17] using historical and modeled data in the Midwest. Earlier monitoring efforts in the CRW [5] resulted in establishing a water quality and quantity baseline and emphasized the need for more perennials on the landscape in order to enhance both environmental and agronomic agro-ecosystem services, promote conservation incentives, and develop viable markets for perennials and integrated crop-livestock production systems. Previous modeling efforts predicted that land resources in the CRW can be multifunctional; they can provide a large number of functions related to social, economic, and environmental prosperity and sustainability [4, 6]. Given the complexity of agroecosystems, and the spatiotemporal variability of the involved biophysical processes within the CRW, modeling is needed for assessing the effect of climate change and alternative management strategies [18]. These strategies should include more permanent plant cover on the landscape and their effect should be assessed on a number of ecosystem services cum biophysical processes (e.g., biomass and grain yield of crops, carbon sequestration in the soil, nitrogen loss, runoff, and soil erosion).

Due to current land-use systems, the CRW faces a number of environmental challenges, including degradation of water quality, threats to biodiversity, and increased flooding, soil erosion, and nutrient loss to both surface and underground water [4, 19]. In addition to these environmental challenges, the effect of spatiotemporal climate variability in the watershed is likely to be intensified by climate change, which is predicted to disrupt many ecosystem functions, altering their capacity to provide goods and services and rendering them more susceptible to degradation [20, 21]. The question facing farmers and researchers alike is how to optimize a diversified agro-ecosystem, maximize production and profits, and minimize negative environmental effect of future agroecosystems [11]. The objectives of this simulated study were to (1) predict the long-term effect of diverse crop rotations on agro-ecosystem services and the relationship between soil types, crop rotations, and their response to historical and future climate change, (2) model the effect of targeted land-use changes on managed agroecosystems, and (3) develop predictive and validation models of future ecosystem services under the A2 climate change scenario if current management practices of high external inputs and conventional tillage are continued.

2. Materials and Methods

2.1. Watershed Physiography

The Chippewa River drains 5,387 km2 of mixed natural and managed ecosystems in West Central MN (Figure 1). Annual crops, predominantly grown in monocultures, occupy between 60 to 94% of the land area in the major agro-ecological regions of the watershed [19]. Recently, however, increased demand for domestic biofuels is expected to influence land-use options [22] and will impact several biophysical processes in managed and natural agroecosystems within the watershed. The CRW has high-value ecological features including seven major lakes, two state parks with prairie, forest and lake areas, numerous wildlife management and waterfowl production areas, and ~3,000 km of intermittent and perennial streams [46].

2.2. Selection of Soil Types

Given the complexity of the agroecosystems and the spatiotemporal variability of the involved biophysical processes in the CRW, we selected 12 representative soil types (series) to cover the major physiographic features of the watershed (Figure 1). General information on these soils is available at http://soils.usda.gov/survey/geography/ssurgo/ detailed characteristics can be found at http://soils.usda.gov/survey/geography/ssurgo/description.html. The soil types (series) were (1) BaB2: (408 meters above sea level, m) Evansville Barnes loam, 2–6% slope, eroded; (2) Fa: (411 m) Kensington Flom silty clay loam, 0–2% slope; (3) BaB: (369 m) Swan Lake Barnes loam, 3% slope; (4) HaA: (361 m) Swan Lake Hamerly clay loam, 0–3% slope; (5) SuA: (383 m) Cyrus Svea loam, 0–2% slope; (6) BbC2: (356 m) Long Beach Barnes-Langhei loams, 6–12% slope, eroded; (7) EvB: (415 m) Glenwood Estherville loam, 2–6% slope; (8) Pf: (365 m) Starbuck Parnell and Flom silty clay loams; (9) J55A: (314 m) Benson Sedgeville loam, channeled, 0–2% slope, occasionally flooded; (10) L33B: (373 m) Sunburg Kandiyohi clay, 2–5% slope; (11) J51A: (312 m) De Graff Bearden-Quam, depression, complex, 0–2% slope; (12) J30A: (347 m) Clara City Tara silt loam, 1–3% slope.

2.3. Ecosystem Services

Data on a number of ecosystem services (or biophysical processes, i.e., total biomass and grain yield of crop rotations, carbon sequestration in the top 1 m soil profile, NO3-N and NH4-N loss, runoff, and soil erosion; Table 1) were compiled from simulated data to quantify and evaluate the long-term effects of climate change, soil types, traditional corn-soybean crop rotation, and more complex crop rotations with corn-soybean-wheat and a forage perennial crop, and their interactions in the watershed.

2.4. Simulation Scenarios and Modeling

Each of one hundred and sixty eight combinations of (1) two global climate change scenarios (Past (P-), based on historical weather local data for the period 1900–1999, and future (F-), based on simulated weather data for the period 2000–2099, climate change scenarios; see below), (2) seven crop rotations (corn-soybean, corn-soybean-wheat- (1–5 years of alfalfa), and continuous alfalfa for eight years), and (3) twelve soil types was simulated for 100 years. We utilized the Agricultural Production Systems Simulator (APSIM 7.3; 23) that integrates modules of cropping, management, and biophysical processes within farming systems for improved decision support. APSIM is a point-based model that runs on a daily time step and is capable of yield prediction across land areas. It estimates crop productivity and potential effects of farming systems at selected locations in the study area based on the site’s soil properties and climate regimes [23]. The soil nitrogen (NO3-N and NH4-N) and carbon transformation within APSIM are largely based on experiences within the CERES model; both variables were estimated to a depth of 1 m in the soil profile. APSIM includes a generic crop module to capture physiological characteristics and processes that are shared between different crop species. The curve number method [24], where soil type, land use, management practices, and antecedent moisture are lumped and represented on average by a curve number, was used to simulate runoff. The conceptual basis of this method is supported by empirical data. Soil erosion was estimated using runoff volume, vegetative cover, land slope, and several soil physical and chemical properties [23].

A database derived from on-station and on-farm research on traditional and alternative cropping systems, including perennial biomass crops, was utilized to calibrate the simulation model. Sensitivity analysis was carried out for each crop in the crop rotation. The database included detailed quantitative measurements on crops, soils, nutrients, and weather variables [25]. To conduct the simulations and subsequent multivariate statistical analyses, we used historical weather data for the past 100 years (compiled from the weather station of West Central Research and Outreach Center, Morris, MN, and will be made available upon request at ∖∖Mw45-ws1Common DataWeather DataWCROC 8AM observations (1886 - present)WCROC Weather1886-2010.xlsx) and will be referred to as the P-climate change scenario, as well as simulated changes in future (21st century) temperature, precipitation, and CO2 concentrations based on the A2 climate change scenario which was generated by the Coupled Global Climate Model (CGCM) 3.1 which was developed by the Canadian Center for Climate Modeling and Analysis, Canada (http://www.ec.gc.ca/ccmac-cccma/default.asp?lang=En) and will be referred to as the F-climate change scenario. The IPCC A2 scenario was selected for simulations (Figure 2) because it assumes “business as usual” in terms of CO2 emissions and is characterized by a high rate of growth in CO2 emissions and most closely reproduced the actual emissions trajectories during the period since the scenario was defined [20].

For each of the 12 locations, representing soil types, conventional management practices, as defined in Jaradat and Weyers [25], and position in the landscape within the watershed (Figure 1), we assessed the effect of current conventional (corn-soybean) and alternative crop rotations (corn-soybean-wheat-alfalfa (with 1 to 5 years of alfalfa following three years of grain crops), and continuous alfalfa for 8 years) on several crop and soil biophysical processes [26]. Conventional tillage and local management practices were used in simulations along with 100-year historical weather data and 100 years of projected weather data based on the A2 climate change scenario as predicted by CGCM 3.1 [20].

2.5. Statistical Analyses

Data derived from each simulation run were subjected to several statistical analyses procedures in order to (1) assess the effect of climate, soil types and their interaction on each ecosystem service using univariate analysis of variance procedures; (2) quantify the effect of, and estimate the variance due, to fixed (crop rotations and crops within crop rotations) and random factors (soils types and their interaction with crop rotation) on these ecosystem services taking the spatial and physiographic (longitude, latitude, and elevation) characteristics of the 12 soil types into consideration; and using the spatial option in the multivariate residual maximum likelihood (REML) procedure; (3) classify the soil types according to their reaction to climate change and identify which ecosystem services significantly contributed to this classification; (4) contrast ecosystem services based on past and future climates and develop predictive models for each ecosystem service; (5) quantify the effect of increased proportion of a perennial crop on biomass and grain yield of crop rotations and develop statistical models to predict its effect on the remaining ecosystem services; (6) develop partial least square regression models to assess the predictability of future status of ecosystem services based on their current status in each soil types. Annual data output was used to carry out statistical procedures described in (1) to (3), whereas normalized aggregate data were used to carry out the statistical procedures in (4) to (6). The largest value of each ecosystem service in all 12 soil types was considered as 1.0 and the remaining 11 values were expressed as a percent of that value. Normalized data is more appropriate to use for comparisons across all ecosystem services due to different units of measurement and to assess multiple agro-ecosystem functions [27]. Several modules in STATISTICA [28] and GenStat software programs [29] were used to carry out statistical analyses and to develop predictive models.

3. Results

The 12 sites selected for the simulation study represented the predominant soil types, cropping systems, and physiographic features in the watershed (Figure 1). In addition, these locations represent different positions in the landscape with slopes ranging from 0.0 to 12%. The North-South and West-East spread of the 12 locations and their proximities/distances from water bodies as well as their soil physical and chemical properties and soil profile characteristics captured the maximum diversity within the watershed and provided necessary data to validate the simulation results.

3.1. Past and Future Climates

Total annual rainfall under P- and F-climate change scenarios were almost equal (611 and 596 mm, resp.); however, differences in minimum and maximum rainfall were large in magnitude and resulted in the respective ranges of 628 and 332 mm, coefficients of variation of 21.1 and 12.4%, skewness of 0.026 and −0.098, and kurtosis of −0.22 and −0.49. Mean annual temperatures and their standard deviations (°C) for P1 and A2 were 5.74 (1.07) and 9.35 (1.58). Although they had almost equal skewness values (0.205 and 0.25, resp.), their respective kurtosis values (0.12 and −0.36) were very different. Carbon dioxide (CO2) is projected to increase (almost linearly from the initial level of 370 to 856 ppm at the end of the 21st century; Figure 2) under the A2 climate change scenario; however, it is expected to display a large level of variation (coefficient of variation = 25.2%) and will be positively skewed (0.47) with a negative kurtosis (−0.95) value.

Deviations of projected mean temperature (d-Temp, °C), total rainfall (d-Rain, mm), and CO2 concentration (d-CO2, ppm) for the next 100 years based on the A2 climate change scenario from the respective mean values derived from historical records for the last 100 years (P-scenario) are presented in Figure 2. Future mean temperature is projected to be higher by 2–8°C, especially towards the second half of the 21st century. Deviations of future from past rainfall fluctuated widely (from −300 to +395 mm), whereas those of CO2 are expected to range from 0.0 (at the start of the simulation period) to ~500 ppm (at the end of the 21st century) above the current level assuming “business as usual” in terms of CO2 emissions [20]. d-Temp was negatively correlated ( 𝑟 = 0 . 6 4 ; 𝑃 < 0 . 0 0 1 ) with past temperature (P1) and positively correlated with projected temperature ( 𝑟 = 0 . 8 5 ; 𝑃 < 0 . 0 0 1 ) and CO2 ( 𝑟 = 0 . 7 1 ; 𝑃 < 0 . 0 0 1 ) under A2 climate change scenario. d-Rain exhibited a parallel trend; it was negatively correlated with past rainfall under P-scenario ( 𝑟 = 0 . 8 8 ; 𝑃 < 0 . 0 0 1 ) and positively correlated with rainfall under A2 ( 𝑟 = 0 . 5 7 ; 𝑃 < 0 . 0 0 1 ). d-CO2, on the other hand, was positively correlated with future temperature under A2 ( 𝑟 = 0 . 8 3 ; 𝑃 < 0 . 0 0 1 ) and with d-Temp ( 𝑟 = 0 . 7 1 ; 𝑃 < 0 . 0 0 1 ).

3.2. Effect of Climate, Soil Types and Their Interaction on Ecosystem Services

Results of univariate analyses of variance of ecosystem services based on simulated response to climate change, soil types, and their interaction (Table 1) indicated that all ecosystem services are likely to be impacted by past and future climates, soil types, and their interaction. The F-values (Table 1) and their probabilities suggested that climate change would be the major factor in determining the future magnitude of biomass and grain yield production, the loss of NH4-N, and runoff, whereas soil types would be the main determinant of the level of carbon sequestration, the loss of NO4-N, and soil erosion. The magnitude of the effect of the interaction between climate and soil types, when compared with the effect of climate and soil types, would be a major contributor to determining future grain yield, runoff, and soil erosion levels.

3.3. Multivariate Analyses of Ecosystem Services

The variance portion attributed to fixed (crop rotations and crops within crop rotations) and random (soils types and their interaction with crop rotation) factors on ecosystem services taking into consideration the spatial and physiographic (longitude, latitude, and elevation) characteristics of the 12 soil types (Table 2) suggested that different factors and their interactions accounted for a wide range of variability in past and future ecosystem services. Crop rotations, whether based on past or future climates, differed significantly for all ecosystem services under study; however, crops within crop rotations differed significantly as to their biomass and grain yield under both climate change scenarios. The soil-related processes (soil carbon, NO3-N, NH4-N, runoff, and erosion) were not significantly impacted by differences between crops within crop rotations. Nevertheless, under future (A2) climate change scenario, all variables except soil carbon differed significantly from those based on past weather.

Variances accounted for by soil types in all ecosystem services were highly significant and ranged from 29% for NH4-N to 74% for soil carbon based on P-scenario and from 32% for NH4-N to 64% for runoff, based on the F-scenario. Variances accounted for by the interaction between crop rotations and soil types were highly significant, except biomass and grain yield based on P-scenario and grain yield based on F-scenario. In general, these variances were smaller in magnitude as compared with those attributed to differences between soil types with three exceptions: larger variance estimates due to the interaction between crop rotations and soil types were predicted for NO3-N and NH4-N based on P-scenario and for NH4-N based on F-scenario.

A number of ecosystem services are expected to be differentially impacted by future climate change scenario depending on how the soil types may react to changes in future weather variables. More variation in NO3-N (48 versus 34%) and runoff (64 versus 54%) will be accounted for by future climate change scenario, whereas slightly less variation in biomass (56 versus 61%), soil carbon (62 versus 74%), and erosion (45 versus 52%) will be accounted for by the same scenario. These projected changes may or may not be associated with changes in variances accounted for by the interaction between crop rotations and soil types (Table 2). Variances accounted for by this interaction are expected to increase for soil carbon (12 versus 23%) and to decrease for NO3-N (54 versus 31%), NH4-N (56 versus 46%), and erosion (40 versus 27%) in response to future climate change scenario. Soil carbon is a variable of interest in studies involving crop rotations and their interaction with soil types and climate change scenarios. The simulation study suggested that there will be a slight increase (4%) in soil carbon sequestration due to increased CO2 as projected by future climate change scenario. Also, soil carbon is one of a few variables where soil types and their interaction with crop rotations accounted for large (86 and 85% based on P- and F-scenarios, resp.) portions of its variation.

3.4. Response of Soil Types to Climate Change

Each soil type was 100% correctly classified based on its biophysical properties and its reaction to past and future climate change scenarios (Figure 3). Two discriminant functions, with significant contributions of all ecosystem services (i.e., biophysical characteristics) based on their Wilks’ λ and R 2 values (except biomass and grain yield), accounted for 0.89 and 0.11 of variation between soil types, respectively, due to past, and for 0.92 and 0.08 of variation, respectively, between soil types due to future climate change scenarios. Soil carbon, runoff, erosion, NO3-N, and NH4-N, in decreasing order, significantly contributed to this discrimination based on their Wilk’s λ and 𝑅 2 values (Figure 3). Similarly, runoff, erosion, soil carbon, NO3-N, and NH4-N, in decreasing order, significantly contributed to this discrimination based on the future climate change scenario.

The scatter of soil types along the first discriminant root (which accounted for most variation) indicated that a few of these soil types (Fa, J30A, L33B, Pf, and SuA) may experience relatively small variation in their biophysical properties due to future climate change. The remaining soil types may undergo medium (BaB, BaB2, EvB, and HaA) or large (BbC2, J51A, and J55A) variations. These future shifts in soils’ biophysical characteristics are largely quantified by their standardized coefficients on both discriminant roots (Figure 3(b)). Soil carbon and runoff had larger standardized coefficients as compared to erosion, NO3-N, and NH4-N, while biomass and grain yield had the smallest standardized coefficient values.

3.5. Spatiotemporal Variation in Ecosystem Services

Correlation coefficients between past and future normalized values, and regression models (Figure 4) describing the nonlinear temporal changes in ecosystem services, indicated that soil types (representing spatial ecosystem services), most likely, will respond to climate change in different manners. The relatively small NH4-N data, as compared to NO3-N data, was not included in Figure 4 for convenience. Temporally, past and future ecosystem services were significantly correlated. The largest 𝑟 -values were found for soil carbon ( 𝑟 = 0 . 9 8 ; 𝑃 < 0 . 0 0 1 ) and erosion ( 𝑟 = 0 . 9 7 ; 𝑃 < 0 . 0 0 1 ); intermediate values were found for biomass ( 𝑟 = 0 . 7 8 ; 𝑃 < 0 . 0 1 ), grain yield ( 𝑟 = 0 . 7 6 ; 𝑃 < 0 . 0 1 ), and NO3-N ( 𝑟 = 0 . 7 3 ; 𝑃 < 0 . 0 1 ), and smallest value was found for runoff ( 𝑟 = 0 . 5 6 ; 𝑃 < 0 . 0 5 ); however, the variations explained by the nonlinear regression models expressing future values of ecosystem services as function of past values were all significant ( 𝑃 < 0 . 0 5 ) and ranged from small ( 𝑅 2 = 0 . 3 1 for runoff) to large ( 𝑅 2 = 0 . 9 6 for soil carbon).

Regression models (expressing future ecosystem services as functions of their respective past ones; Figure 4) for biomass and grain yield mirrored each other as to statistical descriptors, most productive (BaB; P-biomass and P-grain y i e l d = 1 ; F-biomass and F-grain yield = ~0.6 and ~0.42, resp., Figure 4) and least productive (BbC2) soil types, and temporal changes in productivity. Biomass and grain yield under future climate change scenario are predicted to decline, with largest expected values being 0.92 and 0.78 of their respective maximum values. The soil type BaB would be the most productive under past (P-biomass and P-grain yield) but not under future climate change scenario where J51A would be the most productive soil type. If we consider the 0.5 as the agronomically and environmentally acceptable cutoff level for all ecosystem services, then BbC2 would be the only soil type to fall below the biomass cutoff level, and along with EvB and Bab2 soil types would fall below the grain yield cutoff level under future climate change scenario.

Carbon sequestration is expected to increase slightly under future climate change scenario, with two soil types (J55A and J51A) having values above 0.5 while the remaining 10 soil types falling below the 0.5 cutoff level. Runoff and soil erosion are expected to decrease and NO3-N leaching is expected to increase under future climate change scenario. The EvB soil type would experience the largest NO3-N leaching and runoff, and BbC2 would experience the largest soil erosion under past climate change scenario. Two soil types (BbC2 and BaB) would experience the largest NO3-N leaching under future climate change scenario; however, their future values are expected to reach 0.3–0.4 of the maximum.

Most of the nonlinear regression models describing the relationships between each of the six ecosystem services (Figure 4) had negative intercepts and a negative polynomial regression coefficient (except those for runoff and erosion). Notwithstanding the 𝑅 2 values (range from 0.73 to 0.98), this suggests that the future levels of biomass, grain yield, soil carbon, and NO3-N leaching under future climate change scenario would reach a maximum then decrease in relation to their past values, whereas runoff and erosion are expected to decrease with low ( 𝑅 2 = 0 . 3 1 ) and high ( 𝑅 2 = 0 . 9 7 ) probability, respectively, under future climate change scenario. The extreme erosion (and to some extent runoff) values for the BbC2 soil type (slope of 6–12%) may have skewed the distribution of the normalized values of these two ecosystem services.

3.6. Perennial Crop Effect on Ecosystem Services

The effect of increased proportion of a perennial crop (i.e., alfalfa as a nitrogen-fixing forage legume crop) on ecosystem services (Figure 5), whether under past or future climate change scenarios, was nonlinear and positive for biomass, grain yield, and carbon sequestration and nonlinear and negative for NO3-N leaching, runoff, and soil erosion. Maximum expected biomass and grain yield due to the perennial effect would be achieved at 0.45 perennial proportion of the crop rotation under both climate change scenarios. However, the maximum expected biomass and grain yield under future climate change scenario would be 0.97 and 0.8 of their respective maximum values expected under past climate change scenario. In both cases, the conventional corn-soybean crop rotation is expected to produce the least biomass and grain yield under both scenarios. The predictive models estimated biomass under past and future climate change scenarios with 𝑅 2 of 0.58 and 0.74, respectively. The respective 𝑅 2 values for grain yield were 0.62 and 0.64. The effect of percent perennial in the crop rotation on carbon sequestration was almost linear under both climate change scenarios (regression coefficients in the nonlinear portion of the regression models were 0.02 and 0.03 and were associated with large 𝑅 2 values). The continuous perennial is expected to sequester the maximum amount of carbon while the corn-soybean crop rotation is expected to sequester the minimum; however, the differences were minimal as compared with other ecosystem services.

The level of NO3-N leaching is expected to decline under future compared with past climate change scenario as quantified by their respective linear regression coefficients ( 𝑏 1 = 1 . 2 2 and −0.45), nevertheless, the corn-soybean crop rotation would experience the maximum NO3-N leaching level. The regression models predicted a smaller decrease in NO3-N leaching (nonlinear regression coefficient, 𝑏 2 = 0 . 1 7 ) under past, and a larger decrease ( 𝑏 2 = 0 . 8 6 ) under future, climate change scenario when the perennial crop occupies more than 0.8 of the crop rotation. Runoff and soil erosion predictions followed trends similar to that of NO3-N leaching. The corn-soybean crop rotation had the largest runoff and soil erosion predicted values under past and future climate change scenarios. Runoff is predicted to amount to ~0.8 and 0.35 and soil erosion 0.5 and 0.35 of their respective maximum values under continuous perennial due under past and future climate change scenarios, respectively. The nonlinear regression coefficients in both cases ( 𝑏 2 = 0 . 1 1 and −0.15, resp.) were larger than their respective linear regression coefficients ( 𝑏 1 = 0 . 0 4 and −0.05, resp.). Nonetheless, the percent perennial in the crop rotation predicted these ecosystem services with large probabilities ( 𝑅 2 > 0 . 8 9 ).

3.7. Predictability of Future Ecosystem Services under Climate Change

Partial least squares regression modeling (PLS; Figure 6) captured maximum variation in independent (P-) and dependent (F-climate change scenario) ecosystem services at the prediction phase ( 𝑅 2 𝑋 and 𝑅 2 𝑌 , resp.), as well as in the dependent ecosystem services, at the validation phase ( 𝑄 2 𝑌 ) of model development. The PLS regression models differed as to all three statistics (i.e., 𝑅 2 𝑋 , 𝑅 2 𝑌 , and 𝑄 2 𝑌 ). These differences can be attributed to differences between physical and chemical characteristics of different soil types and to the nature and strength of association between different ecosystem services within a particular soil type.

Two partial least squares regression components (PLSC1 and PLSC2) accounted for a maximum ( 𝑄 2 𝑌 = 0 . 8 9 ) to a minimum ( 𝑄 2 𝑌 = 0 . 4 8 ) of the total variation in predicting future (F-) as a function of past (P-) ecosystem services. However, the portion accounted for by each component differed among soil types. Loadings (i.e., correlation coefficients between an ecosystem service and the PLSC) primarily on PLSC1, and to a lesser extent on PLSC2, quantify the power and direction (positive or negative effect) of that ecosystem service in predicting future ecosystem services. Biomass and grain yield loadings on PLSC1 were positive but small as compared with the rest of ecosystem services, except in BbC2 where their loadings on PLSC1 were small and negative. The remaining ecosystem services were mostly clustered in tow groups with positive (e.g., NO3-N, runoff, and erosion) and negative (e.g., soil carbon and NH4-N) loadings on PLSC1, with some variation in these groupings depending on soil type. Loadings on PLSC2 were generally smaller than loadings on PLSC1.

The carbon-nitrogen associations with, and loadings on, PLSC1 and PLSC2 separated 11 of the 12 soil types into three groups. The BbC2, unlike other soil types, displayed positive loadings of soil carbon and NH4-N and negative loadings of NO3-N on PLSC1. Soil carbon and NH4-N had negative and NO3-N had positive loadings on PLSC1 in four soil types (Bab2, BaB, EvB, and HaA) regardless of their loadings on PLSC2. Three soil types (SuA, J55A, and Pf) had negative loadings of soil carbon and NO3-N and positive loading of NH4-N on PLSC1 regardless of their loadings on PLSC2. The remaining four soil types (J51A, Fa, L33B, and J30A) had negative loadings of all three ecosystem services on PLSC1, and mostly negative loadings of NO3-N on PLSC2. The NO3-N was associated with runoff and erosion with positive loadings on PLSC1 in four soil types (Bab2, BaB, EvB, and HaA) and negative loadings in one soil type (BabC2). On the other hand, NH4-N was associated with soil erosion and runoff with positive loadings on PLSC1in three soil types (SuA, J55A, and Pf). Finally, runoff and soil erosion loaded positively, and NO3-N and NH4-N loaded negatively on PLSC1 and were associated with soil carbon in the remaining four soil types (J51A, Fa, L33B, and J30A).

4. Discussion

Agricultural activities influence the environment and many agro-ecosystem services through their effect on water quality, nutrient cycling, microclimate control, biodiversity conservation, and carbon sequestration, at local, regional, and global scales [30]. Pressing challenges of climate change, natural resource depletion, and environmental degradation triggered renewed interest in current and future status and value of ecosystem services by scientists, farmers, and policy makers at a global [31], regional [32], and local levels, including the CRW [6, 19]. Agroecosystems in the CRW are more complex than other resource systems due to the large number of manipulated biophysical processes, inputs, outputs, and their interactions [6, 19, 31]. Productivity of these agroecosystems is heavily constrained by physical ecosystem components, especially climate (e.g., rainfall and temperature) and soils; therefore, the main focus in this analysis was more on relationships among components and processes than on single-factor cause-effect relationships.

The CRW is a source of provisioning (food, feed and biofuel), regulating (climate regulation, flood regulation, and water quality control), cultural (recreational and aesthetic values), and supporting services (nutrient cycling) [6]. Historical land-use changes have had a significant effect on the hydrology of the watershed and therefore on the provision of these ecosystem services. Due to changing agricultural demands, land use in the watershed will continue to change, and we should be able to derive insight into more sustainable management practices for the future [16]. The environmental implications of intensive corn and soybean production in the CRW are serious cause of concern. Runoff and tile drainage from intensively managed fields are well-documented causes of nonpoint source contamination of surface and underground water bodies with sediments and nutrients [26]. Therefore, farmers are facing growing demands to reduce environmental pressure and the greater implementation of diversified agroecosystems may be a productive strategy to build resilience into future agroecosystems.

4.1. Past and Future Climate

Agricultural production in the Midwest is critically dependent on the weather, especially rainfall amount and distribution during the growing season [17]. We selected the A2 scenario because it represents a very heterogeneous world where economic development is regionally oriented and economic growth and technological change are relatively slow [17, 20]. By mid- and late 21st century, A2 emissions are expected to reach 550 ppm and to exceed those of the A1F1 “high” scenario, respectively [17]. However, despite some similarities between weather variables of past and future weather, future weather based on A2 (business as usual; [17]) is expected to radically depart from the long-term distribution of rainfall, temperature, and CO2 concentration as evidenced by its statistical moments. Deviations of future temperature (d-Temp), rainfall (d-Rain), and CO2 (d-CO2) from past records, being positively correlated with future values of their respective variables, gave insight into the expected reaction of crop rotations and soil types to future climate and into their expected effects on ecosystem services. The amount of variation in several ecosystem services accounted for by differences between past and future climate change scenarios (Tables 1 and 2) are indications of the direction and level of impact future weather may have on these services.

The 100-year period was chosen to allow for the dynamics of the ecological effects of land-use change to play out and be captured in simulation results [33]. Empirical results [13] indicated that trends in rainfall and temperature may explain 35–40% of crop yields under the Corn Belt conditions, whereas, according to projections based on future climate change scenario, for each additional degree (°C) of future warming during summer months, corn and soybean yield could potentially decrease by 13 and 16%, respectively. Nevertheless, future increases in CO2 may not be effective at boosting productivity of all crops in the CRW; the expected rise in temperature may offset the positive impact of higher CO2, especially on small grain crops [34, 35].

4.2. Effect of Climate, Soil Type, and Their Interaction on Ecosystem Services

All ecosystem services in the CRW are likely to be impacted by climate, soil types, and their interaction (Table 1). However, proper management of the intensively cropped agroecosystems in CRW under future climate change is expected to improve soil quality and sustain future productivity [32, 36, 37]. In their comprehensive review of the effect of climate change on agriculture, Hatfield et al. [18] concluded that American agriculture, especially in the Midwest, is vulnerable to climate change but has considerable potential for adaptability. In the next few decades, agriculture will have to undergo significant transformation in order to adapt to climate change and protect natural resources and the environment. Therefore, if the current management practices of high external inputs and conventional tillage continued well into the 21st century, ecosystem services, depending on soil type [38], may respond mostly negatively to future climate change, although this response could be mediated by the soil carbon dynamics [38, 39] and modified by appropriate management practices [40, 41]. The latter may include adoption of more complex crop rotations to match soil types [39, 42]. Soil carbon was one of a few variables where soil types and their interaction with crop rotation accounted for large (86 and 85% based on past and future climates, resp. Table 2) portions of its variation. Therefore, matching crop rotations and soil types may constitute the most appropriate mitigation option in the future [43, 44], and adjusting crop sequence proportions may alleviate some of the negative climate impact on ecosystem services [38, 45].

4.3. Multivariate Statistical Analyses of Ecosystem Services

Sources of variation in this simulation study were classified as fixed (crop rotation and crops within crop rotations) and random (soil types and the interaction of crop rotations with soil types) (Table 2). The classification made it possible to quantify the variance components that can be attributed to the random factors [29] assuming that the soil types represented a random sample across the landscape in the watershed and they are expected to respond at random to climate change.

Ecosystem services, especially those related to soil processes, are more dependent on crop rotations rather than on crops within rotations (Table 2), provided that a perennial crop is rotated with annuals [39, 46]. The increasing portion of a perennial crop was accredited with most of the beneficial effects of complex crop rotations by positively impacting soil-related processes [47], especially carbon sequestration [48]. However, apart from a slight increase in carbon sequestration in response to future climate change scenario, the simulation study suggested that the reaction of a soil type to different crop rotations will account for most of the variation in potential carbon sequestration under the A2 climate change scenario.

4.4. Response of Soil Types to Climate Change

The soil resources in the CRW are typical of intensively cropped soils in the Midwest [37] with a wide range of physiographic, physical, chemical, and biochemical characteristics (Figure 1; [25]). However, the predominant cropping system, which is based on corn-soybean crop rotation and external inputs, is supporting a few ecosystem services, mainly grain production. The scatter diagram of soil types along the first discriminant root (which accounted for most variation in Figure 3) illustrates the wide range of response of soil types to past and future climate change scenarios. A few soil types may not drastically react to future climate change (e.g., SuA, Pf, Fa), while others (e.g., J55A and BbC2) are projected to undergo dramatic shifts in their ecosystem services. Intensive tillage and external inputs for sugar beet crops are typical management practices for J55A, while BbC2 has the steepest position (6–12%) in the landscape and is prone to high erosion. The standardized coefficients, derived from discriminant analysis, parallel their Wilks’ lambda and 𝑅 2 values and highlight the greater influence of the soil-related processes (i.e., soil carbon, N leaching, runoff, and erosion) in mediating soil-type response to climate change as compared to the smaller response of provisioning service (i.e., biomass and grain yield). The antagonistic relationship between soil carbon and NO3-N, as quantified by their standardized coefficients on both discriminant functions may indirectly support the generally accepted negative effect of nitrogen fertilizer on net carbon sequestration [49] in some soil types (J55A, BbC2) more than others (e.g., J30A).

The current crop production system, even on a highly productive soil (e.g., BaB and J51A, Figures 3 and 4), will have difficulty attaining consistently a goal of achieving acceptable levels of ecosystem services (e.g., <10 mg NO3-N/L in tile drainage water, even though prudent N fertilization is followed) [50]; therefore, a tradeoff between ecosystem services is inevitable [37]. Nevertheless, long-term experimental results [25, 50] indicated that the type of cropping system and management practices had a greater positive effect on soil performance or ecosystem services than did nitrogen fertilization. Our simulated data, supported by empirical results of a long-term crop rotation experiment within the CRW [25], can be used to prescribe specific management strategies for different soil types or locations in this watershed.

4.5. Spatiotemporal Variation in Ecosystem Services

The spatial variation in ecosystem services (Figure 4) is depicted by the scatter of soil types around the 1 : 1 regression line describing the constant (no-change) relationship between the level of these ecosystem services based on past and future climate change scenarios. Apart from biomass and grain yield, soil types are predicted to undergo major non-linear shifts in their position around the 1 : 1 spatial relationship in response to future climate change. The normalized values allowed for direct comparisons to be made between different ecosystem services and would facilitate the development of simple performance-based index to assess multiple agro-ecosystem functions [27]. Soil carbon exhibited the largest spatial shift (a minimum of 0.1 for EvB to maximum of 1.0 for J55A in response to future climate change), whereas the smallest spatial shift was displayed by runoff (a minimum of ~0.05 for J51A to a maximum of 0.55 for EvB).

The temporal variation is quantified by the non-linear relationships between past and future simulated values of all ecosystem services, except soil carbon, in which the polynomial regression coefficient ( 𝑏 2 = 0 . 0 6 ) was not significant (Figure 4). When solving each non-linear equation assuming that the average simulated value during the past 100 yrs (i.e., the 𝑋 variable) equals 1.0, grain yield (1.44) and soil carbon (1.09) under future climate change are expected to exceed their respective values by 0.44 and 0.09 units, respectively, and with the respective medium ( 𝑅 2 = 0 . 5 8 ) and large ( 𝑅 2 = 0 . 9 6 ) probabilities. Biomass, NO3-N, runoff, and soil erosion under future climate change scenario are expected to be smaller than the maximum value by 0.22, 0.45, 0.41, and 0.24 units, respectively. The expected decrease in simulated biomass can be explained on the basis of increased proportion of the perennial crop in all crop rotations, except the traditional corn-soybean [41]. Increased proportion of the perennial crop in the crop rotation can also explain the expected favorable increase in grain yield, due to the favorable residual nitrogen fertilizer effect [51], and carbon sequestration [41, 47], and decrease in nitrogen leaching, runoff, and soil erosion [52, 53], due to deeper and more vigorous root systems.

The seemingly small increment of carbon sequestered under future climate change scenario can be attributed to the positive feedback in which increases in global temperature (Figure 2) may lower the ability of soils to sequester carbon [54]. In addition, cropping intensification by including more perennials in the crop rotation may result in eventual saturation of the soil profile [48]; hence the small positive increment of carbon sequestered after a 100-year simulation run.

4.6. Perennial Crop Effect on Ecosystem Services

Strategies to optimize crop rotations and deliver ecosystem services need to be developed and tested in the CRW. Current simulation results are expected to contribute to this objective and are supported by empirical results of a long-term crop rotation experiment carried out within the CRW [25]. Those results do indicate some potential rotation strategies for different soil types or locations in the watershed. Simulation and empirical results suggested that adding a perennial to the corn-soybean crop rotation can have substantial positive effects on multiple ecosystem services [26, 47]. Growing alfalfa continuously for 8 years may result in larger biomass and sequestered carbon compared to a more complex crop rotation that includes grain crops (e.g., CSW+4A). Including alfalfa to enhance the complexity of a crop rotation was reported to have increased soil organic carbon sequestration by 20 g C m−2 yr−1 as compared to the corn-soybean rotation [49]. Perennials can mobilize and retain nutrients and water very effectively over extended periods of time as compared with annual crops, thus providing a mechanism to mitigate climate change, and reduce runoff, nutrient loss, and soil erosion [53]. These perennial advantages can be attributed to the fact that more residue diversity and the concomitant increased microbial activity lead to faster decomposition and loss of more carbon from the residue of a complex crop rotation as compared to continuous alfalfa [55].

Regression models suggested that simulated and normalized biomass and grain yield under past and future climate change scenarios are expected to reach their maximum at about 80 and 45% perennial in the crop rotation (Figure 5). Therefore, farmers may have to choose between maximizing biomass or grain yield under a particular management and climate change scenario [32, 40], knowing that a single year of alfalfa is expected to boost grain yield of the succeeding corn crop by ~30% as compared to corn-soybean rotation ([26]; Figure 5). However, this positive effect is expected to gradually diminish as indicated by the negative non-linear regression coefficients for biomass and grain yields ( 𝑏 2 = 1 . 1 2 and −1.82, resp.) under past and future climate change scenarios ( 𝑏 2 = 1 . 8 2 and −1.55, resp.).

Runoff, soil erosion, and NO3-N leaching mirrored each other and are expected to decline under future climate change scenario. The most likely factors, besides the perennial crop effect, are future rainfall patterns and temperature regimes and their interactions (Figure 2). However, NO3-N leaching under future climate change scenario, unlike runoff and soil erosion, was estimated with less certainty ( 𝑅 2 = 0 . 5 7 ). Nitrate leaching from a corn-soybean crop rotation is expected to be the highest and is bound to increase under future climate change scenario (Figure 5). Greater NO3-N concentrations and drain flows from corn-soybean rotations resulted in annual losses of NO3-N that were 35 times greater than NO3-N losses in a perennial system [9]. Other studies [8] estimated annual losses from the traditional corn-soybean rotation at about 30% of applied nitrogen. The perennial crop effect in reducing nitrate leaching is expected to be stronger under future climate change scenario (~50% reduction due to continuous alfalfa) and would level off at about 50–60% perennial in the crop rotation.

4.7. Predictability of Future Ecosystem Services under Climate Change

The partial least squares regression modeling procedure used to predict future ecosystem services is more reliable than ordinary regression because it takes into consideration the variation in the dependent variables, the independent variables, and their joint variation [28, 29]. Furthermore, the validation phase of PLS modeling improves the level of prediction by using a rigorous cross-validation procedure [29]. The level of (un)certainty in the validation models, quantified by the 𝑄 2 𝑌 values, accounted for by PLSC1, ranged from 0.32 (HaA) to 0.61 (Pf) and by PLSC2 from 0.08 (EvB) to 0.38 (HaA). Many factors contributed to these results including patterns and levels of association between dependent and independent variables (i.e., ecosystem services) in different soil types. Also, the uncertainty of the climate model [53], especially the fertilization effect of and response of different crops and (future) crop varieties to CO2, may have contributed to these results.

5. Conclusions

Ecosystem services, whether provisioning, regulating, or supporting services, provided by the CRW are being threatened by anthropogenic as well as climatic factors. Although there will be an antagonistic relationship between ecosystem’s capacity to simultaneously sustain the availability of provisioning and regulating ecosystem services, complementarities among different types may increase and the strength of their tradeoffs may diminish with the spatial complexity of the agricultural landscape. Therefore, an in-depth understanding of agroecosystems is essential for the future of land-use change in the CRW. The simulation results offer strategies to optimize site-specific crop rotations in this watershed under projected climate change. The simulated effect of increasing perennial land use in agroecosystems across the watershed on several biophysical processes suggested that farmers can help improve environmental health through increased perennial cover with sustained carbon sequestration and the concomitant reduction in soil erosion, runoff, and nutrients leaching. Diversifying the corn-soybean crop rotations by including a perennial crop, especially in erosion-prone soil types and locations in the watershed, would provide a strategy to mitigate negative environmental effects caused by corn and soybean production while providing an additional source of income based on new regional markets for food and biomass from perennials and diverse crops. Adaptive management, where stakeholders contribute to optimize resource use and minimize anthropogenic and climatic effect on the production base within the watershed, is expected to help develop multifunctional production systems that can produce standard commodities as well as a wide range of other system services.

Acknowledgments

The authors thank The Walton Family Foundation for financial support, Jon Starr for his technical assistance, and Beth Burmeister for editing the paper. This study was supported by the National Institute of Food and Agriculture, USDA, under Agreement no. 2010-65615-20630. The use of trade, firm, or corporation names in this publication is for the information and convenience of the reader. Such use does not constitute an official endorsement or approval by the United States Department of Agriculture (USDA) or the Agricultural Research Service (ARS) to the exclusion of others that may be suitable. USDA is an equal provider and employer.