*2.1. Study Sites*

This study was conducted in eight pastures within the Georgia piedmont; four at the Animal and Dairy Science Eatonton Beef Research Unit (33.420759◦ N, 83.476555◦ W, Elevation 152–177 m, Eatonton) in Putnam County, GA, USA, and four at J Phil Campbell (JPC) Sr. Research and Education Center (33.887487◦ N, 83.420966◦ W; Elevation 213–259 m, Watkinsville) in Oconee County, GA (Figure 1). Pastures are characterized by moderate and wet winters and long and dry summers. Table 1 describes the study pastures at two locations, their respective area coverage, and number of sampling points in pastures (Matrix). Denser sampling indicates areas of interest (AOIs). The soil types [26] in the study locations, along with other hydrologically important characteristics, are summarized in Table 2.

**Figure 1.** Maps showing; (**a**) Eatonton and (**b**) Watkinsville study sites with concentrated flow paths, pasture boundaries (CHD: continuous grazing with hay distribution, and STR: strategic grazing), runoff collectors, exclusions, and sampling locations (black dots).



**Table 2.** Soil classifications of Eatonton and Watkinsville study pastures.


#### *2.2. Experimental Design*

The study began in 2015 at both study sites with Baseline sampling carried out in spring/summer of 2015 (May–June). The 10-year legacy pasture management was continuous grazing management with free movement of cattle and fixed locations for hay feeding, waterers, and shade. All pastures were grazed continuously with 1.7–2.2 cattle head ha−<sup>1</sup> prior to treatment application.

In May 2016, four replications of each of two grazing management treatments: (1) Continuous grazing with hay distribution (CHD) and (2) strategic rotational grazing (STR) were implemented. Cattle in CHD pastures had continued access to all locations inside pastures including areas with high P transport potential. Hay, however, was fed by distributing it to other locations in the pastures instead of conventional hay feeding at fixed locations in the pastures. This was done to prevent congregation of cattle at vulnerable locations during hay feeding. In this study, areas where cattle tended to congregate and that were also either close to the streams and/or areas in higher elevation, which had high P transport potential with steep slopes (6–27% slope), were designated as areas of interest (AOIs). In STR pastures, exclusions were setup by fencing the areas with high transport potential of P (Figure 1) and were over-seeded to maintain forage productivity and ground cover. Movable farm equipages: shades, waterers (with quick connects at different locations), and hay feeding rings were used to lure cattle to different locations strategically in STR pastures. Cattle density of all pastures during post-treatment was 1.1 cattle head ha<sup>−</sup>1.

In STR pastures, exclusions were over-seeded with pearl millet (*Pennisetum glaucum*), crabgrass (*Digitaria* spp), and cow pea (*Vigna unguiculata*) in the spring, and with crimson clover (*Trifolium incarnatum*), canola (*Brassica napus*), ryegrass (*Lolium*) and cereal rye (*Secale cereale*) in the fall. The exclusions were flash grazed based on the height of forage. Adaptive rotation strategy was used based on forage availability. Shade and waterers were moved to lure animals away from AOIs so that manure deposited by animals was less likely to be transported to the edge of fields or streams.

#### *2.3. Sampling*

Soils were sampled on a 50 m grid (matrix samples) in all pastures, with additional samples taken to ensure coverage with the AOIs. At each sampling point, two soil cores were taken, and each was divided into 0–5 cm, 5–10 cm, and 10–20 cm sections, bagged, placed in a dark cooler, and taken to the lab for processing. Soil sampling was carried out from mid-May to June in 2015 (Baseline) and 2018 (Post-Treatment). Within STR sites, sample points that were inside the excluded areas were termed "exclusions" while all other points outside the exclusions were termed "non-exclusions" to further understand the influence of the exclusions.

Pour-point runoff collectors were set up during Baseline (2015) at the drainage outlets of at least three watersheds within each pasture. Surface runoff samples generated by medium to high intensity rainfall (>20 mm) were collected from pour-point collectors, which had replaceable Nalgene bottles. Runoff samples were collected from June 2015 through April 2016 for Baseline (2015). Post-treatment runoff samples were collected from May 2016 to December 2018 as 2016, 2017, and 2018 sampling years. As 2016 was a drought year, there were no runoff samples for a 10-month period and this year was excluded from the analysis. Samples from each runoff event were brought to the lab and a 50 mL sub-sample was filtered (0.45 μm) within 48 h of the end of each event. The filtered runoff sub-samples and unfiltered runoff samples were stored in acid-washed Nalgene bottles at −10 ◦C for further analyses.

#### *2.4. Analysis of Soil and Water Samples*

Mehlich-1 phosphorus (M1P), which is a measure of plant available P in soil, was calculated using spectrophotometric analysis of soil extracts obtained from a two-component acid mixture [27]. Five grams of each soil sample was mixed with 20 mL of the double acid mixture (HCl and H2SO4) and shaken for 5 min, followed by filtering using a Whatman #42 filter. The extract was then analyzed for M1P employing the molybdenum blue dye procedure with spectrophotometric determination at 882 nm [28]. The analyses were expressed on a dry-weight basis by taking soil moisture into account. The concentrations of M1P were then converted to mass of M1P in kg P ha−<sup>1</sup> using the bulk density of the soil samples calculated for another part of the study.

Dissolved reactive P (DRP) in runoff samples was measured in filtered runoff samples using spectrophotometric analysis on a Tecan Infinite Pro 200 (Tecan Group Ltd., Männedorf, Switzerland) with the Molybdenum-Blue dye [28]. Total Kjeldahl P (TKP) in runoff samples was determined from unfiltered runoff samples. A digestion solution containing HgSO4, K2SO4 and concentrated H2SO4 was used to digest the unfiltered water samples at 114 ◦C for 14 h followed by digestion at 380 ◦C for another 2.5 h [29]. Phosphorus in the digests was determined using the Murphy–Riley procedure [28].

#### *2.5. Conversion of Concentration to Loads of DRP and TKP in Runo*ff

Amount of runoff was estimated based on the amount of precipitation received, hydrologic condition of the soil and cover on the ground surface using the Curve Number (CN) method for estimating runoff [30]. It should be noted that ground cover conditions varied between post treatment calculations and the hydrologic conditions did not change for post treatment calculations. The CN is related to the potential maximum retention after runoff begins (S) as,

$$\mathbf{S} = \text{(1000/CN)} - \text{10} \tag{1}$$

Initial abstraction (Ia), which represents all losses before runoff begins, is related to S as,

$$\mathbf{I}\_a = \mathbf{0}.2\mathbf{S} \tag{2}$$

and was used to estimate the amount of runoff for each rainfall event using Equation (3),

$$\mathbf{Q} = (\mathbf{P} - 0.2\mathbf{S})^2 / (\mathbf{P} + 0.8\mathbf{S}) \tag{3}$$

where Q is the amount of runoff, and P is the amount of precipitation. All values of DRP and TKP concentrations were converted to loads in runoff using the amount of runoff generated during each runoff event.

#### *2.6. Spatial Analysis*

A Trimble R10 GPS unit (Trimble, Sunnyvale, CA) was used to locate sampling points within the study pastures on a 50 m grid and to collect elevation every 2 m (4 cm resolution) for the development of digital elevation models (DEMs). ArcGIS 10.6.1 (ESRI, Redlands, CA, USA) was used to process the data by storing them in geodatabases. The DEMs were used to identify locations for pour-point collectors and delineation of watersheds. Values of the M1P concentrations in mg kg−<sup>1</sup> of each soil depth (0–5, 5–10, and 10–20 cm) at point locations were converted into raster files using "Create TIN" and "TIN to Raster" tools to create a continuous surface (raster file) from the point file (ESRI). The raster files (M1P concentrations) were then used to visualize differences between 2015 and 2018 samples for M1P as 2018 minus 2015 raster for each depth. "Raster calculator" was used to carry out this analysis.

Spatial autocorrelation of M1P was estimated for each pasture under study using the Global Moran's I tool. Clusters of P in pastures were determined from the point values of M1P using the "Hot Spot Analysis" tool in ArcGIS (ESRI). Since sampling was conducted on a 50 m grid, 55 m was used to ensure more than one neighbor for each point, which would also ensure that not all points were neighbors of each other. The Hot Spot tool calculates the Getis-Ord GI\* statistic [31,32], which is a z-score for each point under consideration, along with a *p*-value denoting significance. Higher *z*-values (0 to 1) would mean clustering of either high or low P concentrations.

#### *2.7. Statistical Analysis*

Comparison of the amount of M1P in soil between sampling periods, 2015 and 2018, was carried out for individual treatments (Baseline-CHD-2015 to Post-Treatment-CHD-2018, and Baseline-STR-2015 to Post-Treatment-STR-2018) using one-way analysis of variance (ANOVA) with differences analyzed by the Wilcoxon Test, which compared median M1P concentration. Comparison of the amount of M1P in soil between CHD and STR pastures at each sampling period (2015 or 2018) was also carried out using one-way ANOVA with differences analyzed by the Wilcoxon Test comparing median M1P concentrations. Comparisons between "exclusions" and "non-exclusions" in STR pastures for 2015 and 2018 sampling periods were carried out using one-way ANOVA. Non-exclusions are all other parts of STR pastures outside of fenced exclusions. The differences revealed were analyzed by the Wilcoxon Test. The Wilcoxon Test was used because of the non-normal M1P distribution.

Concentrations and loads of DRP and TKP in runoff were compared between sampling years (2015, 2017, and 2018) for each of the treatments (CHD or STR) and between CHD and STR treatments for each of the sampling years using one-way ANOVA. Differences were analyzed by the Wilcoxon Each Pair Test as DRP and TKP concentrations and loads did not meet the normal distribution assumption for parametric analysis. Simple linear regression model was used to determine the relationships between soil P (0–5 cm layer) and runoff P loads (DRP load, and TKP load) during Baseline and Post-Treatment. Event loads of DRP and TKP of the individual watersheds were fitted with the average soil-P values for each of the individual watersheds for the Baseline and Post-Treatment periods. Harmel et al. [33] did similar analysis on NO3 and PO4 losses from cropland and grazed pastures. The difference in slopes were analyzed by comparing fitted models with runoff P loads as response variables and soil P as the explanatory variable. Significance of interaction between sampling periods (Baseline vs. Post-Treatment) and soil P denotes difference in regression slopes between Baseline and Post-Treatment. Test of significance was conducted at 0.05 level of significance in all cases. All statistical analysis was

carried out using the JMP software package (JMP®, Version 14. SAS Institute Inc., Cary, NC, USA, 1989–2019).

#### **3. Results and Discussion**
