Experimental approach and initial forest response to a simulated ice storm experiment in a northern hardwood forest

Ice storms are a type of extreme winter weather event common to north temperate and boreal forests worldwide. Recent climate modelling studies suggest that these storms may become more frequent and severe under a changing climate. Compared to other types of storm events, relatively little is known about the direct and indirect impacts of these storms on forests, as naturally occurring ice storms are inherently difficult to study. Here we describe a novel experimental approach used to create a suite of ice storms in a mature hardwood forest in New Hampshire, USA. The experiment included five ice storm intensities (0, 6.4, 12.7 and 19.1 mm radial ice accretion) applied in a single year, and one ice storm intensity (12.7 mm) applied in two consecutive years. Results demonstrate the feasibility of this approach for creating experimental ice storms, quantify the increase in fine and coarse woody debris mass and nutrients transferred from the forest canopy to the soil under the different icing conditions, and show an increase in the damage to the forest canopy with increasing icing that evolves over time. In this forest, little damage occurred below 6.4 mm radial ice accretion, moderate damage occurred with up to 12.7 mm of accretion, and significant branch breakage and canopy damage occurred with 19.1 mm of ice. The icing in consecutive years demonstrated an interactive effect of ice storm frequency and severity such that some branches damaged in the first year of icing appeared to remain in the canopy and then fall to the ground in the second year of icing. These results have implications for National Weather Service ice storm warning levels, as they provide a quantitative assessment of ice-load related inputs of forest debris that will be useful to municipalities creating response plans for current and future ice storms.


Introduction
The increase in the frequency, intensity and severity of extreme weather events is emerging as a significant concern of climate change in the Anthropocene [1,2]. These types of episodic events can have a greater impact on natural and managed ecosystems than the more gradual changes in mean temperature and precipitation that are typically associated with climate change [3,4]. Despite the importance of these extreme weather events in shaping ecosystems, relatively little is known about their short-or long-term impacts. They are inherently difficult to study due to infrequent return intervals, high spatial variation in occurrence and intensity, difficulty to forecast reliably, and safety concerns [5,6]. Much of the literature on ecosystem response to extreme weather events is derived from post facto analyses following naturally occurring events such as ice storms, late spring frosts or microbursts [7][8][9]. These studies generally rely on serendipitous pre-event data, lack true controls and rarely have immediate poststorm impact assessments. Further, the historical understanding of ecosystem response to extreme events based on these post facto studies may not be useful in predicting responses to future events, as the underlying template of climate, nutrient reserves, or ecosystem health may have changed over the intervening years [10].
As an alternative to post facto analyses, a growing number of controlled experiments are emerging to help elucidate and model ecosystem responses to changes in the frequency, intensity and severity of extreme events [e.g., 11,12]. These experiments allow for investigation of cause and-effect relationships with rigorous statistical models. They typically include pre-and post-extreme event conditions as well as multiple levels and frequencies of disturbance that allow for determination of response thresholds (i.e., the level of disturbance when an ecosystem irreversibly shifts from one state to another) [4,6,13]. Experiments also have the advantage that they can be implemented under controlled conditions designed to improve safety associated with working in hazardous conditions. Ice storms are a type of extreme winter weather event that are a major driver of forest disturbance in temperate and boreal ecosystems worldwide [14][15][16][17][18]. Climate modeling studies suggest that these types of winter storms will become more frequent and/or more intense under a changing climate [19,20]. Impacts of ice storms on forests depend on a variety of factors including amount of ice accretion, stand age, stand health and species composition [21]. Direct impacts range from twig and branch breakage with light-to-moderate icing, to loss of tree tops and uprooting of entire trees with severe icing [16,[22][23][24]. Based on post facto studies, these changes can lead to a cascade of changes in forest ecosystem structure and function, including increased inputs of coarse and fine woody debris to the forest floor, increases in light below the forest canopy, reductions in photosynthetic capacity and stored plant carbohydrates, alterations of soil temperature and moisture, changes in belowground plant and microbial processes, and altered leaching and gaseous nutrient losses from the soil [5,8,[25][26][27].
Despite this general understanding of the impacts of ice storms on forest ecosystems, manipulation of icing intensity and frequency followed by investigation of ecosystem effects have not been attempted to date. Given the widespread occurrence of these extreme weather events in north temperate and boreal regions, their immediate and lasting impacts, and projections for increased frequency and intensity under a changing climate, there is an urgent need preparation of the manuscript. The specific roles of these authors are articulated in the 'author contributions' section. The funder also provided support in the form of a contract to author FB. FB is an engineer and sole proprietor of Research Designs, and he provided specialized design services for this experiment. The funder had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript relative to this contract. The specific role of FB is articulated in the 'author contributions' section.

PLOS ONE
to better understand the impacts of a range of ice loadings and return intervals on forest ecosystems. The objectives of this paper are to (1) provide a detailed methodological approach for creating experimental ice storms, (2) discuss the efficacy of the design (i.e., our ability to create glaze ice on surfaces), (3) describe possible artifacts of the treatments, and (4) present results on the initial impact of these experimental storms on nutrient transfers in fine and coarse woody debris and tree canopy damage assessments. This work was conducted at the Hubbard Brook Experimental Forest (HBEF) in New Hampshire, builds on a feasibility study at the same site [17], and complements reporting by Fahey et al. [28] on the response of forest canopy structure to these experimental perturbations. Fahey et al. [28] detailed how experimental ice storms increased canopy openness, light transmission and complexity. Here, we describe the initial loss of wood and nutrients from the canopy that contributed to those changes in canopy structure. The detailed description of the methods and initial results on fine and coarse woody debris inputs and tree canopy damage will aid in the design of future studies on ice storm impacts specifically, and extreme weather events more generally, and inform utility companies, emergency management organizations, forest managers, municipalities and the public on how to better prepare for future ice storms.

Site description
The HBEF is a 3,160-ha valley in central New Hampshire, USA (43˚56'46"N, 71˚47'19"W) ( Fig  1A). The climate is humid continental with mean monthly air temperatures ranging from -9 C in January to 18˚C in July [29]. Annual precipitation averages 1400 mm (30% as snow) and is distributed relatively evenly throughout the year. Much of the Hubbard Brook Valley was cut approximately 100 years ago, with subsequent impacts from a hurricane in 1938 and an ice storm in 1998.

Experimental design
Ten 20×30 m plots were established in May 2015 and were blocked by location, with plots 1 through 5 in one block on the eastern side of the study site, and plots 6 through 10 on the western side (Fig 1B). Within the block, plots were randomly assigned to one of five icing treatments, resulting in two replicate plots per treatment (Fig 1B). Treatments were based on target radial ice accretion, which is defined as the radius of ice added to a cylindrical structure such as a wire or branch [30]. Treatments (in mm radial ice accretion) included Control (0 mm), Low (6.4 mm), Mid (12.7 mm), High (19.1 mm), and Mid×2 (12.7 mm applied in two consecutive winters). The Low and Mid treatments were chosen to correspond to National Weather Service Ice Storm Warnings, which are 6.4 mm for mid-Atlantic states and 12.7 mm for New York and New England [31]. The High treatment was chosen to represent an extreme ice storm comparable to the 1998 storm [8].. The minimum distance between plots was 10 m. The rectangular plot shape maximized the ability to ice the interior of plots from outside the plots.
The entire 20×30-m plot was iced during treatments. All response measurements were made within an inner 10×20-m area, allowing for a 5-m buffer zone ( Fig 1C). The inner 10×20-m plot was further subdivided into eight 5×5-m subplots ( Fig 1C). Four of these subplots were used for non-destructive measurements repeated over time; three were used for destructive measurements; and one was undisturbed and reserved for future measurements.

Ice application
To create ice, stream water (~4˚C) was taken from Hubbard Brook by two low pressure/high volume pumps and supplied to two track-equipped utility terrain vehicles (UTV) via moveable fire hoses. Each UTV was equipped with a BB4 1 centrifugal pump that fed high pressure water to a nozzle mounted next to the pump. Pump nozzles were pressurized to flow at approximately 300 liters per minute. All pumps were gasoline powered. One UTV operated along each of the two long sides of the plot, with each spraying water the full length and half the depth of the plot (Fig 1C; S1 Fig). Water was sprayed up through gaps in the tree canopy, with spray nozzles adjusted so that the water would be deposited over the canopy as a fine mist, freezing on contact with vegetation surfaces. The two UTV teams worked simultaneously, spraying layers of ice on the canopy from opposite sides of the plot. The total amount of time of water application was recorded for each UTV ensemble, with approximately equal icing times on both sides. Icing continued until visual inspection of branches within the plots and physical measurements of ice accretion on twigs indicated that the target ice accretion had been achieved. Final measurements of ice accretion were made immediately following the icing treatments on passive ice collectors, as described below.
Weather conditions during icing. Air temperature was measured using two Campbell Scientific 1 HMP 60 thermistors housed in gill solar radiation shields, with one located near plot 9 and one near plot 2. Wind speed was measured with a RM Young 1 model 05103 anemometer at a weather station approximately 5 km east of the experimental plots (250 m altitude). Application dates were chosen to meet criteria assuring effective ice accumulation on canopy surfaces (air temperature < -4˚C) with minimal drift of mist (i.e. wind speed < 4.5 m/ s) as detailed by Rustad and Campbell [17]. Much of the applications took place at night when low temperature requirements were met.
Water volume and chemistry of application. Volume of water applied per plot was calculated as the product of application time and water flow. To construct chemical budgets for the treatments, a 60 ml aliquot of Hubbard Brook stream water was collected prior to each icing, and analyzed for sulfate (SO 4 ), nitrate (NO 3 ) and chloride (Cl) on a Metrohm 1 Ion Chromatograph; calcium (Ca), magnesium (Mg), potassium (K), and sodium (Na) on an Agilent 1 730 ICP optical emission spectrometer; dissolved organic carbon (DOC) and total dissolved nitrogen (TDN) on a Shimadzu 1 TOCV with a TNM-1 nitrogen (N) detector; and ammonium (NH 4 ) by colorimetry on a SEAL Analytical AQ2 discrete analyzer at the Louis C. Wyman Forest Sciences Laboratory, Durham, NH. Dissolved organic N (DON) was calculated as the difference between TDN and inorganic N (NO 3 -N + NH 4 -N). Total mass of elements applied to the canopy as iced 'precipitation' was calculated as the 'precipitation' chemistry from the stream aliquot multiplied by the 'precipitation' volume for each plot.
Ice accretion. Ice accretion was measured using passive collectors constructed from 2.54 cm diameter birch dowels, with six 30-cm long arms radiating out from a central six-way galvanized steel connector in six directions (0˚, 90˚, 180˚, 270˚, up and down; Fig 3). Four collectors were suspended in the forest canopy using parachute cord (to provide a low resistance surface to facilitate release after icing) at a height of~15 m prior to the icing, with one collector in each of the four cardinal corners of the inner plots.
Radial ice accretion on collectors was measured with two methods. First, ice thickness was measured in the field immediately following icing with digital calipers at multiple points along the six dowel arms. Second, after ice thickness measurements were made in the field, each of the six arms of the collectors were cut at the base with a reciprocating saw, being careful to keep the ice on the dowel intact. Ice-coated dowels were placed individually in labeled zip lock bags and transported back to HBEF Headquarters where the ice was melted at room temperature and the volume of meltwater from each dowel measured. The radial ice accretion equivalent (Req) was calculated using the equation from [30]: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi where D = diameter of the dowel (cm), L = length of the dowel (cm); V = volume of melted ice (mL), and ρ i = density of glaze ice (g cm -3 ).

Canopy throughfall volume and chemistry
Canopy throughfall volume was measured with 20 L (28.7 cm diameter) buckets placed on the ground in the center of each subplot (eight per plot). Two of these buckets were acid washed and designated for chemical analyses. After icing, buckets were retrieved, thawed, and the volume of canopy throughfall in each bucket was measured. Aliquots from the two acid-washed buckets per plot were poured off, refrigerated, and analyzed for SO 4 , NO 3 , Cl, Ca, Mg, K, Na, NH 4 , DOC, TDN and DON by the same methods as for the stream water aliquots. Total throughfall flux for each element in each plot was calculated as the mean plot chemistry multiplied by the mean plot volume.

Fine and coarse woody debris
Fine woody debris (i.e., woody material < 2 cm in diameter; hereafter FWD) and foliar litter were collected in plastic baskets (52L×37W×27H cm) that were placed in the center of each of seven subplots for all plots. Foliar litter was collected during the pre-treatment period 5 August, 2015 to 17 January, 2016. Fine woody debris was collected between two to four weeks after icing treatments in 2016 and ten weeks after icing treatments in 2017 to allow all the near-term broken branches to reach the ground. In instances where larger branches lay on the litter baskets, twigs less than 2 cm were clipped around the perimeter of the basket and were included as part of the sample. All foliar litter and FWD samples were oven dried at 60˚C to constant weight, weighed, ground through a Wiley Mill 1 #40 mesh sieve, and analyzed for total carbon (C) and nitrogen (N) on a Thermo Flash 1 EA 1112 Series CN Analyzer 1 and for total Ca, Mg, K, and phosphorus (P) by microwave digest and analysis on an ICP [32]at the Louis C. Wyman Forest Sciences Laboratory, Durham, NH. Total mass of all elements was calculated as total litter mass multiplied by element concentration.
Coarse woody debris (i.e., woody material �2 cm diameter; hereafter CWD) was collected on a 2m×2m quadrat located within each of the three subplots designated for destructive measurements on each of the main plots ( Fig 1C). On November 16, 2015, prior to icing, all CWD was collected from these quadrats. Post-treatment CWD was collected in July and November of 2016 and 2017. Downed branches and tree boles that crossed quadrat boundaries were cut on the border with a hand saw. Coarse woody debris was sorted into five decay classes as described in Table 1 [after 33]. The wet weight of CWD for all categories was measured in the field for each quadrat. A subsample comprised of approximately twenty percent of the total wet weight of CWD for each category was brought back to the laboratory and oven-dried at 60˚C to a constant weight to obtain oven dry mass and a moisture correction factor. Samples were ground and analyzed for nutrients by the same method as for FWD. Total mass of all elements was calculated as total litter mass multiplied by element concentration.

Tree canopy damage assessment
A qualitative tree canopy damage assessment was made pre-and post-treatment, following methods used after the 1998 ice storm and the HBEF pilot ice storm experiment [34]. In this method, all trees greater than 5.0 cm DBH were identified, tagged, measured (DBH, height), and mapped. Basal area was calculated from tree diameter by species. All tagged trees were assigned to an injury class value based on a visual inspection of the branches, using the following crown injury classes: 1 = no visible damage; 2 = 1-49%; 3 = 50-75%; 4 = >75%; 5 = dead. Damage assessments were made on May 12, 2015, August 9, 2016, August 16, 2017, August 7, 2018 and July 1, 2019.

Statistical analysis
Ice, throughfall and precipitation. A series of linear mixed effects models were used to analyze throughfall ice accretion, water volumes and chemistry. Response variables included ice accretion (mm), water volume (mm depth m -2 ) and mass (mg m -2 ) of 10 analytes (Ca 2+ , Mg 2+ , K + , Na + , Cl -, SO 4 2--S, NO 3 --N, NH 4 + -N, DON, DOC); explanatory variables included treatment and a random effect for plot. The lme4 package in R was used to conduct mixed effects analyses [35,36]. As only one sample was collected per plot for precipitation, plot was not used as a random effect for these data. The emmeans package in R was used to conduct multiple comparisons testing and an α of 0.05 was used to assess significance for both throughfall and precipitation [37]. We calculated the amount of ice accretion retained in the canopy as the volume of precipitation minus the volume of throughfall. Fine and coarse woody debris. Linear mixed effect models were used to analyze FWD and CWD total mass, C, N, P, Ca, Mg, and K mass. Response variables included FWD and CWD total mass and nutrients; explanatory variables included treatment (for FWD and CWD), sampling period (for CWD) and treatment × sampling period (for CWD). Plot was incorporated as a random effect. An α of 0.05 was used to compare treatments (for FWD and CWD) and sampling periods (for CWD). Cumulative link mixed models were used to assess decay class as an ordered, categorical variable [38]. CWD mass was used as weights in the decay class models. Differences among treatments were calculated in each sampling period, using treatment as a fixed effect and plot as a random effect. Decay classes from the July 2016 and 2017 samples were compared for Mid×2 plots. Weighted means of decay classes are presented as numbers to illustrate differences among treatments.
Tree canopy damage assessment. Cumulative link mixed models were used to assess damage as an ordered, categorical variable. Fixed effects were treatment × sampling period × species; and tree was used as a random effect. Only trees within the inner 10×20 m buffer were included in the analysis. Eight trees were excluded from analyses from one year of sampling (n 2017 = 2, n 2018 = 6) because they were not found or were not assigned a damage class. Multiple comparisons testing with an α of 0.05 were used to compare treatments and sampling periods. We conducted three rounds of analyses: 1) all species including all treatments, species was not included as a fixed effect in this model; 2) American beech, sugar maple, and yellow birch including all treatments; and 3) American beech, red maple, sugar maple, and yellow birch, including the Control, Mid, Mid×2, and High treatments as there were no red maple trees in the Low treatment. Weighted means of damage classes are presented as numbers to illustrate differences among treatments, sampling periods and species.

Weather conditions during icing
Air temperature and wind speed during the spray applications were generally within the target range of < -4˚C and < 4.5 m/s ( Table 2). Mean temperatures across the plots during spray applications ranged from -12.1 to -5.3˚C, and mean wind speeds across the plots ranged from 0.23 to 1.2 m/s.

Water and chemical inputs
Because stream water used for the applications had greater concentrations of most chemical constituents than ambient precipitation at Hubbard Brook, inputs of many chemicals (especially base cations, Cl and SO 4 ) exceeded those of natural ice storms and comprised inputs comparable in magnitude to annual precipitation inputs at the site (Table 3). Chemical inputs generally increased with increased treatment applications. Throughfall water volume and chemical constituent inputs were generally lower than those for precipitation water and chemical constituent inputs (reflecting water and chemical retention in the canopy as ice) and followed similar patterns to those described for precipitation volume and chemical constituents ( Table 3).

Canopy water retention and ice accretion
During the first winter of treatments (2016), 45 (Plot 3) to 68 (Plot 2) percent of the water sprayed onto plots was retained by the canopy as ice (Table 1). During the second winter of treatments (2017), 9 (Plot 1) and 42 (Plot 8) percent of the water sprayed onto plots was retained by the canopy as ice ( Table 1).
The two methods for measuring ice accretion were in close agreement for smaller treatments (< 10 mm) with few icicles. Specifically, ice accretion measured by calipers was 111 and 99 percent of ice accretion measured by melted water volume for the Low and Mid treatments, respectively (Fig 4). However, when estimated ice accretion was >10 mm, and icicles were more conspicuous, ice accretion measured by calipers was 75, 85 and 82 percent of ice accretion measured by melted water volume for Mid×2-yr1, Mid×2-yr2 and High treatments, respectively (Fig 4).
In general, the Low, Mid×2-yr1 and Mid×2-yr2 treatments were close to target levels (89 to 99 percent, 87 to 115 percent, and 89 to 104 percent of target values based on caliper and volume estimates, respectively), whereas the Mid and the High treatments were below targets (66 to 67 percent and 70 to 86 percent of target values for caliper and volume estimates, respectively; Fig 4). Across all treatments, measured ice accretion showed a stepwise increase in values from Low-to Mid-to High ice accretion targets, resulting in a suite of experimental light to moderate to severe ice storms (Fig 4).

Fine woody debris
Mean pre-treatment inputs of foliar litter across the treatments ranged from 1.67 to 310 g m -2 for total mass, 0.79 to 147 g m -2 for C, 0.01 to 4.61 g m -2 for N, 1.12 to 454 mg m -2 for P, 0.01 to 3.75 g m -2 for Ca, 1.49 to 659 mg m -2 for Mg and 5.20 to 1201 mg m -2 for K. There were no significant pre-treatment differences in foliar litter mass and nutrients among treatments.
In February 2016, the first ice applications resulted in a general stepwise increase in the FWD from the Control to the High treatment plots (Fig 5). Statistically, the amounts of FWD deposited on the Mid and Mid×2 treatment plots were greater than the amount deposited on the Control treatment plots; and the amount of FWD deposited on the High treatment plots was greater than the amount deposited on the Low and Control treatment plots (Fig 5). In 2017, following the second ice applications on the Mid×2 treatment plots, significantly greater FWD was deposited on the Mid×2 treatment plots compared to all other treatments (Fig 5).
Notably, the FWD deposited on the Mid×2 treatment plots was 2.1 times greater after the 2017 treatment compared to the 2016 treatment (P < 0.01).
The mass of FWD nutrients generally followed the same patterns as for total mass for both years. (Table 4).

Coarse woody debris
Total mass and element content of CWD removed from the plots in November 2015 was highly variable, reflecting the localized history of inputs. Across all treatment plots, mean CWD mass ranged from 475.08 to 1564.05 g m -2 , mean C mass ranged from 225.79 to 750.79 g m -2 , mean N mass ranged from 1.48 to 2.72 g m -2 , mean P mass ranged from 0.07 to 0.12 g m -2 , mean Ca mass ranged from 1.75 to 4.19 g m -2 , mean Mg mass ranged from 0.12 to 0.36 g m -2 , and mean K mass ranged from 0.20 to 0.55 g m -2 . No significant pre-treatment differences were observed among the treatment plots for total mass or any of the elemental masses. Weighted mean decay class across all treatment plots was 3.5, with no significant differences among the treatments.
In July 2016, the first ice applications resulted in a general stepwise increase in CWD mass from the Control to the High treatment plots (Fig 6). Significantly greater CWD mass was deposited on the Mid treatment plots compared to the Control plots, and significantly greater CWD mass was deposited on the High treatment plots compared to the other treatment plots   6). Mean decay classes differed significantly among treatments. Decay classes were lowest (i.e., the wood was least decayed) for the High and Mid×2 treatment plots (weighted mean decay classes = 1.3 for both), mid-range for the Mid treatment plots (weighted mean decay class = 1.7), and highest for the Control and Low treatment plots (weighted mean decay classes = 2.3 and 2.6, respectively). Little additional CWD fell between July and November 2016 (Fig 6), and the sample size (only 12 quadrats total with measurable CWD) was too low for meaningful statistical comparisons. Following the second year of ice applications to the Mid×2 plots, a larger mass of CWD was deposited on the Mid×2 plots compared to the other treatment plots (July 2017, Fig 6). This amount was similar to the CWD mass deposited after the 2016 Mid and Mid×2 ice applications. Mean decay classes differed significantly among the treatments. Decay classes were lower in the High, Mid×2 and Mid treatment plots (weighted mean decay classes = 2.6, 2.1, and 2.4, respectively) compared to the Control and Low treatment plots (weighted mean decay classes = 3.7 and 3.6, respectively). For the Mid×2 treatments, decay classes were greater for CWD collected in July 2017 compared to CWD collected in July 2016 (weighted mean decay class = 2.08 vs 1.27, respectively; P < 0.001). Little additional CWD fell between July and November 2017 (Fig 6). Decay classes differed significantly among the treatments, with mean decay class in Mid and Mid×2 treatment plots (weighted mean decay classes = 2.7 and 2.5, respectively) lower than the Control plots (weighted mean decay class = 3.6). The Mid treatment had a lower decay class than the High treatment plots as well (weighted mean decay classes = 2.7 and 3.6, respectively). Thus, nearly all the CWD resulting from the experimental ice applications fell during and immediately after the applications, the CWD was generally less decayed in the Mid, Mid×2 and High treatment plots compared to the Control and Low treatment plots, and the CWD that fell in 2017 had a greater component of older, more decayed wood than the CWD that fell in 2016.
Total mass of CWD nutrient elements generally followed those of total CWD mass ( Table 5).

Tree canopy damage assessments
For the damage assessments for all trees per plot combined, tree canopy damage indices were not significantly different among the treatments prior to icing (May 2015). In August 2016, the combined mean tree canopy damage indices were significantly greater in the Mid and High treatment plots compared to the Control and Low treatment plots, with the Mid×2 at an intermediate value. In August 2017, the combined mean tree canopy damage indices were significantly greater in the Mid×2 treatment plots compared to the Control plots (Fig 7). In August 2018, the combined mean tree canopy damage indices were significantly greater in the Mid, Mid×2, and High treatment plots compared to the Controls, and by July 2019, the combined mean tree canopy damage indices were significantly higher in all four treatments compared to the Controls (Fig 7).
The combined mean canopy tree damage indices increased in all treatments, including the Control, over the five years of the study, but the Mid, Mid×2 and High treatment plots showed a greater increase in damage over time than the Control and Low treatment plots, which resulted in the significantly higher damage indices for the assessments completed in the summers of 2018 and 2019 (Fig 7).
For individual species, American beech had numerically the lowest initial values for mean canopy tree damage indices (1.2 in the Control treatment to 1.6 in the Mid treatment), followed by yellow birch (1.2 in the Mid×2 treatment to 2.3 in the Mid treatment), red maple (1.3 in the Mid×2 treatment to 2.6 in the Control treatment), and then sugar maple (1.6 in the Mid treatment to 2.6 in the Control treatment) (Fig 7).

PLOS ONE
For American beech, no significant differences were observed in mean tree canopy damage indices among treatments prior to icing (May 2015) (Fig 7). By August 2019, the mean tree canopy damage index for American beech trees in the control plots had not changed whereas values for the indices for the treated plots had increased by 46 to 74%. Initial mean tree damage indices were more variable for yellow birch, sugar maple and red maple than for American beech. By August 2019, mean tree canopy damage indices for the trees in the Control plots for these three species showed only minor changes (0 to +18%), whereas values in treated plots had changed by +38 to +133% for yellow birch, -3 to +67% for sugar maple and +50 to +140% for red maple. The largest numerical increases in mean tree damage indices for each species were in the most severe treatments, with mean indices increasing by 1 unit in the High treatment for American beech, 2 units in the High treatment for yellow birch, 1.5 units in the High treatment for sugar maple, and 1.7 units in the Mid×2 treatment for red maple. Overall, the percent change in mean tree canopy damage indices over the four years of the study was numerically greatest for red maple (in Mid×2 treatment), > yellow birch > American beech > sugar maple; the increases in absolute values of mean tree canopy damage index were numerically greater for yellow birch > red maple (in Mid×2 treatment) > sugar maple > American beech.
Due to small sample size, few of the species year × treatment comparisons were statistically significant (Table 6). Of the 24 possible year × treatment pair-wise comparisons prior to icing (May 2015), 18 showed no significant pairwise differences, five showed tree canopy damage indices greater in sugar maple (2), red maple (2) and yellow birch (1) than American beech, and one showed tree canopy damage index in red maple greater than yellow birch. By August 2019, of the 24 possible year × treatment pair-wise comparisons prior to icing, 20 showed no significant pairwise differences, three showed tree canopy damage indices greater in sugar maple (2) and red maple (1) than American beech, and one showed tree canopy damage index in red maple greater than yellow birch.

Discussion
Our methods and results demonstrate the feasibility of (i) creating a gradient of experimental ice storms and (ii) creating experimental ice storms in consecutive years. The treatment levels were chosen to reflect the amount of ice accretion that triggers National Weather Service ice storm warnings from the Mid-Atlantic to the New England and New York states (6.4 and 12.7 mm radial ice accretion), and then to extend this to a more extreme level (19.1 mm radial ice accretion), and to include an extreme frequency (12.7 mm radial ice accretion in two consecutive years). This initial work establishes the framework to evaluate the response to and recovery from this experimental suite of ice storms for decades to come.

Ice accretion
Radial ice accretion was used as a measure of disturbance intensity. Ice accretion is difficult to characterize and quantify as its buildup varies on surfaces of different materials, diameters, shapes and inclinations [39]. Moreover, measurements in larger storms are dangerous. Pleas for better measurements of ice accretion have been made since the early 1960s [40]. Here, a passive ice collector was deployed to create a standardized, uniform and easily replicated collection surface for ice accretion measurements. Ice accretion measured with calipers produced similar values to measurements of water volume when estimated accretions were < 10 mm, with few icicles visibly present. However, caliper ice measurements underestimated ice accretion compared to the water volume measurements when estimated accretion was > 10 mm and abundant icicle formation was observed. This result highlights the difficulty of measuring ice loading with calipers when icicles are present. Combining both methods probably provided the best metric of icing intensity on these experimental plots. More work is needed to develop simple, low cost measurement protocols for ice accretion to standardize collection surfaces and accretion angles and determine how to account for icicle formation.

Fine and coarse woody debris inputs
The most immediate ecological impact of ice storms on forests is creation of canopy gaps as twigs, branches, tree tops and eventually entire trees break and snap under the load of the accumulated ice. In a related paper on the canopy impacts of this ice experiment, Fahey et al. [28] showed significant post-treatment declines in leaf area index of 27 to 37%, increases in Table 6. Total numbers of significant (P<0.05) pair-wise comparisons for pre-treatment (2015) and post-treatment (2019) tree canopy damage assessments.

Results
Total gap light index of up to 200%, and shifts in the vertical structure of the canopy by the summer of 2017 compared to pre-treatment values. They found that changes in canopy properties were significantly related to both ice accretion as well as the annual total mass of FWD coupled with fall foliar litterfall. Results reported here complement the work of Fahey et al. [28] by quantifying the magnitude and timing of inputs of the FWD and CWD creating the changes in canopy structure following icing. Although these inputs were generally commensurate with treatments (Figs 5 and  6), it is noteworthy that the mass of FWD was significantly greater in the 2017 Mid×2-yr2 treatment compared to both the 2016 Mid and Mid×2-yr1, despite the same target icing (13 mm). Similarly, the mass of CWD was significantly greater and more decayed in the 2017 Mid×2-yr2 treatment compared to both the 2016 Mid and Mid×2-yr1. We propose five nonexclusive explanations for this result: (i) the first year of mid-level icing did not break or damage all the susceptible twigs and branches in the canopy; (ii) the first year of icing broke branches that remained suspended in the canopy until they fell in response to the second year of icing, (iii) the first year of icing weakened or damaged branches that were then more susceptible to breaking during the second year of icing; (iv) the crown loss following the first year of ice treatment opened up access for greater water penetration/ice buildup on remaining twigs and branches, and reduced internal crown support that otherwise could have prevented the excessive bending and additional breakage and loss of canopy components; and (v) for FWD, the amount of time between icing and litter collection was greater in 2017 (70 days) than 2016 (13 to 37 days), although mean and maximum wind speeds between icing and collections were comparable. In any case, these results demonstrate an interactive effect of ice storm frequency and severity that would not be identified without our experimental approach.
Carbon inputs in FWD and CWD, in general, followed mass inputs, with step-wise increases in input commensurate with icing in 2016 and a pulse of C input in the Mid×2-yr2 treatment plots in 2017. The FWD C inputs following the single events in this ice storm experiment ranged from 18% of the average annual fine litter C inputs at HBEF in the Low treatment to 82% of the average annual inputs in the Mid×2-yr2 treatment ( Table 7). The CWD C flux inputs following the single events in this experiment ranged from 16% of the average annual CWD C inputs at HBEF in the Low treatment plots, to approximately equal to the average annual CWD C inputs in the Mid, Midx2-yr1 and Midx2-yr 2 treatment plots, to almost five times greater than the average annual CWD C inputs in the high treatment plots (Table 7). These inputs, especially for the High treatment plots, represent a large transfer of C from the canopy to the soil, and accordingly a large amount of organic debris to be managed by foresters and local municipalities. As woody tissues typically have high C/N ratios and high lignin concentrations [41] these inputs may take years to decades to decompose.
Experimental inputs of FWD C and CWD C can also be compared to inputs incurred at HBEF during the natural ice storm of 1998 (Table 7). This storm severely impacted parts of the south facing watersheds at the HBEF [5,23]. The FWD C inputs in the Low treatment plots of this experiment were roughly similar to those recorded at mid elevations at HBEF (31 vs 29 g C/m 2 , respectively) ( Table 7) after the 1998 ice storm, where the amounts of ice accretions were similar (5.7 mm ice for experiment; 5.9 mm ice for natural storm). The FWD C inputs in the Mid and High treatment plots were approximately two to four times the amount of FWD C inputs at high elevations at HBEF after the 1998 ice storm where ice amounts were roughly similar (8.5 to 16.4 mm ice for experiment and 9.4 mm ice for natural storm) ( Table 7). The CWD C inputs in the experimental plots with roughly equivalent icing to that of the natural ice storm (Midx2-yr1, Midx2-yr2 and High) ranged from 27% to 155% the CWD C recorded after the natural ice storm (Table 7). This current experiment was designed after a pilot ice storm experiment conducted at HBEF in 2012 at a similar elevation and forest type to this experiment [17]. Fine woody debris C inputs for the Mid and Mid×2-yr1 treatments (which had the same target ice accretion as the pilot study (12.7 mm)) were approximately half of the FWD C inputs of the pilot study; FWD C inputs for the Mid×2-yr2 treatments (which has the same target ice accretion as the pilot study but was applied in a second consecutive year) were approximately equal to the FWD C inputs of the pilot study (Table 7). Coarse woody debris C inputs for the Mid and Mid×2-yr1 treatments (which had the same target ice accretion as the pilot study (12.7 mm)) were 64 and 53 percent of the CWD C inputs of the pilot study, respectively; CWD C inputs for the Mid×2-yr2 treatments (which had the same target ice accretion as the pilot study but was applied in a second consecutive year) were 84 percent of the pilot study (Table 7).
Differences in FWD C and CWD C inputs among the several experimental and the natural icing events at HBEF reflect the many factors that interact to influence forest response to ice storms. These include natural variations in species composition, stand age, tree and stand health, soils, and antecedent or postcedent event weather conditions, as well as experimental artifacts such as the more rapid icing (~4 hours) during the experiments versus the natural event (~3 days).
In addition to C inputs, fine and coarse woody debris play an important role in recycling nutrients in northern forest ecosystems. Nitrogen, P and Ca can all be limiting in these forests [43]. Additions of these elements in woody materials replenish forest floor nutrients and, as they are slowly mineralized over time, stimulate future growth. Magnesium and K are also critical for tree health and to repair wounds following injury or insect infestation [44]. The simulated ice storms in this study resulted in significant inputs of N, P, Ca, Mg and K to these forest plots. The nutrient inputs reported here for the combined CWD and FWD measured in the Control (2016 and 2017) and Low treatment plots in 2016 were less than or equal to the amounts reported by Gosz et al. [45]for annual litter inputs in perennial tissues for a mixed hardwood site at a comparable elevation at Hubbard Brook ( Table 8). The nutrient inputs reported here for the Mid, Mid×2-yr1 and High treatment plots in 2016 were 137 to1048% higher than those reported by Gosz et al. (1972); and were 318 to 617% higher for the Values are the mean ± SE. 1 FWD and CWD from Fahey et al. [42]. CWD is sum of coarse litterfall (20 g C m -2 yr -1 ) and coarse woody debris (120 g C m -2 yr -1 ).
Mid×2_yr2 treatment plots in 2017. Thus the nutrient transfers in these single simulated storm events were equivalent to up to 10 years of average annual nutrient inputs at this site.

Tree canopy damage
The results from the tree canopy damage assessments demonstrated that tree response to iceinduced crown injury unfolds over a period of years. For all species combined, the indices of canopy damage increased over time, with significant difference between Mid and High treatment plots vs Control and Low treatment plots only emerging after two to three years following icing (Fig 7). Similar trends were evident for the four individual deciduous species which had representation in all the plots. It is still early in the recovery period, and despite the high levels of icing and the increase in tree canopy damage indices, few trees were characterized as class 4 (> 75% damage). Continued monitoring of these trees will be needed to reveal the full extent of the damage over time. These results are consistent with the evolution of disease and mortality following ice storms as observed by Shortle et al. [21] and Deschenes et al. [46]. Delayed, lagged and/or persistent declines in tree health following ice storm-induced crown damage may be due to a combination of C starvation (due to loss of canopy foliage and increased C requirements for wound closure and defense) and secondary attacks by insects, pests and pathogens [21,47]. Four deciduous species were present in all the treatments and provided limited evidence for species-specific differences in response to icing. For example, red maple and yellow birch trees showed the greatest change in tree canopy damage indices over the four years of the study. In contrast, American beech had the lowest initial tree canopy damage indices for all the species, and showed the smallest change in this index in response to the repeated ice applications. It is notable that these responses are for the overstory American beech which were relatively healthy, with little sign of beech bark disease, a major disease that impacts the health of American beech across its range [48]. Observationally, understory American beech trees were considerably bent over by the weight of the ice, and stayed bowed over for the duration of the study, contributing to the canopy affects described by Fahey et al. [49]. The sugar maple in these stands had the highest initial tree canopy damage indices for all the species, and exhibited only moderate increases in response to the treatments. These observations can be contrasted to those reported by Shortle et al. [21]. In their study of 60-120 year old forest stands across New England impacted by the 1998 ice storm, American beech sustained the most damage, followed by yellow birch and then sugar maple. This discrepancy in the susceptibility of Table 8. Combined input of Hubbard Brook ice storm experiment coarse and initial fine woody debris, and historic input of litter in perennial woody tissues reproduced from Gosz et al. [45]. American beech to icing between the two studies may be due to the relative lack of a pre-existing stressor such as beech bark disease, the shorter amount of time between the icing and the sampling (as impacts evolve over time), and the small sample size in this study compared to the more comprehensive review including a greater range of sites impacted by pre-existing conditions, as presented by Shortle et al. [21]. Despite these patterns over time, few significant pair-wise species differences in tree canopy damage emerged among the four deciduous species. This may be due, in part, to the small sample size for most species, or to the fact that not all species were present in all treatments, which precluded additional species in the pairwise comparisons. Notably, sample sizes for the two coniferous species present in the plots-red spruce and balsam fir-were not robust enough for pair-wise comparisons. Coniferous species are known to be more resilient to ice storms than deciduous trees [50]. In order to answer the driving question of which species or age classes of species are more or less susceptible to ice storm damage, more synoptic studies of tree damage, likely following natural icing events, are needed.

Thresholds
Application of an experimental gradient approach to ice storm damage allowed for the identification of several thresholds for this forest stand. Light icing (i.e., target 6.4 mm ice) caused little twig or branch breakage. Although this amount of icing causes slippery conditions on roads and walkways, it had little impact on FWD and CWD in this forest, as reflected in the relatively low amounts of FWD and CWD in the Low treatment plots, which were not significantly different than the Control plots (Figs 5 and 6). The first thresholds in this forest was observed between 6.4-and 12.7-mm target ice accretions, when noticeable amounts of twigs and small branches were observed to break and fall to the ground. This was reflected in the significantly greater accumulations of FWD in the Mid and Mid×2-yr1 treatment plots compared to the Control plots in 2016, and significantly greater accumulation of FWD in the Mid×2-yr2 treatment plots compared to all other plots in 2017. It also contributed to the small increase in CWD in the Mid and Mid×2-yr1 treatment plots compared to the Control and Low treatment plots in 2016, and to the significant increase in CWD in the Mid×2-yr2 treatment plots compared to all other plots in 2017. A second threshold was observed between 12.7-and 19.1-mm target ice accretion, where whole tree tops bent under the weight of accumulated ice and were beginning to snap at the stem. This observation was reflected in the significant pulse of CWD in the High treatment plots compared to all other treatment plots in 2016. Similar tipping points were reported by Lemon [40] in a synthesis of observations from a suite of ice storms affecting forests in New York State, USA. Note that our experiment was conducted under calm conditions (maximum wind speed 1.8 m/s). Wind accentuates forest canopy ice damage, and the tipping points observed here could occur at lower ice accretions under higher wind conditions.

Hydrologic and chemical artifacts of treatments
In conducting field-scale experiments, potential methodological artifacts that may confound results must be considered. For icing experiments, an artifact can be collateral impacts on forest hydrologic and chemical budgets from the water used to generate ice treatments. Here, the water source for icing was the nearby Hubbard Brook. The chemistry of this stream, studied for over 50 years [51], is characterized as dilute [52] and thus was a good candidate as a water source. Results showed that the amount of water applied to experimental plots was within the range of a large precipitation event at Hubbard Brook (i.e., 5-16% of mean annual precipitation), and similar in magnitude to the water input in the natural ice storm of 1998, with water input to the Low treatment plots being less than water inputs in the natural storm (74%) and water inputs in the High treatment plots being more (254%) ( Table 3). Measured throughfall volume was lower than precipitation and ranged from 35 to 140% of the water input measured for the 1998 ice storm.
Chemical inputs in the treatment were measurable but well below commercial forest fertilization levels [53,54]. Here, total N inputs (NO 3 -+ NH 4 + + DON) for the High treatment represented 7% of the average annual precipitation inputs for the site; total base cations input ranged from 63% (K+) to 300% (Ca 2+ ) and 317% (Mg + ) of the average annual precipitation inputs (Table 3). All nutrient inputs were considerably higher than those in the naturally occurring 1998 ice storm (Table 3), representing a possible experimental artifact that needs to be accounted for in future interpretations of biologic and chemical responses to ice storm treatments.

Conclusions
Currently, a knowledge gap exists on the short-to long-term ecological consequences of extreme weather disturbances. These types of events are problematic to study because they are heterogeneous in time and space, difficult to predict, and disrupt scientific and social infrastructure when they occur, making access to experimental instruments and research sites challenging. Extreme event experiments can help fill this gap, but they can be difficult and even risky to undertake. Initial results from this study provide quantitative measurements of FWD and CWD mass and nutrient transfers from forest canopies to soil in response to different intensities and frequencies of canopy ice accretion. These metrics of response can help foresters, emergency managers, town planners, electric power utilities, arborists, winter recreationists and others be better prepared for the aftermath of ice storms. In particular, results suggest that mature northern hardwood forests such as the one studied here will sustain little damage below 6.4 mm radial ice accretion, moderate damage with up to 12.7 mm of accretion, and significant damage to crowns with upwards of 19.1 mm of ice or more. Further, current ice storm warnings that are issued by the National Weather Service in New York and New England for �12.7 mm radial ice accretion may be too high because significant FWD and some CWD will likely have already fallen at that level of icing and damage will have been done. The 6.4 mm ice storm warnings as issued for the Mid-Atlantic states may provide a safer level to alert the public to pending tree damage from ice storms. Results also suggest that first responders need to be prepared for an equal or greater amount of organic debris inputs if large storms (here �12.7 mm radial ice accretion) occur in back-to-back years. Our results suggest that damaged trees, branches or whole trees that do not fall in the first year are more susceptible to further damage in an ensuing year.
This experiment was the first of its kind in creating a suite of experimental ice storms of different intensities and frequencies in a forest ecosystem. Recommendations for future work include: (i) increase intensity of icing to simulate a tipping point that marks a state change characterized by long-lasting alterations of stand composition or structure; (ii) increase frequency of icing within a single year and over a greater number of consecutive years to better understand how resilient trees are to repeat ice storms and identify tipping points characterized by a state change; (iii) conduct a similar gradient of ice storms in different vegetation types and age classes to better understand what vegetation types and age classes are more or less susceptible to ice damage; (iv) conduct the experiment under different wind speeds, as greater damage is likely with higher wind velocities; and (v) increase plot size to encompass the complete canopies and root systems of a larger number of trees to minimize belowground interactions from adjacent non-iced trees. More work is also recommended to develop simple, low cost measurement protocols for ice accretion, including a refined design for the passive collector that optimizes collector surfaces, angles, and approaches to minimize icicle formation (Campbell et al. in review). Such collectors would allow researchers and citizen scientists to gather more standardized information on the intensity and frequency of icing events across the landscape.
The most interesting results from this experiment are likely those that will evolve over time as this forest ecosystem reorganizes and recovers from this suite of experimental ice storms. This study provides an experimental framework to study the short-and long-term impacts of these experimental ice storms on the structure and function of this northern hardwood forest now [28,55,56] and in the decades to come.