Soil P Storage Capacity in Agricultural Treatment Wetlands: Can a System Designed for N Reduction Also Retain P?

Increasing interest focuses on utilizing wetlands to reduce nutrient loading to surface waters. The current study examines soil P storage capacity in three treatment wetlands designed to decrease N loading from an agricultural tile drainage system in Illinois. Adjacent farm field and restored floodplain wetland soils were also evaluated. Results demonstrate that wetland soils sequestered P; however, the magnitude of P retention varied significantly across treatment systems and floodplain wetlands related to differences in soil properties. Soil P storage capacity increased in the direction of water flow, but varied across treatment wetlands; soils ranged from P sinks (5.8 ± 1.5 mg P kg−1) to potential P sources (−17.2 ± 2.0 mg P kg−1). Farm fields displayed the highest water soluble P levels (11.3 ± 1.5 mg P kg−1) and represented the largest potential source of P with a mean storage capacity of −48.6 ± 6.8 mg P kg−1. A phosphorus saturation ratio threshold value of 0.10 differentiated between potential P sources and sinks. Findings suggest wetlands receiving P loadings from tile drainage accumulate soil P over time, but maintaining P removal efficiency in treatment wetlands may require periodic management to decrease soil P concentrations via nutrient removal and/or soil amendments.


Introduction
Phosphorus (P) is required for vegetative and microbial biomass production, providing a key nutrient to maintain both natural systems and profitable agricultural yields (Sims and Sharpley 2005). As a limiting nutrient in many freshwater systems; however, a surplus of P increases aquatic plant and algal growth that can result in water quality impairments and eutrophication in receiving waters (Sharpley et al. 1994). The environmental fate of P is of increasing concern due to the potential for excess P loadings to water bodies that degrade water quality and cause harmful algal blooms across broad spatial scales, including the Great Lakes region of the United States (Conley et al. 2009;Great Lakes Commission 2012). It has been well-established that extensive agricultural tile drainage in the Midwest is a major source of excess nitrate-nitrogen inputs into surface waterways (Goolsby et al. 1999;Sugg 2007); however, these drainage systems have also been identified as an effective conduit for P transport to freshwater river systems (Gentry et al. 2007;Smith et al. 2015). Reducing excess loss of P from non-point sources such as widespread tile drainage systems that support agricultural food and animal production is especially challenging (Daniel et al. 1998).
Wetlands have the capacity to intercept and retain excess nutrients from surface and groundwater sources and are increasingly implemented to reduce nitrogen (N) and P loading to surface water bodies (Vymazal 2007). Specifically, wetlands have been identified as a key strategy for reducing Ninduced hypoxia in the northern Gulf of Mexico (Mitsch et al. 2001a, b), in which microbially mediated denitrification converts excess surface water N to atmospheric N 2 gas (Gardner and White 2010). For example, Berkowitz et al. (2016) Electronic supplementary material The online version of this article (https://doi.org/10.1007/s13157-019-01205-3) contains supplementary material, which is available to authorized users. estimated annual average N removal rates via denitrification of 1073 kg N km −2 in several constructed wetlands in the Atchafalaya River. However, unlike N, the P cycle exhibits various pathways for nutrient sequestration (e.g., soil adsorption, organic matter assimilation and accretion, internal cycling) and outputs (surface or groundwater discharge), yet lacks a mechanism equivalent to denitrification for removing P. As a result, managing P movement on the landscape using wetlands requires an understanding of complex interactions between soil conditions, hydrology, and P sources, forms, transport, and fate. The extent and magnitude of wetland P retention remains largely determined by soil-water interactions (Reddy and DeLaune 2008) and these factors, with the inclusion of land use history, influence the management of created or restored wetlands used for P sequestration (Currie et al. 2017).
Water volume, P loading rates, and water retention time are commonly measured in treatment wetlands and used to assess short term P storage dynamics. However, long term P assimilation and retention is a function of soil P accretion rates and storage capacity . Soil properties that govern P retention and capacity are rarely measured despite finite soil P storage capacities, which determine the potential long term retention of P and associated water quality improvements (White and Reddy 2003). The concepts of soil P sequestration and the capacity of wetlands to improve water quality have been studied extensively in the Everglades (Reddy et al. 2011) but application remains limited in other areas (Dari et al. 2018).
Results of previous work indicate that the P storage capacity of wetland soils varies as a function of clay content and type, organic matter content, redox status, concentrations of Al, Fe, and Ca, and soil pH (Hansen et al. 2002). Because soil properties integrate long term conditions across space and time (Buol et al. 2001), evaluating soil P dynamics accounts for land use history, the presence of legacy P, and remaining P storage capacity (Reddy et al. 2011). Legacy P results from the accumulation of soil P over time from historical land use, primarily related to agricultural runoff (Sharpley et al. 2013;Nair et al. 2015). Phosphorus accumulation largely occurs through P sorption onto mineral soils, which accounts for the majority of long term P retention in wetlands (Richardson and Marshall 1986;Dunne et al. 2006). However, P concentrations can limit or exceed retention capacity of a soil, resulting in the development of potential P sources in wetlands that surpass critical P retention thresholds. As conservation practices integrate wetlands in agricultural landscapes to reduce nutrient loadings to water bodies, it will be imperative to consider historic landuse, legacy P, and remaining soil P storage capacity (Currie et al. 2017).
The concept of soil P storage capacity has not been widely applied in the Midwest region of the United States where constructed wetlands are increasingly being use to capture and treat excess nutrient runoff within agricultural watersheds. Our study examines the soil P storage capacity in constructed wetlands within an agricultural setting in Illinois, which has the highest estimated total area of subsurface drainage of any state in the upper Mississippi River Basin (Goolsby et al. 1999) and is among the highest contributors of total N (16.8%) and P (12.9%) flux to the Gulf of Mexico (Alexander et al. 2008). Illinois is among 12 states mandated by U.S. Environmental Protection Agency to develop nutrient loss reduction strategies to reduce N and P contributions to poor water quality and Gulf hypoxia. Analyses conducted during the development of the Illinois Nutrient Loss Reduction Strategy (NLRS) showed that 82% of N and 48% of P contributing to riverine nutrient export from the state originated from agricultural sources (IL NLRS 2015). Constructed wetlands have been identified in the strategic plan as an effective practice for capturing and treating excess N runoff; however, are currently designated as 0% effective for net reductions in P loads (IL NLRS 2015), likely due the variability of wetland P retention results in the literature (e.g., Kovacic et al. 2000).
Thus, the objectives of this study were to investigate: 1) the soil P storage capacity in constructed wetlands in use for over 10 years, and 2) assess the soil P storage capacity of sequential constructed wetlands designed to reduce N flows to surrounding surface waters. The P dynamics within treatment wetlands designed for N removal is of particular interest, as many studies focus on either reduction of N to address estuarine hypoxia in the northern Gulf of Mexico or freshwater water quality degradation related to excess P including the development of harmful algal blooms (Kovacic et al. 2006). Basic calculations to estimate P removal and P dynamics are also included to aid resource managers in estimating the longevity and effectiveness of treatment wetlands designed for reducing P in impaired watersheds.

Study Area and Methods
The study area is located on the Franklin Research and Demonstration Farm in Lexington, Illinois, where The Nature Conservancy has been granted access to research the effectiveness of conservation practices at reducing nutrient export from agricultural tile drainage systems (The Nature Conservancy 2014). Three experimental treatment wetland systems, East (E), Gully (G) and West (W), were constructed in 2005 (Fig. 1a). Each treatment wetland system consists of three sequential wetland cells where each cell represents 3% of the surrounding tile drainage area, for a total of 9% of the tile drainage area. The wetland systems were designed and built with the intention of reducing agricultural nutrients, specifically N from tile drainage installed at the farm which is maintained under row crops (e.g., corn, soybeans) using conservation tillage techniques. Monitoring of tile drainage inflow, wetland water levels, and nutrient concentrations began in 2007. In addition, a series of restored riparian wetland cells that receive drainage from the surrounding agricultural area and overbank flooding from an adjacent creek have also been established at the farm. The areal extent of Treatment area E, G, and W cells, and the riparian floodplain cells are available in Supplemental Table 1.
Soils within the study area consisted of fine, silty loam to silty clay loam Endoaquolls, Argiudolls, and Hapludolls (Drummer, Flanagan, Saybrook, Sawmill, and Lawson series) with slopes <2% (Soil Survey Staff 2018). Drainage classes ranged from somewhat poorly drained to poorly drained and each wetland sampling area contained ≥50% hydric soil map units; flooding and/or ponding remain frequent within the floodplain areas examined. During construction of the wetland treatment cells, soils were disturbed in order to establish berms and water control structures, limiting the applicability of soil map units for further analysis.
Soils within Treatment areas E, W, and G, Floodplains 1 and 2, and actively farmed soybean fields were sampled for physicochemical properties. Soil properties were evaluated by homogenizing triplicate soil samples located along a transect within the interior of each treatment wetland or floodplain subunit to minimize potential differences in inflow and outflow nutrient concentrations (Fig. 1b). The edge of actively farmed agricultural fields directly adjacent to each treatment wetland were also sampled. The top 7.5 cm of triplicate soil samples were collected within a 0.3 m 2 area, partitioned into 0-2.5 cm and 2.5-7.5 cm sections, representing the zone of highest soil P (Baker et al. 2017). Three locations along the transect were sampled except for cells in Treatment area G and Floodplain 2 cell 2, where cell size limited soil sampling to two locations. Forty four sampling locations and two soil depths resulted in a total of 88 soil samples. Aboveground biomass was removed from the soil surface prior to soil sampling. Soil analyses included moisture content, bulk density, soil organic matter, total carbon (C), total nitrogen (N), total P, water soluble P (WSP), and Mehlich-1 extractable P, Fe, and Al.
Dry weight bulk density was calculated as the total dry weight of soil, corrected for moisture content at 105°C, divided by the volume of the soil core. Soil organic matter was determined using the loss on ignition (LOI) method, where dried, ground soil subsamples were combusted at 550°C for 4 h (Sparks 1996). Total C and N were measured using an Elemental Combustion System (Costech Analytical Technologies, Inc., Valencia, CA) on dried, ground soil subsamples. An additional soil subsample was air dried and ground for WSP and Mehlich-1 extractions.
Water soluble P (WSP) was determined by deionized water extraction for one hour followed by filtration through a 0.45 μm membrane filter (Kuo 1996;Richardson and Reddy 2013) and analyzed using a SEAL AQ2 Automated Discrete Analyzer (SEAL Analyzer, West Sussex, England; AGR-201-A Rev. 1). Method detection limits were 0.015 mg P L −1 for WSP. Mehlich-1 extractable P, Fe, and Al were determined by Mehlich-1 soil solution (0.025 mol L −1 H 2 SO 4 and 0.05 mol L −1 HCl) extraction for five minutes followed by filtration through a Whatman No. 42 filter (Mehlich 1953). Mehlich-1 extractions were sent to the Wetland Biogeochemistry Analytical Services Lab (Louisiana State University, Baton Rouge, LA) for extractable P, Fe, and Al analysis using ICP (Loeppert and Inskeep 1996). Total P was also analyzed at the Wetland Biogeochemistry Analytical Laboratory using the ashingdigestion method (Andersen 1976  The soil P storage capacity factor (CF M1 ) and P saturation ratio (PSR M1 ) was determined using the following formulas (Nair et al. 2013), where the threshold PSR M1 of 0.10 was utilized according to Nair et al. (2004) and Dari et al. (2018): Where the PSR M1 is calculated as the molar ratio of Mehlich-1 extractable P to Mehlich-1 extractable Fe and Al: A univariate generalized linear model evaluated statistical differences in soil properties, in which the dependent variables (e.g., bulk density; LOI; TN; TC; TP; Mehlich-1 P, Fe, Al; WSP; PSR M1 ; CF M1 ) were compared between sample locations (e.g., Treatment area E; Treatment area G; Treatment area W; Floodplain 1; Floodplain 2; Farm field). Treatment cell was included as a random blocking factor to control for the lack of independence in samples collected within sequential treatment cells. The Games-Howell post hoc test was applied to account for uneven sample sizes and associated lack of homogeneity of variance in some parameters.
Differences between sequential treatment cells within each sample location (e.g., Treatment area E cell 1; Treatment area E cell 2; Treatment area E cell 3) were evaluated using analysis of variance (ANOVA) following testing for normality (Shapiro and Wilk 1965) and homogeneity of variance (Levene 1960 test). The Games-Howell post hoc test was applied due to a lack of homogeneity of variance in some parameters. The non-parametric Kruskal-Wallis test (denoted as P′) evaluated critical thresholds in WSP, CF M1 , and PSR M1 due to violation of normality assumptions required for ANOVA. All statistical significance was evaluated at the α = 0.05 level using SPSS version 25.0 (IBM, Inc). In cases where soil parameters did not differ with depth, the depth intervals were pooled for subsequent analysis (Table 1).
Although not significantly different between surface and shallow subsurface depths (P = 0.10), mean soil WSP was slightly higher in surface samples (4.8 mg P kg −1 ) than shallow subsurface soils (3.14 mg P kg −1 ). Water soluble P displayed a range of values in both treatment wetlands and floodplains, remaining within the range reported for Mollisols in the region (1-21 mg P kg −1 ; Sui et al. 1999) and WSP concentrations from a variety of Alfisols, Entisols, Inceptisols, Spodisols, and Ultisols from other regions in the United States (0.2-25 mg P kg −1 ; Dari et al. 2018).
Soils with higher Fe and Al contained higher P due to the associated high P sorption capacity of Fe and Al (Table 1; Hansen et al. 2002). For example, Treatment area G displayed significantly lower Mehlich-1 extractable Fe and Al than wetland Treatment areas E and W (P < 0.001), coinciding with significantly less Mehlich-1 extractable P (P < 0.001). The current data set displays a wider range of Mehlich-1 extractable P values (0.2-63 mg P kg −1 ) than Allen et al. (2006) reported in five common Midwest soils (10-27 mg P kg −1 ) using Mehlich-3 extraction. However, Dari et al. (2018) found similar Mehlich-1 extractable P ranging from 16 to 81 mg P kg −1 in soils from southeast and mid-Atlantic regions where a variety of soil orders were sampled. A wide range of Mehlich-1 extractable Fe and Al in the current study was also observed ( Table 1). The range in Mehlich-1 extractable P, Fe, and Al may be attributed to a variety of factors including soil texture (Hansen et al. 2002;Sekhon et al. 2014), organic matter content (Reddy et al. 1998), legacy P and land management (Sharpley et al. 2013), treatment of agricultural drainage (Braskerud 2002), and soil redox status (Nair et al. 2015).
Soil PSR M1 , and CF M1 did not significantly differ between the surface and shallow subsurface depth intervals examined (P = 0.28-0.66). However, PSR M1 was slightly higher in surface samples than shallow subsurface soils (0.48 vs. 0.41) and CF M1 was lower in surface soils (−17.2 mg P kg −1 ) than in shallow subsurface samples (−13.3 mg P kg −1 ). These results are similar to Baker et al. (2017) in that PSR M1 and CF M1 declined with depth, though the decline was not significantly different. The lack of significant difference was likely an artifact of sampling the upper 7.5 cm of soil in this study, rather than Baker et al. (2017) and others who reported vertical stratification of soil P within the upper 5 cm of the soil surface in the midwestern United States compared with samples collected deeper in the soil profile.
Analysis of soils in the Midwest region often utilize Mehlich-3 or Bray and Kurtz P-1 extractions to determine plant available P (Kovar and Pierzynski 2009); however, previous studies applied Mehlich-1 extractions to investigate P sorption dynamics in various U.S. regions such as the Mid-Atlantic, the Ozark Plateau, the Southeastern Piedmont, Gulf Coastal Plain, and Southern Coastal Plain (Dari et al. 2018;Mukherjee et al. 2009;Nair et al. 2004). As a result, Mehlich-1 extractions were selected here in order to support comparisons from previous studies, which have mainly been conducted in the southeast and mid-Atlantic regions, and to apply established critical threshold and P storage capacity factor methodologies to soils in the Midwest including Mollisols, which have not been extensively evaluated.
Recently, a number of states have moved toward adoption of Mehlich-3 extract solutions due to improved correlations between plant available P and soil test results (Penn et al. 2018). In response, future studies of soil P dynamics should consider using Mehlich-3 extractions since that appears to be the general trend in agronomy and soil science. Relationships between Mehlich-1 and Mehlich-3 test results have been established, with Mehlich-3 extracting more Al, Fe, and P than Mehlich-1 (Mylavarapu et al. 2002). As a result, the data presented herein likely represent an underestimate of P sorption rates compared with studies utilizing Mehlich-3 extracts (Dari et al. 2018); however, the relative differences observed between treatment wetlands and individual treatment cells are not expected to change based upon the extractant selected.

Treatment Area Effects
Examining the wetland Treatment areas (E, G, W), WSP was significantly higher in Treatment area E (4.9 ± 0.7 mg P kg −1 ) than the other two treatment areas (0.40 ± 0.06 mg P kg −1 ; P = 0.002-0.007; Fig. 2). The WSP of Treatment area E was also higher than one of the floodplain wetlands (0.39 ± 0.05 mg P kg −1 ; P < 0.001), but not significantly different than the second floodplain wetland (4.4 ± 0.8 mg P kg −1 ). The farm fields displayed highest WSP within the dataset (11.3 ± 1.5 mg P kg −1 ) remaining significantly higher than all other sample locations (P < 0.001). The WSP in the farm fields was within the range reported in other upland soils under plant and animal production systems, where WSP ranged from 0.2-25 mg P L −1 (see Fig. 1 in Dari et al. 2018).
Treatment area W (0.07 ± 0.01) displayed lower PSR M1 values than the other treatment wetlands, resulting from significantly higher levels of extractable Fe and/or Al (P < 0.002). The restored floodplain areas also exhibited a range of PSR M1 values, with Floodplain 1 significantly lower (0.16 ± 0.01; P = 0.002) than the Floodplain 2 wetland (0.42 ± 0.04). The farm field samples displayed intermittent PSR M1 values (0.42 ± 0.04). Soil PSR M1 in Treatment area G was highly variable (range = 0.22-3.62) as a result of low soil Fe and Al relative to P concentrations at some locations.
Soil P storage capacity factors also varied across the treatment areas, with positive, negative, and near-zero values observed (Fig. 2). Specifically, Treatment area E contained a significantly lower storage capacity factor (P < 0.001) than the other two treatment wetlands, with the negative capacity (−17.3 ± 2.1 mg P kg −1 ) suggesting a potential source of P to receiving water bodies. Soil P storage capacity factors remained near zero at Treatment area G (−0.41 ± 0.11 mg P kg −1 ), suggesting a potential flow through area where limited additional P will be retained in surface soils. Treatment area W displayed positive storage capacity (5.6 ± 1.5 mg P kg −1 ), with additional P assimilation sites likely available to sequester incoming P due to high extractable Al and Fe in conjunction with relatively low P. These findings suggest that treatment wetlands display a gradient of potential P sources, sinks and flow through areas, with P fate related to soil conditions such as Fe, Al, and organic matter content (Hansen et al. 2002) and past management activities related to tillage and fertilizer application (Sharpley et al. 2013). Floodplain wetlands also displayed highly variable P storage capacity factors, ranging from near zero values to −25 mg P kg −1 . Unlike the treatment wetland areas which are subject to directed tile drain inputs, the floodplains receive variable hydrologic and particulate contributions (and potential P sources) from overland flow, groundwater discharge, and overbank flooding from the adjacent creek. As a result, the variability in floodplain soil P dynamics is not unexpected. Farm fields displayed the lowest storage capacity factor (−48.7 ± 6.8 mg P kg −1 ) within the dataset, likely associated with the historical and current surface fertilizer applications common throughout the region.  ; Fig. 3). The critical threshold PSR M1 of 0.10 also differentiated between soil CF M1 , with PSR M1 < 0.1 displaying significantly higher P storage capacity (7.7 ± 2.7 mg P kg −1 ; P′ ≤ 0.001) than soils with PSR M1 > 0.1 (−20.0 ± 12.1 mg P kg −1 ). As reported in other studies WSP was also closely related to soil CF M1 (Fig. 3), where WSP increased significantly when CF M1 was below zero (P′ ≤ 0.011; Dari et al. 2018;Nair et al. 2004). The current data set shares similar critical threshold PSR M1, WSP-CF M1 relationships, and magnitude of CF M1 with other studies. For example, Dari et al. (2018) reports a range of CF M1 values between −25 and 18 mg P kg −1 across two soil types in the southeastern Piedmont (35 soil samples). Nair et al. (2004) evaluated sandy Florida soils impacted by dairy production, which yielded CF M1 values as low as −227 mg P kg −1 , and four unimpacted soils with CF M1 between 4.2 and 11.7 mg P kg −1 . Additionally, Nair and Harris (2014) compared a wetland soil to an adjacent P impacted ditch soil in Florida which yielded CF M1 values of −118 and − 69 mg P kg −1 , respectively, in the top 10 cm of soil. Our data (CF M1 range = −118 -16.8 mg P kg −1 ; 88 soil samples) falls within the general range of values observed in other studies Fig. 2 Mehlich-1 extractable elements (P, Fe, and Al), WSP, PSR M1 , and P storage capacity factors in treatment wetlands, floodplain wetlands, and surrounding farm fields. Lower case letters indicate statistically significant differences; error bars display one standard error of the mean including those within different soil orders and textures (Nair et al. 2004;Dari et al. 2018), landuses (i.e. dairy sprayfields; Nair et al. 2004), and wetland types (Nair and Harris 2014),

Treatment Cell Effects
Mehlich-1 extractable soil P significantly decreased (P < 0.001-0.031; Fig. 4) as water moved through the sequential wetland treatment cells, while soil Fe and Al largely remained constant between cells. This suggests that more P storage sites have been occupied in the first treatment cell relative to the second and third cells. The saturation of initial P treatment cells has been observed elsewhere, including within the inflow areas of Stormwater Treatment Areas and at water control structures utilized during restoration initiatives in the Florida Everglades (Reddy et al. 2011). These trends are also reflected in the WSP and P storage capacity factor data, which show significantly more WSP (P = 0.001-0.037) and significantly lower storage capacity factors (P = 0.001-0.044) in the initial wetland treatment cell than subsequent treatment cells (Fig. 5). For example, Treatment area W displays a near zero storage capacity factor (−0.87 ± 1.6 mg P kg −1 ) within the first cell indicating potential pass through conditions, but exhibits significantly more capacity in the second (7.8 ± 0.5 mg P kg −1 ) and third (10.4 ± 2.1 mg P kg −1 ) treatment cells. Treatment area E also shows increases in P storage capacity across the three wetland cells, but each cell remains below zero (−24 to −8.5 mg P kg −1 ) suggesting limited additional P sequestration may occur without management intervention.

Treatment Cell Performance and Management Implications
The treatment areas have proven effective in reducing dissolved nutrients (e.g., TNC 2014) despite varied P storage capacities. The P storage capacity; however, only accounts for P sorption onto mineral soil sites (i.e. Fe and Al), and does not consider P retention by other mechanisms such as incorporation into soil organic matter or particulate entrapment . Wetlands are effective at capturing particulate matter and can be a significant P retention pathway as P is often associated with particles Kadlec and Wallace 2008). Although we acknowledge particle entrapment as a potentially additional reduction pathway in these treatment areas, we did not quantify this P retention mechanism herein. Capturing particulate matter may be a particularly important P retention pathway in the first cell, where any eroded soil is most likely to accumulate as inflow waters decrease velocity promoting particle settling. A subset of sites suggest that particulate P retention may be occurring in the first cells, with Treatment area G and Floodplain 1 displaying higher total soil P in the first cell than in subsequent cells (P = 0.001-0.004; data not shown). Treatment area W also displayed higher total soil P in the first cell, but results were not statistically significant. Additional research is required to further investigate the relationship between P sorption, particulate, and dissolved P in treatment wetlands receiving water from tile drainage (e.g., Treatment areas E, G, and W) and those linked with surfacewater, groundwater, and overbank flooding inputs (e.g., Floodplains 1 and 2). Fig. 3 Relationship between WSP and PSR M1 identifying the critical threshold value of 0.10 above which increases in WSP are observed (left panel). Note that the WSP increased significantly (P′ = 0.005) above the threshold value; similar to results reported by with Nair et al. (2004) from soils in Florida and Dari et al. (2018) who evaluated soils from multiple regions within the United States. Relationship between P storage capacity factor (CF M1 ) and WSP (r 2 = 0.87) also followed expected patterns, as WSP increased significantly when CF M1 was below zero (P′ = 0.011; right panel) Resource managers continue to express interest in determining nutrient retention efficiencies of treatment wetlands, including developing trajectories designed to inform treatment wetland lifespans (Gale et al. 1994). Additionally, evaluating the current and potential future nutrient dynamics in these systems allows for implementation of management strategies to achieve nutrient reduction goals. For example, Vymazal (2007) indicates that frequent biomass harvesting improves nutrient retention in treatment wetlands and others suggest amending wetland soils to improve water quality outcomes (Ann et al. 1999;Westholm 2006).
In response, the following section estimates the past P retention performance of the Franklin Farm wetlands and suggests P retention trajectories based upon available data. Notably, these calculations make a number of assumptions, as pre-construction soil data at the study site was not available and data only examines the upper 7.5 cm of the soil profile. In order to address the lack of pre-construction data, we estimate the P retention of treatment cells based upon the difference between Mehlich-1 extractable P in the first (inflow) or second cell (intermediate) and the last (outflow) cell. Also, saturation and P sorption likely occurs at depths beyond 7.5 cm, but ground water hydrology data were not available. This approach provides a conservative estimate of P removal, since significant P has likely been accumulated in the third treatment cell and soil saturation (and P sorption) likely occurs below 7.5 cm. Additionally, we utilize differences in Mehlich-1 extractable P for the analysis to represent P likely sorbed following treatment wetland construction, analysis conducted using stronger extraction solutions (e.g., Mehlich-3) would yield higher P sorption estimates. As noted above, the Mechlich-1 extractable P decreases in subsequent treatment cells while extractable Fe and Al remain consistent across cells, suggesting that differences in soil P are due to post construction retention of P.
Soil P accumulation and trajectory calculations were conducted using the following equations: Capacity for future P retention=release kg ð Þ in a Treatment cell ¼ bulk density cellx *soil depth cellx *cell area cellx *CF cellx ð Þ This approach has several limitations, including the lack of pre-construction data, the inability to differentiate between dissolved and particulate P retention, the lack of soil P data within a deeper subsection of the soil profile, and limited knowledge of saturation depth within the treatment wetlands. However, the following allows for an approximation of past near-surface P retention, comparison between treatment wetlands, and promotes the application of treatment wetland lifecycle analysis to improve optimization of P retention in watersheds with poor water quality.
Mehlich-1 extractable P results indicate that each of the treatment wetlands has accumulated P in surface soils. The three treatment wetlands have cumulatively retained a total of 12.9 kg P in the upper 7.5 cm, including 9.9 kg P in the initial treatment cells and 2.9 kg P within the second treatment cells. Examining individual wetlands, estimates yield P removal of 8.0 kg P, 0.04 kg P, and 4.8 kg P in Treatment areas E, G, and W, respectively. The variability in P retention between treatment areas results from a combination of soil characteristics (available Fe and Al discussed above) and the area of the treatment cells. For example, the first two cells of Treatment area E occupy an area of approximately 3240 m 2 , much larger than the first two cells of Treatment area G (930 m 2 ; Supplemental Table 1). Additionally, Treatment area Fig. 4 Mehlich-1 extractable P, Fe, and Al across the direction of flow in wetland treatment cells. Note Mehlich-1 extractable P decreases across sequential treatment cells, whereas Mehlich-1 extractable Fe and Al remain relatively constant between cells. Lower case letters indicate statistically significant differences; error bars display one standard error of the mean E contains significantly more extractable Fe and Al than Treatment area G (Table 1; Figs. 2 and 4).
As noted, additional P removal is likely occurring in the third treatment cells and at depths below 7.5 cm, resulting in lower P removal rates compared to other studies. For example, during the 10 years of operation Treatment area E has removed an estimated 0.49 g P m −2 yr. −1 ; much lower than the 3.0-30 g P m −2 yr. −1 Craft (1996) reported for estuarine constructed wetlands and the 26-71 g P m −2 yr. −1 observed by Braskerud (2002) in freshwater treatment wetlands receiving agricultural runoff in Norway. However, given the lack of pre-construction data soil conditions and saturation depths, these data serve as examples of how P retention calculations can be completed when evaluating P retention dynamics. Results also highlight the need to collect soil samples during the project design and monitoring of hydropatterns to improve nutrient retention accounting. Even if soils can be archived for later analysis, the capacity to determine soil conditions before and after project implementation is essential for determining long term treatment system effectiveness and justifying the costs of treatment wetland construction, management, and monitoring.
Examining potential future P dynamics using the capacity factors presented herein, Treatment areas E and G have the Results from WSP (top panels) and CF M1 (bottom panels) in wetland treatment cells, floodplain wetlands, and adjacent farm fields. Treatment cells (left panels) exhibit WSP decreases and CF M1 increases in the direction of flow. Floodplain wetlands (right panels) displayed a range of values and farm fields contained high levels of WSP suggesting they are a potential sources of P (i.e., negative CF M1 values) compared to other study locations. Note the different scales between panels. Lower case letters indicate statistically significant differences; error bars display one standard error of the mean capacity to release 16.8 and 0.78 kg of P, respectively from the upper 7.5 cm of soil. As noted elsewhere soil CF M1 and potential P release decreases in sequential treatment cells, as seen in Treatment area E cell 1 (7.9 kg P), cell 2 (5.9 kg P), and cell 3 (3.0 kg P). Interestingly, Treatment area W exhibits the capacity to release 0.90 kg P from cell 1 and the potential to adsorb 4.1 and 3.9 kg P from cells 2 and 3 respectively; yielding an overall potential P sink of 7.1 kg. Notably, these values do not directly translate to P retention or release, as the rates of P movement are dictated by a number of factors including P content in inflow water, retention time, redox status, sediment deposition/erosion, plant nutrient uptake, and other properties (Kadlec 2005). Additional research is required to further quantify the fate of retained P in treatment wetlands once storage capacity has been exceeded. This discussion should include the potential use of soil amendments to increase available soil P sorption sites, periodic flushing of treatment wetlands with low nutrient waters to desorb P during select time periods, management of plant biomass to decrease P content, and soil excavation for redistribution onto active agricultural fields to improve P retention while recycling sequestered P and decreasing fertilizer input costs.

Conclusion
Treatment wetlands are increasingly being used for retention of pollutants and excess nutrient reduction, often targeting a specific water quality parameter. Study results suggest that treatment wetlands designed for N retention have sequestered incoming P, and now display a range of potential P sources, sinks, and flow through systems. Notably, the treatment wetlands provide an example of ongoing P saturation in which soil P concentrations decrease and P storage capacities increase in the direction of tile drainage flow, with conditions improving in subsequent treatment cells. Mass balance accounting of P retention and future P trajectories are limited by the lack of pre-construction data, highlighting the need to consider long term monitoring during project design and implementation. The evaluation of soil properties including storage capacities and critical thresholds must be expanded to other regions and soil types as additional research seeks to expand and optimize the use of wetlands for P reduction and water quality improvements. Specifically, evaluation of Mollisols and other soils within the Midwest and Great Lakes regions will be required to achieve the P reduction goals outlined by the Great Lakes Restoration Initiative and other areas in the United States and abroad.