Introduction

In his personal essay titled, “The Median Isn’t the Message,” Stephen Jay Gould notes we often think about medians as “hard realities” and the variation in values around them “a set of transient and imperfect measurements of this hidden essence” (Gould 1991). The same is true for the mean, which is far more commonly reported in the biogeochemical literature than the median. In Gould’s case, the median life expectancy of his mesothelioma diagnosis was 8 months, although he lived for another 20 years: the skewness of the distribution of life expectancies for individuals with that type of cancer mattered to his expected longevity. A similar principle applies to nutrient fluxes in terrestrial ecosystems; when we distill datasets down to their median or mean, we discount the extreme values, the nutrient “cold spots” or “hot spots.” Indeed, in many biogeochemical and other scientific studies, efforts are made to remove these outliers from the dataset. Just like Stephen J. Gould’s prognosis, the median (or the mean) does not represent well the population of nutrient fluxes in soil, nor do they adequately reflect the significance of nutrient flux outliers to the overall biogeochemical structure of the soil.

Areas of elevated biogeochemical activity in soil have long been recognized by the terrestrial ecological community as so-called “hot spots” (HSs; Parkin 1987; McClain et al. 2003; Vidon et al. 2010). Hot spots have been identified quantitatively by a skewed distribution (Kurunc et al. 2011), and as high statistical outliers in nutrient pool sizes, transformation rates, or fluxes (i.e., Johnson et al. 2010; Darrouzet-Nardi and Bowman 2011). Biogeochemical activity can also be “hot” temporally; these so-called “hot moments” (HMs) are defined by a high rate of biogeochemical activity occurring at a specific location for a short period of time relative to the rate observed in the bulk soil over longer time periods (McClain et al. 2003; Vidon et al. 2010). More recently, Bernhardt et al. (2017) recommended using the new term of “ecosystem control points” to better recognize that HSs and HMs are integrated in time and space; however, in this paper, we present and discuss our findings using the HS-HM nomenclature because they occur more commonly in the literature, while still acknowledging spatiotemporal dynamics.

Many investigations of soil HSs assess the dynamics of only one or two nutrients (Bernhardt et al. 2017) over a relatively short timeframe (i.e., within one year; Gibson 1986; Kleb and Wilson 1997; Johnson et al. 2014; Barcellos et al. 2018), and often use destructive sampling techniques when conducted in the field (i.e., Bruckner et al. 1999 Charley and West 1977; Kleb and Wilson 1997; Farley and Fitter 2001; Washburn and Arthur 2003; Mueller et al. 2008; Harms and Grimm 2008; but see Bilyera et al. 2020; Kratz et al. 2022). Without repeatedly sampling the same location over time, it is impossible to know if HSs are elevated for prolonged periods or if they only occur periodically (i.e., as HMs). Additionally, without evaluating several nutrients simultaneously, the relative importance of HS phenomena to overall biogeochemical cycling within an ecosystem cannot be ascertained. Before we can incorporate HSs into terrestrial modeling, we need to better quantify their presence and determine how soil properties and other environmental factors, such as water availability, temperature, and vegetation, regulate their occurrence (Groffman et al. 2009; Chen et al. 2021; Zhao et al. 2021).

There are a number of contributing factors that have been identified to create HSs. Hot spots may occur in localized areas where vegetation concentrates nutrient resources due to enhanced atmospheric dry deposition (via impaction and washout), throughfall and stemflow, litterfall, and decreased erosion (Halvorson et al. 1994; Schlesinger et al. 1996; Weathers et al. 2001; Darrouzet-Nardi and Bowman 2011). Additionally, depressional topographical positions may form HSs relative to upland landscape positions due to the accumulation of organic and inorganic nutrients, and enhanced soil water content promoting microbial activity (Harms and Grimm 2008; Barcellos et al. 2018). On a smaller spatial scale, surface microtopography, such as rocks and preferential flowpaths, can funnel soluble nutrients into the mineral soil to form nutrient HSs (Bundt et al. 2001; Johnson et al. 2011; Frei et al. 2012; Göransson et al. 2014).

In contrast, HMs are regulated by abiotic or biotic controls that are transient at a given location to elevate nutrient pools, processing rates, or fluxes. For example, pulses of nutrients or microbial activity may only be activated immediately following re-wetting events of dry soil (McClain et al. 2003; Harms and Grimm 2008; Barcellos et al. 2018; Barrat et al. 2021), or immediately after catastrophic disturbances such as wildfire (Zhao et al. 2021). In addition, dead and decaying roots or macrofauna can create biopores that accrue nutrients and serve as HMs of microbial activity (Kuzyakov and Blagodatskaya 2015; Hoang et al. 2016). The timeframe and interactions of controls regulate how long the area is “hot,” with HM formation essential for HS occurrence (Bernhardt et al. 2017; Barcellos et al. 2018; Wagner-Riddle et al. 2020), and dictate the respective impact HSs and HMs have on the biogeochemical behavior of the ecosystem. Indeed, the designation of HSs and HMs are defined by individual studies based on the specific temporal and spatial scales considered.

Forests that experience Mediterranean-type climates, like the Sierra Nevada of California, are excellent model systems to explore HS phenomena. Nutrient fluxes in Sierran forests are driven largely by activated hydrologic flowpaths that occur primarily following rains in the fall after a prolonged summer dry-period, and following relatively rapid snowmelt in the spring (Homyak et al. 2014; Williams et al. 1995; Hunsaker and Johnson 2017). Annual precipitation in these ecosystems is highly variable and is becoming more inconsistent with climate change (Zamora-Reyes et al. 2021). In a mixed-conifer forest in the western Sierra, Johnson et al. (2014) reported persistent and elevated magnesium (Mg2+) fluxes over fall and spring at specific soil locations, suggesting HSs; however, other nutrient HSs (ortho-phosphate: PO43−; ammonium: NH4+; nitrate: NO3−; calcium: Ca2+; and sodium: Na+) exhibited more HM phenomena because their fluxes varied spatially with time. Although this and associated studies (Woodward et al. 2013; Johnson et al. 2011; Woodward et al. 2013) improved our understanding of HS controls in Mediterranean forests, nutrient flux measurements were only conducted in the fall and spring of a single year; hence, it is uncertain if the observed patterns persist across multiple years with contrasting precipitation inputs. Understanding the longer-term controls on HS is likely to become increasingly important as the climate continues to warm (Ma et al. 2017).

Using the same plots and sampling design, and building upon the data collected over a one-year period reported by Johnson et al. (2014), we investigated long-term spatial and temporal patterns of nutrient fluxes in a Sierran mixed-conifer forest using ion exchange resin (IER) capsules. Specifically, over a four-year period, we characterized the magnitude and spatial distribution of nutrient fluxes (PO43−, NH4+, NO3−, Ca2+, Mg2+, and Na+; Na+ serves as semi-conservative tracer of water flow because it is not an essential nutrient; Vitousek and Reiners 1975) after the first significant fall rains following the dry summer growing season, and after snowmelt in the spring. We sampled the soil solution at the same spatial locations (horizontally and vertically) repeatedly over time to achieve the following research objectives: (1) characterize the spatiotemporal trends in HS occurrence within and across sampling years, seasons, and soil depths; (2) determine if HSs of different nutrients co-vary with each other over time; and (3) quantify the relative contribution of nutrient HSs to overall biogeochemical fluxes.

Materials and methods

Site descriptions

The study area is located on the western slope of the southern Sierra Nevada, California, USA. This region experiences a Mediterranean-type climate where mean annual air temperature is 9.8 °C and mean annual precipitation is 1325 mm year−1, with 35–60% falling as winter snow (Yang et al. 2021). The vegetation is a mixed-conifer forest with canopy-dominant tree species of sugar pine (Pinus lambertiana Douglas), white fir (Abies concolor (Gordon & Glend.) Lindl. ex Hildebr.), ponderosa pine (Pinus ponderosa P. Lawson & C. Lawson), and incense cedar (Calocedrus deucrrens (Torr.) Florin). Soils are derived from granitic parent material and classified in the Shaver soil series, a coarse-loamy, mixed, superactive, mesic Humic Dystroxerept (Johnson et al. 2014). Additional information on ecosystem productivity, soil biogeochemistry, and stream chemistry of these mixed-conifer forests in Sierran headwater watersheds can be found in Yang et al. (2022) and Hunsaker and Johnson (2017).

Sampling design

Johnson et al. (2014) established a 6- × 6-m grid of resin capsules (UNIBEST PST-1, https://www.unibestinc.com/technology; UNIBEST Corporation, Washington, USA; Yang and Skogley 1992; Dobermann et al. 1994) at two representative plots within a Sierran mixed-conifer forest (Fig. 1). These resin capsules are comprised of a mixed-bed (cation and anion) ion exchange resin with hydrogen and hydroxyl ions as the counterions, enclosed within a water-permeable, inert mesh spherical capsule. The locations used were near the plots where previous field investigations on HS-HM nutrient phenomena took place (Johnson et al. 2011; Woodward et al. 2013). Plot 1 (Lat. 37.061°, Long. − 119.183°) is at an elevation of 1980 m, has an eastern aspect, and is on a 20% slope. Plot 2 (Lat. 37.065°, Long. − 119.203°) is at an elevation of 1910 m, has an eastern aspect, and is on a 5% slope. Plot 1 has less dense vegetation than Plot 2, with 85% and 95% average overstory canopy cover and 45 and 75 m2 ha−1 tree basal area, respectively (Woodward et al. 2013). At each plot, 16 resin capsules were placed carefully at 2-m intervals horizontally within this grid immediately below the O horizon (i.e., 0-cm mineral depth) using a hand trowel to minimize disturbance. Adjacent to these resin capsules, WECSA® Access System units (WECSA® LLC; Montana, USA) were installed at 20-, 40-, and 60-cm depths at an angle that varied according to the microtopography, such that the deeper capsules would all be vertically coincident with their respective 0-cm capsule. The WECSA® system consisted of a thinned-walled PVC casing with a removable cap on one end and a removable resin capsule carrier that fits inside the casing; this allows the delivery and retrieval of a resin capsule so that the same soil location can be sampled repeatedly over time. The WECSA® systems were not installed in mineral soil locations when large roots or rocks prevented vertical placement below the 0-cm resin. This sampling design resulted in a horizontal and vertical matrix consisting of 60 resin capsules in Plot 1 and 63 in Plot 2 that allowed repeated sampling of the same locations with minimal soil disturbance. Only the 0-cm capsules required any subsequent soil disturbance, and this was kept at a minimum by carefully peeling back the O horizon during capsule retrieval and replacement.

Fig. 1
figure 1

General resin capsule layout deployed in Plot 1 and Plot 2. Each gray circle depicts a resin capsule under the O horizon, followed by resin capsules installed using WECSA® systems at 20-, 40-, and 60-cm depths (i.e., 16 locations that each consist of 4 resin capsules placed vertically into the soil)

Full size image

Resin capsules were installed at both plots in August 2011 and retrieved in October 2011 after the first major precipitation event in the fall (defined as ≥ 30 mm of daily precipitation). Then, new resin capsules were installed and left in situ until the completion of snowmelt in June 2012 (hereafter referred to as the Water Year 2012 collection; Johnson et al. 2014). A similar procedure was repeated in Plot 1 for water years (WY) 2013 and 2014 and for Plot 2 for WY 2013, 2014, and 2016 (resin capsules were not installed during WY 2015 in either plot; Table S1). This resulted in a final population of n = 360 at Plot 1 and n = 504 at Plot 2 when considering all spatial locations and sampling periods.

Soil volumetric water content and precipitation data were acquired over the entire study period to evaluate the hydrologic environment experienced by the resin capsule during incubation. Volumetric water content (m3 m−3) was measured in soil with a northern aspect at the Upper Providence Creek meteorological station (Long. 37.075°, Lat. − 119.182°; 1920 m elevation) located ~ 2 km from the resin plots (Critical Zone Observatory; Bales et al. 2022). Values were recorded hourly using Decagon Devices ECHO-TM at 10-, 30-, 60-, and 90-cm depths, and then averaged by day using the dplyr package in R (R Core Team 2008; Hadley et al. 2018). Daily precipitation (mm day−1) was also collected at the Upper Providence meteorological station using a Belfort-TM 50780 shielded weighing rain gauge (Kings River Experimental Watershed; Hunsaker and Safeeq 2018). Soil penetration resistance (kPa; i.e., soil strength) was measured at each plot to identify bulk density differences that may impact water fluxes using a Field Scout SC900 Soil Compaction Meter (Spectrum Technologies, Inc., Aurora, IL, USA) with a 1.27-cm cone size. Soil penetration resistance measurements were made after the last set of resin was collected and taken throughout both plots in undisturbed soil between the locations where WECSA® access systems were installed. Penetration resistance was logged at 106 and 87 points in Plot 1 and Plot 2, respectively, from 0- to 45-cm (maximum instrument depth); these values were then averaged across all recorded points within a plot for each 2.5-cm depth increment.

Laboratory analyses

After the capsules were collected and upon returning to the laboratory at the University of California, Merced, USA, adhering soil was removed with deionized water. The resin capsules were processed following the procedure described in Johnson et al. (2014). Briefly, the resin capsules were extracted with three sequential 20-mL solutions of 2 M KCl on an orbital shaker at low speed (≈ 60 rpm) for 20 min each. The three extracts were then combined (60 mL). This procedure was used to ensure complete exchange of adsorbed ions from the resin (Hart and Binkley 1984). For the WY 2012 samples (data reported in Johnson et al. 2014), extracts were sent to the Oklahoma State University Soil, Water and Forage Laboratory (SWAFL) and analyzed for NH4+-N, NO3−-N, and PO43−-P on a Lachat 8000 Flow-Injection Analyzer (Hach Co., Loveland, CO, USA), and Ca2+, Mg2+, and Na+ were analyzed using a Jarrell Ash Inductively Coupled Plasma Optical Emission Spectrometer (ICP-OES; Thermo Scientific, Waltham, MA, USA). For extracts generated from resin capsules collected after the first water year (i.e., WY 2013, 2014, and 2016), extracts were analyzed for NH4+-N, NO3−-N, and PO43−-P on a Lachat 8000 flow-injection analyzer, and Ca2+, Mg2+, and Na+ on a Perkin-Elmer Optima 5300 DV ICP-OES (Waltham, MA, USA) at the Environmental Analytical Laboratory, University of California, Merced.

Nutrient resin fluxes

Nutrient resin flux (fluxR) describes ion accumulation on the resin for use as an index of soil nutrient status (Skogley and Dobermann 1996; Qian and Schoenau 2002). FluxesR are generally correlated with other measures of nutrient availability (i.e., “buried bags”; Binkley et al. 1986; Hart and Firestone 1989), but are influenced by different processes than these other flux metrics. FluxR is dependent on ion diffusion, mass flow, competition between the resin and other sinks (i.e., plant and microbial uptake), and are impacted by the installation method, making the volume of soil sampled unknown (Binkley 1984; Binkley and Hart 1989; Dobermann et al. 1994). Therefore, fluxR cannot be compared quantitatively to other metrics of nutrient flux.

FluxR was calculated by converting nutrient concentrations in the extracts to a nmol cm−2 day−1 basis; this was achieved by transforming mass to moles for that particular nutrient, multiplying the concentrations by the total extractant volume (60 mL), and then dividing those values by the cross-sectional area of the capsule exposed to the soil (5.7 cm2 for WECSA® deployed capsules and 11.4 cm2 for 0-cm depth capsules; Johnson et al. 2014) and the incubation period (Table S1). Normalization to the number of incubation days allows for a more direct comparison among sampling periods because resin capsules were incubated in the field for different time periods due to the variability in the timing of spring snowmelt and significant fall rain across years. Although the number of days the resin was deployed varied by year (Table S1), high nutrient fluxes generally occur in these ecosystems only during brief periods after snowmelt in the spring and immediately following the first significant rainfall after the summer dry period in the fall (Homyak et al. 2014; Williams et al. 1995; Holloway and Dahlgren 2001; Hunsaker and Johnson 2017).

HS identification

Hot spots were identified as positive moderate outliers using the following equation (Johnson et al. 2014):

$$X>Q3+1.5\left(IQR\right),$$
(1)

where X is the HS value (i.e., outlier), Q3 is the third quartile value (75th percentile) and IQR is the interquartile range (25th-75th percentiles). The 25th percentile, 75th percentile, and IQR for a given nutrient fluxR were determined using all grid locations, depths, seasons, and years, but data from each plot were analyzed separately (Table S1). Plots were not treated as field replicates due to physical and chemical differences that may mask individual plot outliers, and therefore HS patterns. No negative outliers were observed (Johnson et al. 2014).

Statistical analyses

Mixed-effect models were conducted to understand the likelihood of HS formation (binary response) with water year, season, and depth (i.e., the factors hypothesized to co-vary with HS formation). Models were run separately for each nutrient and plot. Water year (categorical; WY 2012 coded as 0), depth (continuous), season (categorical; fall coded as 0, spring as 1), and the interaction between depth and season were used as fixed effects. Random intercepts were used to account for repeated measures on the grid locations (n = 60 and 63 for Plots 1 and 2, respectively) and for the spatial dependency of the vertically associated grid points (n = 16). Models were performed with the lme4 package in R (Bates et al. 2015) and were fit by maximum likelihood with the Laplace Approximation. Simplified models were used for Ca2+ in Plots 1 and 2 and Mg2+ in Plot 1 because the full models would not converge due to the lack of response variables (i.e., no HSs) for WY 2013. In the simplified models, WY 2013 was removed from the dataset, but all other model characteristics remained the same. Multi-collinearity of the fixed variables was inspected by calculating the Variance Inflation Factor (VIF) using the VIF function in the car R package (Fox and Weisberg 2019). All VIF values were less than 4 (except for the simplified HS presence models, which were less than 9) for all independent variables, including the interaction term. Marginal (only fixed effects considered) and conditional (both fixed and random effects considered) R2 values were calculated following Nakagawa et al. (2013) using the R MuMIn package (Barton 2009). When explanatory variables (WY, season, depth) were statistically significant (p < 0.05), Tukey’s Post Hoc tests were performed to evaluate whether means were significantly different from each other using the lsmeans function within the emmeans R package (Lenth et al. 2018). Models were also run with precipitation, season, and depth as explanatory variables; however, large VIF values resulted between precipitation and season due to multicollinearity and therefore these model results are not reported. Correlations between total precipitation and HS count are provided in the Supplemental Information (Fig. S1). We also evaluated the relationship between soil volumetric water content and HS formation; however, because soil moisture measurements were not co-located with resin capsules and moisture is known to be highly heterogenous within the soil profile (Franklin et al. 2021), we use these data to better understand water dynamics of the entire soil system rather than directly linking them to HS formation.

Spearman Rank correlation coefficients were calculated to identify if multiple nutrient HSs commonly occurred together in space using the rcorr function within the R corrplot package (Wei and Simko 2017). Hot spot presence was coded as a binary response variable (no HS was coded as 0; HS as 1) and data from each plot were again analyzed separately.

The contribution of HSs or bulk matrix soil (i.e., not a HS) to the total nutrient fluxR across all depths was also evaluated for a given nutrient in each sampling period and plot separately, using Eq. 2:

$$\frac{{\sum Sampling\,Period\,Nutrient\,Flux}_{HS\,or\,Matrix\,Soil}}{{\sum\,Samping\,Period\,Nutrient\,Fluxes}_{HS\,and\,Matrix\,Soil}} \times 100\%$$
(2)

A one-way analysis of variance was run to identify if the fluxR from HSs was significantly different than the matrix soil. Statistical analyses were carried out using R (version 4.2.3; R Core Team 2021) and all datasets are available for download (Barnes et al. 2023).

Results

Meteorological data

Precipitation varied greatly across sampling periods. During WYs 2012 and 2016, resin capsules were exposed to 121 and 113 mm of precipitation from the first major event in the fall, and 874 and 1103 mm of precipitation after the spring snowmelt, respectively (Fig. S2). This contrasts with WY 2013 and 2014 that had 95 and 72 mm of precipitation in the fall, but only 807 and 529 mm in the spring. Overall, resin capsules were exposed to 1.4 and 1.9 times less total precipitation in WYs 2013 and 2014 than 2016, respectively (WY calculated from October 1 of the previous to September 30 of the named year). Annual precipitation during the study period was an average of 1.7 ± 0.5 (standard deviation) times lower than mean annual precipitation (1325 mm year−1 on average from 2005 to 2017; Yang et al. 2021).

Changing precipitation patterns over the sampling periods corresponded to distinct patterns in soil volumetric water contents over time and with depth (Fig. S2). Volumetric water content followed precipitation events over the course of the year, and was always higher in the spring compared to the fall across all depths. For example, the 60-cm depth was approximately 0.04 m3 m−3 and 0.03 m3 m−3 drier in the fall and spring, respectively, in WY 2013 and 2014 compared to WY 2012 and 2016. Soil water content tended to be similar across depths with notable differences in the shallow 10-cm depth which had the lowest moisture content over the summer.

Nutrient flux distributions

All nutrient distributions were right (i.e., positive) skewed, with the degree of the skewed distribution variable between plots and among nutrients (Table 1; Fig. 2). Sodium exhibited the lowest median fluxesR (Plot 1: 0.06 nmol cm−2 day−1; Plot 2: 0.08 nmol cm−2 day−1), but had a similar skew as the other nutrients (Table 1). Ammonium (Plot 1: 0.15 nmol cm−2 day−1; Plot 2: 0.18 nmol cm−2 day−1) and NO3− (Plot 1: 0.14 nmol cm−2 day−1; Plot 2: 0.14 nmol cm−2 day−1) exhibited similar median fluxesR, but NH4+ exhibited a greater skew in Plot 2 (17.6) than Plot 1 (5.13). Phosphate (Plot 1: 0.46 nmol cm−2 day−1; Plot 2: 0.23 nmol cm−2 day−1) and Mg2+ (Plot 1: 0.93 nmol cm−2 day−1; Plot 2: 0.65 nmol cm−2 day−1) had higher median fluxesR than NH4+, NO3−, and Na+ with moderate skew (4.39 to 8.52). Calcium had the highest median nutrient fluxR (Plot 1: 4.69 nmol cm−2 day−1; Plot 2: 3.48 nmol cm−2 day−1) and tended to have a relatively lower skew than other nutrients within the same plot (Plot 1: 3.18; Plot 2: 7.43).

Table 1 Nutrient resin flux (nmol cm−2 day−1) descriptive information including mean (standard deviation), median, minimum, maximum, interquartile range (25–75th quartiles), and skew for PO43−, NH4+, NO3−, Ca2+, Mg2+, and Na+, calculated for Plot 1 and Plot 2 separately
Full size table
Fig. 2
figure 2

Histogram of log + 1 transformed PO43−, NH4+, NO3−, Ca2+, Mg2+, and Na+ resin fluxes. Fluxes are grouped into 50 bins. The dashed line indicates the median and solid line represents the hot spot (i.e., outlier) cutoff value

Full size image

Nutrient HS occurrence

Hot spots formed most often for NO3− (Plot 1 = 53; Plot 2 = 104) and Na+ (Plot 1 = 54; Plot 2 = 90) across all depths and sampling periods (Plot 1: n = 360, Plot 2: n = 504; Fig. 3). All other nutrients had a similar number of total HS, with 36 to 45 forming in Plot 1 and 56 to 66 in Plot 2.

Fig. 3
figure 3

Distribution of hot spots with water year, season, and depth for each nutrient and plot. When no bar is shown for a specific depth, then no HS occurred during the respective sampling period. No data were collected in Plot 1 for 2016. Total precipitation (mm) during each sampling period is included in parentheses following the water year

Full size image

All nutrients formed HSs repeatedly at the same spatial location in more than one sampling period, but HSs did not occur in all locations (Fig. 4). Nitrate had the greatest number of reoccurring HSs in Plot 1 (15 resin capsule locations), while HS reoccurrence for Na+ was greatest in Plot 2 (23 resin capsule locations). Although HSs did reoccur for Ca2+, Mg2, and Na+, they tended to only happen twice across all sampling periods, whereas some resin capsule locations formed HSs for more than two sampling periods for PO43−, NO3−, and NH4+. For all nutrients, HS reoccurrence was more common at the 0-cm depth than at deeper mineral soil locations. Although HSs often occurred in the same spatial location repeatedly for all nutrients, HSs occasionally persisted over more than two sampling periods for PO43−, NH4+, and NO3−, and were more transient for Ca2+, Mg2+, and Na+.

Fig. 4
figure 4

Number of times spatial locations were designated as a hot spot over all sampling periods in Plots 1 and 2 for PO43−, NH4+, NO3−, Ca2+, Mg2+, and Na.+. Plot 1: n = 360, Plot 2: n = 504. See Fig. 1 for resin location deployment details

Full size image

Relationship of water year, season, and depth with HS formation

Water year was a significant predictor of HS occurrence, particularly for Plot 2 (Tables 2 and 3; Fig. 3). Hot spots were more likely to occur for Ca2+, Mg2+, and Na+ in WY 2012 and WY 2016 than WY 2013 and WY 2014 at both plots (WY 2013 was not included in the model for Ca2+ and Mg2+ at Plot 1 or Ca2+ at Plot 1 and 2 due to singularity; no HSs occurred for these nutrients in WY 2013; Table 3). In contrast, HS occurrence for PO43− and NO3− was not significantly different in 2012, 2013, or 2014, and NH4+ HSs formed more often in 2013 than 2012. Hot spot formation was more likely in 2016 than the other water years for PO43−, NH4+. and NO3−. Overall, Ca2+, Mg2+, and Na+ HSs exhibited more variation among WY than the other nutrients.

Table 2 Results of mixed-effect models for Plot 1 and Plot 2 on hot spot likelihood (binary response) of PO43−, NH4+, NO3−, Ca2+, Mg2+, and Na+ as effected by the fixed effects of water year (categorical, 2012 coded as 0), season (categorical, fall coded as 0), depth (continuous), and season by depth interaction
Full size table
Table 3 Post-hoc results for hot spot mixed-effect models of PO43−, NH4+, NO3−, Ca2+, Mg2+, and Na+
Full size table

Season was occasionally a significant predictor of HS occurrence (Table 2). In Plot 1, PO43−, NH4+, NO3−, and Na+ HSs were greater in the fall after the first major precipitation event compared to following spring snowmelt. Season was only significant in Plot 2 for NH4+; in contrast, HS formation in Plot 1 was more likely in the fall. Calcium and Mg2+ HS formation was not impacted by season.

The likelihood of the occurrence of HSs decreased with depth into the soil profile for all nutrients and in both plots (Table 2). The interaction of season and depth was significant for predicting HS occurrence for Ca2+ in both plots and Na+ at Plot 2 (Tables 2 and 3). Calcium HSs were more likely to occur deeper in the soil profile after snowmelt in the spring (Table 3). In contrast, Na+ HSs persisted deeper into the mineral soil after the first major precipitation event in the fall (Table 3).

Nutrient HS covariance

Many spatial locations exhibited a HS for more than one nutrient within a given location and sampling period (Fig. 5). Indeed, all bi-nutrient correlation coefficients were positive and statistically significant. A HS occurred for two nutrients simultaneously 61 times in Plot 1 and 97 times in Plot 2. It was rare for a HS to form for all six nutrients within the same sampling period (i.e., 7 locations in Plot 1 and 13 in Plot 2). Calcium and Mg2+ had the highest correlation coefficient (Plot 1: ρ = 0.71; Plot 2: ρ = 0.72). Phosphate, NH4+, and NO3− had the next highest correlation coefficients between each other (ρ = 0.59 to 0.68) in both plots.

Fig. 5
figure 5

Spearman rank correlation coefficients (ρ) for hot spot presence calculated for Plot 1 and Plot 2 separately. Circle size is larger when the correlation coefficient is farther from 0. All correlations were positive and statistically significant (p < 0.05)

Full size image

Contribution of HSs to total nutrient flux

Hot spots represented 8 to 17% of the total sampling volume when considering all sampling events. Phosphate HSs occurred in the smallest volume, occupying only 8% of the sampling volume in Plot 1 and 9% in Plot 2. Nitrate HSs only occupied 12% of the sampled volume in Plot 1, but was 16% of the sampling volume in Plot 2. Sodium HSs occurred in a larger volume compared to the other nutrients (Plot 1: 14%, Plot 2: 17%).

The percent of the total nutrient fluxR attributed to HSs varied by nutrient when averaged across all sampling periods (Fig. 6). The contribution of NO3− HSs to the total NO3− fluxR was significantly greater than the contribution from the matrix soil, with HSs contributing 77.9 ± 17.9% (standard error) and 87.9 ± 11.1% of the total NO3− fluxesR in Plot 1 and 2, respectively. Ammonium (Plot 1: 62.2 ± 24.7%; Plot 2: 73.5 ± 21.3%) and PO43− (Plot 1: 61.4 ± 26.8; Plot 2: 56.4 ± 26.1) HS contributions to their respective total nutrient fluxesR followed a similar pattern as NO3−, with HSs representing most of the nutrient fluxR (although this was only statistically significantly in Plot 2 for NH4+). Total Na+ fluxesR were composed of approximately half from HSs (Plot 1: 46.7 ± 43.8; Plot 2: 56.6 ± 45.7) and half from the soil matrix. In contrast, Ca2+ HSs (Plot 1: 26.8 ± 36.8; Plot 2: 33.0 ± 40.0) and Mg2+ HSs (Plot 1: 28.7 ± 34.9; Plot 2: 39.6 ± 40.2) accounted for less than half of the total nutrient fluxR, with HS contribution for Ca2+ in Plot 1 being significantly less. Overall, PO43−, NH4+, NO3−, and Na+ HSs contributed to a relatively larger portion of the total nutrient fluxR compared to Ca2+ and Mg2+.

Fig. 6
figure 6

Percent and ± one standard error (vertical error bars) of the total nutrient resin flux observed from hot spots (orange) and the bulk soil matrix (i.e., non-hot spot; gray) calculated for Plot 1 and Plot 2 separately

Full size image

Discussion

Understanding the temporal and spatial controls that produce locations of elevated nutrients is critical to advance our understanding of HS phenomena. Although HSs are now well-recognized by the biogeochemical community, incorporating HS phenomena into predictive biogeochemical models remains challenging (Landsberg et al. 1991; Groffman et al. 2009; Arora et al. 2022). Most ecological models use nutrient averages for spatial and temporal scaling (Strayer et al. 2003), discounting outliers found in soil microsites (Schimel and Bennett 2004), and are often not spatially explicit (Landsberg 2003) or do not consider multiple nutrient HSs with varying concentrations (i.e., Ryel and Caldwell 1998). Although there has been some recent progress in this area (Arora et al. 2022), we need to use sampling designs that allow for their explicit characterization to advance our mechanistic understanding of HS formation and subsequent inclusion into models. Our study advances knowledge of nutrient HS phenomena by using an experimental design that allowed us to elucidate multi-nutrient HS occurrence at different spatial locations and temporal scales.

Influence of yearly drought on HS formation

Hot spot formation was impacted by WY, especially for Ca2+, Mg2+, and Na+, likely because annual precipitation exhibited large variation across sampling periods (ranged from 558 to 1223 from WY 2012 to 2016; see Table S1 for amount of precipitation resin capsules were exposed to). A multi-year drought started in the fall of 2011 and occurred during our study (Bales et al. 2018), resulting in variable total precipitation that was higher in WYs 2012 (1018 mm year−1) and 2016 (1223 mm year−1) than WYs 2013 (879 mm year−1), 2014 (646 mm year−1), and 2015 (558 mm year−1 and; Fig. S1). Although California has experienced multi-year droughts ten times in the last 100 years, unprecedented higher temperatures with reduced precipitation during this period caused more extreme conditions than previous dry periods (Bales et al. 2018; Griffin and Anchukaitis 2014). Calcium, Mg2+, and Na+ fluxesR and HS formation were greatest during WY 2012 and 2016, and were substantially lower in the extreme drought years (Fig. 3 and S3, Tables 2, 3, S2 and S3). Variable water fluxes across the study period likely impacted not only nutrient transport through the soil, but also altered the abundance of these nutrients by changing the rates of chemical weathering, cation exchange reactions, and the source and quantity of atmospheric deposition (Johnson et al. 1968; Gislason et al. 2009; Aarons et al. 2019).

In contrast, PO43, NH4+, and NO3− exhibited less variability in HS formation (Fig. 3) from WYs 2012 through 2014 (Tables 2 and 3), even though these nutrient fluxesR were also altered by WY (Fig. S3, Table S2 and S3). This suggests that labile C, water, or other limiting resources were accessible to microorganisms, even during drought years, to promote nutrient mineralization and nitrification, and water fluxes were high enough to transport these ions to the resin capsules across all sampling years. There was a notable increase in HSs for all nutrients in WY 2016, which followed the extreme drought WYs (2013–2015). Nutrients likely accumulated during the extreme drought years and were subsequently mobilized to the resin capsules during the following wet year (WY 2016). Rewetting after droughts periods have been found to cause large N and P pulses across multiple humid and temperate forests. Microbes accumulate intracellular solutes to aide in moisture retention during drought, which can result in cell lysis upon rewetting, with additional nutrient contributions from increased rates of litter decomposition and throughfall (Lodge et al. 1994; Fierer and Schimel 2002; Billings and Phillips 2011; Schlesinger et al. 2016). Water fluxes were likely too inconsistent during our study period to generate HSs repeatedly in the same location over time, albeit this occasionally occurred for PO43, NH4+, and NO3− (Fig. 4), making HSs mostly transient rather than persisting in the same location over time.

Influence of seasonal drought on HS formation

Although PO43−, NH4+, and NO3− HSs were relatively insensitive to drought on an annual basis, they were responsive to seasonal drought conditions. We found nutrient fluxesR and HSs for PO43−, NH4+, NO3−, and Na+ to be higher after the first fall rains compared to spring snowmelt (Table 2 and S2, Fig. 3 and S3). Nutrients are known to accumulate during the hot, dry summers that characterize Mediterranean climates (Davidson et al. 1990; Lewis et al. 2006). Previous research in the Sierra Nevada found O horizon interflow had higher nutrient concentrations compared to underlying mineral soil solution, snowmelt, and stream water, although these differences were highly variable in space and for the different nutrients assessed (Woodward 2012; Miller et al. 2005; Johnson et al. 2011). When the first major precipitation event in the fall occurred, O horizon interflow likely became nutrient-rich (Miller et al. 2005, 2006; Johnson et al. 2011) and infiltrated mineral soil via preferential flow paths (Burcar et al. 1994; Bundt et al. 2001) to form HSs. Alternatively, wetter antecedent moisture conditions in the spring compared to the fall may have resulted in rapid matrix flow (i.e., lateral flow through large pores and along preferential flow paths; Burcar et al. 1994), leading to less contact time between the soil solution and resin capsules along with lower diffusion gradients. Therefore, fall rains may have deposited highly concentrated nutrients in a heterogeneous manner, whereas spring snowmelt was more homogeneous.

HS structure with depth

Previous studies have largely focused on characterizing HS occurrence in surficial soil rather than evaluating their formation with depth into the soil. We found HS formation significantly decreased with depth in the soil profile for all nutrients examined (Fig. 3; Table 2). The greater occurrence of nutrient HSs in surficial soil than at greater soil depth is likely due, in part, to higher organic matter content, microbial activity, and water fluxes in the upper soil profile (Jobbágy and Jackson 2001; Sun et al. 2021). However, other biogeochemical differences among these nutrients, such as type and density of their charge, solubility in the soil solution, and the biological demand (plant and microbial), also likely play a role in determining their individual HS structure with soil depth (Wolt 1994). For example, PO43−, NH4+, and NO3− are in high biological demand (i.e., growth-limiting nutrients) such that these nutrients tend to be removed quickly from the soil solution as water moves downward through the soil profile. Hence, competition between the resin capsules and uptake by the soil biota may be responsible for the steep drop in occurrence of these nutrient HSs with soil depth (Binkley 1984). In contrast, Na+ differs from the other nutrients studied because it is not an essential element for plant and microbial growth, and therefore it does not accumulate within their biomass to a significant degree, resulting in Na+ acting like a semi-conservative tracer of water flow (Vitousek and Reiners 1975). This characteristic likely contributed to Na+ having the greatest number of HSs below 0-cm depth out of all nutrients (Fig. 3).

In the previous one-year study at this site using the same sampling devices and locations, PO43−, NH4+, and NO3− HSs were also concentrated in the upper sampling depths, whereas Ca2+, Mg2+, and Na+ HSs were identified at greater depths (Johnson et al. 2014). However, all nutrient HSs tended to persist more at depth in the Johnson et al. (2014) study, possibly because statistical analyses were run without the inclusion of the 0-cm depth resin capsule data. When we ran the same mixed-effect models without the data from the 0-cm depth, all nutrient fluxesR except Na+ still decreased with depth (data not shown). Hot Spot likelihood also decreased with depth at Plot 2 for PO43−, NH4+, NO3−, and Mg2+, although Plot 1 nutrient HSs were no longer significantly impacted by depth. This result highlights that the observed depth trends are not strictly regulated by elevated HS formation at the interface of organic and mineral soil horizons. The dissipation of nutrient HSs with depth (especially the growth-limiting nutrients PO43−, NH4+, and NO3−; Vitousek and Reiners 1975) may be due to plant root uptake because HS occurrence has a negative relationship with rooting patterns in these mixed-conifer forests experiencing a Mediterranean-type climate (Hart and Firestone 1991). Globally, nutrients required by plants in greater amounts are generally found in the upper profile, while others such as Na+ persist deeper in the profile (Jobbágy and Jackson 2001). Therefore, it is not surprising that nutrient HSs follow these broad-scale nutrient patterns.

Influence of HSs on biogeochemical structure

By definition, HSs are found in select spatial compartments, therefore they necessarily occupy a relatively small percentage of the total soil volume. We found that, although PO43−, NH4+, and NO3− HSs only occurred in 8–17% of the total area sampled, on average, HSs were responsible for 56–88% of the total nutrient fluxR (Fig. 6). This finding is similar to that observed previously in surficial soil (3-cm depth) in an alpine ecosystem where HSs only accounted for 14% of the study area but provided more than 50% of the bioavailable inorganic N (Darrouzet-Nardi and Bowman 2011). Additionally, a meta-analysis found rhizosphere HSs accounted for 10–33% of C and N net mineralization while only occupying 8–26% of the soil volume to a 10-cm depth (Finzi et al. 2015). In our study, Ca2+ and Mg2+ HSs accounted for a similar amount of the sampling area as PO43−, NH4+, and NO3−, however the proportion of the total nutrient fluxR attributed to HSs was relatively smaller (an average of 27–40%). Therefore, growth-limiting nutrient HSs, especially PO43−, NH4+, and NO3−, may play a disproportionate role in regulating ecosystem processes and function compared to the volume of soil they represent.

Hot Spot behavior tended to cluster into groups of nutrients that have similar biogeochemical cycles. This is likely because nutrients primarily under biological control (PO43−, NH4+, and NO3−) require high substrate availability and high microbial activity to co-occur with active hydrological flowpaths to form HSs (McClain et al. 2003). In contrast, soil nutrient fluxes that are primarily under geochemical control (i.e., chemical weathering, ion exchange reactions; Ca2+, Mg2+, and Na+) largely require only the co-occurrence of substrate abundance and water flow, and are less dependent on high rates of microbial activity. These drivers likely contributed to biologically versus geochemically controlled nutrient HSs to have: dissimilar responses to yearly and seasonal drought (Tables 2 and 3), contrasting patterns in reoccurrence (Fig. 4), differences in the tendency for nutrient HSs to co-occur within respective groups at the same spatial locations (Fig. 5), and biologically controlled nutrient HSs to contribute to a larger portion of the total fluxR compared to geochemically controlled nutrients (Fig. 6).

Plot characteristics

Hot spot patterns demonstrated variability between plots that is likely due to a variety of environmental differences between them, including soil bulk density (related to soil strength), slope degree, and vegetation density. Generally, Plot 2 had higher soil strength than Plot 1, suggesting that pore connectivity was higher in Plot 2 (Lipiec and Hatano 2003). Additionally, Plot 2 was on a shallower slope than Plot 1 (Plot 1: 20%, Plot 2: 5%), which may have resulted in greater water infiltration in Plot 2 by reducing overland flow. Plot 2 had greater vegetation density than Plot 1, which likely resulted in more homogeneous substrate distribution and greater competition for nutrients between plant roots and resin capsules. Taken together, these factors probably resulted in HS formation to be less seasonally dependent (Table 2) and created more co-occurring HSs across all nutrients in Plot 2 compared to Plot 1 (Fig. 5). Given the inter-plot variability, identifying easily measured variables that co-vary with HS formation is needed for HS incorporation into computer simulation models.

HS ecological significance

The transient nature of HSs likely plays an important role in nutrient competition and utilization between plants and microbes as well as among individual plants. Although HSs can be areas of intense competition for available nutrients between vegetation and soil microorganisms (Jingguo and Bakken 1997), especially N (Kaye and Hart 1997), the extremely dry summers in Mediterranean-type climates limits root development in the surficial O horizon (Hart and Firestone 1991); therefore, soil microorganisms likely obtain nutrients from the organic layer with minimal plant competition until leaching occurs into the mineral soil (Johnson et al. 2011). We found that reoccurrence, or duration, of nutrient HSs in soil tended to be more frequent in surficial soil (i.e., permanent control points; Bernhardt et al. 2017) than deeper in the soil profile (i.e., activated control points; Bernhardt et al. 2017; Fig. 4). Fungal species with vast mycelial networks (i.e., filamentous fungi) can likely exploit these reoccurring HSs by actively foraging for the nutrients contained in these locations; in contrast, more stationary microorganisms (i.e., bacteria and archaea) attached primarily to surfaces or held within soil aggregates are less likely to utilize nutrients contained within these locations as they are more dependent on diffusion and mass flow processes for nutrient uptake (Stark and Hart 1999).

The temporal characteristics of HSs benefits some vegetation species over others. Reoccurring HSs of PO43−, NH4+, or NO3− may stimulate a morphological response of the belowground plant biomass (Robinson 1994; Hutchings et al. 2003; Wang et al. 2018). However, HSs of all nutrients were more commonly transient in nature and thus are more likely to elicit a physiological response in vegetation (Fransen et al. 1999). We hypothesize that plants with low morphological plasticity, but an extensive rooting system that can exhibit physiological plasticity, would likely benefit most from transient HSs (Crick and Grime 1987; Cui and Caldwell 1997; Wang et al. 2018), while reoccurring HSs would be more advantageous for species with high morphological plasticity to exploit these persistent elevated nutrient patches (Jackson and Caldwell 1989; Hutchings and Kroon 1994). Although microorganisms may outcompete vegetation in the short-term for transient nutrient HSs, plants that access reoccurring HS with their roots and mycorrhizal symbionts can acquire more nutrients over the long-term (Cui and Caldwell 1997; Hodge et al. 2000; Schimel and Bennett 2004).

The response of plants to nutrient heterogeneity is impacted by the concentration and chemical composition of the nutrient patch, the nutrient distribution within the soil profile, the duration of the nutrient pulse, and the frequency it occurs (Jackson and Caldwell 1989; Bilbrough and Caldwell 1997; Jingguo and Bakken 1997; Fransen et al. 1999; Wijesinghe et al. 2001). Because HSs have extremely high fluxesR compared to the soil matrix, are composed of co-varying macronutrients, are found at soil depths relevant to roots, persist over time in some locations, and are responsible for a large portion of the total nutrient fluxR, they likely play a disproportionate role in supporting ecosystem structure and function.