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).