Seal Tightly and Store in a Cool Dry Place: Exploring Soil Phosphorus Storage in a Subtropical Treatment Wetland.

Oligotrophic wetlands of the Everglades are often the nal recipients of nutrients from adjacent ecosystems and tend to accumulate phosphorus (P) in their soils. Understanding P source and sink dynamics in wetlands are critical for managing wetland ecosystems and protecting downstream resources. In this study, soil P storage capacity (SPSC) was evaluated within two treatment ow-ways of the Everglades Stormwater Treatment Areas (STAs). This study hypothesized that SPSC will vary between ow-ways, soil depth, and spatially along the inow-to-outow gradient. The P storage capacity in the STAs depend on the proportion of iron, aluminum, calcium, and magnesium (Fe, Al, Ca, and Mg, respectively) to P with oc and recently accreted soils (RAS) being associated more with Ca and Mg and pre-STA soils being associated more with Fe and Al. Phosphorus loss, as indicated from SPSC values would vary between systems and soil depths suggesting a variable condition of P sink and source within and along ow-ways. This result, while limited, demonstrates the applicability of SPSC to wetlands systems and provides information that will aid operational or management decisions associated with improving P retention of the Everglades STAs.


Introduction
Soils and sediments within wetland ecosystems store an appreciable amount of phosphorus (P) relative to other ecosystem components. This dynamic is driven largely by ecosystem position along the terrestrial-aquatic continuum where wetlands exist at the interface between terrestrial and aquatic ecosystems and receive inputs from adjacent upland ecosystems (Reddy and DeLaune 2008;Badiou et al 2018). As such, P accumulated within wetlands from upstream land uses can be both reactive or stable pools of P re ecting historical conditions within the watershed (Nair et al 2015). The storage and retention of P in wetlands are driven by interrelated abiotic and biotic processes. More speci cally, in mineral soils, P storage and retention can be dependent on the pH and concentrations of iron (Fe), aluminum (Al), calcium (Ca), and magnesium (Mg; Reddy et al 1999;Reddy et al 2005;Dierberg et al 2017).
While P storage is a bene t of wetland ecosystems, the potential to release stored P can put downstream ecosystems at risk. Therefore reliable techniques of predicting P storage and release have been developed taking into consideration organic and inorganic P, Fe, and Al as reviewed by Reddy et al (1998), Nair et al (2015), and others Nair and Harris 2004;Nair and Reddy 2013). The Degree of P saturation (DPS) or P saturation ratio (PSR) have been used as indicators of P loss from soils in agricultural applications (Maguire and Sims 2002;Haygarth et al 2005;Kleinman et al 2011). Meanwhile, Nair and Harris (2004) developed the concept of Soil P Storage Capacity (SPSC) where the P sorption capacity and PSR of a given soil are used to determine source/sink conditions relative to the capacity of soil to store P. Later, Nair et al (2015) applied the PSR and SPSC concepts to wetlands, demonstrating that heavily P impacted wetlands could be a potential P source ultimately signifying the value of the PSR-SPSC concept as a valuable tool for evaluating legacy P. Treatment wetlands leverage biogeochemical processes to retain nutrients, sediments, and contaminants as a strategy to protect downstream ecosystems (Kadlec and Wallace 2009). The most notable are the Everglades stormwater treatment areas (STAs); these treatment wetlands in south Florida were constructed to improve the quality of agricultural runoff water originating in the Everglades Agricultural Area (EAA) before entering the downstream Everglades ecosystem (Chen et al 2015). The Everglades ecosystem is a P-growth limited system with ora and fauna adapted to nutrient-poor (i.e. oligotrophic) conditions (Noe et al 2001). To be protective of these sensitive downstream ecosystems the STAs have a regulatory water quality requirement, outlined in the Everglades Forever Act ( § 373 and § 403 Florida Statutes and 62-302.540 Florida Administrative Code;Florida Department of Environmental Protection 2005) to ensure that discharges entering the Everglades do not cause or contribute to the degradation of water quality, ora, and fauna. Therefore understanding the source and sink dynamics of the Everglades STAs is important for the protection of this essential downstream ecosystem. Prior studies have remarked on notable soil P concentrations and storage gradients spatially stretching from in ow to out ow and along the soil pro le (Bhomia et al 2015;Zamorano et al 2018). However, limited information has been synthesized to evaluate SPSC of soils within the STAs. The overall objective of this study is to evaluate soil P sink/source dynamics within the Everglades STAs using the SPSC index. It is expected that as a result of biogeochemical cycling and burial dynamics, SPSC will vary between soil depths and spatially moving from in ow to out ow.

Study Area
Stormwater treatment areas (STAs) of the Everglades are constructed and managed wetlands that provide a natural reduction in water column pollutants and ultimately improve water quality conditions for downstream ecosystems. In South Florida, treatment wetlands were constructed to improve water quality and restore the biological integrity of the Everglades ecosystem by reducing nutrient concentrations of stormwater run-off from upstream land-uses before discharge to more pristine Everglades wetlands. A total of ve STAs with an approximate area of 231 km 2 are located south of Lake Okeechobee on the southern edge of the Everglades Agricultural Area (Fig. 1). The Everglades STAs are established on land formerly dominated by natural wetlands and/or agricultural farmlands under sugarcane production. The primary source of in ow water to the STAs is agricultural runoff originating from approximately 284 km 2 of farmland. Everglades STA treatment cells or ow-ways are comprised of a mixture of emergent aquatic vegetation (EAV) and submerged aquatic vegetation (SAV) communities in several con gurations including EAV and SAV treatment cells arranged in parallel or in-series (Chen et al 2015).
Stormwater Treatment Area-2 (STA-2) has been in operation since June 1999 with an effective treatment area of approximately 63 km 2 divided into eight treatment cells/ ow-ways. This study was conducted in two ow-ways: ow-ways 1 and 3. Both ow-ways receive the same source water from the same in ow canal. The vegetative community of ow-way 1 comprises predominately EAV including Typha domingensis Pers. (cattail) and Cladium jamaicense Crantz (sawgrass), while ow-way 3 is dominated by SAV, including Chara spp. (muskgrass), Potamogeton spp. (pondweed), and Najas guadalupensis Spreng (southern naiad), with approximately a quarter of the ow-way occupied by EAV species. Furthermore, prior to STA-2 construction, ow-way 1 was a historic natural wetland while approximately three-quarters of ow-way 3 was previously farmed and is now managed as an SAV system (Juston and DeBusk 2006).

Soil Sampling
Flocculent ( oc) and soil samples were collected from transect locations in STA-2 ow-ways 1 and 3 extending from water in ow points to the wetland interior. Each transect in STA-2 consisted of 11 stations, and three stations along each transect (in ow, mid ow, and out ow) were identi ed as benchmark stations where triplicate soil samples were collected. Soils were sampled using the push-core method consistent with methods used for prior wetland soil studies (Bruland et al 2007;Osborne et al 2011;Newman et al 2017). A 10-cm diameter polycarbonate core tube was pushed into the soil. Samples were extruded from the soil core tube and partitioned into oc and soil. Floc was characterized as the suspended unconsolidated material present on top of the consolidated soil and was segregated into a secondary sampling container, allowed to settle for 4 hours, and supernatant water was removed via aspiration with the remaining oc material being collected. The consolidated soil underneath the oc layer was sectioned into recently accreted soils (RAS) and pre-STA (sectioned into 0-5 and 5 to a maximum of 15 cm) based upon soil density and visual color differences. Sectioned soil samples were placed in heavy-duty air-tight plastic bags and stored at 4 o C until analysis. For purposes of this study, oc, RAS, and the 0-5 cm pre-STA layer were used for comparison.

Soil metal and P fractionation
Total phosphorus was determined by a combination of ignition at 550° C and acid digestion to convert organic P into inorganic P, followed by analysis for inorganic P in digestates by ascorbic acid techniques (US EPA 1993). Total inorganic P in soil was extracted with 1 M HCl (soil to solution ratio = 1:50; 3 hours extraction time), ltered using 0.45 µm lter, and analyzed for P (US EPA 1993). Water extractable P in soil was extracted with deionized water using soil to water ratio of 1:10, ltered using 0.45 µm lters, and analyzed for P (US EPA 1993). Total and HCl extractable aluminum (Al), iron (Fe), calcium (Ca), and magnesium (Mg) were extracted with 1 M HCl and analyzed using inductively coupled argon plasma spectroscopy (US EPA 1994). Oxalate-extractable Al, Fe, and P were determined by extraction with 0.1 M oxalic acid and 0.175 M ammonium oxalate (McKeague and Day 1966). The suspension was equilibrated for 4 hours in the dark with continuous shaking, centrifuged, ltered through 0.45 µm lters, and analyzed using inductively coupled argon plasma spectroscopy. All chemical analyses were based on wet weight and corrected for moisture content.
The bulk density of the different soil core sections (i.e., oc, RAS, and pre-STA soils) was determined by taking the dry weight of each material divided by the volume of the material (section depth or thickness × area of the core). Soil pH was determined using 1:2 soil and water suspension and loss-on-ignition (LOI) was calculated using percent ash.

Data Analysis
For purposes of this study, the phosphorus saturation ratio (PSR) was broken into two forms, PSR that considers Fe and Al (PSR FeAl ; Eq. 1a) and a PSR that considers Ca and Mg (PSR CaMg ; Eq. 1b). Therefore the total PSR (TPSR) was calculated as the sum of the molar ratio (mmol kg − 1 ) of HCl extracted P to HCl extracted Fe and Al and the ratio of HCl extracted P to HCl extracted Ca and Mg consistent with Eq. 1c (Nair and Reddy 2013;Nair et al. 2015) where HCl extracted P, Fe, Al, Ca and Mg (P HCl , Fe HCl, Al HCl, Ca HCl , and Mg HCl, respectively) are in units of mg kg − 1 . Similarly, soil phosphorus storage capacity (SPSC) expressed as mg kg − 1 was decomposed into its two dominant mineral associations, Eq. 2a, and 2b. Soil phosphorus storage capacity (SPSC x ) was calculated consistent with prior studies where PSR values compared to a PSR threshold (PSR thres,x ) corresponding to the critical P solution concentration, estimated from segmented regression ("segmented" R-package; Muggeo 2008) between PSR x and water extracted total dissolved P consistent prior studies (Nair and Harris 2004;Nair and Reddy 2013). A PSR FeAl threshold value of 0.57 (from segmented regression; Fig. 2) was used to estimate SPSC FeAl . Soil phosphorus storage capacity associated with Ca and Mg (SPSC CaMg ) was calculated similarly to SPSC FeAl (Eq. 2b) where the threshold PSR (PSR thres, CaMg ) with a value of 0.003 was used. The PSR thres, CaMg value was determined based on the segmented regression analysis of water-soluble extracted P and PSR CaMg (Fig. 2). Total SPSC was estimated as the sum of SPSC FeAl and SPSC CaMg for each sample.
While SPSC FeAl relates to both positive and negative storage, the SPSC associated with Ca and Mg will likely be primarily releasable P. In addition to calculating PSR FeAl using HCl extractable fractions, PSR FeAl was also calculated using the oxalate extractable fraction consistent with Eq. 1a to evaluate the inter-compatibility between HCl and oxalate PSR FeAl values. Oxalate extractable SPSC FeAl was also calculated with oxalate extractable PSR FeAl threshold value of 0.77 (from segmented regression; Fig. S1). Oxalate and HCl extractable PSR and SPSC were compared using spearman's rank-sum correlation and an SPSC oxalate-HCl conversion equation was developed using median-based linear model (Komsta 2013). Median-based linear model was used rather than least-squares methods due to the violation of linear model assumptions. All other comparisons were made using HCl extractable PSR FeAl and SPSC FeAl .

PSR
Total SPSC was compared to biogeochemical (oxalate extractable P, total P, N, C, Ca, Mg, Fe, and Al) and physical characteristics (pH, bulk density, % OM and distance downstream) using principal component analysis (PCA, "vegan" R-package; Oksanen and Guillaume 2018). Parameters used to calculate TSPSC (i.e., HCl extracted P, Fe, Al, Ca, and Mg) were excluded from the PCA analysis to reduce co-linearity. Prior to PCA analysis, data were evaluated using Kaiser-Meyer-Olkin measure of sampling adequacy ("REdaS" R-package; Maier 2015) and Bartlett's Test of Sphericity ("REdaS" R-package; Maier 2015).
Soil phosphorus storage capacity was compared across soil depth category and ow-way using Dunn's pairwise comparison. Phosphorus source/sink potential was evaluated by determining if SPSC values were signi cantly different from zero using the One-Sample Wilcoxon rank sign test. To evaluate how SPSC and soil TP change along each STA ow-way, SPSC and TP were compared to fractional distance downstream using spearman's rank-sum correlation. Finally, SPSC was compared to soil TP for each STA ow-way using spearman's rank-sum correlation.
All plots were generated in base R except for the CaMg, FeAl, and P ternary diagram which used the "ggtern" R-package (Hamilton and Ferry 2018). All tables were formatted using " extable" R-package (Gohel 2021). All statistical operations were performed using R (Ver 4.0.4, R Foundation for Statistical Computing, Vienna Austria). Unless otherwise stated, all statistical operations were performed using the base R library. The critical level of signi cance was set at α = 0.05.

Results
( ) ( ) ( ) ( ) As reported above, the segmented regression of PSR and water extracted total dissolved phosphorus (TDP) determined PSR thresholds of 0.003 ± 0.001 and 0.57 ± 0.19 (estimate ± standard error) for HCl extractable PSR thres, CaMg and PSR thres, FeAl , respectively. The upper and lower 95% con dence interval of the threshold for PSR CaMg was 0.001 and 0.004, respectively with an overall model adjusted R 2 value of 0.48. The upper and lower 95% con dence interval of the threshold for PSR FeAl was 0.20 and 0.94, respectively with an overall model adjusted R 2 value of 0.24. It is expected, given the distribution of the data (Fig. 2) that the overall model adjusted R 2 values would be low, however for purposes of this study the break-point where the slope of the line changes is of particular interest. Likewise, oxalate extractable PSR thres, FeAl was 0.77 ± 0.63 with an upper and lower 95% con dence interval of -0.48 and 2.02 and an overall R 2 value of 0.44. Segmented regression method is an iterative approach and for the PSR CaMg threshold, convergence was attained after one iteration (relative change 2.42x10 − 7 ), two iterations (relative change 2.38x10 − 7 ) for the PSR FeAl, HCl threshold, and ve iterations (relative change 1.40x10 − 7 ) for PSR FeAl, Ox .
Oxalate and HCl extractable PSR and SPSC were both positively correlated (r = 0.89, ρ < 0.01 and r = 0.88, ρ < 0.01) for all soil depths (Fig. 3). The relationship between the two calculated PSR indices appears to be non-linear with higher PSR values deviating from the 1:1 line. However, limited data is available in this region. Based on the median-based linear model, the conversion of oxalate and HCl extractable SPSC for the system in this study is nearly a 1:1 relationship with an equation of SPSC HCl = 0.81 × SPSC Ox + 34.00 and an R 2 of 0.79 (Fig. 3).
Overall, Ca HCl and Mg HCl account for a greater proportion by an order of magnitude of mineral-based elements relative to Fe HCl and Al HCl ( Fig. 4 and Table 1). While Ca HCl and Mg HCl are in greater proportion across samples, STA2 ow-way 1 soils contain a relatively greater range of Fe HCl and Al HCl when compared to STA2 ow-way 3. Moreover, STA2 ow-way 3 soils ( oc and RAS) are predominately Ca HCl and Mg HCl ( Fig. 4 and Table 1). Therefore, given the relative proportion of Ca HCl and Mg HCl to P HCl , PSR CaMg values are relatively low ranging from 0.0003 to 0.012 (Table S1). The proportion of Fe HCl and Al HCl was greater with PSR FeAl values range from 0.02 to 5.01 (Table S1). Across the two ow-ways, TSPSC values ranged from − 1090 to 1636 mg kg − 1 (Table 1) with all soil depth categories being signi cantly different than zero except for STA2 ow-way 1 RAS (Table 2).
( )  The KMO criterion of the dataset used in the PCA analysis was 0.75 indicating that the dataset was suitable for factor analysis. Bartlett's Test of Sphericity indicated that the data was signi cantly different from an identity matrix and suitable for factor analysis (χ 2 = 1730, df = 45, ρ < 0.01). The rst two components of the PCA accounted for 66% of the variance in the dataset with components 1 and 2 have eigenvalues greater than one (Fig. S2). In the multi-variate space (Fig. 5), soil depth and ow-ways are generally grouped into two distinct groups. Pre-STA would for both ow-ways group together in the upper right quadrant of the multivariate space (positively loading for PC1 and PC2). Based on the distribution in the ordination space, STA2 ow-way 3 is in uenced (loaded) by the soil TCa and TMg and pH while owway 1 is more in uenced by LOI, TC, and TN. This suggests differences in organic versus mineral (precipitation) biogeochemical pathways/mechanisms. Based on the loading of each variable (direction of arrows; Fig. 5), TSPSC is generally (positively) correlated with BD, TAl, TFe, and to some degree, fractional distance downstream and negatively correlated with P Ox and TP in the multivariate space ( Fig. 5). Meanwhile, TSPSC was generally not correlated with TCa, pH, TN, LOI, or TC. However, these components all correlated with one another to some degree.

Discussion
Wetland soils have a nite capacity to store P with increasing P loading potentially pushing the system from being a P sink to a source (Reddy et al 1999). This shift from a sink to a source can increase the risk of impacting downstream ecosystems via eutrophication (Nair and Reddy 2013). The mobility of P is balanced between P retention capacity and P buffering intensity of the soil (Reddy et al 1999), with soil physicochemical characteristics such as Fe and Al complexes (including amorphous and poorly crystalline forms), Ca and Mg-based mineral and organic matter being important in regulating P sorption and buffering intensity (Syers et al 1973;Chambers and Odum 1990;Reddy and DeLaune 2008). In general, P can be integrated with Fe and Al complexes while P is typically associated with (not bound) Ca and Mg minerals in wetland soils. This complexation versus association is ampli ed by physicochemical dynamics within any given system where uctuations in pH and redox conditions can cause rapid dissociation (Reddy and DeLaune 2008). Moreover, this complexation versus association is apparent in the PSR value observed in this study with PSR CaMg being orders of magnitude lower than that of PSR FeAl across ow-ways and soil depths indicating that P complexation with Fe and Al is much more stable. This relationship of complexation versus association between Fe + Al and Ca + Mg, is especially true for acidic soils however, Ige et al (2005) have demonstrated the use of a Ca and Mg DPS is more appropriate in neutral to alkaline soils to evaluate P sink/source in calcareous soils. Given the different processes of P retention between acidic (Fe + Al) and calcareous (Ca + Mg) soils, and in wetland soils in this study, easily releasable P from this "stored" P, indices of P saturation and retention (i.e. PSR, DPS, or SPSC) are not equitable between the different mineral groups.

PSR Threshold variability
The change point between PSR and water-soluble/water-extractable P (i.e. PSR thres ) is a soil-speci c threshold (Bhadha et al 2010). Several studies have used values ranging from 0.10 to 0.20 as a PSR threshold in predominately sandy soils Nair and Harris 2004;Nair et al 2015 Fe and Al (PSR FeAl ) was determined as 0.57 (Fig. 2) for HCl extractable and 0.77 (Fig. S1)  Differences in PSR thres, FeAl across studies could also be due to variability in soil composition, prior loading, and hydrologic conditions. As reviewed by Bhadha et al (2010) and Dieter et al (2015) repeated drying and wetting of soils can affect the capacity of soils to sequester P through changes to crystalline structure and oxidation of Fe. Meanwhile, P sorption by other metals such as Al, Ca and Mg does not involve oxidation-reduction reaction and is much more in uence by pH and organic matter content (Traina et al 1986;Diaz et al 1994;Staunton and Leprince 1996). It is likely that the hydrologic conditions of the STAs, where dry-out conditions are typically avoided, remain hydrated and contributed to the difference in PSR thres, FeAl relative to other studies discussed above.
Due to the relatively mineral nature of the soils in this study, PSR and SPSC estimates were evaluated with respect to HCl-extractable Ca and Mg. In acidic mineral soils, Fe and Al are more active in P-cycling and soil P buffering capacity tends to be higher relative to circumneutral Ca-rich soils (Richardson 1985;Richardson and Vaithiyanathan 1995). proportion of P to Ca and Mg in this study is relatively low, however, a change point at which watersoluble P starts to increase is apparent (Fig. 2). While PSR FeAl makes up a large percentage of the TPSR, signi cant variability between ow-ways is evident (Table 1 and S1) with Ca being an important component in abiotic retention of P in the STAs (Dierberg et al 2002;Gu 2008;Zamorano et al 2018).
Moreover, the P associated with Ca and Mg can be easily releasable relative to Fe and Al bound P.

SPSC Conversion
Numerous studies have used either oxalate or Mehlich extractable P fractions to estimate PSR and SPSC (Maguire and Sims 2002;Chrysostome et al 2007a;Bhadha et al 2010). However very few, if any, compared the indices from different extractions. This study compared HCl and oxalate extraction PSR and SPSC indices to i) ensure consistency between other studies, and ii) evaluate the potential for a simpler analysis for PSR and SPSC studies in the future. Both extraction procedures have their limitation and preferentially extract particular P fractions, for instance, HCl extracts inorganic P while oxalate P extracts both inorganic and organic P (Reddy and DeLaune 2008). Furthermore, HCl extracted P can be classi ed as the reactive P pool including Fe, Al, Ca and Mg bound P (Reddy et al 2020). In the case of this study, HCl and oxalate extractable PSR and SPSC were comparable (Fig. 3) with the potential of a non-linear relationship between oxalate and HCl PSR. This potential non-linear relationship could be due to differences in soil texture between the two ow-ways, the relative proportions of P (i.e., inorganic versus organic), the extraction e ciencies of the method, or a combination of these methods. While the indices are highly correlated caution should be used applying this across different systems, soil types, and/or nutrient conditions.

Sink-Source Dynamics
While PSR and the PSR thres are important metrics for environmental risk assessments, they do not provide an estimate of P storage capacity or P source/sink dynamics. Originally developed and tested for upland sandy soils Nair and Harris 2004), SPSC has also been applied to wetland soils (Mukherjee et al 2009;Bhadha et al 2010;Nair et al 2015). Nair et al (2015) concluded that soil organic matter in wetland soils did not affect P retention and release below the PSR thres despite the ability of organic matters' ability to complex with most metals and potentially lowering the a nity for P. Moreover, Nair et al (2015) also noted that wetlands that are heavily P impacted could be a P source where vegetation or other ecosystem components may no longer be able to assimilate addition P. These results are congruent with the results presented in this study where TSPSC was not correlated with organic matter (i.e. LOI) and inversely correlated with TP and P Ox (Figs. 5 and 7).
Wetlands soils are not a homogenous column of accreted material but rather layers of reactive biological material, decaying organic matter, and geologic material. The rst layer and possibly the most reactive is the oc layer, a complex matrix of microbes, organic matter, and inorganic material, oc is ecologically important in the aquatic ecosystem and is relevant to several biogeochemical processes (Noe et al 2001;Neto et al 2006). The source/sink dynamics of oc can be variable across systems as noted in this study ( Fig. 5 and Table 1). As observed in this study, EAV systems have a distinct partitioning of TSPSC between soil depths while SAV systems are less notable with respect to surface soils. Generally, oc has a relatively low bulk density (Table S1) and therefore a greater pore space providing the opportunity for more surface area for P to ux onto and from mineral and organic particles within the complex oc matrix.
As demonstrated by prior studies (Nair and Reddy 2013;Nair et al 2015) using SPSC, and this study using TSPSC, SPSC provides a potential indicator of P sink/source ( Fig. 7 and Table 2). This metric provides a relative risk of P mobilization for a particular area or soil compartment. As soils build, P storage capacity increases, suggesting that mineral components are largely responsible for P retention of the soil (Nair et al 2015). As observed in this study, generally TSPSC values for pre-STA soils were positive (Figs. 5 and 6), with SPSC FeAl making up a greater proportion of the overall TSPSC (Fig. S3, S4, and Table S1) suggesting P-retention via mineral association and a relatively stable P-sink. Floc or RAS P-retention is presumably associated with a combination of biological (i.e., microbial biomass) and mineral (to a lesser degree) processes. Additionally, pre-STA soils are relatively deep and effectively removed from active P-cycling. However, SPSC values reported in this study for pre-STA soils are relatively high compared to prior studies and presumably linked to the relatively high estimated PSR FeAl, thres . Regardless of using the estimated PSR FeAl,thres or literature reported values (i.e., 0.1; Fig. S5), pre-STA soils are consistently a sink of P whilst Floc is a source and RAS is variable across the ow-way. Moreover, it is worth noting that TSPSC only provides an index of soil mineral P storage capacity (Nair et al 2015) and does not account for biological storage capacity via microbial biomass, macrophytes, or organic matter.
The formation of P soil gradients in treatment wetlands can have important implications for the cycling of P within wetland ecosystems (Juston and DeBusk 2011;Zamorano et al 2018). Treatment wetlands remove a greater proportion of water column nutrients near the front of a given ow-way with downstream areas receiving comparatively less nutrient loading through biotic and abiotic retention of nutrients (Kadlec and Wallace 2009). In this study oc and RAS TP concentrations declined with distance downstream ( Fig. 7 and Table 3). Concurrent with this declining trend in TP concentration, P sink strength in oc and RAS also generally increased with distance downstream especially in ow-way 3 ( Fig. 7 and Table 3). These results are consistent with the conclusions by Nair et al. (2015), where regions with high TP concentrations (in ow) were also generally P sources (Fig. 6).
Biotic and abiotic processes regulate the retention and release of P from soils to the overlying water column (Reddy et al 2005). As discussed above abiotic processes such as soil Fe, Al, Ca and Mg, pH, and other factors regulate the retention of P, however, biotic factors also play a critical role in the retention capacity of a wetland to sequester P. The observed, at times, contrasting gradients of soil TP concentration and TSPSC highlight the variable P source and sink dynamic along the study ow-ways. In STA-2 ow-way 3, the shift of P source to P sink (Fig. 7) is also re ected in the aggregated internal load reported in a recent study that demonstrated internal loading can be a substantial contributor to water column P concentrations, especially at in ow regions of the ow-way (Jerauld et al 2020).

Conclusions
Wetlands are dynamic ecosystems that transform and retain nutrients through biotic and abiotic processes. For treatment wetlands, these processes are leveraged to improve water quality and protect downstream ecosystems. In the case of the Everglades, understanding and estimating P storage and release (i.e., source/sink dynamics) is vitally important to ongoing and future restoration efforts. Based on the SPSC, the potential for P loss from the oc could be greater in ow-way 1 while RAS in ow-way 1 is at or near storage equilibrium. The storage equilibrium concept is supported by corroborating P accretions rates reported for ow-way 1 (Bhomia et al 2015) and average annual P load retention rates (SFWMD, unpublished data) However, despite signi cant declines in soil TP within ow-way 1, TSPSC was not signi cantly correlated with distance downstream. Based on the results of the principal component analysis, this lack of signi cant TSPSC change within ow-way 1 could be due to biological retention and source/sink dynamics.
The potential for P loss within the STA should not be evaluated using surface soil samples only. Similar to the behavior of upland and wetlands soils elsewhere, TSPSC varied with soil depth across ow-ways.
In addition to changes with soil depth, the percent contribution of Fe, Al, Ca, and Mg change with soil depth as noted by the shift in SPSC FeAl contributing more to the overall storage capacity in Pre-STA soils.
The difference of TSPSC between oc and RAS could vary suggesting P loss from oc in one system ( ow-way 1) versus P retention in both oc and RAS in another ( ow-way 3). While P associated with Ca and Mg is easily releasable, TSPSC only provides an index of temporary soil mineral P storage capacity and does not account for other storage and retention processes. The contribution integrating Fe + Al and Ca + Mg to TSPSC provides some understanding of the storage strength of soil P relative to surface water P concentration at the present time. This is re ective in ow-way 3 where soils are calcitic in nature and surface water out ow total P concentrations are higher and with more variability relative to ow-way 1 (Villapando and King 2017).
Overall, the use of SPSC in a wetlands application can provide a metric of risk of P release (i.e. source/sink status) and provide evidence to inform management decisions. Within treatment wetlands or wetlands that receive high nutrient discharges, eutrophication gradients can develop where soil and other compartments (i.e., vegetation, water column) begin to become enriched with nutrients thereby reducing the overall e ciency of the system to remove and store P shifting from a sink to a source. The e ciency of this storage and retention of P as it relates to SPSC is also dependant on soil composition (i.e., organic versus mineral). Other factors such as soil and P accretion and loading may also play a role. While this was an initial study in the applicability of SPSC to treatment wetlands, additional investigation evaluating soil mineral composition, a ner soil depth pro le, and potentially other metrics of P retention is warranted. Soil sampling locations within the Everglades Stormwater Treatment Area (STA) ow-ways 1 and 3 (FW 1 and 3, respectively). STA-2 ow-way 1 is predominately emergent vegetation and STA-2 ow-way 3 is predominately submerged aquatic vegetation.

Figure 2
Relationship between water extracted total dissolved phosphorus (TDP) and calculated phosphorus saturation ratio (PSR) for iron and aluminum PSR (PSRFeAl; left) and calcium and magnesium (PSRCaMg; right) for stormwater treatment area 2 ow-ways 1 and 3 (STA-2 FW1 and FW3, respectively). Change-point, as determined by segment regression, are identi ed by the vertical dashed line.

Figure 4
Page 23/25 Ternary diagram comparing HCl extracted calcium and magnesium, iron and aluminum, and phosphorous (HCl Ca+Mg, HCl Fe+Al, and HCl P, respectively) for each stormwater treatment area (STA) ow-way and soil depth.

Figure 6
Total phosphorus storage capacity for compared across soil depth category within each stormwater treatment areas (STA) ow-way (FW). Letters above each boxplot indicate statistical similarity as determined by Dunn's multiple pair-wise comparisons within each STA FW (not between).

Figure 7
Mean (± standard error) soil total phosphorus (TP; top), total soil phosphorus storage capacity (TSPSC; middle) by fractional distance downstream and TSPSC compared to soil TP for each soil depth category along stormwater treatment area 2 ow-way 1 (STA2 ow-way 1; bottom left) and STA2 ow-way 3 (bottom right).

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download. SupplementalInformation.docx