Climate and Land-Use Change Effects on Soil Carbon Stocks over 150 Years in Wisconsin, USA

Abstract

Soil organic carbon is a sink for mitigating increased atmospheric carbon. The international initiative “4 per 1000” aims at implementing practical actions on increasing soil carbon storage in soils under agriculture. This requires a fundamental understanding of the soil carbon changes across the globe. Several studies have suggested that the global soil organic carbon stocks (SOCS) have decreased due to global warming and land cover change, while others reported SOCS may increase under climate change and improved soil management. To better understand how a changing climate, land cover, and agricultural activities influence SOCS across large extents and long periods, the spatial and temporal variations of SOCS were estimated using a modified space-for-time substitution method over a 150-year period in the state of Wisconsin, USA. We used legacy soil datasets and environmental factors collected and estimated at different times across the state (169,639 km²) coupled with a machine-learning algorithm. The legacy soil datasets were collected from 1980 to 2002 from 550 soil profiles and harmonized to 0.30 m depth. The environmental factors consisted of 100-m soil property maps, 1-km annual temperature and precipitation maps, 250-m remote-sensing (i.e., Landsat)-derived yearly land cover maps and a 30-m digital elevation model. The model performance was moderate but can provide insights on understanding the impacts of different factors on SOCS changes across a large spatial and temporal extent. SOCS at the 0–0.30 m decreased at a rate of 0.1 ton ha⁻¹ year⁻¹ between 1850 and 1938 and increased at 0.2 ton ha⁻¹ year⁻¹ between 1980 and 2002. The spatial variation in SOCS at 0–0.30 m was mainly affected by land cover and soil types with the largest SOCS found in forest and wetland and Spodosols. The loss between 1850 and 1980 was most likely due to land cover change while the increase between 1980 and 2002 was due to best soil management practices (e.g., decreased erosion, reduced tillage, crop rotation and use of legume and cover crops).

1. Introduction

Soil organic carbon is a potential sink for mitigating increased atmospheric carbon content under climate change. The international initiative “4 per 1000” launched aims at implementing best management practices on increasing soil carbon storage in soils under agricultue [1]. This requires a fundamental understanding of the soil carbon changes across the globe under changing climate and human activities. Previous studies have suggested that global soil organic carbon stocks (SOCS) have decreased due to global warming [2] and extensive land cover change [3]. However, other empirical studies and process-based models suggested that SOCS may increase due to increased carbon production under global warming [4,5,6,7], reduced respiration [8], lower land-use emissions such as afforestation [9], prairie restoration [10], and improved soil management under agriculture [11,12,13,14,15,16,17].

Uncertainty in SOCS estimates from the process-based models comes from two sources: the complex biogeochemical processes involving production and decomposition of soil organic carbon [17,18,19,20,21], and the spatial and temporal variations in SOCS and environmental variables (e.g., climate, topography, land cover, land use, and soil types) [22,23]. Without the accurate representation of all the model input parameters at fine spatial and temporal resolutions, the process-based models may contain large uncertainty [6].

A number of studies have attempted to use space-for-time substitution methods [24,25] to establish empirical models to predict soil properties changes over time and study the driving factors for changes in soil properties across large spatial extents [16,26,27,28]. The assumptions of the methods are that the empirical models established between measured soil properties and spatially varying environmental driving factors at one time can be extrapolated to other times by substituting the environmental factors at that time to predict the spatial variations of soil properties. Compared to traditional process-based models, these empirical models do not require detailed information about the input variables and are easy to implement. However, these models can be location- and time- specific.

To better understand the control of various environmental factors on SOCS and understand how land use and management affect SOCS under climate change, a case study was conducted for the state of Wisconsin over a 150-year period since early European-American settlement. The hypothesis of this study is that SOCS changed rapidly across Wisconsin due to land cover change since early European settlement and the SOCS continued to change under changing climate due to improved agricultural management practices. To achieve this, we adjusted the space-for-time substitution method to estimate the spatial and temporal variations of SOCS using legacy soil datasets and time-varying environmental factors collected at different times coupled with a machine-learning algorithm.

2. Materials and Methods

2.1. Study Area

Wisconsin is located in the north-central USA, with an area of 169,639 km² [20]. It has seven of the 12 soil orders of the United States Department of Agriculture (USDA) Soil Taxonomy: Alfisols (45%), Spodosols (18%), Entisols (12%), Mollisols (11%), Histosols (10%), Inceptisols (4%), and Ultisols (<1%). About 52% of the soils have a frigid soil-temperature regime, and 42% have a mesic soil-temperature regime. Nearly 80% of the soils in Wisconsin were developed from glacial materials (e.g., moraines, outwash, till, loess) with diverse origins, compositions, and thickness [29,30]. The soils in these areas have a gently rolling relief. The Driftless Area in southwestern Wisconsin has not been glaciated in the Quaternary. This area is covered with loess characterized by a somewhat rugged and hilly landscape. Wisconsin lies in the temperate continental climate zone. Northern Wisconsin has warm, wet summers and cold, wet winters with a Dfb climate, while southern Wisconsin has hot, wet summers and cool-to-cold, wet winters with a Dfa climate according to the Köppen–Geiger climate classification system. The mean annual temperature ranges from 4.0 to 9.4 °C and mean annual precipitation ranges from 740 to 960 mm [30].

The lands of Wisconsin have been occupied by humans for thousands of years. The first French immigrants arrived in 1634 and fur trade was the main interest [29]. Many immigrants arrived during the lead mining era in the 1820s and 1830s, but fur trade and lead mining declined in the 1850s. The use of land was intensified in the mid-1800s with a large number of immigrants coming from eastern US and Europe. Land clearing and crop cultivation began at that time [29] with wheat as the primary crop in addition to tobacco and cranberries [29]. Dairy farming took over wheat production in 1915 due to diseases and low wheat prices. Since then, Wisconsin became the leading producer of dairy products [30]. Nowadays, dairy farming is the leading agricultural activity in Wisconsin, followed by beef cattle and corn, greenhouse and nursery products, and soybeans [31].

2.2. Legacy Soil Profile Data

The soil organic carbon (SOC) dataset was collected from the databases of the National Soil Information System (NASIS) that was made of in situ pedon observations over the entire United States. In total, 550 soil profiles (Figure 1) were collected from 1980 to 2002 in Wisconsin that excludes samples with a pH larger than 6.5 that may contain inorganic carbon. Soil samples were collected from soil horizons in each profile to a maximum depth of 2.5 m. The SOC concentration was analyzed using the Walkley–Black method [32]. The horizon-based SOC concentrations were harmonized to a standard depth intervals (0–0.05 m, 0.05–0.15 m, 0.15–0.30 m) specified in the GlobalSoilMap specification [33] using a mass-conservation spline function [34]. We chose the maximum depth to 0.3 m because many soils in Wisconsin have root-restricting layers at that depth [29] and most studies following the Intergovernmental Pannel on Climate Change (IPCC) investigate SOCS within the top 0.3 m.

2.3. Environmental Covariates

Based on the proposed “Scorpan” framework, soil properties can be predicted from other soil and environmental factors, including soil (e.g., prior soil information), climate (e.g., precipitation and temperature), organism (e.g., land cover), relief (e.g., terrain attributes), parent material, time, and space [35]. These environmental factors can be derived from a number of spatial and temporal datasets including remote-sensing products that are available across a large extent and at fine resolutions. Coupled with data from soil profiles at multiple depths, the “Scorpan” framework has been applied to map soil properties and soil types with uncertainty over the past decades, termed “Digital Soil Mapping” [36]. Based on the review [37], SOC could be predicted from vegetation (e.g., historical land-use patterns, NDVI), climatic factors, and terrain attributes (e.g., landscape position, soil erosion processes). In this study, environmental covariates including soil property maps, climate information, land cover maps, and topography were used to develop the model for the period of 1980–2002 (

Table 1. Summary of environmental covariates for developing the random forest model for the period of 1980–2002 and backward prediction to 1850 and 1938.

Category	Variables	Spatial Resolution	Temporal Scale	Source
Soil	Clay	100 m	Stable	[37]
	Silt	100 m	Stable	[37]
	Sand	100 m	Stable	[37]
	Bulk density	100 m	Stable	[37]
	pH	100 m	Stable	[37]
Climate	Annual temperature	1 km	Yearly, 1980-2002	Daymet Version 3
	Annual precipitation	1 km	Yearly, 1980-2002	Daymet Version 3
	Annual temperature	1 km	1938	Same as 1980
	Annual precipitation	1 km	1938	Same as 1980
	Annual temperature	1 km	1850	[41]
	Annual precipitation	1 km	1850	[42,43]
Organism	Land cover	250 m	Yearly, 1938-2002	Resampled from [44]
	Land cover	–	1850	[45,46]
Topography	Elevation (m)	30 m	Stable	Downloaded from the USGS
	Slope	30 m	Stable	Calculated using SAGA GIS [39]
	Aspect	30 m	Stable	Calculated using SAGA GIS [39]
	Hillshade	30 m	Stable	Calculated using SAGA GIS [39]
	Topographic wetness index (TWI)	30 m	Stable	Calculated using SAGA GIS [39]
	Multiresolution index of valley bottom flatness (MRVBF)	30 m	Stable	Calculated using SAGA GIS [39]

). Topography information was assumed stable, whereas climate factors and land cover varied on a yearly basis.

Maps of clay, silt, and sand contents, bulk density, and pH were obtained at 100 m resolution [37] and the same depth intervals of the splined SOC concentrations. Clay, silt, and sand contents, pH, and bulk density were assumed to be stable during 1980–2002. Some properties such as soil pH, texture, and bulk density were only measured for some pedons. These were not used as covariates to construct the model as they were not available for all pedons and depths that had SOC measurements. Instead, we extracted these values from the 100 m soil property maps [37]. This introduced some uncertainty in the SOC model but it would ensure these soil properties were consistently estimated from one source and remained stable over the modeling period.

Annual maximum and minimum temperature and annual precipitation data from Daymet Version 3 (https://daymet.ornl.gov/) were available at 1 km resolution from 1980 to 2002. The climate information was extracted to the soil profiles and the entire state for specific years based on soil sampling dates.

Yearly land cover projections from 1938 to 2002 were obtained from the US Geological Survey (USGS) Modeled Historical Land Use and Land Cover for the Conterminous United States at 250 m resolution [38]. In general, Landsat imagery was used in combination with other ancillary data to project historical land cover types on a yearly basis across the United States [39,40]. Five land cover types, including forest, grassland, cropland, pasture, and wetland, were investigated for their effects on the changes of SOC stock. The information was extracted to the soil profiles and across the state for specific years based on soil sampling dates.

A 30-m digital elevation model (DEM) was downloaded from the USGS (https://data-wi-dnr.opendata.arcgis.com/datasets?q=dem). A number of terrain parameters were calculated using SAGA GIS [41], including slope, aspect, hillshade, topographic wetness index (TWI), and multiresolution index of valley bottom flatness (MRVBF). The terrain attributes were assumed stable during 1980–2002.

2.4. An Empirical Soil Organic Carbon (SOC) Model (1980–2002)

In this study, a random forest algorithm was used to develop the relationship between SOC concentrations at different depths and various stable variables (i.e., clay, silt and sand content, pH, depth, elevation, slope, aspect, hillshade, TWI, and MRVBF), and time-varying variables (i.e., annual maximum and minimum temperature, annual precipitation, and land cover). The samples were randomly split based on soil profiles with 70% of the total profiles for calibration and 30% of the soil profiles for validation. The model fitting was performed in R software (R version 3.4.3). The randomForest function in randomForest package [43] was used along with the train function from the caret package [44] to fit the model using the calibration dataset.

Once the model was established, SOC concentration was predicted onto a 500 × 500 m grid at three depth intervals (0–0.05 m, 0.05–0.15 m, 0.15–0.30 m) for 1980 and 2002. The maps of soil bulk density at the same depth intervals were obtained from the previously published dataset [38] at 100 m resolution and resampled to the 500-m grid using the nearest-neighbor algorithm for prediction. The SOC stock (SOCS, ton ha⁻¹) was calculated using the following equation (Equation (1)):

$$\mathrm{SOCS}\left(\mathrm{ton}h{a}^{-1}\right)=\mathrm{SOC}(\%)\times \mathrm{bulk}\mathrm{density}\left({\mathrm{g}\mathrm{cm}}^{-3}\right)\times \mathrm{depth}\left(\mathrm{cm}\right)$$

(1)

Soil bulk density was considered constant during 1980–2002. The coarse fragments data were not complete in this dataset, so we did not include it in the equation. Given that large areas of Wisconsin have shallow soils [30] and most studies including the Intergovernmental Panel on Climate Change, investigated SOCS within the top 0.3 m, we presented SOCS within 0.3 m for the purpose of comparison.

2.5. Backward Prediction to 1850 and 1938

The years 1850 and 1938 were selected to investigate the impacts of climate and land-use change on historical SOC stock and carbon changes. In the 1850s, a large number of immigrants arrived, and the lands of Wisconsin were first used for agriculture and dairy production [29]. Much of Wisconsin was undisturbed in 1850 and the SOC stock at that time can be considered as a baseline for comparison. USGS had a detailed land cover map as early as 1938. Considerable deforestation occurred between 1850 and 1938.

Annual temperature in 1850 across Wisconsin was predicted based on the changes of temperature from 1850 to 1980 using a coarse scale global CRUTEM4 dataset [45]. Annual precipitation in 1850 was estimated based on the changes of precipitation from 1850 to 1980 using the observed land surface precipitation data from 1850 to 1980 [46,47]. Annual temperature and precipitation in 1938 were set the same as 1980 in this study assuming the change of temperature was minimal across Wisconsin during this period [45]. Land cover data in 1850 of Wisconsin [48] were generated based on the previous work [49]. Land cover data in 1938 were obtained from the USGS (https://landcover-modeling.cr.usgs.gov/projects.php) at a 250-m resolution and resampled to the 500-m grid. The soil properties (clay, sand, and silt) and terrain attributes were considered constant over 150 years.

The SOC concentrations at three depth intervals were predicted from environmental covariates (soil, climate, land cover, and terrain attributes) to 1850 and 1938 using the established random forest model with estimated environmental variables in 1850 and 1938. This is equivalent to the space-for-time substitution method used in other studies to study the driving factors for changes in soil/ecological properties in space and time [26,27,28]. The SOC stock was then calculated using Equation (1).

2.6. Validation and Uncertainty Analysis

The random forest model established in 1980–2002 was validated using the randomly selected validation dataset consisting of 30% of the soil profiles. The coefficient of determination (R²), Lin’s concordance correlation coefficient, root mean squared error (RMSE, %), and mean error (ME, %) were calculated for the calibration and validation processes. Since there is no available measurement on SOC concentration in 1850 and 1938, the backward predicted maps cannot be validated. Bootstrapping was also used to estimate the uncertainty of the SOC stock maps of 1980 and 2002. Bootstrapping is a non-parametric approach that involves repeated random sampling with replacement of the calibration dataset [50]. The bootstrapped samples can be used to build the model and generate maps. The bootstrapping was run 10 times and the standard deviation (SD) of the 10 predicted SOC stocks maps were calculated to represent the prediction uncertainty in this study assuming bulk density was constant over time. The bootstrapping was also implemented in R software (R version 3.4.3).

2.7. Soil Organic Carbon Stock (SOCS) Change and Percent Projected Natural Vegetation Soil Carbon

Rate of change in predicted soil organic carbon stock (SOCS, ton ha⁻¹ year⁻¹) was calculated for three periods (1850−1938, 1938−1980, and 1980−2002) by subtracting the SOC values of different years point by point on the maps. In addition, percent projected natural vegetation soil carbon (PNVSC) was used to indicate carbon change as compared to its original state [28]. This is defined as the percentage comparing the contemporary soil carbon against a hypothetical amount of soil carbon that would be observed in the local landscape today if the areas under managed agroecosystems remained under natural vegetation [28]. In this study, the SOCS in 1850 was considered as the original carbon state under natural vegetation. The SOCS in 1938, 1980, and 2002 were divided by the SOCS in 1850 to obtain the percent PNVSC (%).

3. Results

3.1. Descriptive Statistics of SOC

The soil samples were collected between 1980 to 2002. We split the data into two time periods (1980–1990 and 1991–2002) to summarize SOC for different soil orders and land covers and the changes in SOC concentration with time (

Table 2. Mean soil organic carbon (SOC) concentration (%) at three depth intervals and two time periods for different soil orders and land covers.

Category		No.	1980–1990			No.	1991–2002
			0–0.05 m	0.05–0.15 m	0.15–0.30 m		0–0.05 m	0.05–0.15 m	0.15–0.30 m
Soil order	Alfisols	154	3.2	2.1	0.8	109	2.6	1.8	0.8
	Entisols	20	3.2	2.3	1.3	29	1.7	1.4	0.6
	Inceptisols	12	3.2	1.7	0.9	9	4.7	3.4	0.8
	Mollisols	18	1.9	1.7	1.1	30	1.7	1.6	1.1
	Spodosols	100	6.4	4.0	2.1	52	3.9	2.8	1.1
Land cover	Forest	169	5.4	3.3	1.6	121	3.2	2.3	0.9
	Grassland	3	2.3	1.9	1.0	–	–	–	–
	Cropland	61	2.6	2.0	1.0	57	2.1	1.7	0.9
	Pasture	70	2.6	1.9	0.8	49	1.9	1.5	0.7
	Wetland	12	2.8	1.7	0.7		8	6.7	3.1	0.5

Note: No. refers to the number of soil profiles in each soil order or land cover type at different time intervals.

). Most soils in the legacy soil dataset were Alfisols and Spodosols (Figure 1 and Table 2). Spodosols had higher SOC concentrations than other soils at three depth intervals. Mollisols had the lowest SOC concentrations at 0–0.05 m and 0.05–0.15 m depth, whereas Inceptisols had the lowest SOC concentrations at 0.15–0.30 m depth. SOC concentrations decreased between 1980–1990 to 1991–2002 in most soils except in Inceptisols. Summary statistics of SOC for different soil orders and land covers are given in Appendix A

Table A1. Summary statistics of soil organic carbon (%) at three depth intervals for different soil orders and land covers.

Category		Depth (m)	No.	Min. (%)	Mean (%)	Median (%)	Max. (%)	SD (%)	CV (%)
Soil order	Alfisols	0–0.05	263	0.2	2.9	2.3	10.9	2.0	67.5
	Alfisols	0.05–0.15	263	0.2	2.0	1.8	8.0	1.1	54.2
	Alfisols	0.15–0.30	263	0.0	0.8	0.7	5.4	0.6	73.6
	Entisols	0–0.05	49	0.1	2.3	1.6	19.3	3.0	128.5
	Entisols	0.05–0.15	49	0.1	1.7	1.2	17.4	2.5	140.3
	Entisols	0.15–0.30	49	0.1	0.9	0.6	9.8	1.4	151.0
	Inceptisols	0–0.05	21	0.3	3.9	3.3	16.8	3.6	92.8
	Inceptisols	0.05–0.15	21	0.3	2.4	1.8	12.7	2.6	108.8
	Inceptisols	0.15–0.30	21	0.1	0.8	0.5	3.9	0.9	108.0
	Mollisols	0–0.05	48	0.5	1.8	1.5	5.0	1.2	64.3
	Mollisols	0.05–0.15	48	0.4	1.6	1.4	4.5	1.0	62.0
	Mollisols	0.15–0.30	48	0.4	1.1	0.9	4.1	0.7	64.8
	Spodosols	0–0.05	152	0.0	5.5	3.3	68.4	9.5	171.7
	Spodosols	0.05–0.15	152	0.0	3.6	1.6	57.8	7.6	214.4
	Spodosols	0.15–0.30	152	0.1	1.8	0.9	42.8	4.5	259.1
Land cover	Grassland	0–0.05	3	1.2	2.3	2.8	3.1	1.0	43.1
	Grassland	0.05–0.15	3	1.0	1.9	2.3	2.3	0.7	38.8
	Grassland	0.15–0.30	3	0.6	1.0	1.0	1.5	0.5	45.6
	Cropland	0–0.05	113	0.3	2.4	2.0	9.4	1.7	70.8
	Cropland	0.05–0.15	113	0.3	1.9	1.6	8.0	1.3	66.1
	Cropland	0.15–0.30	113	0.1	0.9	0.8	5.4	0.7	77.7
	Forest	0–0.05	282	0.0	4.5	2.8	68.4	7.3	161.7
	Forest	0.05–0.15	282	0.0	2.9	1.7	57.8	5.8	200.2
	Forest	0.15–0.30	282	0.0	1.3	0.8	42.8	3.4	256.2
	Pasture	0–0.05	116	0.1	2.3	1.9	8.4	1.5	64.9
	Pasture	0.05–0.15	116	0.1	1.7	1.6	4.6	0.8	48.9
	Pasture	0.15–0.30	116	0.0	0.8	0.7	2.7	0.4	50.8
	Wetland	0–0.05	19	0.7	3.9	3.4	10.9	3.1	79.8
	Wetland	0.05–0.15	19	0.5	2.1	1.7	5.3	1.4	66.9
	Wetland	0.15–0.30	19	0.2	0.6	0.5	1.5	0.4	64.8

. Spodosols had the largest coefficient of variation (CV) of SOC at 0–0.05 m (171.7%). Mollisols had the smallest CV of SOC at 0–0.05 m (64.3%).

More than 75% of the soil profiles were sampled under forest and pasture (Table 2). Soils under forest contained the highest amount of SOC. The SOC concentrations decreased between 1980–1990 to 1991–2002 in most land cover types except in the soils of the wetlands. Forest had the largest CV of SOC at 0–0.05 m (161.7%) while pasture had the smallest CV of SOC at 0–0.05 m (64.9%) (Appendix A Table A1, excluding grassland).

3.2. Changing Climate and Land Cover from 1850 to 2002

The climate and land cover changed considerably between 1850 and 2002 (Figure 2). The mean annual temperature was more or less constant between 1850 and 1980 but the maximum temperature increased from 12.1 to 12.5℃ between 1980 and 2002 whereas the minimum temperature increased from 0.1 to 1.4℃ (Figure 2). Annual precipitation increased from 851 to 892 mm between 1850 and 1938 and remained constant during 1938–1980. It then increased to 969 mm in 2002.

Considerable land cover changes occurred between 1850 and 1938 (Figure 2 and Figure 3). In 1850, 79% of the state was covered by forests (Figure 3). Until 1938, the area of forest decreased to 41%. The area of grassland decreased from 7% to 1% which was mainly the result of the conversion of grassland to pasture and cropland. In total, the area of cropland increased by 26% and the area of pasture increased by 25% between 1850 and 1938. The land cover was relatively constant between 1938 and 2002 with slight reforestation in 1980 and 2002 (Figure 2).

3.3. Distribution of SOCS in Wisconsin

The SOC concentration was predicted using a random forest model and it had a coefficient of determination (R²) of 0.89, Lin’s concordance correlation coefficient of 0.90, mean error (ME) of 0.01%, and root mean square error (RMSE) of 2.04% in the calibration dataset and R² of 0.48, Lin’s concordance correlation coefficient of 0.67, ME of 0.62%, and RMSE of 2.22% in the validation dataset (

Table 3. Calibration and validation results of the random forest model to predict SOC concentration (%) in 1980–2002. Note: R2, coefficient of determination; RMSE, root mean square error; ME, mean error.

Dataset	Soil Profiles	SOC Samples	R2	Lin’s Concordance	RMSE (%)	ME (%)
Calibration	412	2,207	0.89	0.90	2.04	0.01
Validation	138	734	0.48	0.67	2.22	0.62

). Although the random forest model was accurate for the calibration data, it performed moderately based on the validation dataset.

SOCS distribution across the state was similar between 1850 and 2002 (Figure 4). The highest SOCS (>100 ton ha⁻¹) were located in the soils of the eastern and northern part, whereas the lowest SOCS (<40 ton ha⁻¹) was located in the soils of the central, northeastern, and northwestern Wisconsin. The estimated uncertainty was slightly larger in 2002 compared to 1980 and the highest uncertainty (>16 ton ha⁻¹) of the SOCS predictions was in the southeast and north where the SOCS was highest and larger than 100 ton ha⁻¹ (Figure 4).

3.4. Importance of Factors Affecting SOC

The importance of factors affecting SOC is listed in Figure 5. The most important factors were depth, followed by pH, aspect, elevation (DEM), annual maximum temperature, hillshade, silt content, slope, annual minimum temperature, MRVBF, TWI, sand and clay content, annual precipitation, and land cover (LC). Note that LC was not as important as other covariates in the model. This was because LC aggregated different soil properties, terrain parameters, and climate variables, while other covariates only represented soil, terrain, and climate variables. The controls of LC on SOC were possibly masked by other covariates that varied in space and time across Wisconsin.

3.5. SOCS Changes over 150 Years

SOCS decreased from 1850 to 1938 in most parts of Wisconsin (Figure 6). The highest decrease occurred in the soils in southeastern Wisconsin with the conversion of forest to cropland and pasture. A slight increase of SOCS was found in the soils of southwestern, central, and northern Wisconsin. The mean rate of decrease in SOCS was less than 0.1 ton ha⁻¹ year⁻¹ (

Table 4. Mean rate of change in soil organic carbon stock (SOCS, ton ha⁻¹ year⁻¹), and a total change in SOCS (Mton) at 0–0.3 m soil depth across Wisconsin for three time periods.

Period	Rate of Change (ton ha−1 year−1)	Total (Mton)	Potential Causes
1850–1938	–0.1–0	–41.4	Land cover changes, mainly deforestation (conversion of forest to cropland and pasture). High erosion rates
1938–1980	0–+0.1	+0.5	Relatively stable climate and a balance between deforestation (1938–1950) and afforestation (1950–1980); erosion controlled in large parts of the state
1980–2002	+0.2	+39.8	Increased carbon production due to increasing temperature, increased irrigation and N fertilization and better soil management (e.g., decreased erosion, reduced tillage, crop rotation and use of legume and cover crops) In some areas, increased carbon loss due to soil erosion, deforestation, and cultivation

). The largest decrease was found in Mollisols (

Table 5. Estimated rate of change in soil organic carbon stocks (SOCS, ton ha^-1 year^-1) at 0–0.3 m soil depth for different soil orders and land cover types in Wisconsin, USA. Note: values in bold indicate an increase in soil organic carbon stock. Values for grassland may be not reliable due to the small area of grassland in Wisconsin.

		1850–1938	1938–1980	1980–2002
Soil order	Alfisols	<–0.1	<+0.1	+0.1
	Entisols	<–0.1	<+0.1	+0.1
	Inceptisols	<–0.1	<–0.1	+0.2
	Mollisols	–0.1	<+0.1	+0.1
	Spodosols	<–0.1	<–0.1	+0.4
Land cover	Forest	<–0.1	<+0.1	+0.1
	Grassland*	<–0.1	<–0.1	+0.3
	Cropland	–	<+0.1	+0.1
	Pasture	<–0.1	<+0.1	+0.1
	Wetland	–	<+0.1	+0.2

) where approximately 80% of the area of forests were converted to cropland and pasture.

Between 1938 and 1980, SOCS increased with a mean rate of less than 0.1 ton ha⁻¹ year⁻¹ across Wisconsin (Table 4 and Figure 6). SOCS increased in all land cover types, except for grassland (Table 5). The increase in SOCS was more significant in the northern part of the state (mainly land cover of the forest) with PNVSC greater than 100% (Figure 6).

Between 1980 and 2002, SOCS increased in most part of the state except for the northern margin and center of Wisconsin at a rate of 0.2 ton ha⁻¹ year⁻¹ (Figure 6 and Table 4). The largest SOCS increase was in Spodosols (0.4 ton ha⁻¹ year⁻¹) which were mostly under forests (82%).

4. Discussion

4.1. Prediction and Distribution of SOCS

Temperature, topography, and soil properties explained the variation and change of the SOCS over time (Figure 5). In general, SOCS was higher in the soils of northern Wisconsin because of higher elevation (Figure A1), lower temperature (Figure A2), and dominance of forest vegetation (Figure 3). The spatial distribution of SOCS in soils under forest (60–80 ton ha⁻¹) was consistent with other studies [15,16] that mapped SOCS in Wisconsin. Central Wisconsin with sandy soils had considerably lower SOCS compared to northern and southern Wisconsin, whereas the southern and eastern part had soils with higher silt and clay contents (Figure A1) and thus higher SOCS. In addition, the spatial distribution of SOCS was related to soil pH at the 0–0.05 m depth (Figure A1), whereby strongly acidic soils (pH < 5.5) in northern Wisconsin were associated with low SOCS and less acidic soils (pH: 5.5–6.5) in southern Wisconsin had high SOCS.

The prediction error is less than 16 ton ha⁻¹ for most parts of the state (Figure 4). The largest errors are in the northern and southeastern parts of the state likely due to the limited data for these regions (Figure 1). This can also be due to the use of covariates that are only available at a coarser resolution (e.g., 100 m for soil property maps) and contain an inherent uncertainty. Furthermore, the uncertainty for the model in the past cannot be evaluated.

4.2. Soil Carbon under Land Cover Change Only

Between 1850 to 1938, climate and in particular annual temperature was relatively stable [45], so the main driver for SOCS change was land cover change. The decrease of SOCS in the southeastern part is likely due to the conversion of forest and grassland to cropland and pasture [51,52]. SOCS changes are also affected by the soil texture and the initial SOC content (refer to

Table A2. Estimated rate of change in soil organic carbon stock (SOCS, ton ha^-1 year^-1) at 0–0.3 m soil depth for different soil orders and land cover in Wisconsin, USA. Note: values in bold indicate an increase in soil organic carbon stock.

Soil Order	Land Cover	1850–1938	1938–1980	1980–2002
Alfisols	Forest	<–0.1	<+0.1	<+0.1
	Grassland	<–0.1	<–0.1	0.2
	Cropland	–	<+0.1	0.1
	Pasture	–	<+0.1	0.1
	Wetland	<–0.1	<+0.1	0.1
Entisols	Forest	<–0.1	<+0.1	0.1
	Grassland	<–0.1	<–0.1	0.1
	Cropland	–	–	–
	Pasture	–	–	–
	Wetland	<–0.1	<+0.1	0.2
Inceptisols	Forest	<–0.1	<–0.1	0.2
	Grassland	<–0.1	<–0.1	0.2
	Cropland	–	–	–
	Pasture	–	–	–
	Wetland	<–0.1	<+0.1	0.2
Mollisols	Forest	–0.1	<+0.1	0.1
	Grassland	<–0.1	<–0.1	0.1
	Cropland	–	–	–
	Pasture	–	–	–
	Wetland	–0.1	<+0.1	0.1
Spodosols	Forest	<–0.1	<–0.1	0.3
	Grassland	<–0.1	<–0.1	0.5
	Cropland	–	<+0.1	0.1
	Pasture	–	–	–
	Wetland	<–0.1	<+0.1	0.3

and [53]). Mollisols had a slightly larger decreasing rate of SOCS change (0.1 ton ha⁻¹ year⁻¹) compared to other soil orders between 1850 and 1938 while Spodosols had larger increasing rates of SOCS changes between 1980 and 2002 compared to other soil orders.

4.3. Soil Carbon under Climate Change and Agricultural Management

There are four possible causes for the increased SOCS. Firstly, after 1938, atmospheric temperature [54] and CO₂ [7,55,56] have slightly increased which increases photosynthesis and net primary production. Secondly, the area of irrigated agriculture has increased (more water more biomass production) with the construction of high capacity wells from 1957 to the present with an annual construction peaked about 300 wells in 1977 [57,58]. Thirdly, the inputs of N fertilizers increased from 31 kg ha⁻¹ in 1964 to 110 kg ha⁻¹ in 1981 and is around 90 kg ha⁻¹ afterward (Figure 7, [59,60]). Lastly, adoption of best management practices such as conservation tillage [61,62,63], crop residue [64], cover crops [65], crop-rotation [66], as well as many new crop cultivars (e.g., corn, soybeans) introduced between 1930 and 1980 that can produce more biomass than the previous cultivars which leads to an increase of crop residues compared to 30 years ago [67]. This could have increased accumulation of organic matter [61].

Decreased SOCS between 1938 to 1980 can be due to increased soil respiration resulting from increasing temperature [45]. The increase in soil respiration per 10℃ rise in temperature (Q10) is about 2.0 and depends on the type of SOC and environmental temperature [61]. Decreased SOCS may also be the result of increased soil erosion, which has been previously reported in the Driftless areas in the southwest of Wisconsin [68,69] and other regions [70,71], and as estimated by others [16,26].

4.4. Implications for Land Use and Management under Climate Change

It was found that SOCS within top 0.3 m decreased at a rate of less than 0.1 ton ha⁻¹ year⁻¹ between 1850 and 1938, increased at a rate of 0.2 ton ha⁻¹ year⁻¹ between 1980 and 2002 (Table 2). Soil is a carbon sink with a potential soil carbon sequestration rate ranging from 0.2 to 0.5 ton ha⁻¹ year⁻¹ in the top 0.3 to 1.0 m [1]. The carbon sequestration values calculated for Wisconsin are slightly lower than reported elsewhere. In Ireland, SOCS increased at rates of 1.0–1.5 tons ha⁻¹ year⁻¹ over 10 years at 0–0.5 m due to grazed swards [72] with N fertilizer applied up to 500 kg N ha^-1 year⁻¹ and at a rate of 2.2–2.5 ton ha⁻¹ year⁻¹ over 16 years at 0–0.3 m due to afforestation of grassland [73].

The current SOC sequestration rate may decrease with climate change. For example, insufficient nutrient supply (e.g., N and P) would lead to a reduction in the net primary production [74]. In addition, diurnal and seasonal variations in temperature [75], rainfall [76], and soil and ecosystem-dependent land management practices [77,78,79] are likely to decrease SOCS. However, researchers have predicted increased SOC sequestration in Wisconsin until 2050 under projected increased air temperature and precipitation and increased conversion of cropland to the forest [16]. This suggests that the goals of the international “4 per 1000” initiative may be feasible under the current projection of climate change and best land management practices.

However, the model used by [16] only included soil samples collected over a short period of time and the projection may fall beyond the range of the environmental covariates. Therefore, it is worth further investigating the location-specific carbon management practices and climate change across a larger spatial extent using this modified space-for-time substitution approach, which requires high-resolution maps of spatial and temporal varying soil properties including soil bulk density, texture, and coarse fragment concentration [61], climate and topography [79], temperature sensitivity [80], as well as vegetation and forest management practices [81], cropland management history [82] and urbanization [83]. Priorities should be given to these soils so that the best land management practices [84,85] and natural climate solutions [86] can be applied accordingly to increase SOCS and sustain the ecosystem services of the soils.

5. Conclusions

To better understand how a changing climate, land cover, and agricultural activities influence soil organic carbon stock (SOCS) across large extents and long periods, a spatial and temporal analysis of SOCS was conducted in Wisconsin over 150 years using a random forest machine learning algorithm. The model performed moderately well and provides some understanding of the impacts of different factors on SOCS changes across a large spatial and temporal extent. From this analysis, the following can be concluded:

• The loss of SOCS between 1850 and 1980 in Wisconsin, USA, was most likely due to land cover change while the increase of SOCS from 1980 to 2002 was the result of improved soil management practices;
• SOCS changes are affected by land cover and differ by soil type;
• Location-specific management should be developed to increase SOCS under climate change;
• This modified space-for-time substitution approach needs to be tested in other regions of the world to further investigate the influences of the location-specific carbon management practices and climate change on SOCS changes.