**1. Introduction**

Phosphorus (P) is an essential, often limiting nutrient in ecosystems across the globe [1,2]. P has been studied in plants, soils, rivers, lakes, and oceans—both modern and ancient—and is at the center of critical questions about Earth's past and future. Today, concerns around P in terrestrial settings focus on soil fertility and crop yields, as well as on agricultural runoff and its effect on eutrophication [3–5]. Climate, vegetation, and weathering and erosion rates are inextricably linked, so in order to improve quantitative constraints on potential P fluxes from soils, each of these factors must be considered. The redox-sensitive metal iron (Fe) is often associated with both P and organic matter/carbon (C) in soils (e.g., [6–9]). The P-Fe oxide pathway is critical for terrestrial P transport; however, it is currently less well-understood than other aspects of terrestrial P cycling, such as P-C associations. Additionally, while much research linking P, weathering, and soil age has been done, such work has typically done on relatively small scales (e.g., chronosequences in Hawaii) relative to the global scale of the questions

that we ask here. Therefore, constraining climate, vegetation, and weathering controls on soil P and Fe is critical for predicting changes to fluxes of P from soils to rivers and beyond, and to improving our understanding of how regional soil fertility may evolve in response to a changing climate.

Prior to the evolution and expansion of land plants, changes in the flux of P from continents to oceans are thought to have been a primary control on marine productivity, and consequently, on the concentration of oxygen in the atmosphere [10,11]. It remains debated whether land plants' colonization of continents would have increased [12] or decreased [13] that P flux. The better we understand the modern distribution of and controls on P in soils, the better we can model global biogeochemical changes in Earth's past.

#### *1.1. P, Fe, and Weathering and Erosion*

P is released from bedrock via the weathering of P-bearing minerals (primarily apatites, which are Ca minerals relatively resistant to weathering [4]). Weathering depends on all soil-forming factors (climate, biology, topography, time, parent material), but is primarily controlled by climate (e.g., precipitation, temperature, seasonality) and time (surface age, length of exposure [14]). Organic matter also serves as an important pool of P in soils ([4,5]), and vegetation, fungi, and microbial life play an important role in both P and Fe liberation and mobilization through symbiotic (mycorrhizal) relationships [15,16]. The presence of plants, then, changes how both P and Fe are distributed in soils—but the links between landscape-scale soil chemistry, plant functional type, climate, and weathering intensity have not been quantified.

P can be depleted from a soil profile within thousands to tens of thousands of years [4,17–20]. In terrestrial vegetation (biomass), P's residence time is ~13 years, and residence in soil is ~600 years without P replenishment by dust [4]. Soil P can also be replenished through dust deposition (e.g., [17,21,22]). The natural progression of a soil, with a slow accumulation of Fe and Al oxides and a loss of mobile cations and nutrients, means that P is typically very low in old and/or intensely weathered soils [19]. Other soil properties, such as clay content, are also linked to weathering and vegetation, and could affect the concentrations of both P and Fe. Clay minerals and clay grain-size fractions can sorb considerable amounts of P and Fe, and clay content typically increases with weathering time and/or intensity [23,24]. Weathering intensity is expected to decrease with latitude due to colder temperatures, drier precipitation regimes, and less vegetation and shallower rooting depths [23,24]. Additionally, because the potential volume of erodible material generally increases as weathering intensity increases (though varies with factors such as slope and uplift rate, e.g., [25–29], weathering is linked to P fluxes not only through direct apatite dissolution, but also through erosion rates, e.g., [30,31]). Understanding how weathering intensity is related to soil P (and Fe) concentrations is important for predicting how their erosional fluxes and transport may change in response to climate change (e.g., [30,31]).

After weathering, dust deposition is an important secondary source of P, particularly for old soils whose P has been depleted through weathering, erosion, and biological use. For example, chronosequence studies in Hawaii have shown that for older soils, dust replenishment of P is critical to maintaining fertility once bedrock sources have been depleted [17,32,33]. Dust is also a critical source of P in the Amazon basin, which receives significant dust fluxes from Africa [18,34–36], compensating for P loss in deeply weathered tropical soils. Dust-sourced P is less significant for areas that receive little dust deposition, for younger soils that still have a regular input of P from bedrock weathering, or for shorter timescales than either for significant dust accumulation or timescales of soils losing their bedrock-sourced P (i.e., anthropogenic timescales). Arvin et al. [21] found that dust deposition is an important P source in montane soils, despite their relatively "fresh" bedrock. The balance between weathering rate and dust accumulation is also relevant for determining how important the latter is in soil P [22]. Dust is also commonly a source of Fe oxides, increasing soil fertility [35,37,38].

### *1.2. P and Fe Kinetics in Soils*

P and Fe are chemically associated in soils in part because of P's affinity to sorb to Fe (oxyhydr)oxides (e.g., [4,9,39–42]). Orthophosphate, the common ion form of free P in soils, is reactive towards those oxides and becomes immobile (not readily bioavailable) when they associate into neoformed minerals or are sorbed onto existing mineral surfaces. Fe-bound P may be a critical mode of P transport among terrestrial ecosystems and from terrestrial/aquatic to marine environments, as oxyhydroxide-bound P is one of the major bioavailable/exchangeable forms of P in fluvial systems [5].

While we know that Fe oxide-adsorption of P is kinetically favorable in soils and other environments where P is being exchanged (e.g., [43–45]), our understanding of first-order controls on Fe concentration and speciation (i.e., Fe3<sup>+</sup> in oxides versus labile Fe2+) in soils could be improved. Examining a full suite of extractable Fe (rather than only dithionite-citrate and/or ammonium oxalate) allows a more nuanced picture of Fe phases in soils, which can ultimately be linked to variable P mobility in different soil redox conditions [43,46,47]. It is possible that regional climate factors (precipitation, temperature) exert control on Fe speciation (i.e., reduced vs. oxidized) in soils via soil moisture. Controls on soil moisture are debated and likely vary from profile to profile [48–53], and pH (related to soil pore space saturation; often microbially-mediated) is a primary control on reduction/oxidation of Fe phases in soils (e.g., [54]). Additionally, soil porewaters can serve as an intermediate step between recalcitrance and fluvial suspension for dissolved constituents (e.g., [54,55]). Exploring these relationships is key for understanding the likelihood of P mobilization under different moisture situations.

#### *1.3. P and Fe in the Geologic Record*

Geologists are interested in constraining many of the processes described above in ancient soils and ecosystems. Weathered P, transported to the oceans via fluvial networks and continental drainage, is typically invoked as a first-order control on marine primary productivity, which is associated with C sedimentation, atmospheric CO2 drawdown, and oxygen production (e.g., [10,11,56]). Similar to P, Fe can stimulate marine productivity on short timescales (e.g., [57–60]), but it is of additional interest because in the past, it served as a P sorption site in marine systems, effectively sequestering P so that it is not bioavailable (the so-called "Fe-P trap") and limiting marine productivity [10].

With these processes in mind, geologists are interested in continent-to-ocean fluxes and mobility of Fe and P through time. However, the continent-to-ocean transport of P, although central to many biogeochemical models, is usually prescribed and not based on actual changes in fossil soil (paleosol) P values or weathering intensity data. Rather, it has been generalized based on crustal P values and fluxes controlled by large-magnitude changes, such as global glaciations or limited modern observations (e.g., [61]). The value in providing robust constraints on the range and distribution of P in soils, as well as climatic and environmental controls on its mobility and bioavailability, will vastly improve the parameters for these models.

In addition to biogeochemical models, geologists are interested in how the evolution of terrestrial vegetation and mycorrhizae (fungal root symbionts) affected P bioavailability, terrestrial mobility, and fluxes to the oceans. Researchers have argued for both an increase [12,56,62,63] and decrease [13,15] in P fluxes to the oceans following the spread of land plants. Constraining how modern vegetation relates to soil P mobility and distributions will help to answer this question.

Finally, the geochemical composition of fossil soils is used to reconstruct a range of climatic and environmental changes; however, these tools are primarily based on small datasets. This work provides ranges of reasonable values for soil chemistries under known environmental and atmospheric conditions, providing critical background information to improve our paleoclimate and paleoenvironment reconstructions.

#### *1.4. This Work*

In the context of these concerns, several big-picture questions remain. Despite much research on P in modern soils, ecosystem (vegetation, anthropogenic impact) and climate (precipitation, weathering) controls on its concentration have not been quantified on landscape to continental scales for generalizable conclusions. Additionally, although P is known to have a strong affinity towards Fe (oxyhydr)oxides in soils, Fe species and soil redox conditions have not been considered as potential first-order controls on landscape-scale P bioavailability. Questions about how climate (precipitation, temperature), soil factors (drainage, soil moisture, clay content), and overall weathering intensity affect soil Fe-P associations needs to be interrogated, to see if Fe species ultimately serves as a mediator for P bioavailability. Moreover, once the distribution of and controls on soil P have been quantified, how do those new constraints affect or constrain our understanding of biogeochemical P models in Earth's past?

We address these questions, along with hypotheses about P in the geologic record, using a large geochemical database of modern soils (n > 5000). The results can inform models of both past and future terrestrial biogeochemical cycles.

#### **2. Materials and Methods**

#### *2.1. Sampling*

Soil samples were collected in the field by the authors, and received from the USDA KSSL soil repository, and from contributors (total n = 404). A majority (n = 247) of these samples are B horizons, and therefore, reflect longer-term, stable pools of P (e.g., P sorbed to clays or Fe (oxyhydr)oxides, which we test for here, or recalcitrant organic matter) and Fe. These soils include Entisols, Inceptisols, Alfisols, Ultisols, Histosols, Mollisols, and Aridisols and provide coverage from low (Hawaii) to high (Iceland) northern latitudes (Supplemental Figure S1). Soils collected in the field or from repositories were screened for anthropogenic influence (i.e., on developed or disturbed land), with anthropogenically-altered soils noted as such and binned separately in our treatments of the data here. To complement these soils, we used a United States Geological Survey (USGS) soil geochemistry interactive report [64] and collated geochemical data for the Top 5 cm, A horizons, and C horizons of all the soils included in that work (n = 4857). See Supplemental Table S1 for summary descriptions of these two datasets, and Supplemental Table S2 for complete sample information for the new dataset presented here. Geochemical data for the USGS dataset [64] are available at https://doi.org/10.3133/sir20175118.

Combining these soil databases yields excellent spatial, climatological, profile depth, drainage, vegetation (as plant functional type/generalized biome), and soil order coverage (Supplemental Figure S1). While the samples are mostly from the United States, variation in soil forming factors across the continent reflect most of the potential global variation. The tropics are somewhat under-sampled, but data from Hawaii are broadly reflective of those latitudes. While there are species-specific relationships between plants and soil geochemistry, generalized biome is acceptable for the continental scale of questions examined in this work, as well as the comparative nature of our inquiries (e.g., [65–68]). Differences between plant functional types typically swamp differences within one plant functional type. Additionally, because biomes on sub-continental scales are typical for global biogeochemical models, and because more specific data are not available for the latitudinal scale of this work, generalizing by biome scale/plant functional type is appropriate.

#### *2.2. Geochemical Analyses*

All analyses were performed on dry soils within the Sheldon lab, as were samples during the USGS analyses. Soils were dried per USDA APHIS regulations upon import and ground to <70 μm for homogeneity during analyses; samples were stored dry in oxic (ambient room) conditions (during the experimental design, tests of pre- and post-drying, with drying temperatures varying between 25, 50, and 100 ◦C, showed no significant geochemical differences). Bulk elemental geochemistry was performed at ALS Laboratories in Vancouver, BC, Canada, using a four-acid digestion and Inductively Coupled Plasma Mass Spectrometer (ICP-MS) analysis. External precision was 0.01% for all major elements and 10 ppm for P2O5, analyzed at ALS. For samples from the USGS dataset, 1σ sample errors were 488 ppm and 1.39 wt% for P2O5 and total Fe, respectively) [64]. The 1σ sample errors for total C and organic C were 4.6% (n.d. for inorganic C [64]).

To analyze Fe speciation in soils, we used a four-step sequential Fe extraction following Poulton and Canfield [69] and Raiswell et al. [70]. This extraction protocol separates operationally-defined pools as ascorbic acid (labile Fe; Feasc), dithionite- (crystalline Fe oxyhydroxides; Fedith), ammonium oxalate- (highly crystalline Fe oxyhydroxides; Feox), and sodium acetate-extractable Fe (Fe in carbonates, analyzed separately from the sequential extraction; Feacet). Initially, we included a chromium-reducible sulfur (CRS [71]) extraction separately to analyze for Fe sulfides (e.g., pyrite), but those concentrations were below detection limit for all non-wetland soils. CRS was limited to wetland soils after the initial exploratory extractions, and raw values were normalized to average soil order densities to account for the large difference in density between Histosols and denser, more mineral-rich soil orders [72]. Sample supernatants were centrifuged and diluted 50×, and Fe was measured via ICP-Optical Emission Spectroscopy and ICP-MS at the University of Michigan.

An internal standard of a local soil, dried at three different temperatures, was included in all runs for inter-run consistency checks and to determine reproducibility. The absolute standard deviation for a sample's reproducibility depends on the concentration of each Fe species, so the relative percent error ((σ/μ) × 100) is used instead. Repeated analyses show that Feasc is reliable within 3.7% of the mean, Fedith, within 1.7% of the mean, and Feox, within 7.0% of the mean. Reproducibility on the Feacet samples was larger, at 15% of the mean, because the means in many of the samples we analyzed were near machine detection limits (~1000 ppm), where a difference of several hundred ppm yields far larger standard deviations than species with higher concentrations. The Feasc pool is considered to be a low constraint because a larger proportion of labile Fe may be in a soil's porewater (rather than solids). Because here we are primarily interested in relative differences between oxide and labile Fe species, and because species' concentrations tend to differ by orders of magnitude, we have deemed this approach acceptable. For further reading on nuances in sample preparation and extraction protocols, see two recent reviews by Raiswell et al. [73] and Algeo and Liu [74].

Fe speciation data were filtered for outliers by removing any samples where Fespec weight percent > total Fe weight percent (i.e., the concentration of a measured Fe pool was higher than the measured total Fe). These discrepancies (n = 107, 95, 139, 106, for Feasc, Feacet, Fedith, and Feox, respectively) likely arose from the 50×-dilution being insufficient, with high Fe concentrations overwhelming the ICP-OES, leading to erroneous measurements. Future work using Fe speciation on soils should consider using stronger dilution factors for high-Fe soils to circumvent instrument limitations.

#### *2.3. Other Data Collated*

Geochemical data for the Top 5 cm, A horizons, and C horizons from 4857 soils around the U.S. are available on an interactive portal [64]. Weathering-relevant cations (Al, Ca, Na, K) and the nutrients of interest (P, Fe) along with latitude and total clay content (clay: A and C horizons only) were compiled. For the USGS dataset, clay content was determined by X-Ray diffraction (for 10 and 14Å clays only). Organic, inorganic, and total carbon data were available for a majority of soils in A and C horizons in the USGS dataset (not all had inorganic carbon). It should be noted that a dry soil sample will have a higher relative weight percent organic carbon due to water loss during drying. Soil orders and drainage capability are not given for soils in this database, but because they span the conterminous U.S., it is a reasonable assumption that they include a variety of soil orders. Clay content as grain-size fraction, parent material, and soil moisture regime were also gathered for the newly-collected dataset from USDA Soil Series reports when available (Supplemental Table S2).

Climate data were gathered from Soil Series pages; when these were not available, the PRISM 30-year averages (1981–2010) were used. If neither Soil Series nor PRISM climate data were available (i.e., for non-US samples), country-specific governmental meteorological data or journal articles were used (see Supplemental Table S2 for details).

The Chemical Index of Alteration (CIA) [75] was calculated to assess the degree of weathering for each soil, which reflects all of the soil-forming factors to some degree, but predominantly soil age and climate [14]. The CIA is thought to reflect the weathering of feldspars to form clay minerals where high values represent more intense weathering. Molar oxide ratios are used (Equation (1)). The CIA of an individual soil may not correspond directly with the CIA values of sediments in rivers (e.g., [76]). However, because our soils' average CIA values (59 ± 15, 42 ± 32, and 57 ± 21 for A, B, and C horizons, respectively) are within 1σ of the North American rivers' average from that study (n = 7; mean CIA = 66 ± 8.5, range 53–73), though with larger ranges (0.5 to 99 for all three horizons), we proceed with its use here.

$$\text{CIA} = (\text{Al}\_2\text{O}\_3 (\text{Al}\_2\text{O}\_3 + \text{CaO} + \text{NaO} + \text{K}\_2\text{O})) \times 100 \tag{1}$$

#### *2.4. Statistical Analyses*

To examine how the total Fe and P correlate with CIA (weathering intensity), a running average with a 10-CIA-unit window was used. Linear least-squares regressions were used to test for significant linear correlations where appropriate, based on the hypotheses that were being tested. Due to the large size of the USGS dataset, *p*-values are extremely low (<10<sup>−</sup>16; essentially 0), meaning it is highly likely that the two variables are not randomly correlated. Strength of correlation should be primarily considered before taking any relationships to be predictive, but for simplicity's sake, coefficients of determination (R<sup>2</sup> values) are reported. Principle Components Analysis (PCA) was performed for USGS database soils by horizon (Top 5 cm, A, and C). These analyses are used to test the expected relationships (e.g., that Fe and P will be positively and linearly related), rather than to build predictive models. All of the statistical analyses were performed in Matlab.

#### **3. Results**

Complete sample information and all geochemical results are found in Supplemental Tables S2–S5 and are available in the online version of this paper. For the larger USGS dataset [64], we have geochemical, drainage, and vegetation (generalized plant functional type) data, but no soil order or Fe speciation information. Results for soil order, drainage, parent material, and Fe speciation analyses are based only on the newly collected dataset (n = 404).
