2.1. Site description

The eddy covariance (EC) tower is located at a Florida Coastal Everglades Long Term Ecological Research site (SRS6; 25.3646°N, 81.0779°W) and adjacent to a US Geological Service monitoring site (SH3) near the mouth of the Shark River, one of the principal drainages of Everglades National Park (ENP) into the Gulf of Mexico (Fig. 1). This subtropical region is characterized by distinct wet (June-October) and dry (November-May) seasons during most years. At least 60% of annual rainfall occurs during the wet season. Water levels in the coastal Everglades are relatively higher (>0.3 m) during the wet season due to increased fresh water discharge, higher rainfall, and because the peak of the annual tidal cycle also occurs during this period (Stumpf and Haines, 1998). The early part of the dry season is characterized by lower and more variable air temperatures (T_A), and periodic rain events associated with the infrequent passage of cold fronts. Annual minimum freshwater discharge rates and annual maximum salinity values occur at the end of the dry season in April and May. Monthly values of environmental variables at the site from 2004 to 2009 are provided in Fig. 2A-D.

At the study site the dominant mangrove species include Rhizophora mangle, Avicennia germinans, and Laguncularia racemosa, and their maximum heights can reach 19 m. Prior to the disturbance caused by Hurricane Wilma, Chen and Twilley (1999) reported a stand density of 7450 trees ha^-1. The forest understory is sparse and comprised of seedlings and juvenile mangroves with an average height of less than 4 m. The region is characterized by mixed semidiurnal tides. The site is flooded 4622 h yr-^-1 (52% of the time) during high tides when surface water depths can reach 0.5 m (Krauss et al., 2006). Peat thickness beneath the forest increases towards the Gulf of Mexico, and at the SRS6 research site a limestone substrate lies between 5 and 6 m beneath the organic soils (Spackman et al., 1966).

Mangrove forests along the Gulf Coast of south Florida are subject to tropical storm disturbance every 3-5 years (Doyle and Girod, 1997; Hopkinson et al., 2008). Hurricane Wilma made landfall as a category 3 storm with an eye over 128-km wide centered over the SRS6 study site. The storm deposited calcitic marl from the Gulf of Mexico and Florida Bay across a 70-km zone along the western border of ENP. Between 5 and 15 cm of these deposits were found inland at a distance from 1.6 to ~5 km from the coast (Smith et al., 2009).

2.2. Meteorological, eddy covariance and energy balance measurements

Fig. 2. (A) Monthly air temperature (°C) at a height of 27 m, (B) pore-water salinity (parts per thousand; ppt), (C) solar irradiance (W m-2 day-1), and (D) total rainfall (mm) before and after Hurricane Wilma. Box plots include 25th, median (50th), and 75th percentiles of half-hourly values during the month, and whiskers indicate 5th and 95th percentiles. [larger image]

Meteorological and EC measurements of NEE, LE (W m^-2), and H (W m^-2) prior to Hurricane Wilma's impact follow the methods described in Barr et al. (2010). Instruments were deployed on 30-m flux tower, which was constructed in June 2003. Environmental variables (incoming and reflected solar irradiance, net radiation, incoming and reflected photosynthetic active radiation (PAR), and wind speed and direction, air temperature, soil temperature, and soil heat fluxes) were measured at 1-s intervals and averaged over 30-min periods. Hydrologic data (water depth and specific conductivity) were recorded every 15 min at a station 30 m south of Shark River and 150 m west of the flux tower. The eddy covariance (EC) system was mounted at 27 m above the surface. The EC consists of a three-dimensional sonic anemometer (model RS-50, Gill Co., Lymington, England) and an open path infrared CO₂ and water vapor (H₂O) gas analyzer (model LI-7500, LI-COR, Inc., Lincoln, Nebraska). High-frequency (10 Hz) measurements were recorded and processed with custom software to derive half-hourly CO₂, latent and sensible heat, and momentum exchanges between the forest and the overlying atmosphere. Spike removal from raw turbulence data was achieved using a two-dimensional coordinate rotation of the wind field (Vickers and Mahrt, 1997). Time lag corrections were made to maximize of CO₂ covariance with vertical wind speed variation. Buoyancy corrections of sonic air temperatures were also introduced (Schotanus et al., 1983). Density variations arising from water vapor were included the estimation of the total constituent flux (Webb et al., 1980). Storage of CO₂ in the air column below the EC system was estimated based on the half-hourly rate of change of CO₂ concentrations at the infrared analyzer level (e.g., Humphreys et al., 2005). The flux data included in this study are of high quality as verified from inter-comparisons with the Ameriflux algorithms (http://public.ornl.gov/ameriflux/standardsgold.shtml). Barr (2005) and Barr et al. (2010) provide further details about field observations and data processing.

We follow the sign convention for land surface-atmosphere scalar fluxes derived from EC data such that negative values refer to fluxes from the atmosphere to the land surface, and positive values refer to fluxes from the land surface to the atmosphere. Estimates of NEE are expressed in µmol (CO₂) m^-2 s^-1, while integrated values over monthly or annual periods are given in g C m^-2 t^-1 and represent -NEE values. Average values ±1 standard deviation are reported separately for daytime and nighttime periods and further partitioned into dry (November-May) and wet (June-October) seasons. Hurricane Wilma destroyed the tower structure, instruments, and electronics. The site was rebuilt and the same measurements described in Barr et al. (2010) made during 2004-2005 were restarted in November 2006. Soil thermocouples were identical in positioning and model to those deployed in 2004-2005. New aspirated and shielded platinum temperature probes (model 41342VC, R.M. Young Co., Traverse City, MI) were used to measure T_A at heights of 27 m, 20 m, and 1.5 m after the storm.

2.3. Gap-filling missing data

Missing or invalid EC fluxes, referred to as 'gaps', occur during power failures, when gas concentrations are out-of-range as occurs during rainfall events (due to raindrops resting on the mirror of the LI-7500 instrument), when turbulence is weak or intermittent (the u* threshold), or when there is insufficient fetch. Short duration gaps (≤4.5 h) in the EC data occurred primarily at night during low turbulence conditions. The u*-threshold below which flux data were excluded (Goulden et al., 1996; Lee et al., 1999) was calculated by first dividing nighttime NEE values into 20 u* classes for each bi-monthly period, and then defining a u* value above which NEE became invariant. For those bi-monthly periods where no clear relationship between NEE and u* existed, u* values were chosen to correspond to an NEE not less than 85% of the maximum bi-monthly NEE (Barr et al., 2010). Median value of all bi-monthly u* threshold quantities calculated before Hurricane Wilma was 0.21 m s^-1, and varied between 0.15 m s^-1 and 0.30 m s^-1. During three bimonthly periods prior to Hurricane Wilma, the u* threshold was >0.25 m s^-1. However, the differences in fluxes calculated during these periods using a u* threshold of 0.21 m s^-1 versus greater values (up to 0.3 m s^-1) were not significant. Therefore, a global u*-threshold of 0.21 m s^-1 was applied to all bi-monthly periods. The same u* threshold analysis was performed on nighttime flux data from 2006 to 2009 collected after the storm. For those periods, flux calculations made during periods when u* values were below the post-storm, bi-monthly average threshold of 0.14 m s^-1 were excluded. Flux data were also excluded when the flux footprint (Schuepp et al., 1990; Schmid, 2002) extended beyond the forest fetch. The footprint exceeded the fetch most frequently during the nighttime with wind direction from northwest to northeast.

Gaps in the dataset exceeding 4.5 h were generally caused by instrument or data acquisition malfunction. Prior to the hurricane, the combined duration of gaps comprised 61%, 28%, and 46% of the total nighttime, daytime, and combined datasets, respectively. Following the storm in 2007, gaps comprised 67%, 49%, and 58% of the total nighttime, daytime, and combined datasets, respectively. In 2008, gaps comprised 56%, 25%, and 41% of total nighttime, daytime, and combined datasets, and in 2009 gaps comprised 51%, 15%, and 34% of the nighttime, daytime, and combined datasets, respectively.

A mean diurnal variation (MDV) method was used to fill short gaps and look-up tables (LUT) were employed to fill longer gaps. The MDV method is used to fill individual gaps using the mean fluxes occurring during the same half-hourly period within a 14-day window centered on the day of the gap. For longer gaps, separate daytime and nighttime LUTs were developed for each 2-month interval beginning on January 1, 2004. Nighttime T_A was better correlated with CO₂ fluxes than T_S and was therefore chosen as the independent variable in the LUT. For each 2-month interval, halfhourly nighttime CO₂ fluxes were partitioned into 20 T_A bins, each containing the same number of values. For daytime half-hourly CO₂ fluxes, a two-dimensional LUT was constructed using 16-PAR and 3-T_A bin categories (Falge et al., 2001). Potential error and bias in the fluxes introduced by the MDV and the LUT were estimated by randomly creating and then re-filling 100 sets of artificial gaps (Moffat et al., 2007) overlapping with valid data periods for each bimonthly period. Barr et al. (2010) demonstrated that the root mean square error (RMSE; 3.56 ± 0.058 µmol (CO₂) m^-2 s^-1) and the bias error (BE; -0.021 ± 0.054 µmol (CO₂) m^-2 s^-1) associated with this technique imputed minimal error in annual NEE estimates. Gapfilled, half-hourly NEE data were integrated to produce monthly and annual total values for the period January 2004 to December 2009. Gap-filled EC data were also integrated during daytime periods only to produce LE and H estimates over daily, monthly, and annual timescales. MDV and LUT methods were used to fill short (≤4.5 h) and long (>4.5 h) duration gaps, respectively in LE and H estimates. Only daytime gaps were filled, and the two-dimensional LUT was constructed using 16 solar irradiance bins and 3 vapor pressure deficit bins. MDV and LUT methods were also used to fill gaps in net radiation during daytime periods, and the LUT was constructed using 16 solar irradiance bins and 3 air temperature bins.

2.4. Remote sensing, stand surveys, and sediment elevations

Values of the Enhanced Vegetation Index (EVI) from the Moderate Resolution Imaging Spectro-radiometer (MODIS) aboard NASA's Terra satellite were used to examine the magnitude and duration of Hurricane Wilma's impact on canopy reflectance and light interception. EVI is particularly sensitive to the vegetation signal, resists saturation in high biomass regions, provides for de-coupling of background and vegetation signals, and reduces signals attributed to atmospheric conditions (Huete et al., 2002). EVI values were derived from the MOD13A1 data product (EOS; http://modis.gsfc.nasa.gov/). The EC study site is included in grid h10v06, with a 500-m spatial resolution. Using GIS (Geographic Information System) software (Matlab Mapping Toolbox, The Mathworks Inc., Natick, MA), the 16-day average EVI values for the pixel corresponding to the location of the study site and the 8 adjacent pixels were extracted for the period from 2000 to 2009. This 9-pixel domain approximates the extent of the EC measurement footprint (see Barr et al., 2010; Fig. 1). Each 16-day composite EVI value consists of up to 64 observations. However, typically fewer than 10 observations were included in a composite value after exclusion of values during cloudy and off-nadir periods. Minimum and maximum EVI values range between 0 and 1.

Surveys of tree growth and mortality, and soil elevations at three sites in the impact zone were initiated several years prior to Hurricane Wilma. The least impacted site on the Lostmans River (LO3) is located 22.1 km from the EC site (Fig. 1). The SH3 site, located 150 m from the EC tower, experienced greater impacts and approximately 100% defoliation. The Big Sable Creek (BSC) site lies 13.7 km from the EC tower and is 400 m from the Gulf of Mexico. This site experienced more exposure to wind and storm surge during Hurricane Wilma compared to the other sites. Stand surveys began at these sites following Hurricane Andrew in September 1992. Circular plots of variable radii were established in which all stems >1.5 m in height were tagged and measured for diameter at breast height (1.4 m). The first surveys following Hurricane Wilma took place in November 2005. Additional surveys of all stems >1.5 m have been carried out intermittently since the disturbance. Sediment elevation table (SET) benchmarks are located at the center of each site. Quarterly surveys of changes in soil surface elevation, with a total error of ±1.4 mm began in March 2002 (Whelan et al., 2005). Relative elevation changes caused by Hurricane Wilma were determined by subtracting the last elevation measurement made before the storm from all subsequent measurements.

2.5. Statistical analyses

A multivariate ridge regression model (Hoerl and Kennard, 1970) was applied to distinguish the direct effects of disturbance on monthly daytime, nighttime, and daytime plus nighttime NEE from the independent, inter-annual variation in these quantities due to non-hurricane environmental drivers. The ridge regression employs the following relationship:

= (A^TA + Γ) A^Tb       (1)

In (1), is the set of regression coefficients that minimizes the sum of squared error of predictions of the response variable b (in this case monthly NEE). The coefficient matrix, A, consists of the explanatory variables indexed by columns and sequential monthly replicates indexed by rows. Explanatory variables were monthly median T_A at 27 m (°C; Fig. 2A), median soil temperature T_S at -5 cm (°C), monthly sum of solar irradiance (MJ m^-2), median surface water salinity (parts per thousand; Fig. 2B), the fraction of time the site was inundated during each month, and monthly rainfall (mm; Fig. 2D). Monthly rainfall was based on the mean value recorded at three nearby monitoring stations located at a distance of ~3 km, ~5 km, and ~13 km from the tower. Ordinary least squares regression (OLSR) minimizes the sum of squared errors of the predicted values. The minimization only applies to data sets used to train the model. However, when applied to new data sets, the resulting coefficients do not necessarily provide a model which minimizes the sum of squared errors. To overcome this problem, regularization (Γ) of the coefficients in the vector is performed (Hoerl and Kennard, 1970). In this procedure, Γ is defined as the product of α and the identity matrix, I (i.e., Γ = αI). Typically, α is determined numerically by selecting the value that minimizes the prediction of MSE from a validation data set.

Three different models were built with the response variables of observed, monthly integrated daytime (-NEE_day,obs), nighttime (NEE_nt,obs), and 24-h total NEE (-NEE_tot,obs). The final value of α in each model was selected to minimize MSE using values from 2004 to 2005. A leave-one-out (LOO) cross validation procedure (Picard and Cook, 1984) was used to compute MSE for each value of α. Data on environmental variables from the years 2006 to 2009 were applied to the NEE models to produce monthly estimates of daytime (-NEE_day,mod), nighttime (NEE_nt,mod), and 24-h total NEE (-NEE_tot,mod for post-storm periods. The differences between modeled and observed monthly NEE values provide estimates of the changes in CO₂ fluxes resulting solely from the disturbance. The average differences in modeled and measured monthly values were compared to zero using a one-sided t-test for the periods 2006-2007, 2008, and 2009. A rejection of the null hypothesis that the mean monthly pre- and post-storm difference in NEE was equal to zero at the 95% confidence level was interpreted as a significant disturbance effect.