Land use change and ecosystem service tradeoffs on California agricultural land

California is the top producer of cow milk and dairy products in the U.S., accounts for two-thirds of nation's total fruit and nut crop, and is second only to Florida in its floriculture production (CDFA 2020). In 2021 the State of California's farms received $51.1 billion for their production output making it the number one agricultural earner in the country. Agricultural land in California includes cropland, livestock farms, agroforestry, pasture and rangeland, and California agriculture employs up to 900 000 workers throughout a typical year (RMN 2022). Despite agricultural importance, the number of farms in California is decreasing while the average farm size increases. Between 2012 and 2017 (the two most recent U.S. Census of Agriculture years), the number of farms in California decreased by 9% while the total acreage of farmland decreased by only 4%, resulting in a 6% increase in average farm size (U.S. Census Bureau 2022). The 2022 estimate for California agriculture was 68 400 farms (−3% from 2017) spanning 24 million acres (−2.1%), with an average farm size of 351 acres (+0.86%) (USDA 2022b). This 'farmland consolidation' is a theme throughout the United States, and more and more is being considered problematic for farmers, particularly Black farmers and new farmers, and for ecosystem health (Union of Concerned Scientists 2021). Much of California also falls within the California Floristic Province, characterized by its high rate of endemism and floral diversity and listed as one of the world's biodiversity hotspots (Myers et al 2000). California is home to nearly 7000 plant species, many of which are found nowhere else on the planet. Due to the amount of land dedicated to farmland, there is inevitable crossover between wildlife and these semi-natural agricultural landscapes.

The California Department of Conservation maintains a detailed map of farmland throughout the state where agricultural land is categorized based on soil quality, irrigation status and cultural importance (Department of Conservation 2018). To maintain a consistent land area for our study, we used the most recent map of farmland in California as the bounds of analysis. All land designated as 'grazing land', 'farmland of local importance', 'prime farmland', 'farmland of statewide importance', or 'unique farmland' (approximately 12.8 million hectares) were included in the results. This approach allowed us to analyze multiple land cover types within farmland areas, rather than strictly studying land designated as 'cropland' in land cover maps. We categorized farmland based on one of five agricultural commission regions according to the California Agricultural Commissioners and Sealers Association (CACSA) (figure 1). County agricultural commissioners are tasked with protecting agricultural interests and the environment, and in 1922 the CACSA formed to facilitate collaboration between commissioners according to their region. The system is complex but 'flexible enough to address local issues', and each region is based on shared agricultural concerns, production, climate, water resources and policy (Huang 2019).

Zoom In Zoom Out Reset image size

We performed analysis of two ES, carbon storage and soil erosion control. Both ES models required a land user/land cover map, for which we used the United States Department of Agriculture's Cropland Data Layer (CDL). This dataset is derived from satellite imagery annually starting in 2008, has a spatial resolution of 30 m and includes individual crop types as well as non-agricultural land cover types (Boryan et al 2011). The California CDL metadata states that overall accuracy for the 2010 and 2020 data layers are 74.65% and 78.8%, respectively. These accuracy values are impacted by land cover classes that are prone to error such as grassland categories (Lark et al 2017). To reduce overall classification error caused by spectrally-similar land cover types (Johnson 2013, Lark et al 2017), we reclassified the CDLs into broader categories (supplementary table 1), with one exception: for the purposes of mapping soil erosion control, we maintained all of the individual non-tree and non-grass crop types from the original CDL, since different crop types result in different erosion potential (Panagos et al 2015).

To measure FP, we compared production yield of five of the top commodities produced in California in 2010 and 2020: almonds, grapes, strawberries, cattle/calves for beef, and milk/cream. State-level commodity production values were obtained from the California Agricultural Statistics reports created by the USDA National Agricultural Statistics Service California Field Office. County-level commodity production statistics for almonds, grapes and strawberries are reported in the County Agricultural Commissioner's Reports for crop years 2010 and 2020 (www.nass.usda.gov/). Statewide milk and beef production statistics are available in California Agricultural Statistics reports, but only the number of dairy or beef cows is reported at the county level. Therefore, we used the number of milk and beef cows in January of years 2010 and 2020 as indicators for milk and beef production.

Soil erosion control is a land parcel's ability to prevent soil erosion and delivery of sediments to nearby water sources (MEA 2005). Lack of soil erosion control has several consequences including loss of soil fertility and productivity, loss of soil carbon, and nutrient loading into water bodies, which in turn negatively impacts water quality (Keesstra et al 2016). Erosion control can be improved on farmland through conservation practices such as cover cropping, hedgerows, contoured slopes, and reduced or no tillage (Mao et al 2020, Du et al 2022). When soil erosion is minimized, farmers benefit economically from increased soil fertility, greater rates of water infiltration which reduces irrigation requirements, and improved crop yields (Alam 2018). Businesses such as hydroelectric dams, municipal water supply and tourism also benefit economically from reduced soil erosion, and communities that share the farmland watershed benefit from better water quality, reduced flood risk, and improved aquatic habitat. The prevention of soil erosion through land management practices on working lands has the potential to maintain or increase carbon storage and biodiversity in soils and boost overall climate change resilience on farmlands (Javadinejad et al 2021).

We modeled soil erosion for 2010 and 2020 using the InVEST Sediment Retention Model (Natural Capital Project 2022). InVEST is a suite of spatially explicit biophysical models that depend largely on land use maps, topography, and regional climate data. The Sediment Retention model calculates annual soil loss (in tons ha⁻¹ yr⁻¹) based on the Universal Soil Loss Equation (USLE), which estimates overland erosion as follows (Renard et al 1991):

$$\begin{align*}{\text{USLE}} = R \cdot K \cdot LS{\text{ }} \cdot C \cdot P\end{align*}$$

where R is the rainfall erosivity factor (MJ mm ha⁻¹ h⁻¹ yr⁻¹), K is the soil erodibility factor, LS is a slope-length and steepness factor (unitless), C is a mitigating cover-management factor (unitless), and P is a support practice factor (unitless). To calculate the rainfall erosivity factor (R) for 2010 and 2020, we downloaded annual precipitation data at 1 km resolution (Thornton 2023). The statewide total precipitation in 2010 (∼25 inches) was more than 50% higher than the total in 2020 (∼16 inches) (Monserrat 2023). The USLE approach is highly sensitive to all model inputs, so to minimize the impacts of extreme values in a single year, we calculated the average precipitation over a five year period for each year (2008–2012 and 2018–2022) (Wang et al 2022). We then applied the formula developed by Renard and Freimund (1994) to calculate rainfall erosivity from the annual averages.

We obtained the soil erodibility data (K) layer from the Soil Survey Geographic Database (SURGGO) (Soil Survey Staff 2005). The InVEST model calculates the LS factor based on a digital elevation model, which is also the reference for the resolution of the model outputs. We used the 30 meter resolution digital elevation model from the Shuttle Radar Topography Mission, available globally (Farr and Kobrick 2000). The practice factor (P) incorporates practices such as contouring and other cross-slope erosion control practices. The only spatial data available regarding contour farming on California farmland is a national inventory of the relative amount of contour farming in each county in 1992, which shows virtually no contour farming in California (U.S. Geological Survey 2023). Additionally, according to a dataset provided to us by the USDA through a Freedom of Information Act request, no payments were made to farmers in California from the Conservation Stewardship Program or the Environmental Quality Incentives Program for either contour farming or strip-cropping for any year from 2010 to 2020. For these reasons, we assume little to no contour farming is occurring in California and thus P values are removed from the model, i.e. set to 1 for all land cover types.

The cover management factor, C, is the critical component for this analysis since it parameterizes the impact of land use and land management on soil erosion. The C factor incorporates erosion potential related to vegetation cover, crop type, tillage practices as well as different land cover types such as agriculture vs. wetlands, forest, or grassland (Panagos et al 2015). Lower C-factor values indicate lower erosion potential due to land cover or management. There is no spatially-explicit data regarding tillage practices in California, so c-factor values are exclusively tied to land cover type. These values were identified from the literature (primarily from Panagos et al 2015, Wartenberg et al 2021) and matched to land cover classes (see supplemental table S1). The resulting avoided soil loss depends entirely on the presence or absence of natural land cover types within agricultural landscapes and estimated C factor values for different agricultural classes.

The InVEST model also calculates the avoided erosion, which is the amount of erosion that is prevented due to the C and P factors:

$$\begin{align*}{\text{AE}} = \left( {R \cdot K \cdot LS} \right) - {\text{USLE}}.\end{align*}$$

AE is highly dependent on R (rainfall's ability to cause erosion). Since the year 2010 was much wetter than 2020, we expected AE to be much higher in 2010 because the overall erosion was much higher. To account for this, we calculated a soil retention index (SR) as the ratio of avoided erosion to the total erosion that would occur without mitigating C or P factors:

$$\begin{align*}{\text{SR}} = {\text{AE}}/\left( {R \cdot K \cdot LS} \right).\end{align*}$$

Model results could not be validated against ground-truth data, so we only assess changes in SR in terms of percent change from 2010 to 2020 rather than absolute values (in tons) (Natural Capital Project 2022). The InVEST model's threshold flow accumulation parameter defines the number of upstream pixels that must flow into a pixel for the pixel to be part of a stream. We determined this parameter by running the Flow Accumulation tool in ArcGIS Pro software, then incrementally modifying the threshold in the display of values to generate a modeled stream layer (ESRI Inc. Redlands, CA). We then overlayed a national streams data layer to qualitatively assess the accuracy of the generated streams layer and determined that 1000 was the best parameter value. Following the completion of the modeling, we limited the results to the boundaries of the Important Farmland data layer, and we summarized results for the entire State of California and for each agricultural commission region.

The removal of CO₂ from the atmosphere is a critical component in mitigating climate change, and in fact most emissions scenarios in integrated assessment models include some form of carbon capture through conservation, restoration, and land management practices (Smith et al 2016a). SOC is the measurable component of soil organic matter (SOM) with estimates of 1500–2400 Pg of carbon stored in the first 1–2 meters of soil globally, more than twice the amount of carbon stored in the atmosphere (Beillouin et al 2022). SOM is also the primary indicator of soil fertility and therefore critical to agricultural sustainability (Schillaci et al 2017, Ma et al 2023). SOC is impacted significantly by land use management and climate. Intensive agricultural practices deplete SOC, limiting the amount of CO₂ that is withdrawn from the atmosphere and entering a feedback cycle in which more industrial agriculture results in more greenhouse gas emissions and less CO₂ removal. Conservation practices on farmland have the capacity to improve SOC storage, increasing agricultural land's ability to regulate climate and provide food security simultaneously.

A range of statistical methods have been used to perform spatial prediction of SOC including simple linear regression, generalized linear regression, kriging, artificial neural network and other machine learning methods (Subburayalu and Slater 2013, Hengl et al 2015, Somarathna et al 2016). RF modeling is a machine learning method that has several advantages including that it is insensitive to both overfitting and noise in the data, and there is no requirement regarding the probability distribution of the predicted variable (Hengl et al 2015, Huettmann et al 2018). Furthermore, machine learning techniques like RF can handle nonlinear relationships between predictors and predicted variables, which is important in measuring ecological processes and ES (Huettmann et al 2018, Manley and Egoh 2022).

To train our RF model, we used soil samples collected by the Rapid Carbon Assessment (RCA), a national effort by the Soil Science Division of the Natural Resource Conservation Service (NRCS) in 2010 to collect and measure baseline SOC stocks (Wills et al 2014). There are 354 sample points in California, each with SOC stock values (units MgC ha⁻¹ yr⁻¹) at various depths up to 100 cm. The quantity and spatial distribution of SOC depends on environmental variables such as climate, topography, soil properties, and land use (Grimm et al 2008, Wiesmeier et al 2011, Pouladi et al 2019). We selected predictive variables based on a review of previous work; the full list is shown in table 1. Average monthly minimum and maximum temperature data and total annual precipitation data were derived from the Daymet 1 km monthly summaries dataset (Thornton 2023). We aggregated the monthly temperature values into yearly averages for each year 2010 and 2020, and we calculated the total precipitation for 2010 and 2020.

Topographic variables included elevation from the SRTM digital elevation dataset, as well as secondary topographic variables such as slope, aspect, curvature, planar curvature, and profile curvature, all of which were derived using ArcGIS Pro software. Soil variables were obtained from the Soil Survey Geographic Database (SSURGO) (Soil Survey Staff 2005) and included soil order, hydrologic group, apparent electric conductivity, erodibility, and texture. For land cover data, it was necessary to distinguish perennial crops with deep roots (tree crops and vineyards) from crops that are typically tilled annually since deep rooted perennials have a greater capacity for SOC storage and erosion control (Jansson et al 2021). Therefore, we used our reclassified CDLs from 2010 and 2020 as input land cover maps (supplementary table 1). Temperature, precipitation, land use, and vegetation greenness (estimated using the Normalized Difference Vegetation Index) were the only explanatory variables that differ between 2010 and 2020. We used data for the year 2010 to train the model with the 2010 RCA SOC values at various depths, then used variable data from 2020 to predict SOC values in 2020.

All data processing was done using ArcGIS Pro, and the RF model training was done using the Forest-based Classification and Regression tool. All NoData values for predictive variables were interpolated by calculating the average value of pixels in a window of 3 × 3 neighboring pixels, a commonly-used data-filling scheme (Fan et al 2020). The assumption is that neighboring cells tend to have similar characteristics. Once the modeled SOC depth was determined, we reran our model after removing the outliers from the original training data (points where SOC stock values at 10 cm depth were >2 standard deviations from the mean) for a final training dataset with 343 points, with 10% of the points reserved for validation. We used the trained model to generate statewide maps of SOC stock at 10 cm depth in the years 2010 and 2020, then limited the results to include only farmland. We summarized data for the entire Important Farmland area in California, as well as at the farming area in each county.

We compared FP to ES delivery for each of California's five agricultural commission regions and at the statewide level on farmland. We summed the SR, SOC, and FP statistics for each region.

The greatest land cover changes in terms of area on California farmland involved the conversions of Shrubland and Grassland (figure 2). Grasslands experienced a net loss of 2.2 million hectares (Mha), representing a 36% loss of area from 2010 (table 2), due primarily to transition to Shrubland (1.5 Mha) but with additional losses to Orchards/Vineyards (568.4 kha), Forest (281.3 kha), and Crops (345.0 kha). Shrublands experienced a net increase of approximately 1.4 Mha (+65% from 2010), due primarily to conversion from Grassland but with additional conversion coming from Forest (226.8 kha). We drilled down to the original Cropland Data Layer land cover categories to identify the primary contributors to changes in Grassland and Shrubland and found that the single largest conversion from 2010 to 2020 was from 'Grassland/Pasture' to 'Shrubland' (1.5 Mha). This number is over eight times larger than the next largest conversion, which was from 'Grassland/Pasture' to 'Evergreen Forest' (179.1 kha).

Zoom In Zoom Out Reset image size

Orchards/Vineyards experienced the largest percent change in area from 2010 to 2020 with a net increase of 900.2 kha (102.2% increase). The change in Orchards/Vineyards was unique compared to all other conversions because it gained more than ten times as much land as it lost (table 2); essentially, there were very few farmers converting from Orchards/Vineyards operations to any other form of land cover. The largest conversions to Orchards/Vineyards came from the following Cropland Data Layer categories: Grassland/Pasture (20% of area increase), Alfalfa (14%), and Fallow/Idle cropland (8%). Further, the conversions of land to Orchards/Vineyards primarily increased the areas of the following Cropland Data Layer categories: Almond (37%), Grapes (18%), Walnuts (10%), and Citrus (9%).

Forest experienced a net increase of 152.8 kha (+9.3% from 2010) due to the conversions of Grassland, Shrubland, and Open land. There were also significant conversions of Forest to the same categories in some regions, but the overall change was expansion of Forest. Crops decreased by 146.3 kha (−11.9% since 2010) with significant conversions to Grassland and Orchards/Vineyards. Impervious land, which included low, medium, and high intensity developed land as well as barren land, experienced a net decrease of 84.8 kha (−36.3% from 2010), mostly lost to Shrubland. Open land lost 78.8 kha (−19.3%), converting mostly to Shrubland, Grassland, and Orchards/Vineyards. Water/Snow experienced very minor changes (2.8 kha), likely due to the addition of the Perennial Snow category in the Cropland Data Layer in 2020. Wetlands experienced a net increase of roughly 190 kha, primarily from Grassland but with additional contributions from Crops, Shrubland, and Orchards/Vineyards. The expansion is a 96.5% increase from 2010, almost doubling the area.

The top ten land cover conversions (in terms of area) for each agricultural region are shown in figure 3, and detailed county results are shown in table 3. The Grassland to Shrubland conversion dominates in the Northern agricultural commission region, especially in Tehama and Mendocino Counties. Further conversion of Grassland to Forest occurred in Mendocino County. The Sacramento Valley, Central Coast, and Southern agricultural commission regions are also dominated by the conversion of Grassland to Shrubland, though to a lesser degree compared to the Northern region. The San Joaquin Valley agricultural commission region, though, is much more diverse in its land cover changes from 2010 to 2020. Large areas of Grassland and Crops were converted to Orchards/Vineyards, with additional exchanges between Grassland and Crops. Conversions to Orchards/Vineyards were concentrated in Fresno, Kern, and Tulare Counties.

Zoom In Zoom Out Reset image size

The statewide production of milk/dairy products, cattle/calves, grapes, and almonds increased between 2010 and 2020, and strawberry production decreased (figure 4). Almond production increased by 89.9%, far surpassing all other commodity changes, from 1.64B lbs. in 2010 to 3.12B lbs. in 2020. All counties that reported almond production reported increased production, with the largest increase being reported in the San Joaquin Valley agricultural commission region (figure 4). Stanislaus and Fresno Counties in San Joaquin Valley increased their production by +356 M lbs. and +342 M lbs., respectively. The Sacramento Valley region also boosted almond production considerably (+118% from 2010). According to our land-use change analysis, this increase in production is associated with a 106.6% Increase (+374 kha) In land taken up by almond groves. There were no reports of almond production in 2010 or 2020 in either the Central Coast or the Southern regions.

Zoom In Zoom Out Reset image size

In California, agricultural production is dominated by the dairy industry in terms of yield. In 2010, 41.9 billion pounds (B lbs.) of milk and milk fat were produced, and in 2020, this rose to 42.9 B lbs., an increase of 2.5%. Despite this increase in yield, 17 out of 29 reporting counties experienced a net decrease in the number of dairy cows in California. In 2020, there were almost 61 k fewer cows producing milk and dairy products according to county reports. Much of the change is a result of decreases in the Southern agricultural commission region (−72 k cows) and the Sacramento Valley agricultural commission region (−24.7 k cows). The largest increase was observed in the San Joaquin Valley agricultural commission region (+30 k). Twenty-two farming counties reported no dairy cows in 2010 or 2020.

Cattle and calf production increased from 2.6 B lbs. in 2010 to 2.7 B lbs. (+1.9% from 2010), with an associated increase in the number of beef cows (+16.5 thousand cows). Forty-four counties reported some number of beef cows (between 400 and 70 thousand) in either year. The largest increase in beef cows was in the San Joaquin Valley agricultural commission region (+31 thousand) due primarily to increases in Tulare County. The Sacramento Valley region also had an increased beef cow count (+21.7 thousand). The largest decrease was observed in the Northern region (−24.5 thousand) due primarily to a drop in Modoc County.

Forty-two counties reported their grape production in 2010 and 2020. According to county reports, there was a 2.04B lbs. increase in grape production overall (+19.4%), with the largest increase occurring in the San Joaquin Valley region. Fresno County and Kern County experienced the largest increases (1.1 B lbs. and 619.8 M lbs., respectively). The Sacramento Valley region also experienced a major increase in grape production. The largest decrease occurred in the Central Coast agricultural commission region. Statewide, increases in grape production correspond to a 132.1% increase in land use by grapes according to our land use change analysis, from 151 kha in 2010–350 kha in 2020.

Despite being one of the top commodities in California, only 14 out of 52 farming counties in California reported strawberry production in 2010 and 2020, and no strawberry production was reported in the Northern region. During this time, counties reported an overall decline in strawberry production of 218.3 M lbs. (−8.4%). The greatest declines occurred in the San Joaquin Valley and Sacramento Valley regions, with the largest contributors coming from declines in Ventura (245.2 M lbs.), Monterey (−110.2 M lbs.) and Santa Cruz (−108.8) Counties. According to our land-use change analysis, the decline in production was associated with a 94% decrease (−12.3 kha) in land used for strawberries.

In terms of production value, the five commodities yielded California farmers $13.97 B in 2010 and $22.3 B in 2020, an increase of almost 60%. This increase is largely due to the 98% increase in total almond production value, from $2.8 B to $5.6 B. Despite large declines in production of strawberries, the production value of strawberries increased by $192 M (+10.7%).

In 2010, our model predicted a larger amount of overall erosion compared to 2020 due to the increased rain in the 2008–2012 period compared to 2018–2022. In 2010 the total bare soil erosion (RKLS) was 124.1 B tons of soil, and in 2020 it was 98.2 B tons (−20.9%). The avoided erosion in 2010 was 122.2 B tons of soil, resulting in a SR index of 0.985. In 2020, the avoided erosion was 96.9 B tons of soil, resulting in an SR of 0.987, a 0.12% increase from 2010 statewide. In both years, the highest erosion occurred in the Northern agricultural commission region where precipitation was greatest, though this was the only region that experienced a decrease in SR (0.14% decline). The largest increase in SR occurred in the San Joaquin Valley agricultural commission region (1.17%), followed by the Southern region (1.05%) and Sacramento Valley (0.9%).

The SOC prediction model that returned the best results predicted SOC at 10 cm depth and had a training R² of 0.902 and a validation R² of 0.317, with a low standard error (0.01) and a significant p-value (0.00). The predictor variables' importance is listed in table 1, which shows that precipitation had the highest impact on SOC results, followed by maximum temperature and soil electrical conductivity. California's annual precipitation in 2010 was ∼25 inches and the average monthly maximum temperature ranged from 7.1 °C to 33.7 °C, and in 2020 the annual precipitation was ∼16 inches with maximum monthly temperatures between 9.4 °C and 33.7 °C (Monserrat 2023). Based on our prediction model, the total SOC stock in the top 10 cm of soil on California farmland in 2010 was 492. million metric tones C (MtC), or approximately 38.5 tC/ha. In 2020, SOC stock increased by 2.4% to 504.1 MtC, or 39.4 tC/ha (figure 3). The largest increase in SOC storage occurred in the Northern agricultural commission region (+7.14%), followed by the Sacramento Valley (+1.93%) and Southern regions (+1.49). The only region that experienced decreased SOC storage was the San Joaquin Valley region (−0.78%), primarily due to decreases in Tulare, Madera, and Mariposa Counties. In both years, the ranking of SOC regions was the same: San Joaquin Valley contains the largest SOC stock, followed closely by the Northern region, then Sacramento Valley, Southern, and finally the Central Coast.

We calculated the average percent change in SR and SOC for each of the top ten land cover transitions in terms of area. The results are listed in table 4. On average, SR decreased for three land cover transitions: grassland to orchard/vineyard (−8.63%), grassland to open land (−7.09%), and shrubland to grassland (−0.4%). The remaining largest transitions generated increased SR on average, with the greatest increases coming from the conversions of crops to grassland (+3.5%), crops to orchard/vineyard (+3.08%), and grassland to shrubland (+1.04%). SOC increased for all of the top ten land cover transitions, with the highest increases coming from the conversions of shrubland to grassland (+10.48%), grassland to crops (+5.03%), and grassland to open land (+4.57%). The conversion of crops to orchard/vineyard also generated a relatively high increase in SOC on average (+4.06%).

Figure 4 shows the increase and decrease in the delivery of ES and FP. In the Northern region, there was a major increase in SOC storage (+7.14%), a minor decrease in SR (−0.14%), increased production of milk and almonds, and large areas of conversion between shrubland, grassland, and forest. In the Sacramento Valley region, there was increased SOC storage (+1.93%), increased SR (+0.9%), and a large increase in almond production, as well as improvements in grape and beef production. In the Central Coast, there was a slightly lower increase in SOC storage (+1.24%), a very slight increase in SR (+0.21%) and declines in all commodity production except for milk. The San Joaquin Valley region is the only region that experienced a decline in SOC storage (−0.78%) with an increase in SR (+1.17%), while increasing the production of all commodities except strawberries and undergoing a number of large land cover conversions including conversion of grassland and crops to orchards/vineyards. The Southern region experienced a moderate increase in SOC storage (+1.49%), slight increase in SR (+1.05%), and decline in the production of all commodities except strawberries.

The encroachment of shrubland on grassland from 2010 to 2020 was a major source of change throughout California. We found that this change resulted in generally increased SR and SOC throughout California, but other studies have shown conflicting results on this. Li et al (2016) found that SOC content generally increases with shrub encroachment on grassland, but that the direction and magnitude of the SOC change varied depending on soil, climate, and plant cover. Another long-term study found that SOC decreased significantly with the conversion of grassland to shrubland, and much more so with the conversion of grassland to crops (Qiu et al 2012). According to the Panagos et al (2015) study of soil erosion cover management factors at a regional scale across Europe, wetlands, forest areas, and shrublands are better at slowing down rainfall compared to other land cover types, allowing water to infiltrate and preventing soil erosion. We used the cover management factors developed by the Panagos et al study, which means the replacement of grassland with shrubland led to enhanced soil erosion control.

Ecological transition from perennial arid or semiarid grasses to woody shrublands is typically interpreted as a form of land degradation and may be the result of the historic California drought from 2011–2017 (Peters et al 2013). The replacement of grasslands with shrublands is often the result of multiple intersecting social and ecological factors including drought, changing fire patterns, overgrazing livestock, soil erosion, non-native plant invasion, and farm abandonment (Queiroz et al 2014, Bestelmeyer et al 2018). For example, in Tehama County, where the conversion from grassland to shrubland is extensive (61.9% of all changed land and 24.9% of all farmland in the county), a long history of fire suppression and heavy grazing has vastly changed the landscape over the past century (Tehama County RCD 2006).

Enhanced SOC and SR may also be attributed to the increase in forest and nut trees in California from 2010 to 2020. Studies have shown afforestation to be an effective land management tool for enhancing carbon sequestration, though the impacts on organic carbon stocks depend on the tree species and litter type (Ma et al 2016, Smith et al 2016b). The increase in forest is largely a result of the conversion of grassland and shrubland to forest (table 4). These conversions may be attributed to recovery and re-growth after wildfires, expansion of timberland, or the encroachment of forest into land that was previously used for grazing (Weatherford 2009, Lutz 2024). The expansion of nut trees is a result of FP, and is discussed in the section below.

There were also large areas where forest were converted back into shrublands, likely the result of wildfires in California. The 2018 Camp fire in Butte County was the most destructive fire in California history, with a burn area of over 62 kha, visible in our land use change analysis as the conversion of forest to shrubland and grassland, indicating the beginning of regrowth in 2019 and 2020. Another large region of forest-to-shrubland conversion is situated in Yuba County, next to Collins Lake, corresponding to the 2017 Cascade Fire, which burned roughly 6.5 kha over several months (Yuba County 2021). Other regions with large areas of forest-to-shrubland conversion occurred in Sikiyou County (20.2 kha), Tehama County (19.1 kha), and Kern County (17.2 kha). California wildfires can significantly diminish soil erosion control and soil carbon storage in post-fire ecosystems, and may even result in ecosystem disservices in the form of flooding and mudslides (Busari et al 2024). Thus, these events likely constrained the changes in SOC and SR in our study.

The increase in wetlands we observed may have resulted in improved SR in our analysis since wetlands help prevent erosion (Panagos et al 2015), but studies have shown that wetlands may decrease SOC stocks unless carefully managed (Xu et al 2019). This land use change may be the result of restoration efforts or sustainable agricultural practices on farmland, or it may result in part from varying accuracy in the 2010 land use map compared to 2020. The largest wetlands expansion was in Merced County and consisted of a 15.2 kha conversion of grassland to 'herbaceous wetlands' (over 3% of the total farmland area in that county), primarily in and around what is known as the Grassland Ecological Area. The Grassland Ecological Area is the largest remaining wetland complex in the western United States according to the Grassland Water District, the managing organization of the wetlands (GWD 2017). However, general manager of the Grassland Water District indicated there has been very little expansion of wetlands in the Grassland Ecological Area over the past several decades, and that these seasonal wetlands predate the 2010 land cover map (personal correspondence with Ric Ortega, 5/30/2024). Future research should attempt to identify the discrepancies between the National Cropland Data Layer datasets from each year, especially with respect to seasonal wetlands.

The most striking change in FP from 2010 to 2020 is the increase in almond production, which is mirrored by the conversion of grassland to orchards, especially in the San Joaquin Valley agricultural commission region. This region was also the only region that experienced a reduction in SOC and the highest increase in SR. The increase in SR is likely the result of the planting of trees, which are perennials with strong roots that help to anchor the soil. We expected increased SOC since perennial trees are strong sequesters of carbon, though the lasting effects differ depending on climate, topography, and specific plant type (Smith et al 2016b, Dass et al 2018). Indeed, our results indicate that the conversion of crops to orchards/vineyard in general result in an increase in SOC. But other studies have shown that grasslands may be a more reliable and resilient carbon sink compared to forests and other trees, and in arid regions, shrublands may be a larger net carbon sink over time compared to grasses (Petrie et al 2015, Dass et al 2018). Conversion of annual crops, which release carbon and disturb soil at each harvest, to perennial almond groves may yield a net increase in SOC storage and erosion control over a period that is longer than our ten-year study. The decline in SOC may in the San Joaquin Valley is also minimal compared to the overall pool of SOC in that region, and may have more to do with the complexity of the landscapes and water flow, rather than the actual land cover type.

The increase in milk production was not mirrored by an increase in the reported number of dairy cows, which may be explained by the increased milk production efficiency. The dairy industry has been marked by several important developments throughout its history in the United States. The introduction of artificial insemination in the 1930s meant that dairy cows could be impregnated without the need for a mate onsite, and improvements in this practice have resulted in higher rates of conception; thus, more cows are producing milk (Smith Thomas 2019). However, the onset of genomic selection in the 21st century transformed the dairy industry during our period of analysis. Official USDA genomic evaluations were released between 2009 and 2016 for different breeds of cows, and genomic selection for production efficiency has meant that cows can now produce up to twice as much milk as they could before this practice was implemented (Wiggans and Carrillo 2022). The increased production in California in 2020, despite a decrease in the number of dairy cows, is likely the result of advancements in this genetic technology.

The loss of grassland and increased number of beef cows indicates an increased concentration of livestock production (high carrying capacity or overstocking) rather than an increase in land area used for grazing. Concentrated animal feeding operations have significant negative impacts on water quality, air quality, and greenhouse gas emissions, and have life-changing health impacts for surrounding communities (Gurian-Sherman 2008, Kravchenko et al 2018). With the 2023 Farm Bill and dozens of federal conservation programs currently under the microscope in the U.S., agricultural incentive programs have recently been analyzed and found to disproportionately benefit these kinds of large animal operations, which environmentalists and scientists argue is a misuse of conservation funding (NSCA 2021, Held 2022, Lenhardt and Egoh 2023).

A 2.38% increase in SOC storage on California farmlands may not, at first, appear significant, but every ton of carbon stored in the soil is equivalent to 3.67 tons of atmospheric carbon dioxide (Romm 2008). Thus, a 11.7 million MgC increase in SOC is the equivalent of removing an additional 47.5 million tons of CO₂ from the atmosphere; the emissions from over ten million passenger cars in one year (EPA 2016). Though some changes in SOC may be attributed to land use change driven by agriculture and other factors, our SOC model indicates that annual precipitation, maximum temperature, soil properties and vegetation greenness are stronger determining factors of SOC in California farmland. Thus, SOC is directly related to SR due to precipitation changes; less erosion means greater capacity to build up SOC from plant matter and less soil displacement (Fissore et al 2017, Hancock et al 2019). Year 2020 was warmer and drier, and we believe the lower overall erosion due to diminished rainfall (particularly in the north) resulted in greater SOC stores. With a warmer and drier future, it will become more critical to understand how land management practices can improve SOC in California farms.

Lal et al (2021) argued that the restoration of SOC on working lands requires 'system-based conservation agriculture', which includes multiple conservation practices that reduce soil disturbance and enhance soil fertility, such as allowing living roots to remain in the ground throughout the year, which facilitates the long-term storage of carbon in the soil (Jackson et al 2017). When SOC storage is improved through these practices, farmers may benefit economically because of increased crop yields and nutrient availability as well as through the sale of carbon credits (Davoudabadi et al 2021). Beyond farms, populations that are particularly vulnerable to climate change may also be considered beneficiaries of carbon storage on agricultural land (Bagstad et al 2014).

The regional differences in how SOC, SR, and FP changed from 2010 to 2020, as well as the ways in which land cover changed regionally, indicate the complex nature of farming in California. California is a vast state that provides a large variety of agricultural products and includes many ecosystems distinguished by their diverse climatological and biological regimes. It is also a diverse state in terms of farm ownership, farm labor practices, and policy, and each region faces different challenges. For example in the Central Coast, farmers face strict food safety regulations, high cost of rent, and supply chain pressures that make the transition to sustainable agriculture extremely difficult (Carlisle et al 2022). In San Joaquin Valley and Sacramento Valley, farmers face strict groundwater regulations and land subsidence, making the shift to sustainable agriculture a seemingly less urgent priority compared to more immediate challenges (Faunt et al 2024). A more regionally targeted approach to conservation incentive programs may help to overcome some of these barriers. Coordination and co-production of policy with agencies such as the agricultural commission regions in California may help facilitate a transition to sustainable farming while maintaining farmer trust and improving access to services.

There are several limitations to our analysis that constrain the accuracy of our results which should be interpreted with caution. The InVEST Sediment Retention model only accounts for overland erosion and not gully or bank erosion, nor does it include mass erosion events like landslides or mudslides (Natural Capital Project 2022). Landslides and mudslides can be significant sources of soil loss, particularly in areas where there is construction or wildfires (Ren et al 2011). Additionally, given the simplicity of the InVEST model, the results are sensitive to most of the input parameters (Natural Capital Project 2022). Both SR and SOC models are constrained by the accuracy of the input data, and though the land cover maps we used had high overall accuracy in California, one potential source of error is the misclassification of crop types or non-agricultural land use types. Further, the accuracy of our sediment retention model results was not validated due to a lack of ground truth data for the years we analyzed. Future work may incorporate a more detailed watershed-level analysis.

Finally, our study compares ESs generated by two different approaches: a biophysical model and a machine learning model. To our knowledge, there have been no studies discussing the merits or drawbacks of this type of comparison. Future work may incorporate updated machine learning methods and higher accuracy data if available. We believe the methods we used will benefit with enhanced technology. Higher resolution, open-access satellite imagery would improve the accuracy and accessibility of ES modeling globally, though multitemporal comparison will be difficult with future iterations of higher-accuracy data in comparison to contemporary or historical data.