Understanding tropical forest abiotic response to hurricanes using experimental manipulations, field observations, and satellite data

With projected increasing intensity of hurricanes and large uncertainty in the path of forest recovery from hurricanes, studies are needed to understand the fundamental response of forests to canopy opening and debris deposition: the response of the abiotic factors underneath the canopy. Through two manipulative experiments and instrumenting prior to Hurricane Maria (2017) in the Luquillo Experimental Forest (LEF) of Puerto Rico, this study found a long recovery time of primary abiotic factors (beneath canopy light, throughfall, and temperature) influenced by the disturbance of canopy opening, as well as complex responses by the secondary abiotic factors (relative humidity, soil moisture, and leaf saturation) influenced by the disturbance of the primary factors. Recovery took 4–5 years for beneath canopy light, while throughfall recovery took 4–9 years and neither had recovered when Hurricane Maria passed 3 years after the second experiment. Air and soil temperature seemingly recovered quickly from each disturbance (<2.5 years in two experiments for $\sim +\mathrm{1}$ ^∘C of change); however, temperature was the most important modulator of secondary factors, which followed the long-term patterns of the throughfall. While the soil remained wetter and relative humidity in the air stayed lower until recovery, leaves in the litter and canopy were wetter and drier, with evidence that leaves dry out faster in low rainfall and saturate faster in high rainfall after disturbance. Comparison of satellite and field data before and after the 2017 hurricanes showed the utility of satellites in expanding the data coverage, but the muted response of the satellite data suggests they measure dense forest as well as thin forest that is not as disturbed by hurricanes. Thus, quick recovery times recorded by satellites should not be assumed representative of all the forest. Data records spanning the multiple manipulative experiments followed by Hurricane Maria in the LEF provide evidence that intermediate hurricane frequency has the most extreme abiotic response (with evidence on almost all abiotic factors tested) versus infrequent or frequent hurricanes.

1 Introduction

Hurricanes are expected to increase in intensity with climate change (Emanuel, 1987; Knutson et al., 2010; Yoshida et al., 2017); thus understanding how tropical forests respond to hurricanes is critical to understanding future forest regimes. Tropical forests are in a cycle of nonequilibrium, a cycle driven by the response to the large step changes of hurricanes (Burslem et al., 2000). Recently, new tools for understanding the nature and duration of the forest hurricane response have become available for use; satellite data can provide landscape-wide qualities of the historical response (Schwartz et al., 2017), and earth system models can provide the long-term forest response given the projections of the increased frequency of hurricanes (Lee et al., 2018). While these tools can provide a large amount of spatially complete, cost-effective, and consistently recorded data, the data need to be placed in context of what is actually happening at the ecosystem level. There is a need for connection between disturbance and recovery at the critical forest scale: for the manner in which landscape-scale data downscale to the more critical forest landscapes and for the measured response of the forest with repeated hurricanes that should be put into a long-term model (Bustamante et al., 2016; Holm et al., 2014). These connections can only be accomplished with the analysis of fine-scale field observational data. Instead of trying to estimate if and when the vegetation has returned to its predisturbance state, insight on ecosystem health can be gained by studying how the abiotic factors respond to the disturbance. Cascading effects due to canopy openness account for most of the shifts in the forest biota and biotic processes (Shiels et al., 2015), and the biotic environment responds to changes in the abiotic environment.

To this end, a manipulative experiment on hurricane disturbance effects was implemented in the Luquillo Experimental Forest (LEF) in northeastern Puerto Rico, with measurements starting January 2003 and continuing through the time of this paper. The LEF represents a tropical wet montane forest, with high rainfall, high productivity, frequent hurricane disturbance, and semifrequent droughts (González et al., 2013; Scatena and Lugo, 1995; Wang et al., 2003). The forest extends from sea level to 1 km peaks. Droughts occur twice as frequently as hurricanes on the island (every 10 and 21 years respectively) and affect the forest often as dry spells or stronger (Scatena, 1995; Waide et al., 2013). The Canopy Trimming Experiment (CTE) was designed to study the key mechanisms behind a tropical forest's response after a major hurricane and to guide how repeated hurricanes might be expected to alter such ecosystems using these key mechanisms (Richardson et al., 2010; Shiels et al., 2015; Shiels and González, 2014). Multiple control and treated plots were created in the forest. In the treated plots, the forest canopy was trimmed and the canopy debris was littered to the forest floor to simulate the canopy changes from a category 3 hurricane (on the Saffir–Simpson scale).

Two large disturbances occurred during the experiment, both of which were measured by satellites as well as the field instruments. In summer 2015 a drought affected the LEF, starting in May 2015. The forest was still experiencing drought conditions until March 2016 (https://droughtmonitor.unl.edu, last access: 18 June 2020), although precipitation increased after September 2015. On 20 September 2017, category 4 Hurricane Maria made a direct hit on the CTE site. A relatively small amount of disturbance was attributed to the offshore passing of Hurricane Irma 2 weeks earlier; the CTE site was on the lee of Hurricane Irma. The drought and the hurricanes provided data beyond the experimental manipulations on how the forest abiotic environment responds to canopy opening and debris modulated with the variance of climate seasonal cycles and irregularities.

A simplified way of thinking about response to canopy opening and debris deposition is to consider three levels of response. Primary factors are only affected by the initial disturbance: more light and throughfall reach the forest floor and temperatures under the canopy increase. Secondary factors are affected by the primary inputs: relative humidity (in the air), soil moisture, and leaf saturation (wetness of canopy and litter leaves) levels change under the canopy. Tertiary factors are biotic, which are affected by primary and secondary factors, the abiotic factors. Research on biotic effects of hurricane disturbance are numerous (for synthesis efforts see Mitchell, 2013; Shiels et al., 2015), but less researched is how the abiotic factors have changed to alter the biotic environment. This study attempts to quantify abiotic response as acute changes from a hurricane disturbance (experimental or otherwise) and recovery from the changes, for primary and secondary factors. Quantifying the responses makes it possible to assess if the experimental trimming data and satellite data are reasonable sources for studying the effect of hurricane disturbance and appear to be measuring the same abiotic system, as well as appreciate if different events cause substantially different responses. This study does not attempt to determine what amount of recovery is considered normal conditions to biotic life, or in other words what would affect tertiary factors, but instead quantifies changes in the abiotic factors that can be used to frame the changes found in biotic factors posthurricane in many previous studies including those of biotic abundance (Shiels et al., 2015), soil biochemistry (Gutiérrez del Arroyo and Silver, 2018), litter decomposition (González et al., 2014; Lodge et al., 2014), and plant reproduction (Zimmerman et al., 2018).

2 Methods

In spring of 2005 (CTE1) and December of 2014 (CTE2), in 0.09 ha square plots near the El Verde Field Station (419 m; 18^∘20^′ N, 65^∘49^′ W), the forest canopy was trimmed in three treatment plots, and the canopy debris was littered to the forest floor. The plot size and trim amounts were based on the patch disturbance after the two most recent hurricanes before 2017, both category 3 hurricanes at the location of El Verde: Hurricane Hugo in September 1989 and Hurricane Georges in September 1998 (Zimmerman et al., 2014). Non-palm trees of substantial size (>15 cm diameter) had their smaller branches (<10 cm diameter) removed. Smaller non-palm trees and all palm trees were trimmed at 3 m height. All the trimming debris was added to the plot from which it was obtained, with the debris in each plot supplemented with outside debris if necessary, to keep the amounts and kinds of debris equal across the plots.

Biotic and abiotic data were collected in the inner 0.04 ha quadrants of the 0.09 ha trimmed plots to minimize edge effects.

Details of the biotic responses to the 2005 experiment have been extensively documented (Richardson et al., 2010; Shiels et al., 2014, 2015; Shiels and González, 2014), but the abiotic responses (after CTE1 or CTE2) were not fully analyzed until now. Primary factor data (beneath canopy light, throughfall, and air and soil temperature) and secondary factor data (air relative humidity, soil moisture, and leaf saturation) were collected in all plots. To account for spatial heterogeneity under the canopy, multiple sensors and measurements in each plot were used, and the results are calculated off of data averaged from all control and treated plots (with quality control).

While the CTE1 and CTE2 data were being collected, abiotic data were also being collected by satellites and by a nearby weather station. The weather station was located on a tower 30 m above the ground, above control (untrimmed) canopy. After Hurricane Maria, comparisons could be drawn between the experiment and the actual hurricane response, as well as an analysis of which aspects of the response were captured by satellite data, MODIS, and AMSR2. It is important to note that Hurricane Maria provided a much larger hurricane trimming effect than the CTEs were designed to simulate.

2.1 Collecting and homogenizing time series data types

Abiotic field data after 2015 were collected subhourly by automated sensors and averaged into daily values. The abiotic field data before 2015 were collected by different sensors or more intermittent methods (soil and litter gravimetric water contents, GWCs, and canopy photos), so the data had to be converted and calibrated from this first period to the post-2015 period in order to make one time series. Satellite data also had to be converted and calibrated to the post-2015 data type. Specific methods of collection, conversion, and calibration of each data type will be detailed in the following subsections.

Many of the data types required calculation of a smoothed data pattern in order to convert and calibrate. In all cases, the smoothing was done using local estimated scatterplot smoothing (LOESS), which fits least-squares polynomials locally to the points. The LOESS degree of smoothing is contingent on the size of the local neighborhood, which here was always chosen to be 1 year of data around each point. The yearly smoothing was done to extract the larger signal from the data and to homogenize the different collection intervals of the data. The automated sensor field data captured larger amounts of background noise than the temporally smoothed rain funnel data and the geographically smoothed satellite data and, to a lesser extent, the geographically smoothed soil and litter GWCs and canopy photo data. The 1-year smoothing neighborhood was chosen to be longer than the longest length of time between repeat measurements across all data types and methods. No smoothing was done across any of the event dates, in CTE treated, control, or satellite data, regardless of whether the data type was affected by each event. These smoothing breaks were used to keep boundary conditions of the LOESS applications more similar.

2.1.1 Beneath canopy light data

Light beneath the canopy was quantified with solar radiation data. Solar radiation data were collected after 2015 by a Campbell Scientific LI200X pyranometer in each plot measuring 400–1100 nm light from sun plus sky radiation. Earlier estimates of solar radiation were made with sets of hemispherical canopy photos, 10 photos in each set in each plot, which were taken approximately every 4 months from 2005 to 2012. Sets of photos were also taken before the first experiment, and once a year from 2015 to 2017. The solar radiation field data were compared to MODIS Aqua and Terra satellite leaf area index (LAI) data at 500 m, 8 d resolution. The Beer–Lambert law (Monsi, 1953) was used to convert the LAI data into solar radiation estimates, calculating the attenuation that the canopy with a specific LAI invokes on the available (above-canopy) light. Annual patterns of photosynthetically active radiation (PAR) extinction coefficients are needed to calculate the attenuation given by the Beer–Lambert law. An annual pattern of these extinction coefficients was solved for by using 2 years of data of the field-measured CTE2 control plot solar radiation, the tower weather station above-canopy solar radiation, and the MODIS LAI data. The three sets of data were interpolated or averaged to daily values, and then the coefficients were calculated on the 2 years of data before the hurricane (so excluding the 2015 drought). These annual patterns were averaged and smoothed into one annual pattern of extinction coefficients that was applied for every year of the MODIS data.

In the reanalysis of the CTE1 data presented here, canopy photos were converted to global solar radiation data with a modified version of the Hemiphot method (ter Steege, 2018) as follows. Images were converted from color to black and white with a threshold, where the threshold is found iteratively for the best separation of background and foreground using the Ridler and Calvard method (Bachelot, 2016); this method requires calibration. Thresholding was later calibrated to have agreement between annual patterns in the photo solar radiation data and annual patterns in the instrument measured solar radiation data measured in the control plots. Next, the black and white images were converted to canopy openness data by calculating openness on concentric rings of the photo representing sky hemisphere with an arc of 1^∘.

Then, PAR was calculated under the canopy for every day of the year before and after each photo, assuming a constant canopy cover for those time periods. The PAR was then made into one daily time series at each photo site by linearly interpolating PAR each day as a fraction of the previous and the next photo's calculated PAR on that day. This roughly interpolated the canopy cover changes due to recovery from the trimming and interpolated seasonal changes in canopy cover as long as the photos were repeated every winter and summer.

The PAR is the sum of direct and indirect light. The direct light was calculated from the path length of the sun's light through the atmosphere to the forest and the atmospheric transmissivity. Atmospheric transmissivity was given variability around the standard tropical value assuming a linear relationship with relative humidity in the air (Winslow et al., 2001) (as measured above canopy). Path length was calculated from the sun's orbital position on each day of the year relative to the forest. Diffuse light was calculated assuming each part of the sky is equally bright and thus diffuse light is a fraction of direct light. Underneath the canopy, PAR can be approximated as the sum of the direct light through all open parts of the canopy and the diffuse light multiplied by 15 % (based on empirical equations; Gates, 1980). Global solar radiation is then approximated as a multiple of PAR (2X; see Escobedo et al., 2009) calibrated to solar radiation measurements from above the canopy at the tower weather station.

2.1.2 Throughfall data

Throughfall data were collected the entire time period with the same method of biweekly recordings of rain funnels. These funnels were 9.2 cm in diameter, with 1000 mm³ volume. Throughfall was also collected in the treated plots subhourly after 2015 with automatic rain gages. The rain gage data that overlapped the rain funnel data were used to calibrate the rain funnel data.

2.1.3 Temperature data

Temperature data were collected after 2015 by a Decagon Devices VP-3 sensor in each plot in the air 2 m up from the ground and a 5TM sensor in each plot in the soil 5 cm down into the ground. Earlier temperature data were collected hourly by a Campbell Scientific 107 sensor in each plot in the air and one in the soil, underneath the canopy. Air temperature above the canopy at the tower weather station was calculated with the same instrument the entire time period, so annual patterns of the ratios of above-canopy air temperature to below-canopy air and soil temperature were used to calibrate the 107 data. First, the ratios were calculated for 2 years of VP-3 data before the hurricane (so excluding the 2015 drought). These annual patterns were averaged and smoothed into one annual ratio pattern for air and one for soil. Then an air and soil annual ratio pattern was calculated for the complete years of the 107 data (so excluding 2005–2007) and the above-canopy data, and the difference between the ratios was used to make one annual correction each for air and soil that was applied for every year of the 107 data. The air temperature field data were compared directly to MODIS Aqua and Terra satellite land surface temperature (LST) data at 1 km, 8 d resolution. MODIS LST measures energy balance at the land surface, so it is not representative of air temperature under the canopy but it will be affected by changes in air temperature. Annual maximums of LST and air temperature are highly correlated across the globe, with correlation strongest in forested areas (Mildrexler et al., 2011), and LST has been shown to respond to forest cover changes in other areas of the tropics (van Leeuwen et al., 2011).

2.1.4 Air relative humidity data

Air relative humidity data were only collected after 2015. They were collected by the same Decagon Devices VP-3 sensors in each plot (2 m up from the ground) that collected air temperature. Because only one kind of instrument collected these data, no conversion was done on these data.

2.1.5 Soil moisture data

Soil volumetric water content (VWC) data were collected after 2015 by reflectometers, with one Decagon Devices 5TM sensor in each plot measuring shallowly at 5 cm deep and three Campbell Scientific CS616 sensors in each plot collecting profiles from the surface to 15 cm deep. The VWC profile data are comparable to measurements of soil moisture collected by drying out soil samples. Such soil samples were collected for GWCs approximately every 3 months from 2003 to 2006 and in 2015, with five in each plot. Some of these soil GWCs have been published previously before this reanalysis (Richardson et al., 2010). Here, soil GWCs were converted to soil VWCs estimates with measurements of soil bulk density recorded at the same time as the GWCs or using average values from each plot if direct measurements were not available. The 2015 overlap period between the smoothed data of the sensors and the soil sample data was used to calibrate the converted data. The shallow soil VWC field data were compared directly to AMSR2 descending and ascending track satellite soil VWC data, at 10 km, 1 d resolution.

2.1.6 Leaf saturation data

Leaf saturation data were collected after 2015 by three Decagon Devices dielectric leaf wetness sensors in the low canopy leaves in each plot 5 m up from the ground and three in the litter leaf layer in each plot. These sensors have similar thermal mass and radiative properties to real leaves, and wetness is measured by the voltage signal output after voltage excitation, which is higher in proportion to the volume of water on the sensor. This voltage output was then assigned 0 % saturation (dry) at the lowest recorded value and 100 % saturation at the highest recorded value. Earlier measurements of litter saturation were made with leaf GWC values from litterbags, 5 in each control plot and 10 in each treated plot. These litterbags were made of air-dried, preweighed leaves, placed in the litter layer immediately after the CTE1 trimming and retrieved for collection approximately every 3 months from 2005 to 2006. These data were published previously (Richardson et al., 2010). The litterbag procedure was repeated for the CTE2 trimming, and four litterbag measurements of GWC were made in 2015. Leaf GWC is proportional to leaf VWC if the assumption of constant leaf bulk density across plots is made. Then the early litterbag data could be converted to saturation percentages using the ratio between the data of the 2015 litterbags and the smoothed data of the dielectric leaf wetness sensors collected at the same time.

2.2 Quantifying abiotic interaction and response

To explore the relationships between primary and secondary abiotic factors, daily means were correlated. All abiotic data after January 2015 were prewhitened by filtering with an autoregressive integrated moving average model (ARIMA; Box et al., 2015) and first-differencing to remove seasonality and trends. The resulting prewhitened data were examined for correlation between primary and secondary factors for periods with daily data (i.e., after CTE2 and continuing after the hurricanes until 2019).

To explore the differences between responses of different abiotic factors, a smooth time series of each factor was computed, as well as annual averages of data starting after each disturbance event and continuing until the next disturbance event. This was done for CTE data as well as satellite and tower weather station data. For the smooth time series, 1-year LOESS neighborhoods were used to reduce the noise in the data and extract the larger signal with seasonality. For the annual time series, averages of every 365 d after an event were computed (e.g., after Hurricane Maria on 20 September 2017, an average was computed from 21 September 2017 to 20 September 2018). Each yearly mean was visualized as a point at the midpoint of each calendar year (1 July), regardless of the starting date of the average, so that the connected annual time series did not change visually in its seasonal relationship to the smooth time series throughout the series of disturbances. Thus, the first point after an event represents the average of day 1 to day 365 (year one), the second point the average of year two, and so on.

Acute change of each factor was quantified by predefined metrics on the time series. The acute change after the hurricane was defined as the change in the control time series or the satellite time series from right before the hurricane to right after the hurricane, on 20 September 2017. The acute change after an experiment disturbance event was defined as the maximum difference between the treated and control time series (in relation to the control time series) on any day between the last day of the canopy trimming (spring 2005, December 2014) and the next 20 September (year 2005 and 2015, respectively), so that the experimental changes could be compared to the hurricane changes. These changes were calculated on the smoothed data, in order to reduce the noise in the data but still account for short-term changes that yearly means would not capture.

Recovery after a CTE experiment was defined as the point in time after the acute change day that the value of the treated data time series is “close” to the value of the control data time series, and afterwards the difference between the treated and control data stayed “small”. While this definition is qualitatively intuitive, there a several choices that need to be made to enact it in a quantitative manner. The recovery length is sensitive to the choice of metric for initial closeness and for postperiod small differences. The recovery length is also sensitive to the summary choice of the raw data and the length of time the data needs to be close to be called “recovered”.

A sensitivity study was performed to explore the effects of the choices on the reported recovery length. Differences in time series were calculated as fractions of the acute change, or

$$\begin{array}{}\text{(1)}& \mathit{\delta }\left(x\right)={\displaystyle \frac{T^{\ast }\left(x\right)-C^{\ast }\left(x\right)}{C^{\ast }\left(x\right)}}/{\displaystyle \frac{T\left(x_{\mathrm{a}}\right)-C\left(x_{\mathrm{a}}\right)}{C\left(x_{\mathrm{a}}\right)}},\end{array}$$

where x is the day the data were measured on, T(x) is the treated data value, C(x) is the control data value, and x_a is the day the acute change was measured on. The time series of T and C are the smooth LOESS time series used to calculate the acute change, but T^∗ and C^∗ were used as the daily (averaged from raw) data time series, the smooth LOESS time series, and a mean time series on the daily data (using 6-month means after each day x). Then, the day of recovery was defined as the day x_r≥x_a, where

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathit{\delta }\left(x_{\mathrm{r}}\right)\le b_{\mathrm{init}},\text{(3)}& {\displaystyle }\max \left(\mathit{\delta }\right(X_{\mathrm{r}}\left)\right)\le b_{\mathrm{post}},\end{array}$$

using postperiod days as $X_{\mathrm{r}}=x_{\mathrm{r}}+\mathrm{1}$, x_r+2, x_r+3, …x_r+n. In order to capture possible seasonal differences in the postperiod, the ending day for this period was chosen as n=182, or until the next event occurred. Lengths up to n=365 were tested, but the longer period did not make much difference and means the results of the CTE2 (only 2.8 years) have the full postperiod for substantially less time. The recovery time will be reported as time since the start of the experiment, x₀, where x₀ is 1 January 2005, for CTE1, and 1 December 2014, for CTE2.

With the definition set in Eqs. (1)–(3), the entire recovery solution space could be studied in the sensitivity study using buffer values 0 to 1. If both buffers are 1, the recovery time will be calculated as x_r=x_a, and if both the buffers are 0, there will only be recovery if the treated and control time series are identical after day x_r. The entire recovery solution space was calculated for three different data summary choices as discussed above: daily, smooth (LOESS), and mean. For Eq. (2), the mean time series is not used as an option, as it would be the same as using Eq. (3) with the mean time series. Thus, the sensitivity equations can be run for six scenarios. While all scenarios were tested, results here focus on the end-members of the least summarized (daily data time series for both Eqs. 2 and 3), the most summarized (smooth and mean time series in Eqs. 2 and 3, respectively), and the intermediate (smooth data time series for both Eqs. 2 and 3). Based on the sensitivity study, buffers in the final recovery calculations were chosen to be small but not too small.

Note that other studies have defined recovery as the year in which the annual maximum value (of the disturbed area) returns to a previous annual maximum value (assumed representative of undisturbed conditions; Lin et al., 2017). While the method used here is dependent on the size of the amount of data smoothing and the size of the buffers, it is able to make use of the parallelly collected control data to calculate more precise recovery lengths than a year. Furthermore, in a frequently disturbed regime such as the LEF, it is difficult to say what year would be representative of undisturbed conditions. Also, this study does not attempt to quantify the biotic perception of disturbance and just focuses on the abiotic effects, as different biotic species will perceive abiotic disturbance at different sizes of buffers.

3 Results

3.1 Homogenized time series

Some of the secondary factors correlated with the primary factors of solar radiation and temperature, but there were no monotonic relationships found with throughfall. The prewhitened (seasonality and trends removed) air relative humidity correlated well with both the prewhitened primary factors of solar radiation and air temperature ($R^{\mathrm{2}}=-\mathrm{0.67}$) across all periods (after CTE2 and after the hurricanes) and all plots (control and treated). The prewhitened leaf saturation (canopy and litter) correlated somewhat with both the prewhitened primary factors of solar radiation ($R^{\mathrm{2}}=-\mathrm{0.35}$) and air temperature ($R^{\mathrm{2}}=-\mathrm{0.49}$) across all periods and plots. The prewhitened soil moisture (shallow and profiles) did not correlate consistently well with any of the primary factors. All significant correlations were highest at zero lags.

The smoothed time series allowed calculation of more detailed responses than if the analysis had been restricted to only calculations on annual averages. The CTE and satellite acute changes after each disturbance event as calculated from the metrics on the smoothed time series (vertical bars on Figs. 1 and 2 as well as Table 1) are much larger than what was seen with the annual averages of the data (differences in dashed lines on Figs. 1 and 2). Throughfall is the only exception to this, and its acute changes could be accurately summarized with the annual average changes (Fig. 1b).

Figure 1 Primary factor time series, or factors that change due to the initial disturbance changes. Orange vertical lines are the periods of Canopy Trimming Experiment (CTE) 1, CTE2, and hurricanes Irma and Maria (appear as one line), sequentially. Daily values of data are represented with fitted smoothed lines, and thick dotted dashed lines connect annual averages of daily data to aid in visualization of the differences between the time series. Red lines are from treated areas and black lines are from control areas (until the hurricanes) beneath the canopy. Green lines are from the tower weather station above the CTE control canopy, and blue lines are from satellite data. Vertical bars show the acute change after an event for CTE data (gray) and satellite data (blue). The time of recovery from each experiment calculated from the smooth data (if seen) is marked with a gray circle. Plots show (a) solar radiation beneath the canopy, (b) throughfall, (c) air temperature, and (d) soil temperature.

Figure 2 Secondary factor time series, or factors that change because of primary factor changes. All markings are the same as in Fig. 1. Plots show (a) air relative humidity, (b) soil moisture shallow, (c) soil moisture profile, (d) low canopy leaf saturation, and (e) litter leaf saturation.

Table 1 Change after each disturbance event as seen by field instruments and satellites, and recovery time from the Canopy Trimming Experiment 1 and 2 (CTE1 and 2).

^a First column is percentage change from control; second column in italics is absolute change from control. ^b First column is percentage change from before Hurricane Maria; second column italics is absolute change from before Hurricane Maria. ^c First column is recovery time based on smooth data; second and third columns are based on smooth-mean (sm-mn) and daily data, respectively, with units in years

3.2 Sensitivity studies on calculated recovery times

The sensitivity studies reflected the trade-offs in recovery time calculation between extracting the seasonal data signal and capturing the daily variation. The more summarized the time series, the smoother the signals of control and treated data and the sooner post-treatment the abiotic factor was calculated as having small differences for an entire 6 months (and thus satisfied the postperiod buffer size requirement for recovery, Eq. 3). However, with the lesser amount of daily variation in the more summarized time series, the control and treated data do not approach each other on any one day (and thus do not satisfy the initial-closeness buffer size requirement for recovery, Eq. 2) until a long period of the two time series is very similar. Conversely, the less summarized the time series, the more daily variation (or noise from an idealized signal) appears in the reported time series and the opposite situation occurred. Figures 3 and 4 show that when the buffers are large, the daily data sensitivity surfaces plot on top (longer recoveries) and the smooth-mean data sensitivity surfaces plot on the bottom (shorter recoveries). When the initial-closeness buffer gets small, the daily data are quite often not calculated as recovering during the time of the experiment. When the postperiod small difference buffer decreases, the choice of summary for the data in Eq. (2) (the initial closeness) matters more, with the smooth-mean data sensitivity surfaces showing longer recovery than the daily data sensitivity surfaces. The third set of sensitivity surfaces, those from the smooth data, show a reasonable compromise in processing amount versus noise amount.

Figure 3 Sensitivity studies on primary factor calculated recovery times. Calculated years until recovery for each factor are plotted on the vertical axis, with variables of buffers on the other axes. The buffers are the required fraction of the acute change the difference in control and treated plots must approach, on the recovery day for initial closeness (init) and for 6 months after for postperiod small differences (post). Plot column one is for after canopy experiment (CTE) 1 and column 2 is for after CTE2. Surface color represents data summary methods as indicated, and absence of values for specified buffer conditions means the factor did not recover during the time period of the experiment with the specified buffers. Plots show (a) solar radiation beneath the canopy, (b) throughfall, (c) air temperature, and (d) soil temperature.

Figure 4 Sensitivity studies on secondary factor calculated recovery times. All markings are the same as in Fig. 3. Plots show (a) air relative humidity, (b) soil moisture shallow, (c) soil moisture profile, (d) low canopy leaf saturation, and (e) litter leaf saturation.

Because the data noise was seen to overwhelm the recovery calculation on abiotic factors with a smaller acute change signal, buffers in Eqs. (1)–(3) had to be weighted by the size of the acute change relative to the annual (undisturbed) range, or for each abiotic factor,

$$\begin{array}{}\text{(4)}& b=b^{\ast }/\sqrt{\left|{\displaystyle \frac{T\left(x_{\mathrm{a}}\right)-C\left(x_{\mathrm{a}}\right)}{\mathrm{annual}\mathrm{\_}\mathrm{range}}}\right|},\end{array}$$

for both b=b_init and b=b_post. Then, $b^{\ast }=\mathrm{0}$ was used for the Eq. (4) calculation of b_init and $b^{\ast }=\mathrm{0.15}$ was used for the Eq. (4) calculation of b_post. Note that if the data were from an idealized system without noise and not a real system, Eqs. (1)–(3) could have been used with the same b_init and b_post for every abiotic factor. Recovery time calculations from the Eq. (4) buffers are in Table 1.

Some of the shortcomings of the data are apparent in the calculations; these affect the accuracy of the resulting recovery times in Table 1. Firstly, abiotic factors with high geospatial variance may have sizable differences even in a recovered state in colocated plots. This is a known issue with throughfall. It is possible that a recovered state for throughfall is observed 4 years after CTE1 (Fig. 1b; see the temporary convergence in the annual averages and the smooth data) instead of later as the defined recovery metrics calculations in Table 1 report. This can also be seen in Fig. 3b, where the sensitivity surfaces have large planes at z≈4 years, and small buffers either cause a lengthy recovery or no recovery to be calculated. Secondly, missing data will affect the calculations. In the air and soil temperature, the missing data in CTE1 (Fig. 1c, d) means the daily data calculate a recovered state (Fig. 3c, d, Table 1) if the initial-closeness buffer is large enough and there is little to no data for 6 months after the recovery day. Thirdly, abiotic data with high amounts of noise will be affected more by the summarization methods. The high daily variance in the canopy saturation data can be seen with the results in Fig. 4d, with the daily data differences never even reaching 100 % of the acute change for the postperiod small difference maximum during the experiment (the daily data sensitivity surface does not exist because Eq. 3 cannot be satisfied).

Despite these shortcomings, the reported recovery times in Table 1 show evidence of being good estimates of recovery times that consider data seasonality instead of just using annual summary data estimates. The CTE annual averages showed convergence in their time series (dashed lines on Figs. 1 and 2) in approximately the same time as the recovery metrics (gray circles on Figs. 1 and 2) but necessarily round up to the next year (or longer if the recovery is late in the year). For example, solar radiation after CTE1 recovered in a calculated 4–5 years (Table 1), or in year five, and the black and red dashed lines in Fig. 1a cross at the fifth point. Throughfall after CTE1 is the most severe exception. The annual averages support recovery in 3–4 years (Fig. 1b), but the calculations for recovery resulted in 9 years (Table 1).

3.3 Comparing experimental manipulations, field observations, and satellite data

The passage of Hurricane Maria, 2.8 years after the second experiment, happened when most of the abiotic factors had not recovered and a few had just recovered. Temperatures after CTE1 and CTE2 recovered in around 2 years, almost half the time it took solar radiation to recover after CTE1 and a less than a third of the time it took throughfall to recover (Table 1). The effect of Hurricane Maria was smaller on the treated plots than the control plots, such that the absolute level of abiotic disturbance on the treated plots was smaller than on the control plots (Figs. 1 and 2). It is expected that the abiotic fluctuations from the hurricane would be smaller in the unrecovered treated plots than in the control plots since there is less vegetation to disturb. The fluctuation is smaller, but furthermore, for most of the abiotic factors, the treated plots are closer to the recovered state after the hurricane than the control plots are. For example, there is more solar radiation reaching the forest floor in the treated plots than in the control plots before Hurricane Maria, but after the hurricane there is less solar radiation reaching the forest floor in the treated plots than in the control plots (Fig. 1a). The same scenario can be seen in the throughfall (Fig. 1b), the temperatures to a lesser extent (Fig. 1c, d), the soil moisture profile (Fig. 2c), and the litter saturation (Fig. 2e). The air relative humidity has the opposite scenario, showing treated plots closer to the recovery state of less humidity in the air after the hurricane (Fig. 2a).

Overall, the patterns of acute changes across the abiotic factors from the experiments and Hurricane Maria are similar (Table 1). With only 0.09 ha plots, edge effects of the nondisturbed forest were expected to lessen the effectiveness of the experiments in simulating hurricane disturbance; yet, the acute changes showed that CTE2 was the most immediately disruptive event across the abiotic factors, more so than the hurricane. The soil moisture increased much more in the treated plots of CTE1 and CTE2 than in the acute change calculated from before to after the hurricane. However, it is impossible to know the true control (no-hurricane scenario) soil moisture level posthurricane.

The satellite data have somewhat similar characteristics to the field data in the control plots (blue vs. black lines in Figs. 1a, c and 2b) in that the magnitude of the acute change is similar (Table 1) and the responses to the summer 2015 drought and Hurricane Maria are in the same direction. Before the hurricane, the (MODIS-LAI-estimated) solar radiation satellite data look very similar to the field data, but they show a smaller change after the hurricane (Table 1) and faster recovery down to previous values than did the field data (Fig. 1a). The LEF lost 51 % of the initial greenness in Hurricane Maria, but the US Caribbean overall lost 31 % of its initial greenness (Van Beusekom et al., 2018), so for the hurricane disturbance, including the area outside the forest would be expected to dampen the measurement of the LAI hurricane disturbance signal. The (MODIS-LST-estimated) temperature satellite data plot between the field air temperature data measured below the canopy and that measured above the control canopy at the tower weather station (black and green lines respectively, Fig. 1c), as might be expected from a LST representative of surface energy balance. These data were strongly affected by the hurricane and quick to recover. The (AMSR2-estimated) shallow soil moisture satellite data have very large spatial smoothing (10 km resolution, containing nonforest and thin-forest areas), showing a drier soil than the CTE. These data were also strongly affected by the hurricane and appear to recover quickly (Fig. 2b).

4 Discussion

The responses in the CTE plots from after Hurricane Maria were very similar to the responses after the two trimming events, which was the aim; nonetheless, it is encouraging how well the experiments worked. However, lacking a control plot for the actual hurricane response, the differences in the seasonal timing of the experiment treatments and Hurricane Maria, as well as sensitivity of the calculations of actual hurricane effects to the data smoothing, make direct comparison of acute changes from the experimental events and actual hurricane disturbances challenging. For these reasons, the quantification of the acute changes in the experimental setup is most useful as a measure of the effect of a hurricane on the abiotic environment, while the quantification of the acute changes from the actual hurricane serves best as a comparison between the field and satellite data, as well as between the relative effects on each abiotic factor for the CTE and the hurricane.

The smaller effect of Hurricane Maria on the treated plots provides evidence that when frequent hurricanes happen, the forests will exhibit abiotic resilience, and thus possibly forests with an intermediate hurricane frequency will have larger abiotic fluctuations due to disturbance than forests with infrequent or frequent hurricanes. Supporting evidence has been found also in the biotic factors of the forest after Hurricane Maria (Hogan et al., 2018), with analysis suggesting tree demographics (the rates of species and stem mortality and growth) were the most dynamic in areas which had the chance to grow some (but not all) trees past the pioneer stage. Intermediate disturbance has long been suggested to keep systems as far from equilibrium states as possible, with the important effect of driving ecosystem diversity (Connell, 1978). Frequent disturbance in the LEF could be regarded as less than a decade (because abiotic factors have not recovered in this time frame), with intermediate frequency longer than a decade but still less than the 60-year long-term return interval for hurricane disturbance in this forest. The 60 year time frame has been estimated to be a long enough period to achieve steady state for time length of biomass turnover (Scatena, 1995), so disturbances, for this system, might be considered infrequent if they happen less than every 60 years.

It is well known that there are issues of scale when comparing outright the values of large-pixel satellite observations to point field observations (Wu and Li, 2009), but the faster recovery seen in the satellite data is interesting. The MODIS LAI may be measuring some low vegetation that grows back rapidly and not recovering canopy, thus decreasing the estimated satellite solar radiation back to undisturbed values more quickly than seen in field observations. The (MODIS-LST-estimated) temperature satellite data and the (AMSR2-estimated) shallow soil moisture satellite data may have had large acute changes and quick recoveries because they are measuring more than just the forest. Above-canopy temperatures are included in the energy balance LST data, and low-permeability areas that flood and dry out are included in the AMSR2 data.

Two of the primary factors, light and water, changed dramatically after the disturbance events (Table 1, Fig. 1a, b). Across the three events, the range of the percentage change in understory solar radiation after disturbance was quite large (214 % to 919 %); it is likely that a sizeable portion of the range is due to the different seasonal timing of the events. The 1998 Hurricane Georges was estimated to have changed the forest light by almost 400 % (Comita et al., 2009), which is within the range seen here. The response of reduced understory light and throughfall (Table 1) was found here to last much longer than the 18 months concluded previously (Richardson et al., 2010). However, it was noted in a related study (Shiels et al., 2010) that the control plot understory solar radiation appeared to be still recovering from the 1998 Hurricane Georges in measurements made in 2003, 5 years after Hurricane Georges. The recovery times calculated here support a continuing recovery from Hurricane Georges in 2003, as the recovery of solar radiation is estimated at 4–5 years (Table 1) and the recovery of throughfall is estimated at 4 least years by the sensitivity study of the quantitative metrics (Fig. 3b). This study had additional information from the second experimental trimming, as well as a longer record of analyzed data from the first trimming and new methods to make a more-continuous record from the intermittent field data. The response may appear in the drier darker season as being recovered (e.g., January 2008, 3 years post-trimming), but it is clear with the longer record that the response is slower to recover. Temperatures of air and soil were much more robust in respect to the changes from the events versus their annual seasonal cycle changes, with approximately 3 % air and 6 % soil acute increases on average, or +0.7 ^∘C air and +1.4 ^∘C soil, recovered a best estimate of 2–2.5 years (Table 1, Fig. 1c, d). But these changes may still be significant to biotic factors. Other studies show that gross primary productivity of the forest is highly sensitive to small increases in air temperature greatly increasing canopy temperature (Pau et al., 2018), so this change that is amplified in the hottest parts of the year (Fig. 1c, d) should not be discounted. As hurricane intensity is expected increase with climate change, there could be a compounding effect of hurricanes and global warming in the future.

Abiotic factors that change because of primary factor changes, or secondary factors, have more complicated recovery paths than the primary factors. Specific timelines for recovery would be expected to be highly influenced by the tree species and soil types, and the rates seen here for all abiotic factors would not necessarily apply to all hurricane-affected tropical forests. Nevertheless, general patterns might be expected to hold. All the secondary factors were clearly affected by the summer 2015 drought and subsequent long-term rainfall levels, as seen by the large magnitude decreases in summer 2015 and the recovery afterwards in air relative humidity, soil moisture, and leaf saturation (Fig. 2). However, daily patterns of the relative humidity in the air and leaf saturation under the canopy were significantly influenced by the temperature and light inputs (based on the results of the residual correlations), while soil moisture may not be influenced much by these inputs. The soil moisture and litter saturation responses from the first trimming present different conclusions when analyzed along with the nearly continuous in situ measurements after the second trimming. Previous studies found very quick recovery of these factors, 3 months and 18 months, respectively (Richardson et al., 2010). However, reanalysis of the data after the first experimental trimming – separating the data into control and treated plots; calculating volume-based percentages of water in the soil and litter instead of mass-based percentages; and, most significantly, looking at the trimonthly collected data from CTE1 in light of the nearly continuously collected data from CTE2 – led this study to draw different conclusions of most likely longer recovery times (Table 1).

The soil moisture increases in all three trimming events (including the hurricane), but the magnitude of the acute change and the time till recovery appears highly dependent on the amount of rainfall (Table 1, Fig. 2b, c). Differences between treated and control sites appear pronounced in dry periods (e.g., spring 2006 and summer 2015), with wet periods obscuring the differences in the sites when the soil may be approaching saturation (e.g., summer 2006). However, the recovery process happens mostly monotonically in the smoothed time series, and the soil profile may be near or at recovery at 2.8 years (Table 1). Soil moisture is higher after disturbance because there is more throughfall and less transpiration (no leaves), but once the leaf area starts to recover the soil moisture recovers quickly.

Conversely, during the dry periods the differences between treated and control sites are obscured for the leaf saturation data. The litter leaves in the second trimming were measured to be wetter and drier following the trim and not uniformly drier as concluded previously (Fig. 2d, e). Data from the second trimming and Hurricane Maria shows that the litter was more saturated immediately following the events, and the low canopy leaves were drier. During periods of low rainfall, the treated plots dry out faster than the control plots, in both litter and low canopy leaves. Sometimes this results in the leaf saturation being lower in the treated plots than the control plots (e.g., summer 2015 and spring 2017). When the rainfall increases after a dry period, during the late-summer rainfall, the treated plot leaf saturation increases much faster than the control plots, suggesting the long-term effect of disturbance on leaf saturation is a more dramatic modulation in saturation by rainfall. Other studies in completely different ecosystems, in the southeastern United States, have seen that litter is able to become more saturated after large storms than before the storms, and they attribute this to the addition of new debris being able to hold on to more water (Van Stan et al., 2017). The litter saturation data from the first experimental trimming (data from 2005 to 2007) do not contradict this conclusion, but due to their record length and collection interval (trimonthly) they are not overly conclusive.

The results do not support a longer or shorter recovery time interval for the second treatment, 10 years after the first (Table 1). The results showed that quantifying recovery times using subannually summarized time series to homogenize data from several sources was a worthwhile effort, in that the abiotic factors can be sorted into quicker and slower recoveries in subannual lengths. However, the definition of the “recovered point” in time will be dependent on what biotic life considers “normal”, which is necessarily different for every organism. The recovery times presented here for different abiotic factors are a starting point for other researchers to frame the changes found in biotic factors posthurricane. However, the trade-offs in recovery time calculation between extracting the seasonal data signal and capturing the daily variation, and the influences of data noise and variability, point to the difficulty of quantifying recovery in an environmental system.

Climate projections predict Puerto Rico air temperature will be +2 ^∘C warmer in the coming century and rainfall will be −20 % to −30 % smaller in the fall and summer wet seasons (Hall et al., 2013; Karmalkar et al., 2013). Effects from future hurricanes on the abiotic factors will be on top of this background change. This means a hurricane could add an acute effect of almost 50 % more to the temperature increase, with a recovery of over 2 years (Table 1). The throughfall after a hurricane was found to increase >100 % with a long recovery of up to 9 years (Table 1). But, given the climate projections of more events like the summer 2015 drought, the more noteworthy effect of future hurricanes may be the litter and low canopy leaves drying out much faster in the drought and saturating faster with rain after the drought. This will create a much more dynamic environment of leaf wetness, which may have implications for biotic factors.

5 Conclusions

The way abiotic characteristics are disturbed and the speed at which they recover will be key to the continued existence of tropical forests under a climate with more intense hurricane activity. Climate projections predict changes that will exacerbate the effects of hurricanes of increasing temperature and dynamically changing leaf wetness. There is evidence here that intermediate hurricane frequency will have the most extreme abiotic response (with evidence on almost all abiotic factors tested) versus infrequent or frequent hurricanes and that satellite data may show a faster recovery than field data looking at canopy response and soil moisture. Caution must be exercised when declaring the recovered point of a forest, as full abiotic canopy closing may take half a decade or longer and not all abiotic factors recover monotonically. Abiotic factor responses to hurricanes are not included in current climate projections. Results from detailed manipulative experiments such as this study are needed in order to begin to quantify abiotic factor responses to hurricanes to add to the climate projections.