Methods
This study took place on BNI Mine Ltd property, located approximately 6 km SE of Center, ND (lat. 47o02’ 52.66” N long. 101°14’31.61” W). This project site falls within the Northern Great Plains ecoregion with temperatures ranging between -11°C and 22°C (NDAWN 2020), and an average of 150 frost-free days per year. Average growing season precipitation is 355 mm (NDAWN 2021). Soils consisted of silt loams, loams, and silty clay loams complexes prior to mining; however, they are then reclassified as “mined-land complex” (USDA-NRCS 2021a) post-mining.
The project site was stripped mined between 2015 and 2016, and subsequent reclamation completed by April/May of 2018. Installation of our research plots occurred during reclamation and comprised of the following factors: two seeding mixtures, deep ripping at one of two phases of reclamation, and incorporation of mulch into the subsoil horizon. All combinations of factors are represented at the designated location.
Half of the treatments were planted to a grass only seed mixture (G) and the other half seeded to a grass and forb seed mixture (G/F) (Table 1). Deep ripping occurred either within the subsoil horizon prior to topsoil replacement (SSR) or across both the subsoil and topsoil horizons after topsoil replacement (TSR). Ripping shanks reached maximum depths of 56 cm. Straw mulch was applied to the subsoil surface at 1130-1360 kg/ha and then incorporated into the surface of the soil via disking. The reference site was seeded with the grass/forb mix and proceeded with standard reclamation practices which did not include ripping or mulch. Each combination was replicated twice and the reference site once, and each unit was approximately 0.6 hectares.
Plant community composition and canopy cover surveys were conducted by randomly establishing three 60-meter transects in each treatment unit. We placed a 0.5 m2 frame every 15-meters along each transect and identified every plant to a species level. We measured abundance by estimated canopy cover of each species and assigned one of the following modified Daubenmire cover classes, 1= trace -1%, 2= 1-2%, 3= 2-5%, 4= 5-10%, 5= 10-20%, 6= 20-30%…13= 90-95%, 14= 95-98%, 15= 98-99%, and 16= 99-100% (Daubenmire, 1959).
Surveys were conducted in 2019 and 2020 during peak production (mid-July) and mid-point values were used for analysis. We referenced the USDA Plant Database (USDA-NRCS 2021b) to classify each species’ seasonality, metabolic pathway, and origin (native versus exotic) after completion of surveys.
We installed three soil moisture access tubes into the soil profile of each treatment unit. Tubes were a minimum of 30-meters in from the plot’s edges, a minimum of five-meters from other units and installed to a depth of one meter. A soil moisture probe was used to take three separate readings at 10, 20, 30, 40, 60, and 100 cm intervals (Delta-T Device Ltd, UK). We rotated the probe 120o before the second and third readings to account for any soil moisture variability within the soil profile. Our analysis used averages from each tube.
Four penetration resistance readings were taken to a depth of one meter with an automated dynamic cone penetrometer (ADCP; Vertek, USA). We obtained each reading by positioning the ADCP approximately three meters away from the associated access tubes in each of the cardinal directions. Readings automatically recorded were later converted to joules per meter. The method for calculating penetration resistance is based on the soil’s capacity to cease work being performed by the penetrometer divided by how far the penetrometer progressed:
\(R_{s}=\frac{W_{s}}{P_{d}}\) (1)
Where: Rs is the soil resistance (N), Wsis the kinetic energy of the (J), and Pd is the depth traveled by the penetrometer through the soil (Equation 1).
We used averaged species cover values calculated from all frames within a transect across each treatment replication. Generalized linear models were used to test the effects of seeding method, ripping technique, and mulching and their interactions on species richness, diversity, and KBG cover. Tukey post-hoc tests followed any significant ANOVA results (p< 0.10). We performed a separate analysis for each year of data for models of species richness, diversity and KBG cover. We ran nonmetric dimensional scaling analyses (NMDS) and permutational multivariate analysis of variance (PERMANOVA) for 2019 and 2020 separately to understand dissimilarities in species composition between each treatment (Oksanen et al. 2019). Ordinations were calculated in three dimensions using Bray-Curtis distance measures. To further assess the differences in community composition we added plant functional trait data as vectors using envfit . We performed all our analyses using R version 4.0.4 (R Development Core Team 2021) and the following packages car (Fox & Weisberg 2019), agricolea (Mendiburu 2020), vegan (Oksanen et al. 2019).
We used linear mixed effect models to analyze each of the soil moisture depth intervals. We included treatment as the fixed effect and year, replication, and observation point random effects. Tukey’s post- hoc tests were performed for those soil moisture depth intervals that returned significant differences between treatments (p< 0.10). We used lme4 (Bates et al. 2015) for our linear mixed effect model and lsmeans (Lenth 2016) for our post- hoc tests. These same packages were used in the following section as well.
We selected three separate depth bins of 0-15, 15-30, and 30-100 cm after a preliminary analysis determined no significant differences for depth intervals falling between 30 and 100 cm. We averaged the four ADCP readings at each observation point to simplify our model random effects structure. We found no interaction between treatment and depth allowing treatment to be a fixed effect in our linear mixed effect model. Volumetric soil moisture readings were included as a covariate to our models to account for the variability of soil moisture at different depths which may influence PR values. We used the 10, 20, and 60 cm soil moisture depth interval averages for the 0-15, 15-30, and 30-100 cm depth bins, respectively. Readings from the ADCP took place in the same location over the course of three years. Therefore, we included year, replication, and observation point (i.e., location associated with the access tube) as nested random effects in our model. We ran this linear mixed effect model separately for all three depths. We performed subsequent Tukey’s post-hoc procedures with 90% confidence intervals when the model returned significant differences between treatments (p< 0.10).