Effects of Riparian Grazing on Distinct Phosphorus Sources

Core ideas

• Biomass and litter are substantial sources of WEP in riparian areas
• Autumn cattle grazing and mowing treatments reduced the areal density of WEP in riparian biomass
• There were no measurable changes in the areal density/concentration of WEP in the litter, organic layer, or Ah horizon post grazing
• Large spatial variability in areal density/concentration of WEP exists in riparian areas

Abbreviations

FTC, freeze-thaw cycle; MBFI, Manitoba Beef and Forage Initiatives; P, phosphorus; WEP, water extractable phosphorus

1 Introduction

The increasing frequency and extent of algal blooms is typically linked to increased nutrient loading into lake and rivers. Phosphorus (P) loading is particularly concerning as this is generally the limiting nutrient in freshwater systems (Schindler et al., 2012). There have been many lab and field studies demonstrating the role and functionality of riparian areas in reducing P loading to surface water in agricultural settings (Yu et al., 2019). Infiltration, absorption, biological uptake, microbial activity, and sedimentation are the key processes that intercept and buffer the delivery of P (Lacas et al., 2005; Owens et al., 2007; McGuire and McDonnell, 2010).Convergence within the landscape coupled with climatic/weather conditions creates variability in hydrologic conditions and pathways, reducing the buffering capacity of riparian areas and ultimately resulting in reduced, inconsistent, and/or unsustainable reductions in P loading relative to many controlled experimental studies (Roberts et al., 2012; Habibiandehkordi et al., 2017).

In cold climates, the reduced infiltration due to frozen ground, limited vegetation uptake, and low microbial activity coupled with a flashy hydrograph during snowmelt creates conditions that further compromise the buffering capacity of riparian areas (Kieta et al., 2018; Nsenga Kumwimba et al., 2023). Additionally, research increasingly shows that riparian areas can contribute P (i.e., net source) from soil and vegetation to the surrounding environment (Roberts et al., 2012). As soil P concentration increases, so does the risk of P loss through leaching and runoff (Habibiandehkordi et al., 2019). Soil P release can be intensified during periods of inundation that often occur during the spring snow melt, due to both to a longer period of soil-water contact and an increased solubility of iron-bound P as soil redox conditions lower (i.e., become anaerobic) (Carlyle and Hill, 2001; Young and Briggs, 2008). Vegetation P can become more mobile through the mineralization of P from decaying vegetation near the soil surface. There is also evidence that the longer vegetation-water contact during periods of inundation will also increase the mass of P leached out of the dead vegetation and contribute to the P available to be lost during runoff (Lozier and Macrae, 2017; Liu et al., 2019b). Both the soil and vegetation P sources can also be affected by freeze-thaw cycles (FTC). Repeated FTCs result in the cell disruption of microbial and plant biomass, releasing inter-cellular P to the surrounding environment (Kieta and Owens, 2019).

Management of riparian areas to maintain or enhance the buffering capacity of P is typically needed in the long term. Unlike nitrogen (N) where N can be significantly lost to the atmosphere through nitrification and denitrification to offset the continued input (Lyu et al., 2021), P is generally only lost through runoff or leaching. Harvesting and removing of biomass from the riparian area for use as forage can be a practice to remove P. Mechanized biomass harvesting may be impractical or unsafe due to steep gradients, wet soil, and other obstacles like trees; however, livestock grazing in riparian areas (riparian pastures) is common in the Canadian Prairies due to the abundance of forage, particularly during drought. Livestock exclusion from riparian areas has been suggested as a best management practice to reduce the direct inputs of P, limit bank erosion, and avoid soil compaction (Krall and Roni, 2023). However, strategies including alternative water sources, rotational grazing, timed-controlled grazing, rest-rotation grazing, and corridor fencing can all reduce those risks (Fitch et al., 2003).

From a surface water quality perspective, understanding the near-surface P distribution, both vertically and longitudinally, will help develop and identify best management practices for reducing P loading from riparian areas. Vertically, there are often four distinctive and identifiable sources of near-surface P: 1) biomass consisting of living standing vegetation; 2) litter consisting of fresh (within the first three years) residues; 3) partially to well-decomposed organic material; and 4) mineral soil (Reid et al., 2018). Longitudinally there often is a strong soil moisture gradient extending from the edge of the waterbody to the field edge. This results in changes in the mass and composition of biomass and litter as well as soil properties including organic matter content and horizon thickness. A better understanding of the spatial variability and relative contributions of the different sources of P is needed to assess the risks and benefits of different management strategies.

Given the timing and processes of P dynamics within riparian areas in cold climates, like the Canadian Prairies, reducing the near-surface concentration of soluble P prior to spring snowmelt could be a strategy to limit the contribution of P from the riparian area to surface water. Therefore, the overall aim of this study is to assess the impacts of short-term autumn cattle grazing and mowing on the sources and distribution of P in riparian areas. The objectives of this study were to quantify 1) the vertical profile of WEP using four distinctive P sources: biomass, litter, organic layer, and Ah horizon; 2) each of the four distinctive P sources in three riparian locations, near the edge of the waterbody (lower), close to the field edge (upper), and in between (middle); and 3) the net change in each of the four sources of WEP in each riparian location in response to grazing, high-density grazing, and mowing (harvesting) of biomass. Understanding how riparian management practices affect the different sources of P can be used to help tailor management strategies in cold climates and ultimately reduce P loss and improve downstream water quality.

2 Methods

2.1 Site description

In [1]:

station <- weathercan::stations_search("Brandon", normals_years = "current")

The stations data frame hasn't been updated in over 4 weeks. Consider running `stations_dl()` to check for updates and make sure you have the most recent stations list available
The stations data frame hasn't been updated in over 4 weeks. Consider running `stations_dl()` to check for updates and make sure you have the most recent stations list available

climate_normal <- weathercan::normals_dl(climate_ids = 5010480, normals_years = "1981-2010") |>
  tidyr::unnest(normals) |>
  dplyr::select(period, temp_daily_average, precip, snow)

The stations data frame hasn't been updated in over 4 weeks. Consider running `stations_dl()` to check for updates and make sure you have the most recent stations list available

A randomized complete block experimental design was used to assess the sources of riparian P and investigate how it changes following cattle grazing or mowing treatments. The four treatments include control, graze, high-density graze, and mowing. Each treatment was replicated in riparian areas surrounding four prairie potholes (wetlands). Samples of biomass, litter, organic layer, and Ah horizon, were collected in three locations both pre- and post-treatment. The three sampling locations aimed to capture the topography of the riparian areas and include near the edge of the waterbody (lower), close to the field edge (upper), and the mid-point (middle). All samples were analyzed for WEP and the net change in each of the four distinctive sources of P following the treatment was evaluated. The study was replicated across three sequential years using the same plots. A workflow diagram showing the experimental setup, field work, sample preparation, and laboratory analysis can be found in Figure S1.

The study was conducted at the Manitoba Beef and Forage Initiatives (MBFI) research farm (50.06\(^\circ\)N, 99.92\(^\circ\)W; 502 AMSL), approximately 25 km north of Brandon, Manitoba, Canada, in the Prairie Pothole region of North America (Figure 1). The normal (1981 – 2010) average daily air temperature was 2.2 \(^\circ\)C, and the cumulative annual precipitation at Brandon was 474.2 mm, with 24.8 % falling as snow (Environment and Climate Change Canada, 2024). The Köppen-Geiger climate classification is cold, without dry season, and with warm summer (Dfb) (Beck et al., 2018). The region is predominantly agricultural land use, including annual crops (grains and oil seeds) and grazing/forage. MBFI is a 260-hectare (ha) research and demonstration farm with a mix of pasture, hay, and forage/silage cropland. Prior to the establishment of MBFI the site was part of the Manitoba Zero Tillage Research Association farm (1993-2014) where annual crops, including oil seeds and grains, were grown. There are also numerous small permanent and ephemeral wetlands (potholes) and associated riparian areas which account for approximately 35% of the total farm land (Manitoba Beef & Forage Initiatives, 2024). The riparian areas surrounding the larger permanent wetlands are fenced off to exclude livestock and are not actively managed. Approximately half the farm has an irregular undulating to hummocky relief (2-5%) with the reminder being nearly level (0-2%). The soils have developed on fine loamy, moderately calcareous glacial till. The drainage class in upper slope positions are well to rapidly draining while lower slope and riparian soils are poorly drained and primarily consist of Humic and Luvic Gleysols. The surface texture class of the riparian soil is a clay loam and pH values range from 7.1 to 8.3 with a mean of 7.6. Generally the surface soil profile can be described by a 1-10 cm organic layer overlying a 10-18 cm Ah horizon (Podolsky and Schindler, 1993). Vegetation was assessed using the foliar cover method for each plot within each of the four riparian areas. There was considerable variability among riparian areas, plots, and sampling locations (upper, midle, and lower). The four most dominant species identified were Sow Thistle (Sonchus arvensis), Smooth Aster (Aster laevis), Kentucky bluegrass (Poa pratensis), and Smooth Brome (Bromus inermis) and the complete assessment can be found in Figure S2. All riparian areas investigated in this study were adjacent to actively grazed pastures.

In [2]:

p6 <- p4 + inset_element(p1, left = 0.6, bottom = 0.6, right = 1, top = 1) + p5
p6

Warning in st_point_on_surface.sfc(sf::st_zm(x)): st_point_on_surface may not
give correct results for longitude/latitude data

Figure 1: Showing a) the location of the study site in southern Manitoba with an inset map of Canada; and b) the locations of the four riparian areas included in this study

2.2 Experimental design

Four riparian areas surrounding permanent wetlands were selected (Figure 1) and subdivided into four approximately 450 \(m^2\) plots. Within each riparian area, each plot was randomly assigned a treatment. The treatments were 1) control, 2) graze, 3) high-density graze, and 4) mow and harvest. The grazing treatments consisted of a five-hour grazing period, with the grazing treatment having 3.1-3.5 animal units per plot and the high-density grazing having 11.75-12 animal units. For the mowing treatment, the vegetation was cut to a height of 10cm, and the vegetation was manually raked out of the plot. The grazed plots were fenced on all four sides, including the edge of the waterbody. The cattle were rotated among the four riparian areas daily over four consecutive days and provided with supplemental water. Treatments were applied early to mid-September, before the first frost, in three consecutive years (2019-2021) (Figure S3). Within each plot three distinctive sampling locations, or topographic positions, were established, adjacent to the edge of the waterbody (Lower), adjacent to the field/pasture (Upper), and at the mid-point (Middle). Samples were collected at each sampling location 1-3 days before and 1-3 days after treatment (including the control) to assess the impact of grazing and mowing. Before and after samples were collected at immediately adjacent locations.

2.3 Sampling and analysis

Four types of samples were collected: 1) biomass, 2) litter, 3) organic layer, and 4) Ah horizon. Using a 0.25 \(m^2\) quadrate, biomass was collected by cutting the standing live vegetation and litter by raking the surface and picking up the previous year’s growth. Both the biomass and litter were dried at 40 \(^\circ C\), weighed, and homogenized using a blade grinder (<1cm). A composite of five soil samples was collected within the same quadrat as the biomass/litter using a 19 mm diameter soil probe and was divided into the organic layer and the top 10 cm of the Ah horizon. The organic layer and Ah soil were air-dried, disaggregated with a mortar and pestle, and passed through a 2-mm sieve. Additional bulk density samples of both the organic layer and Ah and the depth of the organic layer were collected in 2023. Daily air temperature and rainfall data were collected from an onsite station (Figure S3) (Manitoba Agriculture, 2023).

Water Extractable Phosphorus (WEP), an environmental soil and vegetation P test, was used to mimic soil P release into runoff water. Dried and homogenized samples were extracted by shaking (150 RPM) with deionized water for one hour at a mass-to-volume ratio of 1:30 for the biomass and litter samples (1 g) and 1:15 for the organic and Ah samples (2 g). Extractions were gravity filtered through a Whatman 42 filter followed by syringe filtration with a 0.45 \(\mu m\) nylon filter. WEP in the extract was measured spectrophotometrically by the colorimetric molybdate–ascorbic acid method (Murphy and Riley, 1962; Sharpley et al., 2006).

The concentration of WEP (\(mg~kg^{-1}\)) was calculated for all sources of P. In addition, the areal density of WEP was calculated for biomass and litter by combining WEP concentration with the mass of material collected from the quadrat. The vertical profile of WEP within the riparian area assessed from samples collected before treatments were implemented across the 3-year study. For comparison, a rough estimation of areal density WEP in the organic layer and Ah was calculated using the bulk density and depth measurements collected in 2023 (Figure 2 b).

2.4 Statistical analysis

All statistical analysis, plotting, and mapping were undertaken using the R Statistical Software (v4.4.0; R Core Team (2024)), through the RStudio Integrated Development Environment v2023.12.1.402 (RStudio, 2024). All plots and maps were created using the R package ggplot2 (v3.5.1; Wickham (2016)). Country and regional maps were created using data from the rnaturalearth package (Massicotte and South, 2023) and other maps using ESRI imagery and the OpenStreetMap package (Fellows, 2023). Four Linear Mixed Models (R package glmmTMB v1.1.9; Brooks et al. (2017)) were used to investigate the effect of treatment and riparian sampling location (including interaction) on the change in WEP (before — after treatment) of for each of the four distinct sources of P (areal densities for biomass and litter; concentrations for organic matter and Ah). Year and riparian area were included as crossed random factors to control for the variability within years and riparian areas.

Additionally, when investigating the net change in biomass WEP the initial biomass WEP (before applying the treatment) was included in the model as a covariate. This controls for the fact that the magnitude of change in biomass WEP (i.e., before - after) is directly related to the mass of WEP initially available. By controlling for this, the model indicates which treatments resulted in a relatively greater change in biomass, rather than simply absolute change.

The interaction between treatment and riparian sampling location was removed if non-significant (p < 0.05). When a main effect or interaction was significant, post-hoc pairwise comparisons with a Benjamini-Hochberg p-value adjustment were performed (p <0.05). When a main effect or interaction was significant, post-hoc pairwise comparisons with a Benjamini-Hochberg p-value adjustment was used (emmeans v1.10.1; Lenth (2024)). Model assumptions were assessed using DHARMa residual plots (DHARMa v0.4.6; Hartig (2022)), main effects were tested for collinearity (performance v0.12.2; Lüdecke et al. (2021)), and results were presented as type III ANOVA (car v3.1.2; Fox and Weisberg (2019)). For each unique source of WEP, the null hypotheses were: no difference in the net WEP among treatments or riparian sampling locations and no interactions between these two factors.

Pearson correlations were performed to explore relations in WEP concentrations among the four unique P sources for each of the three topographic positions using samples collected before the application of the treatments. These relations were visualized using a scatterplot matrix created using the GGally R package (v2.2.1; Schloerke et al. (2024) )

3 Results and Discussion

3.1 Vertical profiles of P

The biomass, litter, organic layer, and Ah horizon sources of P demonstrated a strong vertical stratification in both the concentration and areal densities of WEP (Figure 2). The median concentrations in the vegetation sources were 82.8 and 39.0 \(mg~kg^{-1}\) for the biomass and litter components, respectively, which is more than an order of magnitude greater than the soil components (0.9 and 3.4 \(mg~kg^{-1}\); Ah and organic, respectively). Considerable variability in the WEP concentration in the biomass and litter sources were observed with interquartile ranges (IQR) of 54.3 and 32.9 \(mg~kg^{-1}\) for the biomass and litter sources, respectively. In contrast, the IQR for the organic and Ah sources both were <2.5 \(mg~kg^{-1}\). Overall, in terms of the areal density of WEP, the top 10 cm of the Ah horizon was the largest source of WEP (42.5 \(mg~m^{-2}\)) followed by the biomass (26.3 \(mg~m^{-2}\)), organic layer (14.3 \(mg~m^{-2}\)), and lastly the litter (13.7 \(mg~m^{-2}\)). Although it should be noted that these are only rough estimates for the organic layer and Ah horizon. Nevertheless, the vertical profile of WEP in riparian areas (Figure 2) observed in this study supports the concept that a measure of P in soil alone is likely missing a large proportion of the near-surface P that can be potentially lost during the spring snowmelt (Liu et al., 2019a; b; Cober et al., 2019). The substantial proportion of WEP above the soil surface provides evidence that managing the biomass in riparian areas in autumn may reduce the contribution of P lost directly from this area during spring. Specifically, the harvesting of this biomass results in an export of P which can maintain or enhance the buffering or storage capacity of P derived from upslope sources further improving downstream water quality (Kelly et al., 2007; Hille et al., 2019).

In [3]:

p3 <- p1+p2 + plot_layout(guides = 'collect') & theme(legend.position = 'bottom', legend.title = element_blank())  

p3

Warning: Removed 2 rows containing non-finite outside the scale range
(`stat_boxplot()`).
Removed 2 rows containing non-finite outside the scale range
(`stat_boxplot()`).

Figure 2: Vertical and longitudinal profiles of a) WEP concentration and b) WEP content in the riparian areas prior to grazing and mowing treatments.

3.2 Longitudinal profiles of P

Prior to grazing and mowing treatments, the median WEP concentrations were similar among the upper, mid, and lower positions for the biomass samples. There was a small topographic trend in the WEP concentration for both the Ah and organic litter P sources where the concentration decreased from the upper through to the lower sampling locations. The WEP concentrations in the Ah and organic layer were found to be significantly (p < 0.001) and positively correlated (r² = 0.40) (Figure 2 and Figure S5). This topographic pattern is consistent with other studies and is likely due to the rapid physical and geochemical retention of upslope derived P within the first 5 m of the riparian area (Syversen and Borch, 2005).The litter showed the opposite topographic trend with higher WEP concentrations in the lower sampling locations. There was a significant (p < 0.001) positive correlation (r² = 0.34)between the WEP concentration in the biomass and litter samples suggesting that biomass with a high WEP concentration produces litter with a high WEP concentration (Figure S5). There was no correlation (p > 0.05) between the Ah and biomass WEP concentrations suggesting that higher soil WEP concentration does not result in biomass with elevated WEP concentrations at this study site (Figure S5). There was no correlation (p > 0.05) between the Ah and biomass WEP concentrations suggesting that higher soil WEP concentration does not result in biomass with elevated WEP concentrations at this study site (Figure S5). The variability is greatest in the Ah (IQR = 32.0 \(mg~kg^{-1}\)) and biomass (IQR = 23.3 \(mg~kg^{-1}\)) sources. The variability of the other two sources were similar with IQRs of 15.6 and 14.3 \(mg~kg^{-1}\) for the litter and organic layer, respectively. Although there is some evidence that plants in P-rich environments will also be enriched in P (e.g., Kröger et al., 2007)). For the biomass and litter sources the lower riparian locations had greater areal densities of WEP whereas the organic and Ah sources had greater areal densities of WEP in the upper riparian locations. The longitudinal gradient of WEP showed an inverted symmetry where the biomass WEP was largest near the lower sampling location and the Ah soil WEP was largest in the upper sampling location adjacent to the fields (Figure 2 b). The high soil water content in the lower location created conditions that favor high biomass production (Figure S4) ) and higher WEP concentration (Figure 2 a). The higher bulk density was most likely due to the lower soil organic matter content and the higher WEP concentration may be related to the interception of P-rich runoff from upslope areas (Tomer et al., 2007). Understanding and quantifying the sources and patterns of P within riparian areas is a key part of assessing the risk of P loss as it helps to inform management decisions and target the largest sources of P (Reid et al., 2018).

3.3 Impacts of grazing and mowing on P sources

There was considerable variation across all treatments and riparian locations in all four P sources. This high variability in WEP areal density/concentration is best reflected in the control treatment where the expected difference was 0 (Figures 3-6), but WEP losses and gains were still observed despite no treatment being applied. However, despite this variability, several patterns demonstrating relationships among treatments and vertical and longitudinal P emerged.

Results of the linear mixed model of areal density of biomass WEP show a significant effect of treatment (X² = 24.8, df = 3, p < 0.001) and riparian location (X² = 15.7, df = 2, p < 0.001). Post-hoc comparisons showed that the net biomass WEP for the high-density grazing and mowing treatments were similar (p>0.05) but significantly (p<0.05) different from the control and graze treatments (Figure 3 a and Table 1).The mowing and high-density grazing reduced the average WEP areal density by 7.4 and 4.2 \(mg~m^{-2}\) relative to the control, respectively. The reduction in biomass WEP was significantly (p<0.05) greater in the lower sampling locations as compared to the upper and mid locations (Figure 3 b and Table 1) with a difference in average WEP of 10.2 \(mg~m^{-2}\) between the lower and upper locations of the riparian area.

In [4]:

p3 <- p1 + p2 + plot_layout(guides = 'collect') & theme(legend.position = 'bottom', legend.title = element_blank())
p3

Warning: Removed 2 rows containing non-finite outside the scale range
(`stat_boxplot()`).
Removed 2 rows containing non-finite outside the scale range
(`stat_boxplot()`).

Figure 3: Change in riparian biomass WEP following grazing or mowing in each riparian location. Within each plot significant differences (p<0.05) between treatments or riparian locations are denoted with different letters. Lower sampling locations are adjacent to the edge of the waterbody and Upper locations are adjacent to the field.

In [5]:

In [6]:

pairs |>
  mutate(Contrast = fct_recode(Contrast, "Control - Graze" = "Control - Regular Graze", "High Graze - Graze" = "High Graze - Regular Graze", "Mow - Graze" = "Mow - Regular Graze")) |>
  group_by(type) |>
  gt() |>
  fmt_number(columns = c("Estimate", "SE", "t ratio"), decimal = 2)|>
  fmt_number(columns = c("p value"), decimal = 3)|>
  sub_small_vals(threshold = 0.001) |>
  tab_style(style =  cell_text(weight = "bold", align = "center"), locations =  cells_row_groups()) |>
  tab_options(column_labels.font.weight = "bold")

Table 1: Results of the post-hoc pairwise comparisons with a Benjamini-Hochberg p value adjustment for differences in the net biomass WEP (\(mg~m^{-2}\)) between the four treatments and three riparian sampling locations.

Contrast	Estimate	SE	df	t ratio	p value
Treatment
Control - High Graze	−4.83	2.42	132	−2.00	0.072
Control - Mow	−8.52	2.42	132	−3.52	0.002
Control - Graze	2.47	2.40	132	1.03	0.306
High Graze - Mow	−3.69	2.43	132	−1.51	0.159
High Graze - Graze	7.30	2.42	132	3.02	0.006
Mow - Graze	10.99	2.42	132	4.55	<0.001
Location
Lower - Middle	−7.94	2.43	132	−3.26	0.002
Lower - Upper	−9.82	2.57	132	−3.83	<0.001
Middle - Upper	−1.87	2.11	132	−0.89	0.377

The model looking at areal density of litter WEP showed no significant impacts of either treatment (X² = 1.15, df = 3, p = 0.23) or riparian location (X² = 4.30, df = 2, p = 0.56) (Figure 4). ). In contrast, the model exploring WEP concentration in the organic layer detected no significant difference among riparian locations (X² = 0.57, df = 2, p = 0.75) but did find a significant effect of treatment (X² = 8.24, df = 3, p = 0.04). However, the post-hoc pairwise comparisons (Table 2) found no significant differences (p <0.05) among the treatments. Finally, there was no significant effect of treatment (X² = 2.59, df = 3, p = 0.46) or riparian position (X² = 1.17, df = 2, p = 0.56) on the concentration of WEP in the Ah horizon (Figure 6).

In [7]:

p1 <- ggplot(data = filter(plot_data, measure == "p_total")) +
  theme_bw(base_size = 12) +
  geom_rect(data = df, aes(xmin = x1, xmax = x2, ymin = y1, ymax = y2, fill = difference), alpha = 0.15) +
  #scale_fill_manual(values = c("white", "black")) +
  #ggnewscale::new_scale_fill() +
  geom_boxplot(aes(x = treatment, y = diff, fill = location)) +
  labs(y = expression(paste("Net WEP Difference (", mg~m^{-2}, ")")), x = "Treatment") +
  #scale_fill_viridis_d(name = "Location", begin = 0.3, end = 1) +
  scale_fill_manual(values = c("white", "black", "#35608DFF", "#2FB47CFF", "#FDE725FF")) +
  guides(fill = guide_legend(override.aes = list(colour = "black", size = 1))) +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 0, vjust = 1, hjust = 0.5),
        legend.position = "bottom",
        legend.title = element_blank())
p1

Warning: Removed 11 rows containing non-finite outside the scale range
(`stat_boxplot()`).

Figure 4: Change in riparian litter WEP following grazing or mowing in each of the riparian locations. No significant effect of treatment or riparian location on the litter WEP content was detected. Lower sampling locations are adjacent to the edge of the waterbody and Upper locations are adjacent to the field.

In [8]:

p1 <- ggplot(data = organic_diff) +
  theme_bw(base_size = 12) + 
  geom_rect(data = df, aes(xmin = x1, xmax = x2, ymin = y1, ymax = y2, fill = difference), alpha = 0.15) +
  #scale_fill_manual(values = c("white", "black")) +
  #ggnewscale::new_scale_fill() +
  geom_boxplot(aes(x = treatment, y = diff, fill = location)) +
  #scale_fill_viridis_d(name = "Location", begin = 0.3, end = 1) +
  scale_fill_manual(values = c("white", "black", "#35608DFF", "#2FB47CFF", "#FDE725FF")) +
  guides(fill = guide_legend(override.aes = list(colour = "black", size = 1))) +
  labs(y = expression(paste("Net WEP Difference (", mg~kg^{-1}, ")")), x = "Treatment") +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 0, vjust = 1, hjust = 0.5),
        legend.position = 'bottom', legend.title = element_blank())
p1
#ggsave(plot = p1, filename = "Figures/Organics_WEP.png", width = 150, height = 100, units = "mm", dpi = 600)

Figure 5: Change in riparian organic layer WEP concentration following grazing or mowing in each of the riparian locations. A significant effect of treatment was detected; however, the post-hoc analysis was not able to detect any significant (p < 0.05) pairwise contrasts. Lower sampling locations are adjacent to the edge of the waterbody and Upper locations are adjacent to the field.

In [9]:

In [10]:

pairs |>
  gt() |>
  fmt_number(columns = c("Estimate", "SE", "t ratio"), decimal = 2)|>
  fmt_number(columns = c("p value"), decimal = 3)|>
  sub_small_vals(threshold = 0.001) |>
  tab_style(style =  cell_text(weight = "bold", align = "center"), locations =  cells_row_groups()) |>
  tab_options(column_labels.font.weight = "bold")

Table 2: Results of the post-hoc pairwise comparisons with a Benjamini-Hochberg p value adjustment for differences in the net organic layer WEP (\(mg~kg^{-1}\)) between the four treatments.

Contrast	Estimate	SE	df	t ratio	p value
Control - Graze	−1.49	0.59	135	−2.50	0.066
Control - High Graze	−0.63	0.59	135	−1.05	0.353
Control - Mow	−1.38	0.59	135	−2.32	0.066
Graze - High Graze	0.86	0.59	135	1.45	0.299
Graze - Mow	0.11	0.59	135	0.18	0.856
High Graze - Mow	−0.75	0.59	135	−1.27	0.311

In [11]:

p2 <- ggplot(data = soil_diff) +
  theme_bw(base_size = 12) + 
  geom_rect(data = df, aes(xmin = x1, xmax = x2, ymin = y1, ymax = y2, fill = difference), alpha = 0.15) +
  #scale_fill_manual(values = c("white", "black")) +
  #ggnewscale::new_scale_fill() +
  geom_boxplot(aes(x = treatment, y = diff, fill = location)) +
  #scale_fill_viridis_d(name = "Location", begin = 0.3, end = 1) +
  scale_fill_manual(values = c("white", "black", "#35608DFF", "#2FB47CFF", "#FDE725FF")) +
  guides(fill = guide_legend(override.aes = list(colour = "black", size = 1))) +
  labs(y = expression(paste("Net WEP Difference (", mg~kg^{-1}, ")")), x = "Treatment") +
  theme(axis.title.x = element_blank(),
        axis.text.x = element_text(angle = 0, vjust = 1, hjust = 0.5),
        legend.position = 'bottom', legend.title = element_blank())
p2
#ggsave(plot = p2, filename = "Figures/Soil_WEP.png", width = 150, height = 100, units = "mm", dpi = 600)

Figure 6: Change in riparian Ah layer (0-10cm) WEP concentration following grazing or mowing in each of the riparian locations. No significant effect of treatment or location was detected. Lower sampling locations are adjacent to the edge of the waterbody and Upper locations are adjacent to the field.

Taken together, these results suggest that short-term autumn high-density grazing may be a potential management tool that can reduce the mass of P lost directly from the riparian area (Figure 3 a). In addition to managing P loss, grazing riparian areas can also provide an essential source of forage, particularly during drought. Mechanized harvesting of biomass could also achieve this reduction in P loss (Figure 3 a) if the landscape and soil conditions are favorable. Despite the cycling of nutrients by the removal of P through grazing of biomass (Figure 3) and the deposition through excretion, no differences were detected in the litter and Ah sources of P (Figure 4, and 6). The models did detect a significant effect of treatment on the organic layer WEP; however, the pairwise comparisons were not able to detect any significant differences and the exact nature of the impact of the treatments remains unclear. The ability to detect changes in the WEP sources in riparian areas is difficult due to spatial variability in both the pre- and post-grazing treatments. Even within the control plots, both net addition and removal of WEP were detected and in many cases the variability was similar to that of the other treatments. This inherent variability (i.e., pre-grazing) likely results from of a combination of hydrological factors like ground water fluctuations, soil attributes such as texture, ecological dynamics involving plant community composition, and anthropogenic influences like historical land management practices (McClain et al., 2003; Vidon et al., 2010). In particular, the species cover information (Figure S2) demonstrates a wide range in species composition and abundance, this coupled with the variation in P release with different vegetation species may explain some of the observed variability (Cober et al., 2018).

3.4 Sources of variability and uncertainty in P sources

The Prairie pothole wetlands regularly experience high water levels in the early spring. Automated observations made with a water level logger adjacent to one plot between October 2020 and May 2021 showed that the lower, middle, and upper sampling points experienced inundation for approximately 21, 11, and zero days, respectively (Noyes et al., 2024). The annual weather conditions and topography of riparian areas surrounding the wetlands will impact the length and extent of flooding. Prolonged contact with water has been shown to increase the mass of WEP lost in both soil (Young and Briggs, 2008) and vegetation (Lozier and Macrae, 2017) and may also explain some of the observed variability. As reported by Podolsky and Schindler (1993), the soils surrounding these potholes are typically low in CaCo₃ and have a neutral to slight alkaline pH. In this pH range (6.5 to 7.5) P availability is typically at its highest and not expected to precipitate with Ca. A more detailed soil chemical analysis, particularly Fe and Mn, along with soil saturation duration information (i.e., redox) would be needed to fully assess the potential for P loss during the spring (Walton et al., 2020). The WEP protocol used for both soil and vegetation samples are not likely to capture mobilize redox-sensitive P from the soil (Walton et al., 2020) or enhanced P leaching from vegetation (Lozier and Macrae, 2017). Similarly, the WEP protocol also does not capture the enhanced P release from soil and vegetation that results repeated freeze-thaw cycles (Liu et al., 2013; Lozier and Macrae, 2017). However, temperature sensors placed at the soil surface adjacent to one plot recorded four freeze-thaw cycles between Oct 2020 and May 2021 and found that surface temperatures fluctuations are moderated in this region by the relatively persistent snowpack (Noyes et al., 2024).), reducing the potential effects of freeze-thaw cycles on P release. However, both the prolonged contact with water and freeze-thaw cycles are not captured in the WEP protocols and may result in an underestimation of the potential for P loss from the each of the four distinctive sources of P in riparian areas.

In addition to climatic effects, there may be variability in P as a side effect of the study design. One source of variability could be from added urine and manure in grazed areas which likely created additional hotspots of P that may carry forward to subsequent years (Subedi et al., 2020; Donohoe et al., 2021). However, there was no indication of P accumulation due to grazing in any of the four distinctive P sources over the 3-year study period. The highest concentrations of WEP were typically found in the second year of the study (Figure S6). This suggests that other biophysical processes regulated by weather conditions (Figure S3) were of greater importance in controlling the WEP concentrations than P additions from cattle urine and manure. Another source of variability may have been from sampling. As there was significant variability among plots, the single 0.25 \(m^2\) sampling quadrat within each riparian location may have been insufficient to capture the spatial variability. Therefore, larger composite and/or several sampling locations within each upper, middle and lower locations are recommended. Appropriate sampling design becomes critical as the scale of observation of similar research increases to the farm scale, and so will the range and sources of variability. As the scope of research is expanded to the farm level, the importance of using an appropriate sampling design becomes increasingly critical (Hale et al., 2014).

The single 0.25 \(m^2\) sampling quadrate within each riparian location may have been insufficient to capture the spatial variability. Therefore, larger composite and/or several sampling locations within each upper, middle and lower locations are recommended. There was no indication of P accumulation due to grazing in any of the four distinctive P sources over the 3-year study period. The highest concentrations of WEP were typically found in the second year of the study (Figure S6). This suggests that other biophysical processes regulated by weather conditions (Figure S3) were of greater importance in controlling the WEP concentration than any possible P additions from cattle urine and manure. Appropriate sampling design becomes critical as the scale of observation of similar research increases to the farm scale, and so will the range and sources of variability. As the scope of research is expanded to the farm level, the importance of using an appropriate sampling design becomes increasingly critical (Hale et al., 2014).

3.5 Management implications

Autumn was selected for the mowing and grazing treatments for three reasons. The first was to reduce the mass of biomass P available that can contribute to the P loss during the spring snowmelt. Second, drier soil conditions reduce the extent of pugging and soil compaction, which limits the disruption of soil structure and damage to plants (Batey, 2009). ). Lastly, the Prairie potholes and associated riparian areas are important breeding habitats for migratory birds. Grazing can negatively affect these species, but late-season grazing may reduce this potential ecological impact (Stanley and Knopf, 2002). However, the type of grazing system (timing, stocking rate, and density, etc.) may impact habitat quality and breeding success (Carnochan et al., 2018; Hansen et al., 2019; Kraft et al., 2021). Corridor fencing at the edge of the waterbody and alternative water sources were used in this study to limit livestock access in order to prevent bank erosion and protect water quality (e.g., direct deposition) (Dauwalter et al., 2018). Scaling this to the farm level would require virtual fencing or infrastructure (Aarons et al., 2013) and time (to conduct short-term grazing), especially in Prairie pothole regions where there are numerous and small riparian areas (Sovell et al., 2000; Hubbard et al., 2004; Hulvey et al., 2021; Manitoba Agriculture, 2024). The long-term impacts of repeated grazing of riparian areas also need to be considered. From a nutrient loss reduction perspective, a shift in the magnitude of P sources could be expected as less biomass is available to be added to the litter source, affecting the organic layer and Ah sources of P. The regular inclusion of cattle will also introduce a new manure source of P, which can spatially redistribute P and initially be more water soluble and readily transported (Franzluebbers et al., 2019). Grazing can also reduce the litter layer through trampling increasing the soil-vegetation contact and speeding up the decomposition process. These changes in biomass and litter quantities may result in changes to habitat structure. Although this study generally considers environmental implications, forage management practices also have an agronomic effect which should be taken into consideration when developing best management practices (Subedi et al., 2020).

4 Conclusion

Biomass and litter are significant sources of near-surface WEP in riparian areas that have been historically disregarded in studies. Management of the biomass prior to the onset of winter conditions in cold climates has the potential to reduce the mass of P directly lost during the spring snowmelt and maintain or enhance the nutrient buffering capacity. The results from this experiment demonstrated that short-term high-density cattle grazing and mowing both resulted in a reduction in the mass of biomass WEP, particularly in the lower riparian locations. The grazing and mowing treatments had no detectable effect on the other three near-surface sources of WEP. However, detecting changes in the near-surface sources of WEP is challenging due to high spatial variability.

Additional work on riparian management strategies is needed to address the specific challenges posed by cold climates. In these regions, the runoff and nutrient losses occur predominately during the spring snowmelt period when the ability of riparian areas to trap and retain nutrients is diminished. Further, the repeated FTC of the vegetation and soils increases the potential P losses during this key time. Continued research to identify, quantify, and manage these sources of P to improve water quality remains a priority. In addition to improving water quality, the development of riparian management strategies should prioritize the protection other ecological goods and services and recognize these areas as an integral part of the farm.

Acknowledgements

This project was undertaken with the financial support of the Government of Canada through the federal Department of Environment and Climate Change and a Lake Winnipeg Basin Program grant awarded to the Manitoba Association of Watersheds. Additional research funding was provided through a Brandon University Research Committee grant awarded to AK. Thank you to A. Avila, M. Luna, C Sobchuk, and A. Tan for all the help with lab and field work. Special thanks to the Manitoba Beef and Forage Initiatives research farm staff for the use of their facilities and managing the cattle grazing and mowing treatments. Lastly, thank you to R. Canart and M. Elsinger for helping to develop the experimental design.

Data availability

Data and source code for analysis and manuscript available on GitHub: https://github.com/alex-koiter/riparian-grazing-manuscript

Conflict of interest statement

The authors have no competing interests to declare that are relevant to the content of this article.

Author contributions

The authors confirm contribution to the paper as follows: study conception and design: A. Koiter; data collection: T. Malone; analysis and interpretation of results: A. Koiter; draft manuscript preparation: A. Koiter and T. Malone. All authors reviewed the results and approved the final version of the manuscript.

References

Aarons, S.R., A.R. Melland, and L. Dorling. 2013. Dairy farm impacts of fencing riparian land: Pasture production and farm productivity. Journal of Environmental Management 130: 255–266. doi: 10.1016/j.jenvman.2013.08.060.

Batey, Tom. 2009. Soil compaction and soil management a review. Soil Use and Management 25(4): 335–345. doi: 10.1111/j.1475-2743.2009.00236.x.

Beck, H.E., N.E. Zimmermann, T.R. McVicar, N. Vergopolan, A. Berg, et al. 2018. Present and future Köppen-Geiger climate classification maps at 1-km resolution. Scientific Data 5(1): 180214. doi: 10.1038/sdata.2018.214.

Brooks, M.E., K. Kristensen, K.J. van Benthem, A. Magnusson, C.W. Berg, et al. 2017. glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. The R Journal 9(2): 378–400. doi: 10.32614/RJ-2017-066.

Carlyle, G.C., and A.R. Hill. 2001. Groundwater phosphate dynamics in a river riparian zone: Effects of hydrologic flowpaths, lithology and redox chemistry. Journal of Hydrology 247(3): 151–168. doi: 10.1016/S0022-1694(01)00375-4.

Carnochan, S.J., C.C. De Ruyck, and N. Koper. 2018. Effects of twice-over rotational grazing on songbird nesting success in years with and without flooding. Rangeland Ecology & Management 71(6): 776–782. doi: 10.1016/j.rama.2018.04.013.

Cober, J.R., M.L. Macrae, and L.L. Van Eerd. 2018. Nutrient release from living and terminated cover crops under variable freezethaw cycles. Agronomy Journal 110(3): 1036–1045. doi: 10.2134/agronj2017.08.0449.

Cober, J.R., M.L. Macrae, and L.L. Van Eerd. 2019. Winter phosphorus release from cover crops and linkages with runoff chemistry. Journal of Environmental Quality 48(4): 907–914. doi: 10.2134/jeq2018.08.0307.

Dauwalter, D.C., K.A. Fesenmyer, S.W. Miller, and T. Porter. 2018. Response of riparian vegetation, instream habitat, and aquatic biota to riparian grazing exclosures. North American Journal of Fisheries Management 38(5): 1187–1200. doi: 10.1002/nafm.10224.

Donohoe, G., D. Flaten, F. Omonijo, and K. Ominski. 2021. Short-term impacts of winter bale grazing beef cows on forage production and soil nutrient status in the eastern canadian prairies. Canadian Journal of Soil Science 101(4): 717–733. doi: 10.1139/cjss-2021-0028.

Environment, and Climate Change Canada. 2024. Canadian Climate Normals. https://climate.weather.gc.ca/climate_normals/index_e.html.

Fellows, I. 2023. OpenStreetMap: Access to open street map raster images.

Fitch, L., B. Adams, and K. O’Shaughnessy. 2003. Caring for the green zone: Riparian areas and grazing management. 3rd ed. Cows; Fish Program, Lethbridge, Alberta.

Fox, J., and S. Weisberg. 2019. An R companion to applied regression. Third. Sage, Thousand Oaks CA.

Franzluebbers, A.J., M.H. Poore, S.R. Freeman, and J.R. Rogers. 2019. Soil-surface nutrient distributions in grazed pastures of North Carolina. Journal of Soil and Water Conservation 74(6): 571–583. doi: 10.2489/jswc.74.6.571.

Habibiandehkordi, R., D.A. Lobb, P.N. Owens, and D.N. Flaten. 2019. Effectiveness of vegetated buffer strips in controlling legacy phosphorus exports from agricultural land. Journal of Environmental Quality 48(2): 314–321. doi: 10.2134/jeq2018.04.0129.

Habibiandehkordi, R., D.A. Lobb, S.C. Sheppard, D.N. Flaten, and P.N. Owens. 2017. Uncertainties in vegetated buffer strip function in controlling phosphorus export from agricultural land in the canadian prairies. Environmental Science and Pollution Research 24(22): 18372–18382. doi: 10.1007/s11356-017-9406-6.

Hale, R., P. Reich, T. Daniel, P.S. Lake, and T.R. Cavagnaro. 2014. Scales that matter: Guiding effective monitoring of soil properties in restored riparian zones. Geoderma 228-229: 173–181. doi: 10.1016/j.geoderma.2013.09.019.

Hansen, B.D., H.S. Fraser, and C.S. Jones. 2019. Livestock grazing effects on riparian bird breeding behaviour in agricultural landscapes. Agriculture, Ecosystems & Environment 270-271: 93–102. doi: 10.1016/j.agee.2018.10.016.

Hartig, F. 2022. DHARMa: Residual diagnostics for hierarchical (multi-level / mixed) regression models.

Hille, S., D. Graeber, B. Kronvang, G.H. Rubæk, N. Onnen, et al. 2019. Management options to reduce phosphorus leaching from vegetated buffer strips. Journal of Environmental Quality 48(2): 322–329. doi: 10.2134/jeq2018.01.0042.

Hubbard, R.K., G.L. Newton, and G.M. Hill. 2004. Water quality and the grazing animal1. Journal of Animal Science 82(suppl_13): E255–E263. doi: 10.2527/2004.8213_supplE255x.

Hulvey, K.B., C.D. Mellon, and A.R. Kleinhesselink. 2021. Rotational grazing can mitigate ecosystem service trade-offs between livestock production and water quality in semi-arid rangelands. Journal of Applied Ecology 58(10): 2113–2123. doi: 10.1111/1365-2664.13954.

Kelly, J.M., J.L. Kovar, R. Sokolowsky, and T.B. Moorman. 2007. Phosphorus uptake during four years by different vegetative cover types in a riparian buffer. Nutrient Cycling in Agroecosystems 78(3): 239–251. doi: 10.1007/s10705-007-9088-4.

Kieta, K.A., and P.N. Owens. 2019. Phosphorus release from shoots of phleum pretense l. After repeated freeze-thaw cycles and harvests. Ecological Engineering 127: 204–211. doi: 10.1016/j.ecoleng.2018.11.024.

Kieta, K.A., P.N. Owens, D.A. Lobb, J.A. Vanrobaeys, and D.N. Flaten. 2018. Phosphorus dynamics in vegetated buffer strips in cold climates: A review. Environmental Reviews 26(3): 255–272. doi: 10.1139/er-2017-0077.

Kraft, J.D., D.A. Haukos, M.R. Bain, M.B. Rice, S. Robinson, et al. 2021. Using grazing to manage herbaceous structure for a heterogeneity-dependent bird. The Journal of Wildlife Management 85(2): 354–368. doi: 10.1002/jwmg.21984.

Krall, M., and P. Roni. 2023. Effects of livestock exclusion on stream habitat and aquatic biota: A review and recommendations for implementation and monitoring. North American Journal of Fisheries Management 43(2): 476–504. doi: 10.1002/nafm.10863.

Kröger, R., M.M. Holland, M.T. Moore, and C.M. Cooper. 2007. Plant senescence: A mechanism for nutrient release in temperate agricultural wetlands. Environmental Pollution 146(1): 114–119. doi: 10.1016/j.envpol.2006.06.005.

Lacas, J.-G., M. Voltz, V. Gouy, N. Carluer, and J.-J. Gril. 2005. Using grassed strips to limit pesticide transfer to surface water: A review. Agronomy for Sustainable Development 25(2): 253–266. doi: 10.1051/agro:2005001.

Lenth, R.V. 2024. Emmeans: Estimated marginal means, aka least-squares means.

Liu, J., H.M. Baulch, M.L. Macrae, H.F. Wilson, J.A. Elliott, et al. 2019a. Agricultural water quality in cold climates: processes, drivers, management options, and research needs. Journal of Environmental Quality 48(4): 792–802. doi: 10.2134/jeq2019.05.0220.

Liu, J., R. Khalaf, B. Ulén, and G. Bergkvist. 2013. Potential phosphorus release from catch crop shoots and roots after freezing-thawing. Plant and Soil 371(1): 543–557. doi: 10.1007/s11104-013-1716-y.

Liu, J., M.L. Macrae, J.A. Elliott, H.M. Baulch, H.F. Wilson, et al. 2019b. Impacts of cover crops and crop residues on phosphorus losses in cold climates: a review. Journal of Environmental Quality 48(4): 850–868. doi: 10.2134/jeq2019.03.0119.

Lozier, T.M., and M.L. Macrae. 2017. Potential phosphorus mobilization from above-soil winter vegetation assessed from laboratory water extractions following freezethaw cycles. Canadian Water Resources Journal 42(3): 276–288. doi: 10.1080/07011784.2017.1331140.

Lüdecke, D., M.S. Ben-Shachar, I. Patil, P. Waggoner, and D. Makowski. 2021. performance: An R package for assessment, comparison and testing of statistical models. Journal of Open Source Software 6(60): 3139. doi: 10.21105/joss.03139.

Lyu, C., X. Li, P. Yuan, Y. Song, H. Gao, et al. 2021. Nitrogen retention effect of riparian zones in agricultural areas: A meta-analysis. Journal of Cleaner Production 315: 128143. doi: 10.1016/j.jclepro.2021.128143.

Manitoba Agriculture. 2023. Manitoba agriculture weather program. https://www.gov.mb.ca/agriculture/weather/manitoba-ag-weather.html.

Manitoba Agriculture. 2024. Agriculture rotational grazing. https://www.gov.mb.ca/agriculture/.

Manitoba Beef & Forage Initiatives. 2024. MBFI farm stations. https://www.mbfi.ca/farm-station.

Massicotte, P., and A. South. 2023. Rnaturalearth: World map data from natural earth.

McClain, M.E., E.W. Boyer, C.L. Dent, S.E. Gergel, N.B. Grimm, et al. 2003. Biogeochemical hot spots and hot moments at the interface of terrestrial and aquatic ecosystems. Ecosystems 6(4): 301312. doi: 10.1007/s10021-003-0161-9.

McGuire, K.J., and J.J. McDonnell. 2010. Hydrological connectivity of hillslopes and streams: Characteristic time scales and nonlinearities. Water Resources Research 46(10): W10543. doi: 10.1029/2010WR009341.

Murphy, J., and J.P. Riley. 1962. A modified single solution method for the determination of phosphate in natural waters. Analytica Chimica Acta 27: 31–36. doi: 10.1016/S0003-2670(00)88444-5.

Noyes, I., A.J. Koiter, H. Jarvie, J.M. Plach, D.A. Lobb, et al. 2024. Potential phosphorus mobilization from riparian vegetation following freezing. Journal of Environmental Management.

Nsenga Kumwimba, M., J. Huang, M. Dzakpasu, K. De Silva, O.E. Ohore, et al. 2023. An updated review of the efficacy of buffer zones in warm/temperate and cold climates: Insights into processes and drivers of nutrient retention. Journal of Environmental Management 336: 117646. doi: 10.1016/j.jenvman.2023.117646.

Owens, P.N., J.H. Duzant, L.K. Deeks, G.A. Wood, R.P.C. Morgan, et al. 2007. Evaluation of contrasting buffer features within an agricultural landscape for reducing sediment and sediment-associated phosphorus delivery to surface waters. Soil Use and Management 23(s1): 165–175. doi: 10.1111/j.1475-2743.2007.00121.x.

Podolsky, G.P., and D. Schindler. 1993. Soils of the manitoba zero tillage research association farm. Winnipeg, Manitoba.

R Core Team. 2024. R: A language and environment for statistical computing. http://www.R-project.org.

Reid, K., K. Schneider, and B. McConkey. 2018. Components of phosphorus loss from agricultural landscapes, and how to incorporate them into risk assessment tools. Frontiers in Earth Science 6. https://www.frontiersin.org/article/10.3389/feart.2018.00135.

Roberts, W.M., M.I. Stutter, and P.M. Haygarth. 2012. Phosphorus retention and remobilization in vegetated buffer strips: A review. Journal of Environmental Quality 41(2): 389–399. doi: 10.2134/jeq2010.0543.

RStudio. 2024. RStudio: Integrated development environment for r. http://www.rstudio.org/.

Schindler, D.W., R.E. Hecky, and G.K. McCullough. 2012. The rapid eutrophication of lake winnipeg: Greening under global change. Journal of Great Lakes Research 38, Supplement 3: 6–13. doi: 10.1016/j.jglr.2012.04.003.

Schloerke, B., D. Cook, J. Larmarange, F. Briatte, M. Marbach, et al. 2024. GGally: Extension to ’ggplot2’.

Sharpley, A.N., P.J.A. Kleinman, and J.L. Weld. 2006. Environmental soil phosphorus indices. In: Carter, M.R. and Gregorich, E.G., editors. 2nd ed. CRC Press, Boca Raton, FL, U.S.A

Sovell, L.A., B. Vondracek, J.A. Frost, and K.G. Mumford. 2000. Impacts of Rotational Grazing and Riparian Buffers on Physicochemical and Biological Characteristicsof Southeastern Minnesota, USA, Streams. Environmental Management 26(6): 629–641. doi: 10.1007/s002670010121.

Stanley, T.R., and F.L. Knopf. 2002. Avian responses to late-season grazing in a shrub-willow floodplain. Conservation Biology 16(1): 225–231. doi: 10.1046/j.1523-1739.2002.00269.x.

Subedi, A., D. Franklin, M. Cabrera, A. McPherson, and S. Dahal. 2020. Grazing systems to retain and redistribute soil phosphorus and to reduce phosphorus losses in runoff. Soil Systems 4(4): 66. doi: 10.3390/soilsystems4040066.

Syversen, N., and H. Borch. 2005. Retention of soil particle fractions and phosphorus in cold-climate buffer zones. Ecological Engineering 25(4): 382–394. http://www.sciencedirect.com/science/article/pii/S092585740500145X.

Tomer, M.D., T.B. Moorman, J.L. Kovar, D.E. James, and M.R. Burkart. 2007. Spatial patterns of sediment and phosphorus in a riparian buffer in western iowa. Journal of Soil and Water Conservation 62(5): 329–338.

Vidon, P., C. Allan, D. Burns, T.P. Duval, N. Gurwick, et al. 2010. Hot spots and hot moments in riparian zones: Potential for improved water quality management. Journal of the American Water Resources Association 46(2): 278–298. doi: 10.1111/j.1752-1688.2010.00420.x.

Walton, C.R., D. Zak, J. Audet, R.J. Petersen, J. Lange, et al. 2020. Wetland buffer zones for nitrogen and phosphorus retention: Impacts of soil type, hydrology and vegetation. Science of The Total Environment 727: 138709. doi: 10.1016/j.scitotenv.2020.138709.

Wickham, H. 2016. ggplot2: Elegant graphics for data analysis. Springer-Verlag, New York NY U.S.A.

Young, E.O., and R.D. Briggs. 2008. Phosphorus concentrations in soil and subsurface water: A field study among cropland and riparian buffers. Journal of Environmental Quality 37(1): 69–78. doi: 10.2134/jeq2006.0422.

Yu, C., P. Duan, Z. Yu, and B. Gao. 2019. Experimental and model investigations of vegetative filter strips for contaminant removal: A review. Ecological Engineering 126: 25–36. doi: 10.1016/j.ecoleng.2018.10.020.

Supplemental materials

Figure S1: Workflow diagram showing the experimental setup (yellow), field work (green), sample preparation (brown), and laboratory analysis (blue).

Figure S2: Initial year (2019) cover assessment using the foliar cover method for each plot within the four riparian locations

Figure S3: Average daily air temperature and cumulative rainfall over the growing season over the three year study. Red bars indicate sampling dates

Figure S4: a) Mass of biomass and litter before grazing and mowing (2019-2021) and b) the bulk density of the organic layer and 10 cm Ah horizon (2023)

Figure S5: Generalized pairs plot showing the data and relationships between WEP concentration between the different sources of Phosphorus at the lower (purple), middle (blue), and lower (green) topographic positions. Data set only includes samples collected before grazing and mowing treatments were applied. Corr indicates the pearson correlation coefficient. *** p-value < 0.001, ** p-value < 0.01, ** p-value < 0.05.

Figure S6: Mean and standard deviation WEP concentration for each of the different sources of Phosphorus at each topographic position over the three year period of observations. Data set only includes samples collected before grazing and mowing treatments were applied.