Dynamics of dissolved reactive phosphorus loss from phosphorus source and sink soils in tile-drained systems

Results and Discussion

Annual Dissolved Reactive Phosphorus Loads and Flow-Weighted Dissolved Reactive Phosphorus Concentration. Across plots and years, annual DRP loads measured in tile discharge ranged from 0 to 27 g ha^–1 (table 4). Also, in most instances, annual DRP loads from P source soils (plots 10, 11, and 12) were an order of magnitude higher than annual DRP loads from P sink soils (remaining plots) but, in some years, plots 26 and 43 had DRP loads as high as P source soils (table 4). In general, these loads were low compared to average annual DRP loads of 15 to 99 g ha^–1 reported in the same study site for water years 2003 and 2005 (Hernandez-Ramirez et al. 2011). Nevertheless, nutrient loads lost from crop production systems are dependent on discharge (Williams et al. 2014). Therefore, a possible contributing factor for the differences in DRP loads for the treatments at the WQFS are the higher annual tile drain discharges of 25 to 81 m³ reported in Hernandez-Ramirez et al. (2011) compared to 10 to 63 m³ reported in Welikhe et al. (2020a).

Annual FDRP concentrations from the study plots ranged from 2 to 23 ug L^–1 (table 4) and were lower than those observed in the same study site, i.e., 4 to 67 ug L^–1, by Hernandez-Ramirez et al. (2011). Also, these FDRP concentrations were low compared to FDRP concentrations of 80 to 160 ug L^–1 and 58 to 231 ug L^–1, reported for maize–soybean (Glycine max L.) rotations on a Bennington silt loam and a Pewamo clay loam in Ohio, and on mollisols in Illinois, respectively (King et al. 2015; Gentry et al. 2007). In the three previous study sites (WQFS, Ohio, and Illinois), P fertilizer was applied every other year at planting with maize to maintain optimum Bray P1/Mehlich 3P soil test levels for crop growth (King et al. 2015; Hernandez-Ramirez et al. 2011; Gentry et al. 2007). In previous, related work, Welikhe et al. (2020a) showed that soils with negative SPSC values, which correspond to soils with P levels above the critical soil test P (STP) level of 22 mg P kg^–1 (table 2), were more prone to desorbing P and had a greater risk of losing DRP to the tile drains. Furthermore, the authors showed that P recommendations based on a build-up and maintenance approach lowered a soil’s SPSC value below the threshold SPSC of 0. Given neither King et al. (2015), Hernandez-Ramirez et al. (2011), nor Gentry et al. (2007) reported SPSC and STP levels, direct comparisons between P status of soils in these studies versus the present study could not be made. Nevertheless, it is possible that the approach to P fertilizer management used by King et al. (2015), Hernandez-Ramirez et al. (2011), and Gentry et al. (2007) turned the soils into P sources as highlighted in Welikhe et al. (2020a). Soil type is another possible reason for the differences in FDRP concentrations lost to tile discharge among the studies. Pewamo clay loam is the major soil (~60%) in the tile drain contributing areas examined by King et al. 2015. Soils with high clay content have a greater risk of losing DRP to tile drains due to their tendency to develop preferential flow pathways (Simard et al. 2000), which could explain the higher FDRP losses in those soils. In the water years 2012 and 2013, FDRP concentrations >20 ug L^–1 were measured in tile discharge from plot 11 (table 4). Correll (1999) and Sharpley et al. (2003) reported that 20 ug P L^–1 was the concentration above which eutrophication was accelerated. However, all annual FDRP concentrations (table 4) were below USEPA’s 50 ug P L^–1 acceptable water quality limit (USEPA 2002).

Discharge Events: General Description, Dissolved Reactive Phosphorus Loads, and Flow-Weighted Dissolved Reactive Phosphorus Concentrations. Based on our discharge event criteria, a total of 75 events (30 on P source soils [plots 10, 11, and 12] and 45 on P sink soils [all other plots]) were selected (figure 2). Most of the events took place during fall (September 22 to December 20), winter (December 21 to March 19), and spring (March 20 to June 20). During summer (June 21 to September 21), there was no discharge in most tiles except for tiles 30, 43, and 44 in the summer of 2011. This discharge pattern is common in the region as a result of the interactions between precipitation, runoff, infiltration, and evapotranspiration (Gentry et al. 2007). Event discharge lasted an average of three or four days across tiles. Eight events lasted longer than the average duration (>4 days per tile). The Q and Q_max varied between 145 and 7,993 L d^–1 plot^–1, and 77 and 4,048 L d^–1 plot^–1. Across events, Q_ave ranged between 36 and 1,995 L d^–1 plot^–1. Forty-eight percent of the discharge events had D_rel values of 0%, indicating an abrupt rising limb. Previous studies in the region reported similar findings where tile discharge peaked early and abruptly during storm events (Gentry et al. 2000; Schilling and Helmers 2008). Data for all general characteristics, DRP, and FDRP concentrations for each event can be found in supplementary materials (table S1).

• Download figure
• Open in new tab
• Download powerpoint

Figure 2

A bubble plot of flow-weighted dissolved reactive phosphorus (FWDRP) concentrations (ug L^–1), total discharge (L d^–1 plot^–1), and dissolved reactive phosphorus (DRP) loads (g ha^–1) for the selected discharge events between October of 2010 and September of 2013 at the Water Quality Field Station (WQFS). The dashed line marks the concentration (20 ug L^–1) above which eutrophication is accelerated (Correll 1999).

Across all tiles, DRP loads and FDRP concentrations measured during selected discharge events ranged between 0 and 8.9 g P ha^–1 and 2 and 45 ug P L^–1, respectively (figure 2). There was no distinct trend of loss for both DRP load and FDRP in successive events in either P source or sink soils. However, when compared to P source soils, P sink soils had very low (some 0 g ha^–1) event DRP loads. Also, average FDRP concentrations from P source soils (with the exception of plot 10) were an order of magnitude higher than average FDRP concentrations from P sink soils. The event FDRP concentrations were below USEPA’s 50 ug L^–1 acceptable limit (USEPA 2002), with a few events (six events) having concentrations above the eutrophication acceleration limit 20 ug P L^–1 (Correll 1999; Sharpley et al. 2003). It is important to note that tile drainage can directly transport DRP to surface water bodies (e.g., stream, river, and lakes), bypassing the soil matrix where DRP could be sequestered (Smith et al. 2015; Reid et al. 2012). However, since most fields are at some distance from surface water bodies, discharge at the edge-of-field (drain outlets) must travel through drainage ditches, neighboring fields, or buffer areas where attenuation or amplification of the DRP losses could occur (Reid et al. 2018; Haygarth et al. 2005; McDowell et al. 2004).

Concentration-Discharge Relationships. The daily DRP flux had a linear relationship with daily discharge for both P source (R² = 0.61) and sink (R² = 0.73) soils when data were log-transformed (figure 3). The slope (b = 0.93) of the daily DRP load-discharge relationship for source soils had a 95% confidence limit between 0.84 and 1.01, indicating that b was not significantly different from 1 and that DRP concentrations vary little compared to discharge (figure 3a). According to Godsey et al. (2009), this observation of chemostasis implies that there is a mechanism at play in the catchment that is buffering variations in solute concentrations over a large range of discharge. Previous studies have shown that possible mechanisms include continuous weathering of bedrock for solutes such as magnesium (Mg), sodium (Na), and calcium (Ca) (Godsey et al. 2009), or the continuous release of legacy nutrients (e.g., NO₃^– and P) in agricultural soils (Basu et al. 2010, 2011; Thompson et al. 2011). In Menezes-Blackburn et al. (2016), P desorption rates were positively correlated (r ≥ 0.5; p ≤ 0.01) with different STP concentrations in different tests, i.e., FeO strips, Olsen, oxalate, and NaOH-EDTA, suggesting that increasing STP level was the main driver of P resupply from solid phase P into soil solution, and hence P bioavailability and lability. Further, Welikhe et al. (2020a) showed that once STP levels increase above an environmental STP threshold of 22 mg kg^–1 (which corresponds to an SPSC threshold of 0), the soils turn into P source soils with an eight-fold greater risk of solid-phase P release into the soil solution. Sharpley et al. (1994) showed that these chronic transfers of solid phase P to runoff is symptomatic of long-term, repeated fertilizer and manure applications, which can lead to increases in solid phase P and ensuing desorption of P to runoff water. Thus, it is well-established that accumulated P in soils and DRP in surface and subsurface runoff are strongly tied (Vadas et al. 2005; Welikhe et al. 2020a). In this study, soils classified as P source soils were those with negative SPSC values. According to Nair et al. (2015), the magnitude of negative SPSC values are estimations of the loosely held legacy P stock in soils most prone to desorption into soil solution. Therefore, similar to other catchments with legacy nutrient stores, our P source soils acted as mass stores that continually buffered DRP concentrations against changing discharge. The chemostatic behavior in P source soils provides evidence to support the hypothesis of Williams et al. (2016) that P source soils could act as a source of labile P similar to surface-applied fertilizers, resulting in continuous delivery of DRP to tile drains throughout the year.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3

Daily dissolved reactive phosphorus (DRP) flux (kg ha^–1) to discharge (mm) relationship presented on a log-log scale for (a) phosphorus source soils and (b) phosphorus sink soils. Exponent (b) and coefficient (a) values are defined by equation 1, and R² values were determined by linear regression on log transformed dissolved reactive phosphorus flux and discharge values.

The slope (b = 0.92) of the daily DRP concentrations-discharge relationship for P sink soils had a 95% confidence limit between 0.87 and 0.96, indicating that b was significantly <1, and that lower DRP concentrations occurred in conjunction with high discharge totals (figure 3b). Although this slope for sink soils was not significantly different from that of source soils (b = 0.93), the C-Q relationship for sink soils suggests solute dilution with high discharge, which has been attributed to source limited or reaction rate-limited conditions (Godsey et al. 2009; Bieroza et al. 2018; Duncan et al. 2017a, 2017b). In P sink soils, the magnitude of positive SPSC values is an estimate of legacy P held in stable forms in these soils (Nair et al. 2015). Since desorption rates are much slower than sorption rates (Menezes-Blackburn et al. 2016), the dilution behavior in sink soils is primarily controlled by this hysteretic behavior of P between the solid and liquid phase in soils.

At the event scale, a similar linear relationship between DRP flux from both P source and sink soils and discharge when log-transformed was observed (R² = 0.64 and R² = 0.78, respectively) (figure 4). However, the slopes (b) of the log-log relationships highlighted a contrasting DRP dynamic at the event scale when compared to the daily scale. The slope (b = 0.88) of the event DRP-discharge relationship for source soils had a 95% confidence limit between 0.75 and 0.99, indicating that b was significantly <1, i.e., a dilution effect (figure 4a). On the other hand, the slope (b = 0.94) of the event DRP-discharge relationship for sink soils had a 95% confidence limit between 0.86 and 1.01, indicating that b was not significantly different from 1; i.e., chemostatsis (figure 4b).

• Download figure
• Open in new tab
• Download powerpoint

Figure 4

Event dissolved reactive phosphorus (DRP) flux (kg ha^–1) to discharge (mm) relationship presented on a log-log scale for (a) phosphorus source soils and (b) phosphorus sink soils. Exponent (b) and coefficient (a) values are defined by equation 1, and R² values were determined by linear regression on log transformed dissolved reactive phosphorus flux and discharge values.

Previous authors (Ruark et al. 2009; Duncan et al. 2017a, 2017b) have also reported contrasting C-Q patterns when data are analyzed at different time scales, e.g., daily, monthly, seasonal, or event time scales. These studies showed that the variation of C-Q relationships across different time scales results from the different processes influencing C-Q behavior at different temporal scales. Although it is possible that the event-scale dilution (b < 1) observed in the P source soils could result from a single contribution of a high-source P pool that becomes increasingly exhausted throughout the event, the pervasive and persistent nature of legacy P in agricultural soils (Kleinman 2017) makes it unlikely that legacy stores become depleted on the timescale of an individual discharge event. On the other hand, the b < 1 may be the result of a dilution of the DRP-rich water from P-rich surface soils by P-poor groundwater. Subsurface discharge water is composed of both shallow groundwater and precipitation water that has infiltrated through the soil and moved downward either by matrix or preferential flow (Stamm et al. 2002). Previous studies that used stable water isotopes and tracers to monitor movement of water into tile drains showed that the initial contributions into tile discharge most likely originated from surface soils and moved through active preferential flow paths (Greve et al. 2012; Williams et al. 2016). The latter studies further showed that as the event progressed and the soil matrix neared saturation, tile discharge transitioned from predominantly preferential flow (event water) to a mixture of preferential flow, matrix flow, and shallow groundwater (Greve et al. 2012; Williams et al. 2016). Even though water movement and mixing was not monitored in this study, the occurrence of lower DRP concentrations in conjunction with high discharge events suggests DRP concentrations may have been diluted by matrix and shallow groundwater that are poor in P. This interpretation is supported by Williams et al. (2016), who observed that elevated DRP concentrations in tile discharge coincided with peak event water contribution from preferential flow pathways before the waters, i.e., preferential and matrix flow, mixed along the flow pathway. Additionally, the effect of water table rise and antecedent soil moisture (both not measured in this study) on hysteresis patterns was documented by Macrae et al. (2010) and Wagner et al. (2008) in their NO₃-N study where, instead of a solute dilution, they observed solute flushing (mobilization). With NO₃-N and unlike with P, which is not highly mobile and easily vertically distributed in soils (Baker et al. 2017), rising water tables mobilized NO₃-N that is well-distributed in the soil matrix and is also found in high concentrations in shallow groundwater. In contrast, the chemostatic behavior (b = 1) observed at the event scale in sink soils suggests that these surface soils are just as P poor as matrix and shallow groundwater. Since the stable forms of P in P sink soils are not readily desorbed (Nair et al. 2015), water leaching through a P sink surface soil may have a similar and low P concentration as that of matrix and groundwater. Based on these results, the proposed conceptual models illustrating DRP loss in P source and P sink soils are presented in figure 5.

• Download figure
• Open in new tab
• Download powerpoint

Figure 5

(a) A proposed conceptual model illustrating dissolved reactive phosphorus (DRP) loss in tile drained fields from phosphorus (P) source soils where (i) discharge at the start of the event when P rich preferential flow (event water) from the surface rapidly flows to tile drains, (ii) discharge during the event when P rich event water is mixed with P poor water from the soil matrix and shallow groundwater, resulting in an overall dilution, and (iii) when the discharge recedes, contributions from soil matrix and shallow groundwater are absent, event water with high P concentrations continues. (b) A proposed conceptual model illustrating DRP loss in tile drained fields from P sink soils where (i) discharge at the start of the event when P poor preferential flow (event water) from the surface rapidly flows to tile drains, (ii) discharge during the event when P poor waters (preferential, matrix, and shallow groundwater) mix, and (iii) when soil matrix and shallow groundwater recede, event water with low P concentrations continues.

Our results from both timescales highlight the importance of P-rich surface soils and hydrological connectivity for DRP loss to tile drains. The continuous release of DRP from P source soils highlights how legacy P continues to undermine P mitigation efforts. Unfortunately, phytomining (growing and harvesting crops without P fertilization) takes decades or longer to draw down P (McCollum 1991; Schärer et al. 2007; Sharpley and Rekolainen 1997). On the other hand, there are a few timelier and cost-effective materials arising as waste streams from various processes, which can be used to sequester P from P-rich soils. These materials include fluidized gas desulfurization (FGD) gypsum, red muds, crushed concrete, Fe gels, etc. (Kleinman et al. 2019; Murphy and Stevens 2010; Egemose et al. 2012; Chardon et al. 2012; Weng et al. 2012). Among these amendments, FGD gypsum has gained traction due to its added benefit of acting as a source of sulfur (S), and because it can be used for the treatment of sea-salt impacted soils (Murphy and Stevens 2010). However, cost and possible adverse effects from the use of the remaining P amendments has hindered their adoption (Kleinman et al. 2019). Also, given macropores are generally believed to be the major pathway through which DRP is transported to tile drains, some studies recommend tillage to disrupt macropore flow pathways and decrease hydrological connectivity (Williams et al. 2016). Although tillage could potentially reduce DRP losses through tile discharge, it could also increase particulate P losses (Verbree et al. 2010), therefore it may not be a suitable conservation practice.

Hysteresis Patterns. Pearson’s correlation analysis (combined data for P source and P sink soils) of discharge characteristics descriptors of C-Q hysteresis identified significant relationships among ΔR, Q, and D_rel (table 5). ΔR was positively correlated to Q (p < 0.05) and D_rel (p < 0.01), suggesting that an increase in Q and D_rel would lead to a shift in loop trajectories from anticlockwise to clockwise. This observation has been reported in previous studies in which DRP anticlockwise loops were mostly associated with low discharge events (Williams et al. 2018; Bowes et al. 2005; Chow et al. 2017; Bieroza and Heathwaite 2015). The ΔC_new was not significantly correlated with any of the discharge event characteristics. Figure 6 shows a plot of ΔC_new (%) and ΔR (%) that summarizes C-Q hysteresis loop types of the DRP from sink and source soils during the selected discharge events. Potential nutrient sources and delivery pathways to surface water are usually inferred from C-Q hysteresis loops where clockwise hysteresis infers proximal sources with exhaustible solute supply or intense discharge, and anticlockwise hysteresis infers distal sources or transport-limited systems (Bowes et al. 2015). However, as noted by Williams et al. (2018), interpretation of hysteresis loops in tile-drained systems is more complex as even distant sources are directly connected to streams by tile drains.

• Download figure
• Open in new tab
• Download powerpoint

Figure 6

Plot of relative change in solute concentrations (ΔC_new [%]) versus change in rotational patterns (ΔR [%]) for the C-Q hysteresis loops of dissolved reactive phosphorus (DRP). The i, j, and k terms in the plot labels correspond to ith season (winter [W], spring [Sp], summer [Su], and fall [F]) in a water year (2011, 2012, and 2013), the jth plot at the Water Quality Field Station (P10, P11, P12, P26, P30, P32, P43, and P44), and the kth discharge event for the specified plot. Supplementary table S1 provides detailed information on discharge events. Illustrations of the typical C-Q relationships (C, dashed blue line; Q, continuous brown line) are presented for each of the regions A through D of the ΔC_new (%) versus ΔR plot. A few events in phosphorus source soils showed clockwise hysteresis.

The use of low resolution data (daily sampling) did not accurately capture the variation in DRP concentrations (ΔC_new [%]). For example, for some events, initial Q was also Q_max in a discharge event, therefore, C_i = C_s, resulting in discharge events whose solute trend (flushing or dilution) could not be determined given they plotted on the y = 0 line (figure 6). These discharge events included 42% (30% anticlockwise, 12% clockwise) and 25% (18% anticklockwise, 7% clockwise) of events in P source and sink soils, respectively.

Of the remaining discharge events, our results show that hysteresis behavior of events in P source soils was variable with all of the following behaviors exhibited: anticlockwise with dilution (21%), anticlockwise with flushing (7%), clockwise with dilution (6%), clockwise with flushing (9%), and no hysteresis behavior (15%), respectively. According to Williams et al. (2018) this variability in hysteresis behavior suggests that multiple flow pathways and transport mechanisms are involved in DRP loss to tile drains. The predominant anticlockwise and dilution pattern of events in P source soils (region C of figure 6) has two implications. First, tile drains facilitate the rapid transport of DRP from P rich surface soils and the contribution of event water to tile discharge (Vidon and Cuadra 2011; Williams et al. 2016). Second, the mixing between preferential water, matrix water, and shallow groundwater may delay the peak in DRP concentration; thus a dilution of high DRP concentrations when waters mix results in anticlockwise hysteresis (high DRP concetrations on the falling limb versus the rising limb of the hydrograph) (Williams et al. 2018). The majority dilution pattern of events in P source soils is consistent with our results (b < 1), from the C-Q slope analysis at the event scale (b < 1; figure 4a). It reinforces our interpretation of the slopes of the C-Q relationships where the initial, high DRP concentrations in event water progressively mixes with lower DRP concentrations in shallow groundwater, resulting in a dilution as discharge approaches its peak. The anticlockwise hysteresis (illustration in region C of figure 6) further shows that the dilution pattern reverses as discharge recedes, with soil matrix water and shallow groundwater progressively becoming disconnected from the tile drains, whereas event water from P-rich surface soils remains hydrologically connected through preferential pathways. In contrast, during other events, insufficient rise of shallow groundwater around tile-drains may prevent the mixing of P-rich and P-poor waters. This minimal mixing coupled with rapid transport of DRP and the contribution of event water to tile discharge, may result in clockwise hysteresis (high DRP concetrations on the rising limb versus the falling limb of the hydrograph) (Williams et al. 2018, 2016; Vidon and Cuadra 2011). Also, the positive correlations of ΔR (%) with Q and D_rel (table 5) show that the shift from anticlockwise to clockwise hysteresis may be the result of an increase in Q and D_rel. Finally, the variation in DRP response could be the result of rapid exchanges between P pools due to sorption/desorption and biological processes under varying antecedent conditions (Williams et al. 2018). This variable hysteresis behavior and solute trend underscores the challenge of interpreting the contributing mechanisms to better manage DRP loss from P source soils in tile drained systems. However, we note that our results may be an artifact of low temporal resolution sampling. Events in P sink soils seemed to generally have no hysteresis or solute trend, i.e., 67% of the events were plotted at the origin. These results support the chemostatic (b = 1) finding in the event C-Q slope analyses of sink soils. Among the remaining events whose solute trends were identified, 2% and 6% were anticlockwise with flushing and clockwise with dilution, respectively.

Even though tile drains are a major pathway for DRP loss as shown in this study and others (Welikhe et al. 2020a; Gentry et al. 2007; Macrae et al. 2007; Gelbrecht et al. 2005), our results suggest that the water table dynamics during discharge events, soil P status, soil P stratification, the amount of discharge, and the number of days it takes to reach peak discharge, also control nutrient delivery and hysteresis patterns. Also, there are other factors not accounted for in our study that may influence DRP-Q dynamics and are worthy of future study. For example, detailed assessments of the influence of crop residues in no-till systems on DRP concentrations in tile effluent are needed. Numerous studies have reported that DRP loss can be greater in runoff (surface and subsurface) from no-till or conservation till systems than from conventional till systems (Bundy et al. 2001; Daverede et al. 2003; Tiessen et al. 2010; Liu et al. 2014). Earlier studies attributed the increased DRP losses to nutrient stratification at the soil surface in conservation and no-till systems compared to conventional systems (Bundy et al. 2001; Daverede et al. 2003). However, recent literature also suggests that the crop residues remaining on the soil surface in conservation and no-till systems release P when exposed to multiple freeze-thaw cycles in cold climates from late fall to early spring, thus acting as an additional source of DRP (Liu et al. 2019, 2014; Lozier et al. 2017).