Land Use and Land Cover Affect the Depth Distribution of Soil Carbon: Insights From a Large Database of Soil Profiles

Soils contain a large and dynamic fraction of global terrestrial carbon stocks. The distribution of soil carbon (SC) with depth varies among ecosystems and land uses and is an important factor in calculating SC stocks and their vulnerabilities. Systematic analysis of SC depth distributions across databases of SC profiles has been challenging due to the heterogeneity of soil profile measurements, which vary in depth sampling. Here, we fit over 40,000 SC depth profiles to an exponential decline relationship with depth to determine SC concentration at the top of the mineral soil, minimum SC concentration at depth, and the characteristic “length” of SC concentration decline with depth. Fitting these parameters allowed profile characteristics to be analyzed across a large and heterogeneous dataset. We then assessed the differences in these depth parameters across soil orders and land cover types and between soil profiles with or without a history of tillage, as represented by the presence of an Ap horizon. We found that historically tilled soils had more gradual decreases of SC with depth (greater e-folding depth or Z∗), deeper SC profiles, lower SC concentrations at the top of the mineral soil, and lower total SC stocks integrated to 30 cm. The large database of profiles allowed these results to be confirmed across different land cover types and spatial areas within the Continental United States, providing robust evidence for systematic impacts of historical tillage on SC stocks and depth distributions.


INTRODUCTION
Soils contain a substantial fraction of global terrestrial carbon (C) stocks, accounting for more than the C content of vegetation and the atmosphere combined (Jobbágy and Jackson, 2000;Ciais et al., 2013). While the spatial distributions of soil types, soil characteristics, and soil C (SC) stocks have been mapped at regional to global scales (Nachtergaele et al., 2010;Hengl et al., 2014), the role of the third dimension, depth, is increasingly recognized as an important (if indirect) control on SC stability (Schmidt et al., 2011). The vertical distribution of soil C with depth may hold critical information about the vulnerability of SC to land cover change, land management practices, and other disturbances (Rumpel and Kögel-Knabner, 2011).
Historical changes in SC stocks due to land management practices such as deforestation and agricultural tillage have been responsible for a significant fraction of historical anthropogenic carbon dioxide (CO 2 ) emissions to the atmosphere (Eglin et al., 2010;Houghton et al., 2012). Conventionally tilled soils can lose from 20 to 50% of SC relative to unmanaged or no-till practices (Mann, 1985;Senthilkumar et al., 2009;Bruun et al., 2015). Tillage reduces SC by disrupting aggregates (Six et al., 2000;Wright and Hons, 2005) and by increasing susceptibility to erosion (Harden et al., 1999;Doetterl et al., 2016). However, tillage also redistributes SC through the soil profile, potentially offsetting SC losses in upper layers by increasing storage in deeper soil layers (Baker et al., 2007;Angers and Eriksen-Hamel, 2008).
While comparisons among sites in agricultural experiments have established broad patterns of SC losses under different land management scenarios, understanding how SC sensitivity to disturbance varies with depth across climates, soil types, and land cover types has been challenging. Individual agricultural experiments have been sparsely distributed and can require long-term monitoring, limiting the availability of results that are applicable across large spatial scales. Two-dimensional scaling approaches involving geostatistical modeling such as SoilGrids (Hengl et al., 2014) and the Harmonized World Soil Database (HWSD; Nachtergaele et al., 2010) typically lack detailed information about the depth distribution of SC, yet it is clear that subsoils can store as much as or more C than topsoils (Jobbágy and Jackson, 2000).
Previous studies of SC depth distributions have provided important mechanistic insights (Sun et al., 2020), but generalization of findings across soil types, ecosystems, and spatial scales can be limited by the number and distribution of soil profiles. Emerging databases of soil profiles can now facilitate more comprehensive studies of variations in soil properties across depth, particularly when coupled with sitespecific and geospatial data on land use and agricultural history. One such database, the International Soil Carbon Network (ISCN) database (Nave et al., 2016), contains more than 70,000 soil profiles and includes a wide variety of depth-defined soil parameters, such as organic and total C content, total nitrogen, bulk density, particle size distribution, horizon designations, and various chemical factors all specified to depth layers. A key advantage of this dataset compared to gridded spatial products such as SoilGrids and HWSD is the site-specific nature of the individual measurements. The ISCN database contains data at the individual soil profile level, reported for the specific layers, depths, and horizons that were measured in the field. As a result, analyses using ISCN database profiles build on individual measured depth profiles and can be matched directly to specific locations without the use of geostatistical scaling methods.
We conducted a statistical analysis of SC depth distributions using the ISCN database to determine how tillage affected SC stocks and their vertical distributions across a range of climates and environmental conditions. We harmonized the heterogeneous depth sampling of profiles in the ISCN database by fitting observed C densities to a standard relationship (exponential decline with depth), facilitating statistical analysis of key depth distribution parameters across the database. We applied this technique to the ISCN database with three key objectives: (1) to determine the feasibility of applying a common depth function across variations in climate and soil types in a large soil profile database; (2) to evaluate the influence of past agricultural tillage on SC stocks and depth distributions across climate, soil order, and environmental conditions; and (3) to publicly release a dataset of SC depth distribution parameters across the ISCN database to facilitate future studies of threedimensional SC stocks and vulnerabilities. We hypothesized that a history of tillage would be associated with lower surface SC concentrations, slower declines in SC concentrations with depth, and lower topsoil SC stocks.

International Soil Carbon Network Database
For this analysis, we used the 3rd generation ISCN database (Nave et al., 2016). The database has layer-and profile-level attributes spanning a range of geographic, physical, chemical, and ecological variables, synthesized across 39 individual dataset contributions. The available data for any given profile or layer depends on the measurements reported by each contributor. Some profiles in the ISCN dataset were sampled by distinct horizons and others were sampled by regular depth increments, depending on the approach of the data originators. Thus, we define layers in this manuscript as sampled depth increments, which may or may not correspond to distinct soil horizons. In our study, we used this depth-explicit data to make statistical comparisons of C concentration and C density. Carbon concentrations included analytical determination of organic C (37,261 profiles) and/or total C (16,401 profiles); for calculating the carbon densities we used the reported layer thickness, C concentration (preference given to organic C), and bulk density. Bulk densities supplied by ISCN data contributors were computed by a range of techniques; in our C density computations we preferentially used the fine earth (<2 mm) bulk density (bd_samp in the database); if that was not available we used either the whole soil bulk density or whichever variant was reported by the data contributor. The database also included soil taxonomic information for each profile and for each layer; we used the available taxonomic information to group soil observation points by USDA taxonomic order. In any given profile, we used the presence of an Ap horizon to identify that it had a history of tillage, following Nave et al. (2013Nave et al. ( , 2018Nave et al. ( , 2019. Also following methods of Nave et al. (2013Nave et al. ( , 2018Nave et al. ( , 2019, for each geo-located soil observation within the continental United States, we derived land cover classification from the satellite-derived, 30-m resolution National Land Cover Dataset (NLCD), including products from 1992 (Vogelmann et al., 2001), 2001(Homer et al., 2004(Fry et al., 2011, and 2011 (Homer et al., 2015). We chose the year of NLCD data based on the reported date of the soil observation so that land cover classifications would be consistent with when the soil was collected, in case land cover in that location had changed between NLCD collections. Soil observations measured prior to 1996 were matched with the 1992 NLCD data. Profiles measured from 1996 to 2003 were matched with 2001 NLCD data. Profiles measured from 2004 to 2009 were matched with 2006 NLCD data, and profiles measured after 2009 were matched with 2011 NLCD data. Because the majority of fitted profiles were located in the Continental United States, and the NLCD dataset is limited to the Continental United States, we focused our analysis on that region.

Fitting Depth Attenuation of Soil C Content
Based on numerous studies in which SC was reported to decline exponentially with depth under various conditions (Davidson et al., 2006;Meersmans et al., 2009;Sun et al., 2020), we fit C content data in each profile to a continuous piecewise function in which C declines exponentially with depth until reaching a depth below which carbon remains constant, following Rosenbloom et al. (2006): where C is C content [either C concentration (%) or C density (gC cm −3 )], z is depth (defined as increasing relative to the surface of the mineral soil), C s is C at the mineral soil surface, Z * is an empirical depth scaling parameter, Z min is the depth at which C stops decreasing and stabilizes to a constant value, and C deep is equal to the C at depth Z min , i.e., C deep = C s e −Z min /Z * . Figure 1 shows example soil profiles, fits, and parameters to illustrate the approach. Eq. 1 was fit to each soil profile using non-linear least squares regression (the curve_fit function in the python scipy package version 1.1.0), including C s , Z * , and Z min as free parameters. Fits were attempted on all profiles with C data for at least three mineral soil layers in the ISCN database. First, organic horizons were determined using a cutoff of %C > 20% based on United States Department of Agriculture Natural Resources Conservation Service (NRCS) criteria for distinguishing organic and mineral horizons (Soil Survey Staff, 2014), and profile depths were adjusted so that the top of the first mineral horizon was defined to z = 0. Eq. (1) was then fit to each profile, using %C and/or C density (gC m −3 ) depending on which metrics were available for each profile. C density was previously calculated across the ISCN database by multiplying %C by bulk density in each layer, for profiles where measurements of bulk density were available from data contributors. Organic C and total C were reported separately in the ISCN database. We used organic C when it was available, and total C otherwise. Note that total C may include significant amounts of inorganic C in carbonates in some soils, particularly those with higher pH (Batjes, 1996;Lal and Kimble, 2000). The analysis was repeated using only profiles with organic C measurements, which yielded similar results (see Supplementary Figures). After fitting, we filtered the resulting database of depth attenuation parameters to limit the statistical analysis to high-quality fits. See Table 1 for a summary of profiles omitted from the analysis. Confidently interpreting a profile fit required at least five depth layers to be included in the regression, so profiles with less than five mineral soil layers with reported SC data were excluded. The primary analysis excluded profiles with R 2 < 0.9 between the fit and observed C in the mineral horizon profile to focus the analysis on the parameters determined with the best confidence. To evaluate whether eliminating these profiles affected the results, we repeated the analysis using an R 2 threshold of 0.75, which yielded similar results (see Supplementary Figures).
The fitting approach condensed each profile's set of soil layers to a common set of three parameters (C s , Z * , and Z min ) along with C deep , which was calculated from the other three parameters. This facilitated a harmonized synthesis and statistical analysis of patterns of SC attenuation with depth across the large ISCN database of soil profiles. For fits using C density (which incorporated bulk density), Eq. 1 was integrated vertically to calculate total SC stock in the profile to any specified depth (assuming no changes in C deep with depth): For this analysis, we calculated SC stocks to 30-cm depth to examine landscape variations in topsoils. However, future analyses using this dataset could extend to deeper depths as long as shallower profiles were excluded to avoid extrapolation beyond measured profile depths.

Statistical Analyses
We focused our statistical analyses on the distribution of soil profile parameters (from the fitting procedure described above) as they varied across soil order, spatial areas, and landcover types. We used a Chi 2 test (R command: chisq.test) to determine differences in the distributions of the binned response variables (Z * , Z min , C surf , C deep , C stock ). For these non-normally distributed response variables, we determined medians and 95% confidence intervals around the medians using the mean of 10,000 replicates generated using bootstrap resampling (R command: sample). Differences in mean response variables among the land cover classes by Ap horizon were tested using analysis of variance (R command: aov) and post hoc differences were determined using Tukey's Honest Significance Difference test (R command: TukeyHSD); we tested for interactions between land cover class and Ap horizon. To test for differences between forested and unforested land cover classes, we re-grouped the samples into forested and unforested types. In addition to examining distributions pooled across the whole dataset, we controlled for spatial variations in climate and other drivers by dividing the Continental United States into a 0.5 • (latitude/longitude) spatial grid. Profiles within each grid cell were further divided based on soil order, and depth distribution parameters were compared between the sets of profiles with or without an Ap horizon within each group sharing a grid cell and soil order using a paired t-test (R command: t.test, paired = T). Statistical analyses were done using R statistical software (R Core Team, 2019), and  "Rejected" rows show the number of profiles that were removed from the analysis for the reason in the column. "Rejected (% of total)" rows show rejected profiles as a percentage of the whole dataset. Rejected (% of remaining) rows show rejected profiles as a percentage of the profiles remaining in the previous step. "Remaining" rows show the number of profiles remaining after cumulative rejections.

RESULTS
The ISCN database included 71,134 total profiles, but not all were suitable for the depth profile fitting approach ( Table 1).
16% of profiles lacked any sampled mineral soil layers with C concentration measurements, and 56% lacked any sampled mineral soil layers with both C concentration and bulk density information necessary for determining C density. Confidently interpreting a depth profile required at least five measured layers, removing an additional 27% of profiles for %C and 14% of profiles for C density. Some of the remaining profiles were excluded because they did not yield successful fits due to nonconvergence of the fitting algorithm or yielded unrealistic results Frontiers in Environmental Science | www.frontiersin.org (e.g., Z * < 0). These steps left 38,612 C concentration profiles (54% of all profiles) and 20,151 C density profiles (28%) on which further analysis was possible. Removing profiles with R 2 < 0.9 further reduced the number of included C concentration profiles by 23% of successfully fitted profiles, and C density profiles by 28% of successfully fitted profiles. We repeated the analysis using an R 2 cutoff of 0.75, which did not substantially affect the results (see Supplementary Material). We also repeated the analysis using only profiles that reported organic C, which yielded similar results. A large fraction of successful fits had high R 2 values (Figures 2A,B). Z * values had an approximately log-normal distribution across the database, peaking between 10 and 30 cm for both C density and C concentration but with a long tail indicating significant numbers of soil profiles with more gradual declines in C concentration with depth (Figures 2B,C). C s values had similar log-normal distributions centered around 1-2% C with a small number of profiles having surface C concentrations greater than 10% for C concentration, and centered around 0.01-0.03 g C cm −3 with a small number of profiles having surface C densities greater than 0.15 g C cm −3 .
The vast majority of fitted profiles in the ISCN database were located in the Continental United States (Figure 3). This was largely due to the prevalence of soil profiles from NRCS surveys, which made up approximately 87% of the analyzed data. Thus, we focused our analysis on profiles within the Continental United States (Figure 4). C stocks calculated to 30 cm ranged up to 15 kg C m −2 , with the highest stocks occurring in the Northwest and central plains regions. The deepest SC profiles, as indicated by Z * and Z min values, were located in the central plains region. C s and C deep did not have clear continental-scale spatial patterns.
A key hypothesis of this analysis was that historical disturbance of the soil profile through tillage would influence SC stocks and their vertical distributions. Thus, we analyzed how soil depth distribution parameters and integrated profile C stocks varied between soil profiles with and without an Ap horizon indicating a history of tillage ( Figure 5). The distribution of Z * differed significantly between soils with and without an Ap horizon (Chi 2 = 100,899, d.f. = 187, P < 0.0001). Pooled over all soil orders, the distribution of profiles with an Ap horizon was skewed toward higher Z * values (median: 26.0 ± 0.9 cm) compared to the distribution of profiles without an Ap horizon (median: 18.7 ± 0.5 cm). This indicates that SC concentrations decreased more slowly with depth in profiles with a history of tillage. Because deep-rooted plants and thick SC-rich horizons are a key characteristic of Mollisols, we separately investigated depth parameter distributions within Mollisol soil profiles. These distributions yielded similar results to all soils (Chi 2 = 17,074, d.f. = 187, P < 0.0001), with the distribution of Z * in Mollisol profiles with Ap horizons skewed toward higher Z * than non-Ap profiles (Medians: Has Ap = 42 ± 1 cm; No Ap = 31 ± 2 cm). The Z * distribution for Mollisols was centered at deeper values (median: 37 ± 2 cm) than the Z * distribution for all soils (median: 22 ± 1 cm; Chi 2 = 161380, d.f. = 187, P < 0.0001).
Differences in soil profile parameter distributions with tillage history were also reflected in C s , Z min , C deep , and total profile C stock. The distribution of C s was narrower and had a higher fraction of low C s values for profiles with Ap horizons compared to profiles without Ap horizons (Chi 2 = 168,000, d.f. = 199, P < 0.0001), indicating a lower incidence of high surface C concentrations, which was consistent with physical mixing of a SC-rich surface layer into deeper layers. Z min distributions were skewed toward deeper Z min for profiles with Ap horizons (Medians: Has Ap = 58 ± 2 cm, No Ap = 43 ± 1 cm; Chi 2 = 42551, d.f. = 99, P < 0.0001) and were much broader and deeper for Mollisols than for all soils (Chi 2 = 70,900, d.f. = 99, P < 0.0001). C deep distributions were similar between Mollisols and all soils, and distributions with Ap horizons had lower median C density than soils without Ap horizons (Medians: Has Ap = 0.0024 ± 0.0000 g cm −3 , No Ap = 0.0033 ± 0.0000 g cm −3 ), while distributions of C deep for profiles without Ap horizons had a longer tail of higher values (Chi 2 = 276628, d.f. = 345, P < 0.0001). The distributions of SC stocks to 30 cm were centered around higher stocks for Mollisols than for all soils (Medians: Mollisol = 6.3 ± 0.3 kg C m −2 , All soil = 4.7 ± 0.1 kg C m −2 , Chi 2 = 164062, d.f. = 335, P < 0.0001), and soils with Ap horizons were somewhat skewed toward lower SC stocks and contained fewer profiles with high SC stocks (Medians: Has Ap = 4.0 ± 0.1 kg C m −2 , No Ap = 5.2 ± 0.1 kg C m −2 ; Chi 2 = 94397, d.f. = 335, P < 0.0001).
Land cover affects SC stocks and their depth distributions and may interact with land use history (for example, previously tilled soils in some areas have been reforested). We used NLCD data to identify land cover classes for each profile and calculated mean profile parameters for profiles with or without Ap horizons within the most common land cover types (Figures 6A-C). Generally, profiles with Ap horizons had lower C s [F (1 , 9044) = 1160, P < 0.0001] and lower SC stocks over 0 to 30 cm depth [F (1 , 9044) = 457, P < 0.0001]. Grasslands had the highest Z * values and deciduous forests had the lowest. Forests had lower Z * [F (1 , 9056) = 663, P < 0.0001] and higher C s [F (1 , 9054) = 358, P < 0.0001] than non-forest landcover types. However, profiles with Ap horizons had similar mean C s for all landcover types, with the exceptions of mixed and evergreen forests. Note that non-crop profiles with Ap horizons could represent previously tilled soils that were subsequently converted to other land cover types (e.g., reforested), but may also include agricultural sites that were not detected in the satellite-based landcover products.
Ap horizon differences were also consistent when soils were divided by soil order (Figures 6D-F). Across all soil orders, the presence of an Ap horizon increased Z * by 9-37%. This trend was statistically significant (P < 0.01) except in Aridisols, Entisols, Spodosols, and Ultisols ( Figure 6D). The presence of an Ap horizon decreased C s by 12-59%, which was statistically significant except in Aridisols and Vertisols (Figure 6E). The presence of an Ap horizon significantly decreased SC stock to 30 cm by 15-46% (P < 0.01) except in Aridisols, Entisols, Spodosols, and Vertisols, where there was no significant effect of the Ap horizon (P > 0.10, Figure 6F).
To control for climatic and regional factors, we divided the Continental United States into a 0.5 • (latitude and longitude) spatial grid and calculated pooled differences in integrated SC stocks to 30 cm between profiles with or without Ap horizons within the same soil order for each grid cell (Figure 7). This grid  spacing was likely too coarse to control for many differences in soil forming factors but did allow for the control of continentalscale climate and spatial biases. Note that the purpose of this stratification was not to conduct an explicit spatial or regional analysis, but only to facilitate comparisons between profiles while mitigating confounding effects due to spatial variations in climate and soil order. When combining soil orders within each grid cell, profiles with Ap horizons had significantly lower SC stocks than profiles without Ap horizons (t = −8.23, d.f. = 852, P < 0.0001). The distribution of differences between profiles from the same soil order and grid cell with and without Ap horizons was skewed toward negative values, with a long tail of cells where Ap horizons were associated with more than 5 kg C m −2 lower SC stocks ( Figure 7C). When expressed as a percentage difference in SC stock between profiles with and without Ap horizons, the skew was more apparent, with a median percent difference of 17% lower SC stocks in profiles with Ap horizons compared to profiles without Ap horizons within a grid cell and soil order (Figure 7D). A quarter of grid cells had more than 40% lower SC stocks in profiles with Ap horizons than profiles of the same soil order without Ap horizons. SC stocks in soils with Ap horizons had a significant linear relationship [Intercept: 2.77 (95% CI: 2.4, 3.2) kg C m −2 ; Slope: 0.35 (95% CI: 0.30, 0.40); P < 0.0001, R 2 = 0.18, RMSE = 3.3 kg C m −2 ] with SC stocks in soils without Ap horizons within the same grid cell and soil order (Figure 8). The slope of the relationship being below 1 indicates that profiles with Ap horizons had systematically lower SC stocks than profiles without Ap horizons and that the potential for topsoil SC loss under tillage increases with the magnitude of untilled SC stock.

DISCUSSION
Large datasets of soil profiles offer opportunities to examine spatial variations in biogeochemistry using a number of statistical or mathematical approaches, and we found that applying a common analysis to SC data elucidated patterns of C density and concentration with depth. The high R 2 values for the vast majority of soil profiles with enough depth points to fit (R 2 ≥ 0.9 for more than 29,000 SC concentration profiles and more than 14,000 SC density profiles; Table 1 and Figure 2) indicated that this method was successful in determining profile parameters across the dataset. This harmonized approach addressed difficulties in comparing vertical patterns across a dataset in which the depths and layers measured were variable, and facilitated analysis of total C stocks, surface C densities, and the depth scales of the profiles in an efficient statistical framework. This approach provides advantages over gridded or spatially modeled soil datasets such as the HWSD and SoilGrids. In particular, the parameters in our analysis are traceable back to individual measured profiles rather than being averaged or statistically modeled using spatial relationships or environmental covariates. This allowed depth patterns to be investigated for individual sites in the context of their profile properties and land use histories.
Depth characterization of SC by the formulation of Eq. 1 elucidates some important attributes that may be generalizable and useful in further analyses. Parameters Z min , C deep , C s , and Z * likely inform different aspects of C cycling. For example, Z min likely reflects a threshold depth above which carbon retention is dominated by biotic processes such as root input, microbial processing, and stabilization of C in the biotic medium (White et al., 2012;Lawrence et al., 2015). Below Z min , C concentration is more likely dominated by hydrology and mineralogy. In other words, Z min represents a "biotic" vs. "quasi-biotic" transition. Such thresholds for C storage and turnover were discussed by Lawrence et al. (2015) as being closely coupled to surface area of clays and oxides, reaching a saturation level that changes with mineral transformations occurring in the B horizon subsoil. White et al. (2012) also reported that such a threshold in soil depth exists for soil mineral transformations because of the abrupt changes in hydrology, porosity, and redox reactions. Significantly, the parameter C deep and its upper depth limit Z min reflect a metric for detecting a fraction of C stocks that has been insulated from surface processes under historical conditions, and that may be transformed under new conditions that alter Z min , such as an increase in rooting depth. Thus, upon shifts in ecosystem and/or management factors that change Z min , deep carbon may be more readily decomposed (Hicks Pries et al., 2018). Meanwhile it is clear that the C s parameter is sensitive to land cover and varies spatially. C s was systematically higher in forests than in grasslands, crops, or shrublands, most likely due to the higher aboveground fraction of litter deposition in forests compared to the larger fraction of NPP allocated to fine roots in non-forests ( Figure 6B). The Z * parameter, or "e-folding" for the depth curve, was generally higher in grasslands and croplands than in forests, potentially due to deeper fine root distributions ( Figure 6A). However, the variability in Z * and its dependence on a range of biotic and abiotic processes provides ample opportunities for further analysis.
Land cover and land use histories were tied to significant variations in SC depth parameters. Profiles with and without an Ap horizon (indicating a history of tillage) showed legacies in agricultural disturbance that affected SC stocks and their vertical distributions. Tillage was associated with lower surface C, a more gradual decline (higher Z * ) in C density with depth, and slightly lower C stocks integrated to 30 cm (Figures 5, 6). These effects are consistent with tillage mixing SC from the surface through the plow layer of the soil profile. Differences in historically tilled soils were evident even in profiles with non-agricultural land covers, suggesting that subsequent changes in land cover (e.g., reforestation) did not erase the impact of tillage history. These results are consistent with previous studies FIGURE 6 | Mean Z*, C s , and SC-stocks-to-30 cm for profiles with and without Ap horizons grouped according to NLCD land cover type (A-C) or profile soil order (D-F). Error bars show the standard error of the mean. (Baker et al., 2007;Angers and Eriksen-Hamel, 2008;Sanderman et al., 2017;Nave et al., 2018).
Comparison of profiles within a grid of constrained spatial areas ( Figure 7A) was consistent with the overall finding that tillage history was associated with lower C stocks for topsoils (to a depth of 30 cm). While previous studies have relied on spatial modeling (e.g., Sanderman et al., 2017) or extrapolated results from a limited number of observations (e.g., Angers and Eriksen-Hamel, 2008), our analysis recovers these tillage-SC relationships from a large, broadly distributed dataset of tens of thousands of soil profile measurements traceable back to individual profiles.
In addition to the current analysis, the database of soil depth distribution parameters produced by our analysis has broad potential applications for constraining depth patterns of SC in both empirical and modeling analyses. Current land surface models are now incorporating depth distributions of SC stocks by modeling SC dynamics across multiple layers (e.g., Koven et al., 2013;Burke et al., 2017), though such models have not yet generally included agricultural tillage effects on SC cycling. Our database provides a potential tool for evaluating different model formulations and parameterizations based on the simulated depth distributions of SC stocks in addition to databases of total C stocks (to specified depth) that are commonly used in model evaluation. The use of individual profiles rather than a spatially averaged, gridded database allows variations in SC profiles due to variations in land use history, ecosystem type, and topography to be compared. This could allow land surface models that include sub-grid heterogeneity such as agricultural versus forest lands (e.g., Shevliakova et al., 2009) to be evaluated by comparing subgrid fractions directly with appropriately matched measurements rather than averaging over sub-grid variability. Our database can also facilitate large-scale analysis of controls on SC depth distributions such as topography, slope and aspect, climate, and ecosystem type in addition to the land cover and land use history factors that were the focus of this analysis. While our analysis was primarily limited to the continental United States, an expanded  Points below the 1:1 line indicate that profiles with an Ap horizon had lower SC stocks than profiles without an Ap horizon within the same soil order and spatial area.
dataset of SC depth distribution parameters could be generated from broader datasets with global distributions, such as the WoSIS database (Batjes et al., 2020).

CONCLUSION
We fit a harmonized depth distribution function of SC to a database of over 40,000 soil profiles across the continental United States, yielding a new database of depth distribution parameters. We used the database to compare depth parameters across soil orders, land cover classes, and land use histories. Several depth parameters were sensitive to soil order and land cover; in most cases, tillage indicated by the presence of Ap horizons reduced surface SC concentrations, increased the length scale (Z * ) of SC decline with depth, and decreased SC stocks in the top 30 cm. The database of depth profile parameters is applicable to statistical analyses and land surface model evaluation across different ecosystem types and large spatial gradients.

DATA AVAILABILITY STATEMENT
All fitting and analysis code, along with calculated depth distribution parameters for each profile, are archived at doi: 10. 5281/zenodo.3924496. ISCN data are available through the ISCN data portal: http://iscn.fluxdata.org/data/access-data/databasereports/.

AUTHOR CONTRIBUTIONS
BS analyzed the data and led writing of the manuscript. CT conducted statistical analysis. JH, LN, CK, JO'D, and BS conceived the study. YH, CT, and BS calculated depth distribution fits. UM contributed to land cover data. All authors contributed to writing and editing of the manuscript.

ACKNOWLEDGMENTS
This work was facilitated by the ISCN and benefited from extensive data contributions by the USDA-Natural Resources Conservation Service (National Cooperative Soil Survey) and the United States Geological Survey.