On the difference in the net ecosystem exchange of CO2 between deciduous and evergreen forests in the southeastern United States Authors: Kimberly A. Novick, A. Christopher Oishi, Eric J. Ward, Mario Siqueira, Jehn-Yih Juang, and Paul Stoy This is the peer reviewed version of the following article: see full citation below, which has been published in Global Change Biology which can be found in its final form at http:// dx.doi.org/10.1111/gcb.12723. This article may be used for non-commercial purposes in accordance with Wiley Terms and Conditions for Self-Archiving. [If there is a citation, include it on a line by itself, if grey literature, include the full date of this version (July 4, 2014). include a doi if possible] Made available through Montana State University’s ScholarWorks scholarworks.montana.edu http://scholarworks.montana.edu/ http://scholarworks.montana.edu/ http://dx.doi.org/10.1111/gcb.12723 http://dx.doi.org/10.1111/gcb.12723 https://onlinelibrary.wiley.com/journal/13652486 c o n s u m pti o n ( B ol st a d & S w a n k, 1 9 9 7; Oi s hi et al., 2 0 1 0; F o r d et al., 2 0 1 1). T h e y a r e al s o i m p o rt a nt si n k s f o r a nt h r o p o g e ni c e mi s si o n s of C O 2 ( B r o w n & S c h r o e d e r, 1 9 9 9; Al b a ni et al., 2 0 0 6). W h e n s u r v e yi n g p r e vi o u sl y p u bli s h e d e sti m at e s of t h e n et e c o s y st e m e x c h a n g e of C O 2 ( N E E) f r o m c o nti n e nt al U S e c o s y st e m s, it i s cl e a r t h at s o ut h e a st e r n U S f o r e st s a r e a m o n g t h e st r o n g e st c a r b o n si n k s ( Fi g. 1, T a bl e S 1). F o r i n st a n c e, i n s e v e r al s o ut h e a st e r n U S f o r e st s, t h e mi ni m u m o b s e r v e d N E E ( r e p r e s e nti n g t h e st r o n g e st o b s e r v e d a n n u al c a r b o n si n k) w a s < 7 0 0 g C m 2 , w hi c h pl a c e s t h e s e f o r e st s wit hi n t h e t o p 1 0 % of o b s e r v e d a n n u al c a r b o n si n k s n ot o nl y i n t h e c o nti n e nt al U nit e d St at e s (i. e. Fi g. 1) b ut I ntr o d u cti o n T h e f o r e st s of t h e w a r m a n d m e si c s o ut h e a st e r n U nit e d St at e s c o nf e r a n u m b e r of e c o n o mi c b e n e fit s a n d e c o- l o gi c al s e r vi c e s t h at a r e di r e ctl y li n k e d t o p att e r n s of e c o s y st e m c a r b o n a n d w at e r c y cli n g. S p e ci fi c all y, t h e y p r o d u c e m o r e ti m b e r t h a n a n y ot h e r f o r e st r y r e gi o n i n t h e U nit e d St at e s ( W e a r & G r ei s, 2 0 1 2), a n d t h e y pl a y k e y r ol e s i n c o nt r olli n g t h e q u a ntit y a n d q u alit y of st r e a m w at e r a n d w at e r a v ail a bl e f o r h u m a n O n t h e diff er e n c e i n t h e n et e c o s y st e m e x c h a n g e of C O 2 b et w e e n d e ci d u o u s a n d e v er gr e e n f or e st s i n t h e s o ut h e a st er n U nit e d St at e s K I M B E R L Y A . N O V I C K 1 , 2 , A . C H R I S T O P H E R O I S H I 2 , 3 , E R I C J . W A R D 2 , 4 , M A R I O B . S . S I Q U E I R A 2 , 5 , J E H N - Y I H J U A N G 2 , 6 a n d P A U L C. S T O Y 2 , 7 1 S c h o ol of P u bli c a n d E n vir o n me nt al Aff airs, I n di a n a U ni versit y – Bl o o mi n gt o n, 7 0 2 N. W al n ut Gr o ve A ve n ue, Bl o o mi n gt o n, I N 4 7 4 0 5, U S A, 2 Ni c h ol as S c h o ol of t he E n vir o n me nt, D u ke U ni versit y, B o x 9 0 3 2 8, D ur h a m, N C 2 7 7 0 8, U S A, 3 C o weet a H y dr ol o gi c L a b or at or y, S o ut her n Rese ar c h St ati o n, U S D A F orest Ser vi ce, 3 1 6 0 C o weet a L a b R o a d, Ott o, N C 2 8 7 6 3, U S A, 4 De p art me nt of F orestr y a n d E n vir o n me nt al Res o ur ces, N ort h C ar oli n a St ate U ni versit y, 9 2 0 M ai n C a m p us Dri ve S uite 3 0 0, R alei g h, N C 2 7 6 9 5, U S A, 5 De p art me nt of Me c h a ni c al E n gi neeri n g, U ni versit y of Br asili a, Br asili a, Br a zil, 6 De p art me nt of Ge o gr a p h y, N ati o n al T ai w a n U ni versit y, N o. 1 Se c. 4 R o ose velt R o a d, T ai pei 1 0 6 1 7, T ai w a n, 7 De p art me nt of L a n d Res o ur ces a n d E n vir o n me nt al St u dies, M o nt a n a St ate U ni versit y, 3 3 4 Le o n J o h ns o n H all, B o ze m a n, M T 5 9 7 1 7, U S A A b str a ct T h e s o ut h e a st e r n U nit e d St at e s i s e x p e ri e n ci n g a r a pi d r e gi o n al i n c r e a s e i n t h e r ati o of pi n e t o d e ci d u o u s f o r e st e c o- s y st e m s at t h e s a m e ti m e it i s e x p e ri e n ci n g c h a n g e s i n cli m at e. T hi s st u d y i s f o c u s e d o n e x pl o ri n g h o w t h e s e s hift s will aff e ct t h e c a r b o n si n k c a p a cit y of s o ut h e a st e r n U S f o r e st s, w hi c h w e s h o w h e r e a r e a m o n g t h e st r o n g e st c a r b o n si n k s i n t h e c o nti n e nt al U nit e d St at e s. U si n g ei g ht- y e a r-l o n g e d d y c o v a ri a n c e r e c o r d s c oll e ct e d a b o v e a h a r d w o o d d e ci d u o u s f o r e st ( H W) a n d a pi n e pl a nt ati o n ( P P) c o-l o c at e d i n N o rt h C a r oli n a, U S A, w e s h o w t h at t h e n et e c o s y st e m e x c h a n g e of C O 2 ( N E E) w a s m o r e v a ri a bl e i n P P, c o nt ri b uti n g t o v a ri a bilit y i n t h e diff e r e n c e i n N E E b et w e e n t h e t w o sit e s (D N E E) at a r a n g e of ti m e s c al e s, i n cl u di n g t h e i nt e r a n n u al ti m e s c al e. B e c a u s e t h e v a ri a bilit y i n e v a p ot r a n s pi r a-ti o n ( E T) w a s n e a rl y i d e nti c al a c r o s s t h e t w o sit e s o v e r a r a n g e of ti m e s c al e s, t h e f a ct o r s t h at d et e r mi n e d t h e v a ri a bil-it y i n D N E E w e r e d o mi n at e d b y t h o s e t h at t e n d t o d e c o u pl e N E E f r o m E T. O n e s u c h f a ct o r w a s w at e r u s e ef fi ci e n c y, w hi c h c h a n g e d d r a m ati c all y i n r e s p o n s e t o d r o u g ht a n d al s o t e n d e d t o i n c r e a s e m o n ot o ni c all y i n n o n d r o u g ht y e a r s (P < 0. 0 0 1 i n P P). F a ct o r s t h at v a r y o v e r s e a s o n al ti m e s c al e s w e r e st r o n g d et e r mi n a nt s of t h e N E E i n t h e H W sit e; h o w e v e r, s e a s o n alit y w a s l e s s i m p o rt a nt i n t h e P P sit e, w h e r e si g ni fi c a nt a m o u nt s of c a r b o n w e r e a s si mil at e d o ut si d e of t h e a cti v e s e a s o n, r e p r e s e nti n g a n i m p o rt a nt a d v a nt a g e of e v e r g r e e n t r e e s i n w a r m, t e m p e r at e cli m at e s. A d diti o n al v a ri a bilit y i n t h e fl u x e s at l o n g-ti m e s c al e s m a y b e att ri b ut a bl e t o sl o wl y e v ol vi n g f a ct o r s, i n cl u di n g c a n o p y st r u ct u r e a n d i n c r e a s e s i n d o r m a nt s e a s o n ai r t e m p e r at u r e. T a k e n t o g et h e r, st u d y r e s ult s s u g g e st t h at t h e c a r b o n si n k i n t h e s o ut h e a st e r n U nit e d St at e s m a y b e c o m e m o r e v a ri a bl e i n t h e f ut u r e, o wi n g t o a p r e di ct e d i n c r e a s e i n d r o u g ht f r e q u e n c y a n d a n i n c r e a s e i n t h e f r a cti o n al c o v e r of s o ut h e r n pi n e s. experience similar climate or are located at approxi- mately the same latitude (Fig. 1, Stoy et al., 2008; Powell et al., 2008). Complicating efforts to predict how these patterns of land-use change will affect regional carbon and water cycling is the fact that NEE in both southeastern US for- ests types is highly variable from year to year. For instance, the range of annual NEE exceeds 300 g C m2 in six of the eight southeastern US forests ecosystems surveyed in Fig 1. This variability may reflect patterns of canopy development, but it also may reflect sensitiv- ity to meteorological drivers, which themselves are nonstationary. Indeed, over the past century, changes in annual temperature, atmospheric CO2 concentration, and hydrologic regime have been observed in the east- ern United States (Ford et al., 2011; Laseter et al., 2012; Brzostek et al., 2014) and are predicted to continue in the future (Huntington, 2006; Dai, 2011). Thus, efforts to predict the magnitude and variability in regional NEE in the southeastern United States must consider the extent to which southeastern US hardwood and pine forests differentially respond to climate forcings. A number of independent investigations have char- acterized the relation between the carbon and/or water vapor exchange of southeastern US forests and their meteorological drivers (Wilson & Baldocchi, 2001; Lai et al., 2002; Sch€afer et al., 2002, 2003; Pataki & Oren, 2003; Clark et al., 2004; Wullschleger et al., 2004; Palm- roth et al., 2005; Mccarthy et al., 2006; Novick et al., 2009a; Noormets et al., 2010; Sun et al., 2010; Ward et al., 2013; Whelan et al., 2013). In general, these stud- ies reveal that the magnitude of NEE tends to be highly sensitive to hydrologic conditions (Baldocchi, 1997; Wilson & Baldocchi, 2001; Stoy et al., 2005; Noormets et al., 2010; Whelan et al., 2013) and to patterns of Fig. 1 The mean (symbols) and range (lines) of NEE derived from eddy covariance flux towers in the continental United States. The sites were selected by identifying all sites in the Ameriflux network with at least 4 years of available data, and then extracting previ- ously published estimates of NEE from the literature (see supplementary information for more details). Grasslands and shrublands are identified with squares, agricultural sites with upward pointing triangles, evergreen forests with circles, deciduous forests with dia- monds, and wetland sites with right pointing triangles. Southeastern US forests are shown in black. Negative NEE indicates a carbon sink. also globally (see Baldocchi, 2008, who reported a mean NEE of 183 � 270 g C m2 from 506 site years of data representing a range of global biomes). Further, the mean value of annual NEE in most southeastern US forests is < �400 g C m2 , which is outside of the range of observed NEE in two-thirds of the ecosystems surveyed. However, in the southeastern United States, the mag- nitude of NEE and related ecosystem benefits are sensi- tive to ongoing shifts in the land cover regime affecting the relative distribution of deciduous and evergreen forests. Much of the eastern United States, which was historically dominated by late-successional hardwood forests and fire-dominated pine woodlands (Stanturf et al., 2002; Abrams, 2003), was cleared and converted to farmland between 1800 and 1930 (Hall et al., 2002; McEwan et al., 2011). Since then, agricultural abandon- ment beginning in the mid-20th century has resulted in an increase in the fractional forest cover of the south- eastern United States (Ramankutty et al., 2011; Wear & Greis, 2012). Following well-established patterns of eco- logical succession (Oosting, 1942), many of these aban- doned farm fields have developed into mixed deciduous hardwood forests. However, due to a dou- bling of softwood timber production in the southeast- ern United States over the last 40 years (Wear & Greis, 2012), pine plantations now comprise over 160 000 km2 in the region. The abundance of southern pine forests is expected to increase over coming decades, with an important shift toward establishment on cleared hard- wood forest lands (Wear & Greis, 2012). While both southern pine and hardwood forests are capable of assimilating large amounts of CO2 from the atmo- sphere, the range of annual NEE can differ consider- ably between these two ecosystems, even when they canopy development related to disturbance (Clark et al., 2004; Mccarthy et al., 2006; Stoy et al., 2008). However, only a few of these studies have attempted such an investigation in colocated pine and hardwood forests, which permits a disentangling of how the two ecosystem types differentially respond to meteorologi- cal forcings. Fewer still attempt such an investigation over sufficiently long time periods (i.e. near-decadal) to evaluate the effect of drivers that vary over time scales commensurate with climate change and canopy devel- opment. The principle objective of this study, which leverages long-term flux records collected in colocated pine and hardwood forests in central NC (USA), was to explore how carbon and water fluxes in these important forest types differentially respond to climatic forcings that operate over timescales ranging from hours to years. In meeting this objective, we will specifically address the following research questions: ● To what extent do the magnitude and variability in ecosystem-scale carbon and water fluxes differ between the two principle southeastern US forest types? ● What timescales of variability are associated with the largest differences in fluxes between the two sites? ● What biophysical mechanisms are most likely responsible for observed differences? The answers to these questions will reflect the rela- (Katul et al., 2010; Manzoni et al., 2011; Keenan et al., 2013). Methods The study sites The pine plantation (PP) and hardwood forest (HW) study sites were located in the Blackwood Division of the Duke Forest near Durham, North Carolina (35° 580 41″N, 79° 050 59″W, 163 m a.s.l.), and the PP site was situated within an ambient CO2 ring of the Duke Free-Air Carbon Enrichment (FACE) experiment. These ecosystems vary in vegetation cover and canopy structure, but experience nearly identical climatic and edaphic conditions. The pine plantation was established in 1983 following a clear cut and a burn. Pinus taeda L. (Loblolly pine) seedlings were planted at a 2.0 m by 2.4 m spacing with pine density reduced to approximately 1100 trees ha�1 after 15 years. Canopy height increased from 16 meters in 2001 to over 20 meters in 2008. The hardwood forest is classified as an uneven-aged (90–110 year old) oak (Quercus)–hickory (Carya) forest. The stand is dominated by hickories [Carya tomentosa (Poir.) Nutt., C. glabra (P. Mill). Sweet.], yellow poplar (Liriodendron tulipifera L.), sweetgum (Liquidambar styraciflua L.), and oaks (Quercus alba L., Q. michauxi Nutt., Q. phellos L.). Coniferous species including Pinus taeda L. and Juniperus virginiana L. make up a minor component of the over- and understory, respectively (Pataki & Oren, 2003; Oishi et al., 2008). The forest was not man- aged after establishment, and mean canopy height is 25 m. Both ecosystems have little topographic variation and lie on Enon silt loam, and a clay pan underlies the research sites at a depth of ~ 35–50 cm, imposing similar constraints on root-water access for both ecosystems. Long-term mean annual temperature and precipitation are 15.5 °C and 1146 mm. More details about the study sites are available in Stoy et al. (2006a,b) and Novick et al. (2009a). Eddy covariance, meteorological, and leaf area measurements NEE and the latent heat flux (LE) were measured above the canopies from 2001 to 2008 using eddy covariance systems comprised of triaxial sonic anemometers (CSAT3, Campbell Scientific, Logan, UT, USA) and open-path infrared gas analyz- ers (LI-7500, Li-Cor, Lincoln, NE, USA), with the flux towers decommissioned in 2009. Wind and scalar concentration mea- surements were collected at 10 Hz, and for most of the study period, fluxes were processed in real time. The gas analyzers were calibrated twice each year (at the beginning and end of the growing season), with little variation in the calibration coef- ficients observed over time. Half-hourly meteorological and edaphic measurements were also collected in each site, includ- ing air temperature (Ta), vapor pressure deficit (D), photosyn- thetically active radiation (PAR), soil moisture content (h) and precipitation (PPT). Friction velocity (u*) was calculated from the high-frequency wind velocity data. In each site for a portion of the study period, a multilayer concentration monitoring tive advantages for ecosystem carbon uptake of decidu- ousness vs. evergreenness in the warm, mesic southeastern US. Predictions drawn from classic stud- ies on the topic include: (i) greater magnitude of NEE in the evergreen forest, reflecting the relatively large difference between gross ecosystem productivity and ecosystem respiration predicted for mid-successional forests (Odum, 1969), (ii) greater drought tolerance in the deciduous site (Givnish, 2002), and (iii) relatively less seasonal contrast in evergreen forest NEE (Hollin- ger, 1992; Givnish, 2002). With this study, we will build upon these previous efforts and extend the scope of analysis to characterize the variability in NEE driven by other factors that vary over time scales commensurate with climate change and canopy development (i.e. near-decadal). At these timescales, we expect that biophysical mechanisms responsible for these differences, if they exist, may include a range of processes already shown to impact long-term forest carbon and water cycling in other eco- systems, including canopy development and recovery from disturbance (Law et al., 2001; Amiro et al., 2010), long-term trends in temperature and the length of the growing season (Dragoni et al., 2011; Keenan et al., 2014), and long-term variation in water use efficiency (a) (b) Fig. 2 The average NEE (panel a) and ET (panel b) in PP as a function of the maximum permissible fraction of the flux foot- print that overlies the fertilized portion of the study area (FFERT). Gray lines represent 95% confidence intervals determined from a nonparametric bootstrap. model, Stoy et al. (2006b) determined that the effect of this clear cut was evident in the flux records. Thus, the HW flux data were filtered to remove data collected when more than 30% of the flux footprint originated from the clear cut. These data were further adjusted as described in Stoy et al. (2006b) to account for the fact that the application of the footprint filter in 2003 also selectively filters ideal climate conditions. Eddy covariance data processing – NEE fluxes The Webb–Pearman–Leuning correction (Webb et al., 1980) for air density fluctuations was applied to the half-hourly NEE and LE fluxes. The NEE data were initially screened to remove any flux measurements whose absolute magnitude exceeded 100 lmol m�2 s�1. The data were then screened to remove measurements collected during suboptimal meteorological conditions using the online eddy covariance processing tool (http://www.bgc-jena.mpg.de/~MDIwork/eddyproc), which employs the u* filtering method of Reichstein et al. (2005). These filtering and quality control procedures create a num- ber of gaps in the data records that must be gapfilled to derive annual flux estimates (Falge et al., 2001). Furthermore, to understand the processes that control NEE, it is desirable to partition measured NEE fluxes into the primary components GEP and RE, where NEE = �GEP + RE, and a negative NEE flux denotes carbon assimilation by the ecosystem. We adopted the gapfilling and partitioning methods of Van Gorsel et al. (2009), hereafter VG2009, which is designed to exclude data collected when nocturnal vertical and horizontal advec- tion fluxes prevent a straightforward interpretation of mea- sured turbulent fluxes as representative of the integrated ecosystem scalar sources and sinks. The approach assumes system was installed to measure the mean water vapor and CO2 concentration at 10 different levels throughout the canopy volume using a Li-Cor 6262 CO2/H2O infrared gas analyzer. These profiling systems included a multi-port gas sampling manifold to sample each level for 1 min (45 s sampling and 15 s purging) with a repeated cycle of 10 min. Profile data were recorded as 30 min averages. Evapotranspiration (ET, mm) was derived from the latent heat flux time series using the tempera- ture-dependent latent heat of vaporization. Ground-based leaf area index (LAI) measurements were made in PP during the entire eight year study period (Mccar- thy et al., 2007), and from 2001 to 2005 in HW (Oishi et al., 2008) using litterfall and light-penetration measurements, and then interpolated to a daily time step. To generate a complete eight year record of leaf area index in HW, MODIS LAI data (Myneni et al., 2002) were obtained for the pixel that includes the hardwood site. The MODIS data were interpolated and scaled so that the mean minimum and maximum MODIS LAI from 2001 to 2005 match the mean minimum and maximum ground-based LAI from the same time period. The interpo- lated and scaled MODIS time series from 2006 to 2008 was then appended to the ground-based measurements to gener- ate the eight year record for the HW site. Estimates of the monthly Palmer Drought Severity Index (PDSI) at the study sites were obtained from the National Oceanic and Atmo- spheric Administration (NOAA, http://www.ncdc.noaa.gov/ temp-and-precip/drought/historical-palmers.php). Flux footprint considerations Fluxes were screened to remove measurements collected when the footprint was not representative of the study ecosys- tem. In PP, the flux records may be contaminated by fluxes originating from a small portion of the study site receiving annual fertilization from 2005 to 2008 as part of FACE proto- col. We used the footprint model of Hsieh et al. (2000), extended to two dimensions by Detto et al. (2006), to estimate the flux footprint for every half-hour averaging period, and then quantified the fraction of the footprint that overlaid the fertilized area (FFERT). We then iteratively estimated the total daytime NEE and ET measured from 2005 to 2008 (i.e. the fer- tilization period) by progressively increasing the maximum permissible value of FFERT. We found that the total NEE remained relatively stable as long as FFERT was less than 0.10 (Fig. 2a). A similar relation was not observed in the ET fluxes, though the change in ET as the FFERT filter was increased from 0.0 to 0.10 was small (<5%). As FFERT increased further, NEE decreased (i.e. became larger sink for CO2), and ET continued to increase (Fig. 2b). We note that FFERT values greater than 0.5 comprised less than 0.5% of the data record. On the basis of these results, we excluded all PP flux measurements associ- ated with FFERT > 0.10. While the results are not shown here, we note that the study results are largely invariant for FFERT thresholds as low as 0.03, and also do not change if the fertil- ization sector filter is applied to data collected before the nutri- ent amendments began (i.e. years 2001 to 2004). A clear cut on private land 200 m south of the HW tower occurred in late November 2002. Using the same footprint that the contribution of advection to the total ecosystem flux is minimized during early evening hours, when the sum of mea- sured eddy covariance fluxes and canopy storage fluxes should be representative of total ecosystem respiration fluxes. During periods when the gas profile systems were active, stor- age fluxes of CO2 were estimated by integrating temporal changes in CO2 concentration within the canopy (Juang et al., 2006; Novick et al., 2014); otherwise, they were estimated using site-specific relationships between the ratio of storage to turbulent CO2 fluxes and friction velocity (u*). Data representing the sum of turbulent and storage fluxes collected within two hours of the time of the observed monthly peak in measured RE were used to parameterize a temperature-dependent model of the form RE = e1exp(e2Ta), where e1 is a parameter related to the base respiration rate, and e2 describes the temperature sensitivity of RE. The latter was estimated on an annual time step to avoid situations in which the temperature range of the data used to parameterize the model was much smaller than the temperature range of the data needing to be gapfilled. The base respiration rate parameter was then estimated on a weekly time step using all acceptable, early evening data collected within 3 weeks (i.e. by using a moving 6 week window). This continuous RE record was used to gapfill missing nocturnal data. Missing daytime data, which represent <5% of all missing data, were gapfilled with the approach of Lasslop et al. (2010), hereafter L2010, which is supported by the online flux tool. The L2010 approach links GEP to PAR via: GEP ¼ abPAR aPARþ b ð1Þ where a is the mean apparent ecosystem quantum yield, and b is the maximum assimilation rate which is modified to account for humidity effects according to: b ¼ bo when D� 1 kPa ð2Þ b ¼ bo exp ½�kðD�DoÞ� when D[ 1 kPa where bo is a reference assimilation rate at low D, k is the humidity response parameter, and Do is set to 1 kPa. Follow- ing L2010, daytime respiration is modeled as a function of air temperature according to: RE ¼ ro exp Eo 1 Tref � To � 1 Ta � To � �� � ð3Þ Eddy covariance data processing – LE fluxes The same u* filter described in the previous section was also applied to the LE flux records. In addition, measured LE fluxes outside of a conservative expectation window of �200 to 800 W m2 were removed from the time series. For consis- tency with a previous study, LE records were gapfilled using the marginal distribution methodology (Reichstein et al., 2005) described in Novick et al. (2009a). ET flux magnitudes pre- sented here differ slightly (<5%) from those presented Novick et al. (2009a) owing to different filtering procedures, and the fact that Novick et al. (2009a) did not include the footprint fil- ters, but were nonetheless within the uncertainty bounds of annual ET flux in the study forests (Stoy et al., 2006a; Novick et al., 2009a). Spectral analysis of flux records To explore the timescales of variability in NEE and ET, we relied on a wavelet spectral decomposition. We refer readers to Stoy et al. (2005, 2009), Torrence & Compo (1998), and Ka- tul et al. (1998) for a thorough discussion of the utility and application of wavelet transforms, which are an appropriate tool for spectral analysis of nonstationary and discontinuous time series (Katul et al., 1998; Stoy et al., 2005; Dietze et al., 2011). Briefly, the purpose of this spectral analysis was to quantify the energetic frequencies/timescales in a signal (i.e. NEE or ET), and to relate those energetic timescales to the timescales of variability in meteorological drivers. Determin- ing the wavelet spectra of a single time series (i.e. NEE or ET) permits an identification of the fraction of the variability in the signal attributable to processes operating at specific timescales/frequencies. The determination of a wavelet cospectra describing interactions between two time series permits an identification of the frequencies at which they most significantly coresonate. We employed an orthonormal wavelet transformation with the Haar basis function to characterize the energetic frequen- cies of NEE and ET by relying only on measured data, and not gapfilled data (which were set = 0), so that the spectral decomposition would not be strongly affected by the choice of gapfilling approach. The filtered time series were normalized to have zero mean and unit variance, and were zero-padded to a length equal to the power of 2 greater than the length of the data record (in this case, 218, noting that the 8 year data records contain 140 256 half-hourly data points). The wavelet coefficients were determined and summed within dyadic time scales ranging from 21 (hourly) to 216 (nearly 4 years) as the wavelet coefficients associated with longer (i.e. ~ 8 year, or 217 h) timescales are generally unreliable (Stoy et al., 2005). We also relied on the wavelet cospectra to explore the extent to which the difference in NEE between the two sites (DNEE) coresonates with meteorological factors, including PAR, Ta, D, h, and atmospheric CO2 concentration. Specifically, DNEE was determined at the half-hourly timestep and then normalized to have zero mean and unit variance, as were the meteorologi- cal drivers. The cross-spectra were determined after Katul et al. (1998). All wavelet analyses were performed in Matlab using the Wavelet Toolbox (Mathworks, Natick, MA, USA). The parameter ro is a base respiration rate and the parame- ter Eo is a temperature response parameter. The reference tem- perature parameter Tref was set to 15 °C and the parameter To was set to �46 °C. The parameters were fit using 4–12 day moving windows. The VG2009 flux records were compared to annual fluxes estimated in an earlier effort by Stoy et al. (2006b), hereafter S2006. The S2006 approach does not rely on the Reichstein et al. (2005) u* filter; rather, data quality control is largely gov-erned by the magnitude of the atmospheric stability parameter (see Novick et al., 2004). The S2006 fluxes were gapfilled and partitioned using a light-response model for NEE that is simi- lar to the L2010 model, but is parameterized at a different time step and does not include a humidity response parameter. WUEi ¼ GEP ET D ð4Þ The second metric, here termed the ‘ecosystem marginal water use efficiency’ (kE), is derived from applications of sto- matal optimization theory to explore the effect of variations in atmospheric CO2, D, and hydrologic conditions on leaf gas exchange (Katul et al., 2009, 2010; Manzoni et al., 2011). If sto- mata are assumed to function so as to maximize carbon gain while minimizing water loss, then the Hamiltonian of the sys- tem can be expressed as H=A-kT and optimality is achievedwith respect to the control variable (stomatal conductance g) via. @H @g ¼ 0 ð5Þ where A and T are leaf-level assimilation and transpiration rates, respectively, g is stomatal conductance to CO2, and k (=@A=@E) is the leaf-level marginal water use efficiency, repre- senting the carbon cost of a unit of transpired water. When Eqn (5) is coupled to a linearized variant of the Farquhar assimilation model (Farquhar et al., 1980; Katul et al., 2010) and to Fickian-diffusion type expressions for A and E, it may be shown, following Katul et al. (2010), that k ¼ @A @E ¼ 1 aca A E � �2 D; ð6Þ where a is the ratio of molecular diffusivities between H2O and CO2, and ca is the atmospheric CO2 concentration. We note that this parameter k is an important input into a relatively new suite of stomatal conductance models based in optimization theory (Katul et al., 2009; Katul et al., 2010; Manzoni et al., 2011), though here we are principally focused on exploring the utility of a big-leaf (i.e. ecosystem-scale) analog to Eqn (6) (i.e. kE) to explain cross-site differences in the dynamics of GEP and ET. The parameter kE is determined from: kE ¼ 1 aca GEP ET � �2 D ð7Þ The units of WUEi are expressed as g C kg�1 H2O/hPa for consistency with Keenan et al. (2013), and the units of kE are expressed as lmol m�2 s�1 for consistency with Katul et al. (2009). Interpretation of both WUEi and kE requires that ET be dominated by transpiration. During the growing season, LAI at the two study sites typically exceeds 5 m2 m�2, such that little radiation reaches the canopy floor, and soil evaporation is assumed to be small (Stoy et al., 2006a; Oishi et al., 2008). Thus, we limit the analysis of water use efficiency to data col- lected when leaf area index is relatively stationary (May–Sep- tember), and when at least 2 days have elapsed since a rain event, after which it is assumed that the contribution of can- opy interception evaporation to total ET is negligible, and soil evaporation is likely limited not only by radiation but also by soil moisture. Thus, for continued discussions of these water use indices, we assume ET and E are indistinguishable. The quantities WUEi and kE are related, but with some key differences. Mathematically, given Fickian-diffusion type rela- tions for GEP [GEP= g/(ca – ci)] and ET (ET = agD), where g and ci are here interpreted as integrated canopy averages, WUEi reduces to WUEi ¼ GEP ET D ¼ gðca � ciÞ agD D ¼ ca � ci a ð8Þ Or in other words, variations in WUEi are due solely to varia- tions in the quantity (ca–ci). The marginal water use efficiency may then be related to WUEi as: kE ¼ 1 aca GEP ET � �2 D ¼ 1 aca gðca � ciÞ agD � �2 D ¼ 1 acaD WUE2 i ¼ 1 a3caD ðca � ciÞWUEi; ð9Þ Thus, both kE and WUEi are sensitive to variations in the quantity (ca – ci), which could be due to increases in ca or increases in D that induce stomatal closure and thereby reduce ci. However, they will become increasing decoupled with increasingly large changes in (ca – ci) and D, reflecting the fact that these two driving variables may affect the carbon cost of a unit of transpired water beyond their effect on WUEi. Both quantities are useful for interpreting how gas exchange dynamics respond to varying climatic conditions, though the parameter kE has the advantage of being derived from a theo- retical framework grounded in the assumption that stomata function to minimize water loss while maximizing carbon gain. Results Meteorological and hydrological variability during the study period The study period was marked by variability in meteo- rological and hydrological forcing factors operating at a Phenological and seasonality indices We relied on an approach for active season delineation that first identifies the days of each year when the time derivative of GEP (dGEP/dt) reached a minimum and maximum, respectively, representing the days when GEP was changing most rapidly in spring and in fall. Before identifying the mini- mum and maximum dGEP/dt, a smoothed spline curve was fitted to the daily time series of GEP. We choose to use GEP, and not LAI, to delineate the active season for two reasons: (i) ground-based measurements of LAI were not performed at a sufficiently fine temporal frequency in the HW site (i.e. weekly or every 2 weeks) for the duration of the study period, and (ii) recent work has shown that the relation between assimilation capacity and leaf area index varies considerably over the course of the growing season (Bauerle et al., 2012). Ecosystem-scale water use efficiency We considered two approaches for characterizing water use effi- ciency at each site. The first is the ‘inherent water use efficiency’ (WUEi) metric proposed by Beer et al. (2009) and more recently used by Keenan et al. (2013). The inherent water use efficiency recognizes the potential for variation in D to modify the total water use efficiency (i.e. GEP/ET, presuming ET is dominated by transpiration, Jasechko et al., 2013), and is formulated as: range of time scales, and included two severe drought events. The first began late in 2001 and persisted until the latter part of the growing season in 2002, and the second drought period began early in the 3rd quarter of 2007 and persisted through the 1st quarter of 2008. During both of these drought events, the regional Pal- mer Drought Severity Index (PDSI) dipped below �4 (Fig. 3), which is considered the threshold for extreme drought. Moderate drought events were experienced in 2001 and 2005 (Fig. 3), when soil moisture, precipita- tion, and the PDSI were below average, but excursions were not as severe as in 2002 or 2007. The study sites experienced a severe ice storm event in December 2002. The storm reduced leaf area index in PP by 12–14% in 2003, though the canopy quickly recovered (Mccarthy et al., 2006). The ice storm had lit- tle effect on canopy structure in HW, where the absence of leaves during the dormant season reduced branch ice accumulation and thus damage. The study period was characterized by a long-term increasing trend in active season D during non-drought years (slope = 0.05 kPa yr�1, P = 0.04), and a long-term trend in dormant season air temperature after exclud- ing drought years (slope = 0.03 degrees C yr�1, P = 0.03). There was no observed trend in mean annual temperature. Annual flux magnitudes and trends within and across sites The choice of QC/gapfilling methodology (e.g. VG2009 vs. S2006) affects the magnitudes of the fluxes (Table 1, Fig. 4). We note that while S2006 estimated fluxes for year 2005, we omitted them from this comparison as that study did not employ a footprint filter to remove measurements originating from the fertilized sector of PP. During the years for which fluxes estimated with the VG2009 and S2006 approach were available (2001– (a) (b) (c) (d) (e) (f) Fig. 3 Climatic variables measured during the study period. Panel (a) shows the interpolated Palmer Drought Severity Index (PDSI). Other panels show quarterly deviations from the study period mean air temperature (Ta, b), precipitation (PPT, c), vapor pressure defi- cit (D, d), soil moisture content (h, e), and photosynthetically active radiation (PAR, f), as measured in the pine site. The two drought events are highlighted in gray, and the date of the ice storm is indicated with the asterisk. 2005 for HW, 2001 – 2004 for PP), VG2009 and S2006 estimates of mean annual NEE differed by 60 g C m�2 yr�1 in the PP and 88 g C m�2 yr�1 in HW, on average, representing 14% and 20% of the mean annual NEE, respectively. Over the course of the comparison period, the VG2009 and S2006 estimates of PP GEP and RE agreed to within 200 and 130 g C m�2 yr�1 , or 9% and 7%, respectively. The VG2009 and S2006 estimates of HW GEP agreed to within 70 g C m�2 yr�1 (or 4%), and to within 160 g C m�2 yr�1 (or 11%) for HW RE. Particularly large differences between PP GEP and RE fluxes estimated by the two methods were observed in year 2001, which is not surprising as much of the first 3 months of PP data were missing in 2001 due to a tech- nical problem. In general, variations in the annual NEE fluxes between the two approaches are well correlated in the PP site (r = 0.93), and moderately correlated in the HW site (r = 0.66, Fig. 4), which experienced less annual variability as discussed below. For the remain- ing analysis, we use the VG2009 flux records unless otherwise noted. Over the entire study period, both ecosystems were strong carbon sinks. Mean annual NEE in PP Year VG2009 Approach S2006 Approach NEE RE GPP ET NEE RE GPP Pine 2001 �564 (�492) 1760 (1078) �2324 (�1570) 740 (508) �610 (�443) 1340 (907) �1950 (�1350) 2002 �307 (�142) 1745 (1262) �2042 (�1406) 678 (416) �270 (�138) 1610 (972) �1860 (�1111) 2003 �173 (�165) 1665 (1138) �1838 (�1303) 832 (526) �220 (�298) 1730 (888) �1950 (�1186) 2004 �523 (�302) 1733 (1257) �2256 (�1558) 860 (551) �420 (�219) 1760 (1193) �2190 (�1412) 2005 �635 (�391) 1769 (1233) �2404 (�1624) 888 (586) 2006 �712 (�380) 2121 (1557) �2833 (�1946) 916 (577) 2007 �547 (�244) 1824 (1159) �2371 (�1403) 793 (478) 2008 �561 (�140) 2136 (1490) �2678 (�1630) 876 (549) Mean �506 (�282) 1844 (1272) �2344 (�1555) 823 (524) �381 (�274) 1607 (990) �1989 (�1265) SD 173 (131) 181 (169) 317 (197) 81 (56) 173 (130) 189 (140) 130 (140) Hardwood 2001 �367 (�455) 1472 (1060) �1839 (�1516) 622 (452) �510 (�598) 1200 (650) �1710 (�1248) 2002 �341 (�411) 1310 (888) �1651 (�1300) 608 (436) �390 (�445) 1320 (797) �1710 (�1232) 2003 �267 (�332) 1430 (1068) �1696 (�1401) 747 (510) �400 (�531) 1250 (724) �1650 (�1254) 2004 �384 (�453) 1415 (985) �1799 (�1439) 702 (489) �420 (�573) 1313 (705) �1750 (�1278) 2005 �422 (�502) 1460 (962) �1882 (�1464) 814 (599) �490 (�570) 1230 (700) �1720 (�1269) 2006 �602 (�666) 1742 (1148) �2344 (�1815) 783 (553) 2007 �403 (�494) 1852 (1194) �2255 (�1688) 691 (490) 2008 �432 (�634) 1960 (1199) �2393 (�1832) 794 (562) Mean �402 (�494) 1585 (1063) �1989 (�1561) 720 (508) �444 (�543) 1261 (715) �1705 (�1258) SD 96 (110) 243 (113) 305 (192) 78 (56) 52 (60) 53 (54) 35 (15) ( = �506 g C m�2 yr�1) represented a larger carbon sink than mean annual NEE in HW ( = �402 g C m�2 yr�1), though interannual variation was higher in PP (r = 173 g C m�2 yr�1) than HW (r = 96 g C m�2 yr�1). The magnitude of both NEE and RE, and thus GEP, was higher in the pine plantation as compared to the hardwood forest (Table 1). The interannual variability in GEP fluxes was similar in the two sites (r = 317 and 305 C m�2 yr�1 in PP and HW, respectively), and the variability in RE fluxes was higher in the HW site (r = 181 and 234 C m�2 yr�1 in PP and HW, respectively). A stronger coupling between GEP and RE was observed in the hardwood site (r2 = 0.92, P = 0.0001) when compared to the pine site (r2 = 0.81, P = 0.003). The magnitude of ET was higher in PP (Table 1), and in both sites interan- nual variability in ET was small relative to the variabil- ity in carbon dioxide fluxes. Spectral characteristics of NEE and ET in the two sites The normalized wavelet spectra of PP and HW ET were nearly identical (Fig. 5b); however, significant differ- Table 1 Annual fluxes of the net ecosystem exchange of CO2 (NEE), gross ecosystem productivity (GEP), ecosystem respiration (RE) (in g m�2), and evapotranspiration (ET, in mm) at the pine plantation and hardwood forest ecosystems in the Duke Forest, NC. The data show the annual fluxes, and the growing season fluxes in parentheses. Carbon fluxes are calculated using two differ- ent approaches to gapfilling and partitioning: the Van Gorsel et al. (2009) approach (VG2009) and the Stoy et al. (2006b) approach (S2006) 3rd quarters of the two severe drought years (i.e. years 2002 and 2007). We also parameterized the model for the 2nd and 3rd quarters of the years preceding the drought years to provide a basis for comparison. In this applica- tion, we did not include the D response function [Eqn (2)], so that differences in the model parameters (i.e. b) reflect variations in D and h. During the growing season of 2002 and 2007, the assimilation rate in PP was reduced for all PAR, and was reduced at high PAR by 24% during the 2nd quarter of 2002, 5% during the 3rd quarter of 2002, 14% during the 2nd quarter of 2007, and 48% during the 3rd quarter of 2007 (Fig. 6a–d). These reductions are relative to the modeled assimilation rate during the same quarter of the year preceding each drought year (i.e. 2001 and 2006). The respective reductions in HW assimilation were in general not as large (10–21% reduction, Fig. 6a– d). The effect of the drought events on respiration rate at a given temperature was also not as pronounced (Fig. 6e–h), and only in the 2nd quarter of 2007 was res- piration at high air temperatures (i.e. Ta > 20 °C) reduced in both PP and HW (by 0.04 to 0.06 g C m2 s�1, or by 34% and 24%, respectively, Fig. 6f). Relative to years preceding the drought events, grow- ing season ET was reduced in both PP and HW (Table 1), though in general, the reductions were not as great as those observed for GEP (relative decreases of (a) (b) (c) (d) (e) (f) (g) (h) Fig. 4 A comparison between the annual net ecosystem exchange (NEE, a & b), gross ecosystem productivity (GEP, c & d), ecosystem respiration (RE, e & f), and evapotranspiration (ET, g & h) estimated using the VG2009 approach (closed symbols) and those estimated with the S2006 approach (open symbols). ences emerged in the spectra for NEE (Fig. 5a). In particular, the PP NEE was more sensitive to factors that vary over daily-to-weekly timescales, and to fac- tors that vary over long (i.e. interannual timescales). Indeed, in PP, the fraction of variability attributable to long timescales is similar to that attributable to factors that vary seasonally. In contrast, HW NEE was driven principally by factors varying over seasonal timescales, with relatively little interannual variability. At daily-to-weekly timescales, the difference in NEE between PP and HW (i.e. DNEE) was most strongly dri- ven by variables that tend to change during the course of a drought (i.e. D, h, and Ta, Fig. 5c). At seasonal timescales, DNEE resonates most closely with variation in h and Ta (Fig. 5c). At interannual timescales, varia- tion in DNEE resonates most with h and atmospheric CO2 (Fig. 5c). Drought effects on carbon and water fluxes The strong linkages between DNEE and the meteorologi- cal variables h, D, and Ta at a range of timescales (e.g. Fig. 5c) motivate a closer examination of the drought effects on NEE in both sites. Toward that end, a nonlin- ear least-squares regression was used to parameterize the GEP model of Eqn (1) and the RE model of Eqn (3) using observed (i.e. not gapfilled) data from the 2nd and 8–13% in PP and 3–12% in HW for the drought years as compared to the preceding years). Ice storm effects Damage from the December, 2002 ice storm was pri- marily confined to PP. Annual GEP in 2003 was reduced by >400 g C m2 relative to year 2001, and by >200 g C m2 relative to year 2002 (the drought year). However, GEP recovered quickly, and in 2004 and 2005, annual GEP was comparable to that observed in 2001 (Table 1). Neither the S2006 nor the VG2009 (a) (b) (c) Fig. 5 The wavelet spectra of (a) NEE and (b) ET in the two for- ests, illustrating that the timescales of variability in NEE dif- fered considerably between the two ecosystems at weekly to interannual timescales, whereas the timescales of variability in ET are nearly identical. Panel (c) shows the wavelet cospectra between the difference in NEE between PP and HW (DNEE) and a range of meteorological drivers, including photosyntheti- cally active radiation (PAR), air temperature (Ta), soil moisture content (h), vapor pressure deficit (D), atmospheric CO2 concen- tration (CO2), and also monotonically increasing time. The abbreviation ‘IA’ indicates ‘interannual’ timescales. estimates of RE showed a large change in 2003 as com- pared to 2002 or 2004. This is consistent with previous work focused on the effect of ice storms in managed pine plantations (Mccarthy et al., 2006), which suggests that the carbon content of downed biomass during an ice storm event is ~30% of annual net ecosystem pro- ductivity, but will take >25 years to decompose, which means the contribution of the decomposition of downed biomass to total RE in any given year is rela- tively small. It is worth noting that the smallest carbon sink was observed in 2003 in both sites, suggesting that the effects of the ice storm, moderate and severe droughts in 2001–2002, and low PAR in 2003 (Fig. 3) have combined effects among forest types. Seasonality of the fluxes Considerable variability was observed in the timing of the start and the end of the active season in the PP site (Fig. 7). The start of the active season in PP (hereafter SAS) occurred at DOY 108 on average, with the earliest SAS observed in 2007 (DOY 96), and the latest SAS observed in 2003 (DOY 125). The end of the active sea- son (hereafter EAS) in the pine site occurred on DOY 286 on average, and was similarly variable (ranging from DOY 258 to DOY 306). In HW, the timing of the SAS and EAS were less variable than in PP (DOY 102 to 114 for SAS, and DOY 268 to 288for EAS). The length of the active season (hereafter LAS) was 178 days on aver- age in the pine site, and 167 days on average in the hardwood site. The LAS alone was not a significant determinant of annual NEE/GEP/RE fluxes in either site, though the interaction between LAS and drought status was a significant predictor for annual GEP and RE in the HW site (P < 0.10 from an ANOVA considering LAS, a binary variable for drought status equal to 1 in 2002 and 2007, and 0 otherwise, and their interaction in a linear model). In PP, ~25% of annual GEP occurred outside of the active season (Fig. 8), and ~13% of annual GEP occurred during winter months (January–February, November– December, Fig. 8). Dormant season PP GEP was signifi- cantly related to dormant season temperature (r2 = 0.65, P = 0.06). PP and HW GEP were relatively similar dur- ing the growing season (Table 1, Fig. 8). Trends in water use efficiency The annual growing season WUEi and kE were strongly correlated (r2 > 0.94, P < 0.001 in both sites, Fig. 9a,b), and in many years, the magnitude of WUEi and kE were similar in both sites (Fig. 9). Local maxima in both WUEi and kE were observed in both sites during the two extreme drought years (Fig. 9 a,b), reflecting the explored the relationship between the water use effi- ciency metrics and the allometric index, which is defined as the growing season average sapwood area divided by the product of leaf area and canopy height (AS : AL/h, Ward et al., 2013), and which should be related to stomatal conductance rates at reference envi- ronmental conditions (Novick et al., 2009b). We did not find a significant relationship between AS : AL/h and either WUEi or kE when all years, or even nondrought years, were considered. When the analysis was limited to nondrought years after the ice storm disturbance (2003–2006, and 2008), a significant relationship between AS : AL/h and both WUEi and kE emerged (P = 0.01 and 0.01, respectively). However, the evolu- tion of AS : AL/h and CO2 were also strongly correlated during those years (r2 = 0.96). (a) (b) (c) (d) (e) (f) (g) (h) Fig. 6 (a–d) The difference in GEP at a given flux of photosynthetically active radiation (PAR) in the 2nd and 3rd quarters of the two severe drought years (2002 and 2007) relative to the GEP in the same quarter of the preceding year. PP data are represented with dark black lines, and HW data with gray lines. The bottom four panels show difference in RE as a function Ta for the same time periods. The trends shown here were determined from site-specific, quarterly parameterizations of the GEP and RE models of Eqns (1) & (3). The data show the 95% confidence intervals estimated from a nonparametric bootstrap. direct dependence of both metrics on D. When the drought years were excluded, both WUEi and kE tended to increase with atmospheric CO2 concentration in HW (Fig. 9c,d, P = 0.08 and 0.14, respectively) and in PP (P = 0.0011 and 0.001, respectively). However, in HW at least, the marginal water use efficiency tended to be less sensitive to ca than WUEi. Specifically, when both were scaled to standard normal variables, the slope of the relationship with ca was 0.76 for WUEi as compared to 0.68 for kE. Following Eqns. 8 & 9, the trends in the water use efficiency variables may also be affected by the long-term trend in D over the course of the study period, which is significant even when drought-years were excluded (see Fig. 2d, P<0.05). In PP, which experienced significant changes in can- opy structure during the study period, we also Discussion The timescales of variability in NEE differed consider- ably between the two forests (Fig. 5a), with relatively larger variability in HW NEE at seasonal timescales, and relatively larger variability in PP NEE over interannual timescales. The modes of variability in ET are nearly identical between the two sites (Fig. 5b), which would indicate that ET in both sites is strongly driven by meteorological forcings, and the factors driv- ing variability in DNEE are limited to eco-physiological processes that tend to decouple the dynamics of NEE and ET. In the following subsections, we will illustrate that those factors are likely: (i) variations in GEP and water use efficiency during drought periods, (ii) cross- site differences in the seasonality of NEE, and (iii) eco- system-responses to slowly evolving bio-meteorological drivers. Together, the results highlight some of the rela- tive advantages of evergreenness vs. deciduousness in determining NEE in the southeastern United States. In particular, we will discuss how the relative invariability Fig. 7 The day of year (doy) of the start and end of the active seasons (panels a and b), and the resulting length of the active season (panel c) for each year in the study period. Fig. 8 Mean monthly average NEE (top row), GEP (middle row), and RE (bottom row) in the two study sites. of PP NEE at seasonal timescales represents a distinct advantage for evergreen forests growing in warm cli- mates; however, increased drought sensitivity of GEP in PP relative to HW represents a disadvantage that may become more important if future climates are characterized by increased frequency or severity of drought (Huntington, 2006; Dai, 2011). Differences in the magnitude of carbon fluxes When aggregated over the entire 8 year study period, PP tended to capture and store more carbon than HW (mean annual DNEE = 100 g C m2, Table 1). This result refines earlier work based on shorter time series from these sites (Stoy et al., 2008), and generally agrees with theoretical predictions for the carbon storage capacity of late-successional forests (i.e. Odum, 1969; Tang et al., 2014). However, this result is sensitive to the high degree of interannual variability in NEE in both sites (Table 1), and the related fact that the annual DNEE ● Drought affects GEP more than RE, and affects PP GEP more than HW GEP: Drought decreased GEP capacity during the severe droughts in both sites, whereas the drought effects on RE were mixed and generally small (Fig. 6), consistent with previous work (Ciais et al., 2005; Schwalm et al., 2010). The relative reduction to growing season PP GEP during drought years (~23%) exceeded that observed in the HW (<14%, Fig. 6). ● Ecosystem marginal water use efficiency mediates drought effects: The ecosystem-scale water use effi- ciency increased during drought years in both sites, but not identically so (Fig. 9), Thus, water use effi- ciency is an important variable driving the decou- pling between DNEE and ET. Flux seasonality Net ecosystem exchange resonated strongly with fac- tors that vary seasonally in the HW, but not the PP (Fig. 5). Further, the length of the active season was not a significant predictor for any of the annual PP carbon fluxes, though its interaction with drought status was related to HW GEP and RE (P < 0.10). Thus, with respect to HW, our results confirm previous work showing that LAS is significant control on interannual Fig. 9 Trends in annual intrinsic (WUEi) and ecosystem marginal water use efficiency (kE). Panels a & b show the time series for the two variables over the study period. Panels c & d show the relation between mean growing season CO2 concentration and the water use efficiency variables. Drought years are indicated with gray circles. Regression lines for the relation between WUEi/kE and ca for nondrought years are shown with solid lines for HW and dashed lines for PP. They are significant at the 90% confidence levels for all regressions except HW kE (P = 0.14). Triangles show PP data, and squares show HW data. was not significantly different from zero (P = 0.18). The long-term magnitudes of NEE are also sensitive to the gapfilling and partitioning approach (Table 1, and see also Stoy et al., 2006b). We selected the van Gorsel et al. (2009) approach because the annual RE estimated by other approaches in general did not exceed estimates of annual soil respiration alone (determined from cham- ber-based measurements to be 1383 � 152 g C m�2 yr�1 in PP, Oishi et al., 2013; Palmroth et al., 2005; and 1209 � 99 g C m�2 yr�1 in HW, Oishi et al., 2013). Still, in both PP and HW, the ratio of soil respira- tion to total respiration estimated by the Van Gorsel approach is near unity in the dormant season (~1.06 � 0.21 and 1.05 � 0.25 for PP and HW, respec- tively), and thus tower-derived RE fluxes are still likely underestimated in the winter months. Drought effects The DNEE resonated most strongly with variables that change most during drought events (i.e. D, h, Ta) at daily-to-weekly timescales, and at interannual time- scales (Fig. 5), suggesting that variation in hydrologic condition is an important, if not the predominant, fac- tor driving DNEE. Our results allow us to draw the following conclusions about drought effects in the study sites: (~6 lmol/mol/kPa) is less than the observed increases in kE, which may reflect the influence of an unidenti- fied process, or uncertainty in the data themselves. The observed agreement between DNEE and CO2 at long time scales may be confounded by the effect of other variables like Ta, which may be slowly evolving in concert with increases in CO2. Here, Ta was relatively poorly related with DNEE at long-time scales (Fig. 5c, d), and no long-term trend was observed in mean annual Ta (P = 0.20) or the LAS (Fig. 7). However, there did exist a significant trend in dormant season air tem- perature (slope = 0.03 degrees C yr�1, P = 0.03 after excluding drought years), which is related to the vari- ability in dormant season PP GEP (r2 = 0.65, P = 0.06). Factors related to canopy development may also be important in determining DNEE, even if they do not affect water use efficiency. A direct effect of canopy development on DNEE would likely be realized as a long-term flux trend that is monotonic from 2003 to 2008 (i.e. during the period after the ice storm distur- bance). Such a trend was suggested but not significant for PP NEE (P = 0.20 for years 2003–2008, Fig. 4). Long- term trends in PP GEP and RE were observed, but they were similar to long-term trends observed in the slowly growing HW site (Fig. 4). Thus, while we cannot rule out canopy development as an significant driver of DNEE, the data do not indicate that it is the predomi- nant controlling variable. Finally, while we were able to show here that DNEE is strongly sensitive to factors that vary over interannu- al timescales, we recognize that our ability to confi- dently identify these factors would be facilitated by longer data sets. Thus, efforts should be made to lengthen the time series at long-running flux sites, even if such an effort is not possible for PP and HW (which were decommissioned in 2009). Implications for the variability in future regional carbon and water fluxes Within the context of ongoing climate change and regional land-cover change, this study highlights sev- eral processes important to future regional carbon and water cycling. First, the variability in ET is similar in both sites (Fig. 5b) and relatively small, consistent with results from a previous study investigating the relative temporal invariance of ET (Oishi et al., 2010). Second, predicted increases in drought frequency (Dai, 2011), if realized, will likely promote a more variable and gener- ally smaller regional carbon sink, due to the sensitivity of GEP, and in particular PP GEP, to drought (Fig. 6). Third, future increases in atmospheric CO2 will likely promote increases in water use efficiency, though the extent to which those increases will lead to decreased flux variability in deciduous forests (Dragoni et al., 2011), but highlight that the delineation of an active sea- son may not be appropriate in PP. Indeed, given that PP and HW GEP were relatively similar during the growing season (Table 1, Fig. 8), much of the difference in annual PP and HW GEP (which exceeds 300 g C m2 on average, Table 1) may be attributed to pine physio- logical functioning when the HW forest is inactive. These results represent an important advantage of ever- green ecosystems growing in warm climates, and refine earlier work proposing that the advantages of an extended photosynthetic period in evergreen trees tend to be offset by relatively low photosynthetic rates dur- ing the active season (Hollinger, 1992; Givnish, 2002). Relationship between the fluxes and factors that are changing over timescales commensurate with climate change At annual to interannual timescales, large variation in DNEE was observed (Fig. 5) that resonated most closely with variations in soil moisture content (reflecting drought effects that have already been discussed), and with variation in atmosphere CO2. It has long been hypothesized that rising CO2 may increase ecosystem water use efficiency (Field et al., 1995), representing a physiological adaptation to prevent excessive water loss when CO2 is less limiting (Katul et al., 2010; Med- lyn et al., 2011). Consistent with the results of Keenan et al. (2013), the WUEi and kE were correlated with CO2 in both PP (P < 0.0011 for both metrics) and HW (P = 0.08 and .14 for WUEi and kE) when drought years were excluded, with a more rapid increase in PP WUEi and kE for a given rise in CO2 (Fig. 9). We considered that changes in canopy structure could also affect the evolution of water use efficiency in PP. No significant trend between the allometric relationship (i.e. AS : AL/ h) and WUEi or kE was observed during nondrought years, though a significant relationship did emerge when only drought-free years following the ice storm disturbance were considered (P = 0.01 and 0.02 for WUEi and kE). Thus, variation in canopy structure may be a relatively more important determinant of water use efficiency in years following a disturbance event, though we caution that AS : AL/h and CO2 were also highly correlated during this time period. The trends in water use efficiency are likely also affected by the observed long-term increase in growing season D in non-drought years (P = 0.04, slope = 0.05 kPa/year). From Eqn. 9, it can be shown that the increases ca and D during the course of the study period should each promote an increase in kE on the order of 3 lmol/mol/ kPa, excluding drought years and assuming constant ci. Finally, we note that the combined rate of increase transpiration, increased GEP, or both under elevated CO2 remains an open question (Medlyn et al., 2011; Bar- ton et al., 2012; Ward et al., 2013). And finally, the capacity of PP to assimilate CO2 year-round (Fig. 8) represents an important carbon uptake advantage of evergreen trees in warm temperature climates, and should motivate future efforts to link regional pine for- est functioning to meteorological trends occurring when deciduous forests are dormant. Acknowledgements This research was sponsored by the Office of Science (BER), US Department of Energy (FG02-95SER62083), by the NSF Graduate Research Fellowship Program (DGE1106401), and the James B. Duke Fellowship Program. We thank Gaby Katul and Ram Oren for their long-term oversight of the Duke Forest flux towers, and for helpful feedback on an earlier version of this manuscript. The authors confirm that they have no interest or relationship, financial, or otherwise that might be perceived as influencing objectivity with respect to the work described in this study. References Abrams MD (2003) Where has all the white oak gone? BioScience, 53, 927–939. Albani M, Medvigy D, Hurtt GC, Moorcroft PR (2006) The contributions of land-use change, CO2 fertilization, and climate variability to the Eastern US carbon sink. Global Change Biology, 12, 2370–2390. Amiro BD, Barr AG, Barr JG et al. (2010) Ecosystem carbon dioxide fluxes after distur- bance in forests of North America. Journal of Geophyscial Research – Biogeosciences, 115, G00K02. Baldocchi DD (1997) Measuring and modeling carbon dioxide and water vapor exchange over a temperate broad-leaved forest during the 1995 summer drought. Plant, Cell and Environment, 20, 1108–1122. Baldocchi D (2008) Breathing of the terrestrial biosphere: lessons learned from a glo- bal network of carbon dioxide flux measurement systems. Australian Journal of Bot- any, 56, 1–26. Barton CVM, Duursma RA, Medlyn BE et al. (2012) Effects of elevated atmospheric [CO2] on instantaneous transpiration efficiency at leaf and canopy scales in Euca- lyptus saligna. Global Change Biology, 18, 585–595. Bauerle WL, Oren R, Way DA et al. (2012) Photoperiodic regulation of the seasonal pattern of photosynthetic capacity and the implications for carbon cycling. Proceed- ings of the National Academy of Sciences of the United States of America, 109, 8612– 8617. Beer C, Ciais P, Reichstein M et al. (2009) Temporal and among-site variability of inherent water use efficiency at the ecosystem level. Global Biogeochemical Cycles, 23, 1–13. Bolstad PV, Swank WT (1997) Cumulative impacts of landuse on water quality in a southern Appalachian watershed. Journal of the American Water Resources Associa- tion, 33, 519–533. Brown SL, Schroeder PE (1999) Spatial patterns of aboveground production and mor- tality of woody biomass for eastern US forests. Ecological Applications, 9, 968–980. Brzostek ER, Dragoni D, Schmid HP et al. (2014) Chronic water stress reduces tree growth and the carbon sink of deciduous hardwood forests. Global Change Biology, 20, 2531–2539. Ciais P, Reichstein M, Viovy N et al. (2005) Europe-wide reduction in primary pro- ductivity caused by the heat and drought in 2003. Nature, 437, 529–533. Clark KL, Gholz HL, Castro MS (2004) Carbon dynamics along a chronosequence of slash pine plantations in north Florida. Ecological Applications, 14, 1154–1171. Dai AG (2011) Characteristics and trends in various forms of the palmer drought severity index during 1900–2008. Journal of Geophysical Research-Atmospheres, 116, 1–26. Detto M, Montaldo N, Albertson JD, Mancini M, Katul G (2006) Soil Moisture and Veg- etation Controls on Evapotranspiration in a Heterogeneous Mediterranean Ecosystem on Sardinia. Water Resources Research, Italy. pp. 42. Dietze M, Vargos R, Richardson AD et al. (2011) Characterizing the performance of ecosystem models across time scales: a spectral analysis of the North American Carbon Program site-level synthesis. Jounral of Geophysical Research, 116, G04029. Dragoni D, Schmid HP, Wayson CA, Potter H, Grimmond CSB, Randolph JC (2011) Evidence of increased net ecosystem productivity associated with a longer vege- tated season in a deciduous forest in south-central Indiana, USA. Global Change Biology, 17, 886–897. Falge E, Baldocchi D, Olson R et al. (2001) Gap filling strategies for defensible annual sums of net ecosystem exchange. Agricultural and Forest Meteorology, 107, 43–69. Farquhar GD, Caemmerer SV, Berry JA (1980) A biochemical model of photosynthetic CO2 assimilation in leaves of C-3 species. Planta, 149, 78–90. Field CB, Jackson RB, Mooney HA (1995) Stomatal responses to increased CO2 - implications from the plant to the global scale. Plant Cell and Environment, 18, 1214–1225. Ford CR, Laseter SH, Swank WT, Vose JM (2011) Can forest management be used to sustain water-based ecosystem services in the face of climate change? Ecological Applications, 21, 2049–2067. Givnish TJ (2002) Adaptive significance of evergreen vs. deciduous leaves: solving the triple paradox. Silva Fennica, 36, 703–743. Hall B, Motzkin G, Foster DR, Syfert M, Burk J (2002) Three hundred years of forest and land-use change in Massachusetts, USA. Journal of Biogeography, 29, 1319–1335. Hollinger DY (1992) Leaf and simulated whole-canopy photosynthesis in 2 cooccur- ring tree species. Ecology, 73, 1–14. Hsieh CI, Katul GG, Chi TW (2000) An approximate analytical model for footprint estimation of scalar fluxes in thermally stratified atmospheric flows. Advances in Water Resources, 23, 765–772. Huntington TG (2006) Evidence for intensification of the global water cycle: review and synthesis. Journal of Hydrology, 319, 83–95. Jasechko S, Sharp ZD, Gibson JJ, Birks SJ, Yi Y, Fawcett PJ (2013) Terrestrial water fluxes dominated by transpiration. Nature, 496, 347–+. Juang JY, Katul GG, Siqueira MBS et al. (2006) Modeling nighttime ecosystem respira- tion from measured CO2 concentration and air temperature profiles using inverse methods. Journal of Geophysical Research-Atmospheres, 111, 1–16. Katul GG, Schieldge J, Hsieh C-I, Vidakovic B (1998) Skin temperature perturbations by surface layer turbulence above a grass surface. Water Resources Research, 34, 1265–1274. Katul GG, Vidakovic B, Albertson JD (2001) Estimating global and local scaling expo- nents in turbulent flows using wavelet transformations. Physics of Fluids, 13, 241–250. Katul GG, Palmroth S, Oren R (2009) Leaf stomatal responses to vapour pressure defi- cit under current and CO2-enriched atmosphere explained by the economics of gas exchange. Plant Cell and Environment, 32, 968–979. Katul G, Manzoni S, Palmroth S, Oren R (2010) A stomatal optimization theory to describe the effects of atmospheric CO2 on leaf photosynthesis and transpiration. Annals of Botany, 105, 431–442. Keenan TF, Hollinger DY, Bohrer G, Dragoni D, Munger JW, Schmid HP, Richardson AD (2013) Increase in forest water-use efficiency as atmospheric carbon dioxide concentrations rise. Nature, 499, 324–327. Keenan TF, Gray AJ, Friedl MA et al. (2014) Net carbon update has increased through warming-induced changes in temperatre forest phenology. Nature Climate Change, 4, 598–604. Lai CT, Katul G, Butnor J et al. (2002) Modelling the limits on the response of net car- bon exchange to fertilization in a south-eastern pine forest. Plant Cell and Environ- ment, 25, 1095–1119. Laseter SH, Ford CR, Vose JM, Swift LW (2012) Long-term temperature and precipita- tion trends at the Coweeta Hydrologic Laboratory, Otto, North Carolina, USA. Hydrology Research, 43, 890–901. Lasslop G, Reichstein M, Papale D et al. (2010) Separation of net ecosystem exchange into assimilation and respiration using a light response curve approach: critical issues and global evaluation. Global Change Biology, 16, 187–208. Law BE, Thornton PE, Irvine J, Anthoni PM, Van Tuyl S (2001) Carbon storage and fluxes in ponderonsa pine forests at different developmental stages. Global Change Biology, 7, 755–777. Manzoni S, Vico G, Katul G, Fay PA, Polley W, Palmroth S, Porporato A (2011) Opti- mizing stomatal conductance for maximum carbon gain under water stress: a meta-analysis across plant functional types and climates. Functional Ecology, 25, 456–467. Mccarthy HR, Oren R, Kim HS, Johnsen KH, Maier C, Pritchard SG, Davis MA (2006) Interaction of ice storms and management practices on current carbon sequestra- tion in forests with potential mitigation under future CO2 atmosphere. Journal of Geophysical Research-Atmospheres, 111, 1–10. Mccarthy HR, Oren R, Finzi AC, Ellsworth DS, Kim HS, Johnsen KH, Millar B (2007) Temporal dynamics and spatial variability in the enhancement of canopy leaf area under elevated atmospheric CO2. Global Change Biology, 13, 2479–2497. McEwan RW, Dyer JM, Pederson N (2011) Multiple interacting ecosystem drivers: toward an encompassing hypothesis of oak forest dynamics across eastern North America. Ecography, 34, 244–256. Medlyn BE, Duursma RA, Eamus D et al. (2011) Reconciling the optimal and empiri- cal approaches to modelling stomatal conductance. Global Change Biology, 17, 2134–2144. Myneni RB, Hoffman S, Knyazikhin Y et al. (2002) Global products of vegetation leaf area and fraction absorbed PAR from year one of MODIS data. Remote Sensing of Environment, 83, 214–231. Noormets A, Gavazzi MJ, Mcnulty SG, Domec JC, Sun G, King JS, Chen JQ (2010) Response of carbon fluxes to drought in a coastal plain loblolly pine forest. Global Change Biology, 16, 272–287. Novick KA, Stoy PC, Katul GG, Ellsworth DS, Siqueira MBS, Juang J, Oren R (2004) Carbon dioxide and water vapor exchange in a warm temperate grassland. Oecolo- gia, 138, 259–274. Novick KA, Oren R, Stoy PC, Siqueira MS, Katul G (2009a) Nocturnal evapotranspi- ration in eddy covariance records from three co-located ecosystems in the Southeastern US: the effect of gapfilling methods on estimates of annual fluxes. Agricultural and Forest Meteorology, 149, 1491–1504. Novick K, Oren R, Stoy P, Juang J-Y, Siqueira M, Katul GG (2009b) The relationship between reference canopy conductance and simplified hydraulic architecture. Advances in Water Resources, 32, 808–819. Novick KA, Brantley ST, Miniat CF, Walker J, Vose J (2014) Inferring the contribution of advection to total ecosystem scalar fluxes over a tall forest in complex terrain. Agricultural and Forest Meteorology, 185, 1–13. Odum EP (1969) The strategy of ecosystem development. Science, 164, 262–270. Oishi AC, Oren R, Stoy PC (2008) Estimating components of forest evapotranspira- tion: a footprint approach for scaling sap flux measurements. Agricultural and For- est Meteorology, 148, 1719–1732. Oishi AC, Oren R, Novick KA, Palmroth S, Katul GG (2010) Interannual invariability of forest evapotranspiration and its consequence to water flow downstream. Eco- systems, 13, 421–436. Oishi AC, Palmroth S, Butnor JR, Johnsen KH, Oren R (2013) Spatial and temporal variability of soil CO2 efflux in three proximate temperate forest ecosystems. Agri- cultural and Forest Meteorology, 171, 256–269. Oosting HJ (1942) An ecological analysis of the plant communities of Piedmont, North Carolina. American Midland Naturalist, 28, 1–126. Palmroth S, Maier CA, Mccarthy HR et al. (2005) Contrasting responses to drought of forest floor CO2 efflux in a Loblolly pine plantation and a nearby Oak-Hickory for- est. Global Change Biology, 11, 421–434. Pataki DE, Oren R (2003) Species differences in stomatal control of water loss at the canopy scale in a mature bottomland deciduous forest. Advances in Water Resources, 26, 1267–1278. Powell TL, Gholz HL, Clark KL, Starr G, Cropper WP, Martin TA (2008) Carbon exchange of a mature, naturally regenerated pine forest in north Flordia. Global Change Biology, 14, 2523–2538. Ramankutty N, Heller E, Rhemtulla J (2011) Prevailing myths about agricultural abandonment and forest regrowth in the United States. Annals of the Association of American Geographers, 100, 502–512. Reichstein M, Falge E, Baldocchi D et al. (2005) On the separation of net ecosystem exchange into assimilation and ecosystem respiration: review and improved algo- rithm. Global Change Biology, 11, 1424–1439. Sch€afer KVR, Oren R, Lai CT, Katul GG (2002) Hydrologic balance in an intact tem- perate forest ecosystem under ambient and elevated atmospheric CO2 concentra- tion. Global Change Biology, 8, 895–911. Sch€afer KVR, Oren R, Ellsworth DS et al. (2003) Exposure to an enriched CO2 atmo- sphere alters carbon assimilation and allocation in a pine forest ecosystem. Global Change Biology, 9, 1378–1400. Schwalm CR, Williams CA, Schaefer K et al. (2010) Assimilation exceeds respiration sensitivity to drought: a FLUXNET synthesis. Global Change Biology, 16, 657–670. Stanturf JA, Wade DD, Waldrop TA, Kennard DK, Achtemeier GL (2002) Fire in southern forest landscapes. In: Southern Forest Resource Assessment (eds Wear DN, Greis JG), pp. 607–630. USDA Forest Service - Southern Research Station, Ashe- ville, NC, USA. Stoy P, Katul GG, Siqueira M et al. (2005) Variability in net ecosystem exchange from hourly to inter-annual time scales at adjacent pine and hardwood forests: a wave- let analysis. Tree Physiology, 25, 887–890. Stoy PC, Katul GG, Siqueira MBS et al. (2006a) Separating the effects of climate and vegetation on evapotranspiration along a successional chronosequence in the southeastern US. Global Change Biology, 12, 2115–2135. Stoy PC, Katul GG, Siqueira MBS, Juang JY, Novick KA, Uebelherr JM, Oren R (2006b) An evaluation of models for partitioning eddy covariance-measured net ecosystem exchange into photosynthesis and respiration. Agricultural and Forest Meteorology, 141, 2–18. Stoy PC, Katul GG, Siqueira MBS et al. (2008) Role of vegetation in determining car- bon sequestration along ecological succession in the southeastern United States. Global Change Biology, 14, 1409–1427. Stoy PC, Richardson AD, Baldocchi DD et al. (2009) Biosphere-atmosphere exchange of CO2 in relation to climate: a cross-biome analysis across multiple time scales. Biogeosciences, 6, 2297–2312. Sun G, Noormets A, Gavazzi MJ et al. (2010) Energy and water balance of two con- trasting loblolly pine plantations on the lower coastal plain of North Carolina, USA. Forest Ecology and Management, 259, 1299–1310. Tang J, Luyssaert S, Richardson AD, Kutsch W, Jannssens IA (2014) Steeper declines in forest photosynthesis than respiration explain age-driven decreases in forest growth. Proceedings of the National Academy of Sciences of the United States of America, 111, 8856–8860. Torrence C, Compo GP (1998) A practical guide to wavelet analysis. Bulletin of the American Meteorological Society, 79, 61–78. Van Gorsel E, Delpierre N, Leuning R et al. (2009) Estimating nocturnal ecosystem respiration from the vertical turbulent flux and change in storage of CO2. Agricul- tural and Forest Meteorology, 149, 1919–1930. Ward EJ, Oren R, Bell DM, Clark JS, McCarthy HR, Kim HS, Domec JC (2013) The effects of elevated CO2 and nitrogen fertilization on stomatal conductance esti- mated from 11 years of scaled sap flux measurements at Duke FACE. Tree Physiol- ogy, 33, 135–151. Wear DN, Greis JG (2012) The Southern Forest Futures Project: Summary Report. (ed Report GT) USDA Forest Service, Asheville, NC, USA. Webb EK, Pearman GI, Leuning R (1980) Correction of flux measurements for density effects due to heat and water-vapor transfer. Quarterly Journal of the Royal Meteoro- logical Society, 106, 85–100. Whelan A, Mitchell R, Staudhammer CL, Starr G (2013) The cyclic occurrence of fire and its role in carbon dynamics along an edaphic moisture gradient in longleaf pine ecosystems. PLoS ONE, 8, e54045. Wilson KB, Baldocchi DD (2001) Comparing independent estimates of carbon dioxide exchange over five years at a deciduous forest in the southern United States. Jour- nal of Geophysical Research, 106, 34167–34178. Wullschleger SD, Mclaughlin SB, Ayres MP (2004) High-resolution analysis of stem increment and sap flow for loblolly pine trees attacked by southern pine beetle. Canadian Journal of Forest Research-Revue Canadienne De Recherche Forestiere, 34, 2387–2393. Supporting Information Additional Supporting Information may be found in the online version of this article: Table S1. The previously published estimates of the annual Net Ecosystem Exchange of CO2(NEE) used to create Fig. 1 in the main text are shown in Table S1. These sites were ini- tially selected from among Ameriflux sites with at least five site years of available data, though sites were included in the table even if less than 5 years of reported estimates of NEE were recovered in the literature.