Managing for structural complexity to enhance forest ecosystem health and resiliency is increasingly incorporated in silvicultural treatments. High spatial variability in stands managed for structural complexity could obscure the effects of forest management on surface soils. Yet few studies have assessed how within-stand variation in forest structure and other local controls influence surface soil organic matter dynamics over time following timber harvests. We used a stratified random sampling design to capture variation in stand age, legacy structure, soil type, and topography in a second-growth, oak-hardwood forest in the northeastern U.S. We compared surface soil carbon and nitrogen content and availability in 15 harvested stands managed to promote tree regeneration (n = 144 plots) and five unharvested controls (n = 48 plots). We also examined changes over time since harvest in just the harvested stands using a 22-year chronosequence. Forest management strongly influenced surface soil carbon and nitrogen dynamics. The timber harvests had lower soil carbon and nitrogen, microbial biomass, and carbon mineralization but higher nitrogen mineralization. These differences were more pronounced in the drier, less fertile soil type than in more moist, fertile soils. Across the 22-year chronosequence, topography, soil type, and downed woody material density dictated the direction of changes in surface soil carbon and nitrogen over time. Soil carbon and nitrogen accrued over time at drier, higher elevation (~300 m) sites and under higher densities of fine woody material but declined at lower elevations (~180 m) and under lower fine woody material. Proximity to legacy trees was associated with higher soil carbon and nitrogen concentrations and availability. Our findings underscore the importance of silvicultural practices that retain structural legacies and downed woody material in shaping surface soil carbon and nitrogen dynamics over time. Our results also highlight how accounting for spatial variation in local controls on soil carbon and nitrogen, such as topography, can improve detection of changes from forest management practices that increase spatial heterogeneity within stands, such as irregular shelterwood and seed tree regeneration methods.

Sampling design 

We conducted our work at Yale-Myers Forest (41° 57' N, 72°07’ W), a 3,213-ha mixed hardwood, second-growth forest in Connecticut, USA. We stratified our sampling by two soil types, nine timber harvest years, four topographic classes, and two distances relative to legacy overstory trees, which yielded 144 plots in 15 harvested stands managed through the regeneration treatments. We also established 48 plots in five unharvested stands following a similar sampling design.

To stratify sampling by soil type, we first identified all U.S. soil series at Yale-Myers Forest. We limited our sampling to oak-hardwood stands, which are most commonly managed through shelterwood and seed tree regeneration methods. Three soil series collectively accounted for 86% of oak-hardwood stands at the Forest: Nipmuck-Brookfield complex (39%), Paxton-Montauk fine sandy loams (26%), and Woodbridge fine sandy loams (21%). We excluded other soil types since they were not adequately replicated. We aggregated the Paxton-Montauk and Woodbridge series since they share the same geological features and parent material and differ only in drainage class (well-drained versus moderately well-drained, respectively). 

For each of these soil types, we identified areas of stands from nine harvest years that spanned the 1994 to 2016 chronosequence (2 soil types ´ 9 timber harvest years = 18 study areas). We preferentially selected stands that included areas with each of the two soil types to capture within-stand heterogeneity. Within these 18 areas, we used the “cluster” function in R (version 4.0.3) with the CLARA K-medoids method to run an unsupervised classification on slope, aspect, and elevation raster bands from publicly available Digital Elevation Model data to create four classes with the most distinct topographic features. We randomly located points at least 20 m from the harvest edge for each of the resulting 72 polygons (2 soil types ´ 9 timber harvest years ´ 4 terrain classes). Finally, we established plots 2 m and 10 m from the closest legacy tree (≥20 cm diameter at breast height [DBH] at 1.37 m) to stratify sampling by legacy overstory structure (Fig. 1). This sampling design yielded 144 plots in the harvested stands managed to promote tree regeneration (2 soil types ´ 9 timber harvest years ´ 4 terrain classes ´ 2 distances relative to legacy trees).

For the unharvested stands, we chose five stands with the same soil series (Nipmuck-Brookfield complex and Paxton-Montauk/Woodbridge) that had mean elevation (236–272 m, mean = 251 m) and slope (slope: 1.24–21.1°, mean = 10.7°) values within the range of the plots in the irregular shelterwood/seed tree stands (179–305 m, mean = 246 m; 0.64–28.9°, mean = 9.9°). Across the five stands, we identified three areas with no recorded timber harvests since Yale University acquired the property in 1930, two areas that were crown thinned for firewood 40+ years ago (1973 and 1979), and two areas that were crown thinned for firewood 35-40 years ago (1980 and 1984). We used a similar stratified random sampling design as the regeneration treatments to locate points within each of the unharvested areas (2 soil types ´ 3 time-since-thinning classes ´ 4 terrain classes ´ 2 tree distances = 48 plots).

Vegetation and soil sampling

For vegetation sampling, we identified the species and measured DBH of all living trees and shrubs ≥5-cm DBH in a 5-m radius fixed area plot and ≥1-cm DBH in a nested 1-m radius plot (Fig. 1). We calculated coarse woody material (CWM) volume for the portion of each piece of downed wood >10 cm in diameter that fell within the 1-m radius plot and tallied all fine woody material (FWM) pieces ≥2-cm diameter but less than 10 cm.  

We collected and pooled eight 2-cm diameter, 10-cm deep soil cores in the 1-m radius plots. In plots with low soil bulk density, we collected an additional 1–4 cores to ensure we had sufficient sample for the lab assays. We recorded the exact depth of each core to calculate bulk density. We kept field collected soils in coolers and stored them at 4°C for less than 4 weeks prior to analysis. We measured volumetric water content (VWC) to a depth of 12 cm at four locations using a HS2 HydroSense Soil Moisture Sensor (Campbell Scientific) and measured the depth of the Oa horizon at three locations within each 1-m radius plot. We also measured soil temperature to a depth of 30 cm using an HI 145 T-Shaped Soil Thermometer (Hanna Instruments, Inc.) and recorded the date and time of day of the measurement for temperature corrections.

Laboratory assays

We weighed field-moist samples for use in calculating soil bulk density and homogenized and passed each sample through a 4-mm sieve. We mixed samples with water in a 1-to-1 volumetric ratio and measured the pH of the supernatant after 10 min using a benchtop meter (VWR sympHony Sb70p). We measured gravimetric water content (GWC) by oven-drying field-moist samples to constant mass at 105°C. We estimated water holding capacity (WHC) by saturating each soil sample and allowing it to drain freely for 2-h. We estimated SOM content by calculating mass loss on ignition of soils heated at 550°C for 12-h in a muffle furnace. We used a ball mill to grind air-dried subsamples to a fine powder and analyzed %C and %N with an elemental analyzer (Costech ESC 4010, Costech Analytical Technologies Inc.Valencia, California, USA).

We used a series of laboratory assays to assess microbially available C and N, which we expected be more responsive to management than C and N concentrations. We estimated active microbial biomass with a modified substrate-induced respiration (SIR) method that uses autolyzed yeast extract as a labile C substrate. We then measured rates of CO₂ efflux over a 4-h incubation period with an Infra-Red Gas Analyzer (Li-COR model Li-7000, Li-Cor Biosciences, Lincoln, Nebraska, USA). To estimate potential C and N mineralization rates, we incubated samples for 30-d at 20°C under a humid atmosphere and maintained 65% WHC, which is within the optimal range for microbial activity (Langenheder & Prosser, 2008; Paul et al., 2001). For the C mineralization assay, we calculated cumulative C-CO₂ production rates over the 30-d incubation period by integrating CO₂ efflux measurements sampled over 24-h periods at days 1, 4, 11, 24, and 30. For the 30-d net potential N mineralization assay, we shook 6 g dry-weight-equivalent of each soil sample with 25 mL 2M KCl, refrigerated the samples for 12 h, filtered the samples with Whatman Grade 42 filters, and analyzed [NH₄⁺] and [NO₃^-] using a flow analyzer (Astoria 2, Astoria‐Pacific, Clackamas, Oregon, USA). We repeated this procedure for both the initial measurements and day 30 measurements. Net potential N mineralization was calculated as the difference between the initial and final sum of [NH₄⁺] and [NO₃^-] over the 30-d incubation period.