Multi‐year carbon budget of a mature commercial short rotation coppice willow plantation

Energy derived from second generation perennial energy crops is projected to play an increasingly important role in the decarbonization of the energy sector. Such energy crops are expected to deliver net greenhouse gas emissions reductions through fossil fuel displacement and have potential for increasing soil carbon (C) storage. Despite this, few empirical studies have quantified the ecosystem‐level C balance of energy crops and the evidence base to inform energy policy remains limited. Here, the temporal dynamics and magnitude of net ecosystem carbon dioxide (CO2) exchange (NEE) were quantified at a mature short rotation coppice (SRC) willow plantation in Lincolnshire, United Kingdom, under commercial growing conditions. Eddy covariance flux observations of NEE were performed over a four‐year production cycle and combined with biomass yield data to estimate the net ecosystem carbon balance (NECB) of the SRC. The magnitude of annual NEE ranged from −147 ± 70 to −502 ± 84 g CO2‐C m−2 year−1 with the magnitude of annual CO2 capture increasing over the production cycle. Defoliation during an unexpected outbreak of willow leaf beetle impacted gross ecosystem production, ecosystem respiration, and net ecosystem exchange during the second growth season. The NECB was −87 ± 303 g CO2‐C m−2 for the complete production cycle after accounting for C export at harvest (1,183 g C m−2), and was approximately CO2‐C neutral (−21 g CO2‐C m−2 year−1) when annualized. The results of this study are consistent with studies of soil organic C which have shown limited changes following conversion to SRC willow. In the context of global decarbonization, the study indicates that the primary benefit of SRC willow production at the site is through displacement of fossil fuel emissions.

lower greenhouse gas (GHG) intensity relative to fossil fuels (Creutzig et al., 2015). In the future, energy crops with C capture and storage (BECCS) might also prove viable for delivering negative carbon dioxide (CO 2 ) emissions required to meet national and international GHG emissions reductions targets (Energy Technologies Institute, 2015;Rosen, 2018). Whilst delivering genuine net GHG emissions reductions require a full assessment of all biogenic GHG emissions from cultivation to energy production (Thornley, Gilbert, Shackley, & Hammond, 2015), for low input perennial (second generation) bioenergy crops, it is generally accepted that the major determinant of the site-level GHG balance relates to changes in soil organic carbon (SOC) following land use change to bioenergy (Richards et al., 2017), which in turn is modified by site-specific management practices (Qin et al., 2018).
Short rotation coppice (SRC) willow (Salix spp.) is a second generation (e.g. non-food) energy crop that has been planted over a wide geographical range in the UK and northern areas of continental Europe (Don et al., 2012;Karp & Shield, 2008;Rowe et al., 2016). SRC exploits fast growing woody perennial species that regrow rapidly after harvesting at short (e.g. 2-to 4-year) intervals (Dimitriou, Mola-Yudego, Aronsson, & Eriksson, 2012). Managing willow as SRC aims to maximize biomass yields with limited input requirements, resulting in high energetic returns (Rowe et al., 2011). Compared to annual crops, there is potential to enhance soil C storage under perennial SRC as a function of high productivity, regular C inputs from litter and fine root turnover, and minimal tillage requirements (Karp & Shield, 2008). However, it remains uncertain whether perennial cropping systems capture enough atmospheric CO 2 to totally offset GHG emissions produced over the complete bioenergy life cycle (Robertson et al., 2017;Rowe et al., 2016).
Detecting changes in SOC under SRC is challenging. The rate of change in SOC is slow and many decades may be required to reach a new equilibrium following land use change (Nemo et al., 2017;Smith, 2004). In contrast to methods employing direct measurements of changes in soil C (Rowe et al., 2016), tracking year on year C fluxes from photosynthetic uptake versus plant and soil respiration can provide a means of understanding the overall trajectory of a system before changes in SOC become detectable. This requires multiple years of data to derive consensus, as the net ecosystem exchange (NEE) of CO 2 is highly dynamic across multiple timescales (e.g. Stoy et al., 2005), responding to time since establishment and harvest, management interventions and climatic forcing. The latter drives plant and phenology and productivity, influences site-level management decisions (Caslin, Finnan, Johnston, McCracken, & Walsh, 2010) and controls decomposition processes in litter and soils (Davidson & Janssens, 2006).
The most effective way of measuring the short term dynamics in net ecosystem CO 2 exchange is the eddy covariance (EC) technique which quantifies turbulent surface-atmosphere exchanges of energy and mass (Baldocchi, ) at ecosystem scale. In contrast to longer term perspectives provided by stock-based C accounting (Agostini, Gregory, & Richter, 2015;Berhongaray, Verlinden, Broeckx, Janssens, & Ceulemans, 2017;Dimitriou et al., 2012;Rowe et al., 2016), EC provides direct, high frequency, and quasi-continuous observations of NEE from a flux footprint representative of tens to hundreds of meters (Arriga et al., 2017;Baldocchi, 2003). To date, a limited number of EC studies have focused on the early stages of land use change to SRC willow (Grelle, Aronsson, Weslien, Klemedtsson, & Lindroth, 2007;Harris et al., 2017;Harris, Spake, & Taylor, 2015); however, commercial SRC plantations have suggested economic life spans of 20-30 years (Don et al., 2012) and empirical data obtained at young SRC sites may not reflect the dynamics of mature plantations (Walter, Don, & Flessa, 2015). In order to inform bioenergy policy, observational data from across a wider range of SRC willow production systems are required for model development and testing (e.g. Dondini et al., 2016), and as inputs to bioenergy life cycle analyses (e.g. Agostini et al., 2015) and geographical assessments (e.g. Hammar, Ericsson, Sundberg, & Hansson, 2014;Pogson et al., 2016).
This study reports multi-year flux observations (January 2014 to November 2017) of NEE obtained at a mature commercial SRC willow plantation in Lincolnshire, UK. The SRC is representative of an age class of SRC willow plantations established during the early 2000s. Flux measurements are reported for a near-complete, 4-year production cycle. An unexpected outbreak of willow beetle (Phratora vulgatissima) and leaf defoliation during the second and to a lesser degree the third year of the observation period provided an opportunity to follow the effects of insect defoliation on ecosystem CO 2 fluxes. Outbreaks of willow leaf beetle and associated reductions in biomass yield are widely recognized as a major economic pest for commercial SRC production (Bjorkman & Eklund, 2006;Bjorkman, Hoglund, Eklund, & Larsson, 2000). Despite such recognition, it is unclear how frequently such outbreaks occur in commercial SRC production systems, or how such events might impact the sustainability and profitability of bioenergy production if large scale plantings required to meet bioenergy targets are realized.
The overall aim of the study was to characterize the magnitude and dynamics of CO 2 exchange at the mature SRC willow plantation. The specific objectives of this study were: (i) to characterize seasonal and inter-annual variability in net ecosystem CO 2 exchange (NEE) and its component fluxes and (ii) to quantify the net ecosystem CO 2 -C budget (NECB) of the mature SRC plantation.

| Site description
The flux measurement site is a 9.4 ha commercial SRC willow plantation located approximately 10 km north of Lincoln, England, United Kingdom. The climate is temperate maritime (Cfb, Köppen, Volken, & Brönnimann, 2011), characterized by cool summers, mild winters, and a thermal growing season from April to October. The closest UK Met Office station with complete long-term records (RAF Waddington, 53° 10' 30'' N; 0° 31' 15.6'' W 68 m absl) is located approximately 15 km south of the observation site. Mean (±standard deviation, SD) annual air temperature was 9.8° ± 0.7°C for the 1981-2010 period. July (16.9 ± 1.2°C) is the warmest month and January (4.0° ± 1.6°C) is the coldest. The mean annual (1981 to 2010) precipitation was 614 ± 93.5 mm/ year. Precipitation is distributed approximately evenly over the annual cycle. South westerly winds predominate at this location.

| Land management
The SRC willow plantation was established on former arable land in 2000. The regional crop rotation prior to conversion was winter wheat (Triticum aestivum) with oilseed rape (Brasica napus) and spring barley (Hordeum vulgare L.) as break crops. During conversion, the field was cultivated (plough, power harrow, and flexi-tyne) and sprayed with nonselective weedkiller (Glythosate). A mixture of five commercial varieties of Salix viminalis were randomly planted as an integrated mix (block planted) across the field with the aim of minimizing the impacts of pests and diseases. Approximately 15,000 willows ha −1 were planted in double rows with within and between row spacing of 0.75 m and 1.5 m, correspondingly. The plantation was cut back (coppiced) in 2001 to encourage branching and is typically managed using three-year rotation cycles. Stool survival in 2018 based on counts from 20 random 4 by 5 m plots was 80%. During this study, an outbreak of willow leaf beetle (Phratora vulgatissima) resulted in severe defoliation and subsequently low biomass accrual during the second year of the rotation cycle (see below). In response to this, the land manager delayed the harvest by one year in order to maximize biomass yield and economic return. Biomass was harvested on 2017-09-27 using specialized equipment cutting the willow directly to chip. Other than harvesting, the plantation has been minimally managed with the most recent input in 2011, when micronutrients (Fibrophos 0.5 Mg/ha and lime) and compost wood waste (15 Mg/ha) where applied following harvest. The site is rain fed with no requirement for irrigation.

| Crop development
The SRC plantation crop was characterized by changes in aboveground biomass over the production cycle. At the start of the observation period in January 2014, willow stools were dormant and areas of bare soil and litter remaining from the previous harvest were present across the site. In the first year of the rotation (2014), bud burst occurred during April with notable new shoot extension and leaf area by May. Full canopy closure was observed by July 2014. The SRC stand attained a mean height of 2.0 ± 0.02 m by autumn 2014. In subsequent years, new shoots developed at the apex of the preceding season's woody growth resulting in earlier leaf displays compared to 2014. Mean stand height had increased to 2.9 m ± 0.2 m by June 2015 and was 3.5 ± 1.2 m at the time of harvest in 2017. The active growth phase lasted from April to October in all years, but bud burst and leaf out occurred notably earlier during warmer than average conditions in spring 2017 compared to other years of the rotation. Higher than average air temperatures during autumn 2016 also resulted in a later onset of leaf senescence and fall relative to other years.

| Willow leaf beetle
A P. vulgatissima population was observed in the plantation during all years of the observation period. Leaf damage caused by P. vulgatissima herbivory was most severe during the second year (2015) of the rotation cycle. As damage from a severe pest outbreak was not expected no quantitative measures of beetle density were made; however, in 2015 the population of P. vulgatissima reached a density at which progressive, widespread and severe leaf damage was observed on all willow trees between July and September ( Figure 1). The majority of willow leaves were brown and skeletonized by July ( Figure 1b) and did not recover before the end of the growing season ( Figure 1c). In 2014, 2016, and 2017, herbivory by P. vulgatissima was more selective and localized. Herbivory in 2016 and 2017 mainly affected young leaves (between July-August) on certain individual willow stools, possibly caused by preferential feeding on certain willow varieties. Observed leaf damage in 2016 and 2017 was mostly replaced by compensatory growth by September.

| Flux instrumentation
Sensible and latent heat fluxes (LE and H, respectively) and NEE were monitored using an open-path EC system.
Flux observations commenced on 2014-01-09 during the winter dormant period following the previous harvest (on 2013-10-31) and ended on 2017-11-26. A Solent R3-50 sonic anemometer (Gill Instruments Ltd., Lymington, UK) was used to measure the three components of atmospheric turbulence (u, v, w; m/s) and sonic temperature (T sonic ; °C). An LI7500 infrared H 2 O/CO 2 gas analyser (IRGA; LI-COR Biosciences, Lincoln, Nebraska) was used to measure the atmospheric mass density of vapor (g H 2 O/m 3 ) and CO 2 (mg CO 2 /m 3 ) and barometric pressure (P air ; kPa). Air temperature (T a ; °C) and relative humidity (RH; %) were measured at a height of 2 m using a HC2A-S probe (Rotronic, Bassersdorf, Switzerland). EC sensors were sampled and logged at 20 Hz using a CR3000 Micrologger (Campbell Scientific Inc., Logan Utah). The EC system was installed at the center of the plantation to maximize the available fetch under prevailing wind conditions. The available fetch was limited to 110 m to the north of the EC tower, extending to 210 m to the southwest. EC sensors were mounted on an extendible pneumatic mast (Clarke Masts Ltd., Binstead, UK), the height of which was increased several times per growing season, to maintain the EC measurement height at a minimum of 2 m above the willow canopy.

| Ancillary measurements
A range of micrometeorological measurements was made from a separate scaffold tower located at the northern end of the plantation approximately 110 m from the EC mast. The net radiation (R net ; W/m 2 ) and its incoming and outgoing short-and long-wave components (SW in , SW out , LW in , and LW out , respectively; W/m 2 ) were measured above the canopy using a CNR1 net radiometer (Kipp & Zonen BV, Delft, The Netherlands). Secondary measurements of air temperature (T air , °C) and relative humidity (RH, %) were made at 7 m above the soil surface using a HC2A-S probe (Rotronic, Bassersdorf, Switzerland).
Soil physics were measured close to the scaffold tower. Soil heat fluxes (G; W/m 2 ) were monitored using two HFP01 soil heat flux (Hukseflux BV, BV, Delft, The Netherlands) plates installed 0.03 m below the soil surface. Soil temperature (T soil ; °C) was measured at a depth of 0.05 m using two PT107 soil thermocouples (Campbell Scientific Inc.). Soil volumetric water content (VWC; m 3 /m 3 ) was measured using two CS616 time domain reflectometers installed vertically to measure the VWC of the top 0.3 m of the soil profile (Campbell Scientific Inc.). Precipitation was measured using a tipping bucket rain gauge (P; 0.5 mm sensitivity; Observator Instruments BV, The Netherlands) installed in an open area of low grass approximately 10 m north of the plantation. These sensors were scanned at 0.1 Hz and logged as 30-minute means (sums for precipitation) using a CR1000 Datalogger (Campbell Scientific Inc.).

| Flux data handing
Thirty-minute flux densities (hereafter fluxes) of sensible and latent heat (LE and H) and net ecosystem CO 2 exchange (NEE) were computed from the raw EC data using EddyPRO ® Flux Calculation Software (LI-COR Biosciences, Lincoln, Nebraska; Fratini & Mauder, 2014). Raw EC data were screened for statistical outliers  and other physically implausible values (Vickers & Mahrt, 1997). Sonic anemometer data were rotated using a two-dimensional coordinate rotation procedure (Wilczak, Oncley, & Stage, 2001) and corrected for imperfect cosine response (Nakai, Van Der Molen, Gash, & Kodama, 2006). Time lags between the vertical wind speed and concentration measurements were removed using a cross-correlation procedure. Uncorrected fluxes were calculated as the mean covariance between the vertical wind speed (w) and the respective atmospheric scalar using 30minute block averages (Baldocchi, 2003). Fluxes were corrected for high (Moncrieff et al., 1997) and low frequency cospectral attenuation (Moncrieff, Clement, Finnigan, & Meyers, 2004). H fluxes were corrected for the influence of atmospheric humidity (Schotanus, Nieuwstadt, & Bruin, 1983). LE then CO 2 fluxes were adjusted for fluctuations in atmospheric density (Webb, Pearman, & Leuning, 1980). Random uncertainties for 30-minute flux observations related to sampling error were estimated as standard deviations derived from a variance of covariance approach (Finkelstein & Sims, 2001). No CO 2 profile measurements were available to estimate the CO 2 storage term. CO 2 storage was assumed negligible at the low observation height and NEE was assumed equal to the turbulent CO 2 flux. The micrometeorological sign convention is adopted where positive values represent fluxes from ecosystem to atmosphere and negatives describe the reverse.

| Quality control
Quality control (QC) procedures were applied to ensure only high quality turbulent flux data were retained for analysis. Thirty-minute flux data were screened for statistical outliers using the median absolute deviation approach (Sachs, 2013) following recommendations in Papale et al. (2006). Fluxes were also excluded when: the results of the stationarity (steady-state) test result deviated by more than 100% (Foken et al., 2004); when the automatic gain control (AGC) of the LI7500 was >10% above its baseline (Ruppert, Mauder, Thomas, & Lüers, 2006); and when fluxes were outside the range −200 < H > 500 W/ m 2 , −100 < LE > 600 W/m 2 , and −50 < NEE > 30 μmol CO 2 m −2 s −1 . Periods of low turbulent mixing were identified using a friction velocity (u*) threshold approach (Papale et al., 2006;Reichstein, Moffat, Maria, Wutzler, & Sickel, 2016), and CO 2 fluxes were excluded when u* < 0.14 m/s. As the site has a limited homogeneous fetch, the spatial representativeness of measured fluxes was assessed using a two-dimensional implementation of the Kormann and Meixner (2001) flux footprint model (Neftel, Spirig, & Ammann, 2008). Fluxes were considered representative and retained for analysis when the footprint model indicated >80% of the flux originated within the boundaries of the plantation. EC data retention after system downtime and the application of QC are summarized in Table 1. Overall energy balance closure (EBC; Figure 2) for the EC system at this location (slope of 0.96) was towards the high end of the 0.55-0.99 range attained for EC sites, globally (Leuning, Gorsel, Massman, & Isaac, 2012;Stoy et al., 2013;Wilson et al., 2002).  2014  78  68  72  71  48  58  66  50  57   2015  88  73  79  72  52  61  65  52  58   2016  88  71  79  76  49  61  70  50  59   2017  81  71  75  71  49  59  66  50  58   All years  84  70  76  72  49  60  67  51  58 Values represent the percentage of all potentially available eddy covariance observations that were retained for analysis after system outages and data rejection during quality control (see main text for details). 2.9 | Light response of CO 2 exchange An empirical modeling approach was used to quantify and compare seasonal and inter-seasonal variations in photosynthetic activity and respiration rates during the months of the growing seasons (April to October). All observations of NEE that passed QC during each calendar month were used to parameterize a modified Michaelis-Menton equation as a function of SW in , using: (Carrara et al., 2003;Falge et al., 2001) where α (μmol CO 2 J −1 ) is the ecosystem quantum yield, GPP 900 (μmol CO 2 m −2 s −1 ) is the rate of photosynthesis when SW in is 900 W/m 2 and R m (μmol CO 2 m −2 s −1 ) is the mean respiration rate. The model was selected over other forms as it allows the maximum rate of photosynthesis to be compared at a standardized level of SW in .

| Gap-filling and flux partitioning
Gap-filling of flux data and the partitioning of NEE into estimates of gross ecosystem production (GEP) and total ecosystem respiration were performed using the REddyProc Package version 0.8-2/r14 (Reichstein et al., 2016) for the R statistical Language (R Core Team, 2017). Data gap-filling of H, LE, and NEE was performed using the marginal distribution sampling (MDS) approach (Reichstein et al., 2005). To enable an annual sum to be computed for 2014, missing data at the start of January was gap-filled using EC data collected from 2014-01-09 onwards. Details of the MDS and flux partitioning algorithms have been described in detail and evaluated elsewhere (Desai et al., 2008;Moffat et al., 2007;Papale et al., 2006;Reichstein et al., 2005) and are not repeated here. Gaps in prognostic micrometeorological variables (SW in , T air , and vapor pressure deficit, VPD) required for MDS gap-filling were filled using observations obtained at weather stations located in adjacent fields. T air was used as the driving temperature for flux partitioning as a number of data gaps were present in the T soil record. Uncertainty for individual gap-filled fluxes was estimated as the standard deviation of the observations averaged to fill data gaps (Reichstein et al., 2016). No uncertainties were estimated for GEP and TER as the partitioned CO 2 fluxes represent modeled quantities.

| Biomass yield and net ecosystem C balance
In agricultural systems, time-integrated NEE must be combined with data on harvested yields to determine whether a site is a net CO 2 -C source or sink. C exported as harvested biomass (C export , g C m −2 ) was estimated from farm records. Exported biomass was converted to units of C using a factor of 0.49 (Harris et al., 2017), and assuming a biomass moisture content of 50% at harvest. It was assumed that all exported C would be converted back to atmospheric CO 2 during combustion for bioenergy generation. On the basis of past research at the observation site, it was assumed that emissions/removals of methane (CH 4 ) and nitrous oxide (N 2 O) Finch et al., 2014) were negligible, and that the net ecosystem C balance (NECB, g CO 2 -C m −2 year −1 ) could be approximated as: where all variables have been defined above.

| Environmental conditions
The SRC willow plantation was characterized by seasonal and inter-seasonal variation in environmental conditions during the study period ( Figure 3). Mean annual T air was higher than the 30-year average (of 9.8 ± 0.7°C) during all years of the study ( Table 2). The majority of months (36 out of 48) experienced mean temperatures that were warmer than long-term climate averages ( Fig S1). The spring and early summer months of 2014 and 2017, and the late summer period (August and September) of 2016 were over one standard deviation (SD) warmer than long-term monthly means. December 2015 experienced a mean T air that was more than 5° C higher than normal. Annual precipitation sums (Table 2) were within one SD of the 1981 to 2010 mean rainfall average (of 614 ± 93.5 mm/year). The year 2015 was notably drier than other years due to lower rainfall during the first half of the year. The annual cycle of volumetric soil water content (VWC) showed broadly similar seasonal patterns across years. Absolute VWC ranged from 26% (September 2017) to 51% (June 2014).

| Diurnal and seasonal pattern
Fingerprint plots showing diurnal and seasonal variations in observed and gap-filled NEE are presented in Figure 4. The general pattern was positive NEE during nocturnal periods and outside the growing season, interspersed by periods of negative NEE during the daytime of the main growth phases. In 2014, daytime NEE started to become more negative from late May onwards as the photosynthetic activity of the vegetation increased as new stems developed and leaf area increased. The highest net uptake rates of 2014 were observed after full canopy closure in July. By contrast, earlier leaf out on the stems of preceding seasons' (1) (2) NECB = NEE − C export , growth in later years was associated with more negative CO 2 uptake during spring and earlier seasonal maxima. In 2015, daytime NEE became progressively more positive as defoliation by P. vulgatissima progressed from July onwards. The most negative NEE of the study period was observed during June 2016 when mean (±standard error of the mean) midday (10:00 to 14:00 UTC) fluxes were −24.6 ± 0.7 μmol CO 2 m −2 s −1 . These maximum rates of net CO 2 uptake occurred before any notable leaf damage by willow beetle had been observed in the plantation during that season. The most negative early-and late-season values were measured during the warm conditions of spring and autumn of 2017 and 2016, respectively. The net CO 2 uptake period ended abruptly with the biomass harvest in autumn 2017. Monthly mean nocturnal NEE showed clear seasonal trends, attaining maximum values (mean ± SD) of 7.3 ± 2 and 7.2 ± 3 μmol CO 2 m −2 s −1 during warm conditions of July in 2014 and 2016, respectively.

| NEE-light responses
Parameters derived from NEE-light response curves (Equation ; see example light response in Figure 5) for each month of the growing season revealed seasonal and between-year differences in ecosystem photosynthesis and respiration rates (Figure 6). SW in explained between 50% and 93% of observed variation in monthly NEE ( Figure  6d). Monthly quantum yield (α) values (Figure 6a) were between 0.025 ± 0.003 μmol CO 2 J −1 (October 2014) and 0.18 ± 0.009 μmol CO 2 J −1 (June 2015). Maximum rates of photosynthesis at SW in of 900 W/m 2 (GEP 900 ) ranged from 4.74 ± 0.14 μmol CO 2 m −2 s −1 in April 2014, to a maximum of 38.4 ± 0.6 μmol CO 2 m −2 s −1 in June 2016 (Figure 6b). Monthly mean respiration rates (R m ) were positively correlated with GEP 900 (R m = 0.15 * GEP 900 + 1.7 μmol CO 2 m −2 s −1 , r 2 = 0.79, p < 0.05) and ranged from 2.31 ± 0.1 in April 2014 to 7.58 ± 0.2 μmol CO 2 m −2 s −1 in July 2014 (Figure 6c). All three parameters were lower during the early growing season of 2014 compared to later years, clearly reflecting the early developmental stage of the plantation. In 2015, parameter estimates and net CO 2 uptake rates (see example in Figure 5) were notably lower than for other years during the main P. vulgatissima outbreak between July and September. Other notable interannual variations in parameter estimates reflected the influence of warm meteorological conditions on photosynthesis and respiration rates, such as the higher parameter values estimated during warm spring conditions during April 2017, as well as for the warm autumn period during September and October in 2016 ( Figure 6).

| Ecosystem photosynthesis and respiration
The seasonal growth and decay of flux partitioned estimates of GEP and TER and gap-filled NEE are presented as daily sums in Figure 7. The same data are presented as accumulated monthly totals in Figure S2. The magnitude of peak season GEP (monthly mean ± SD) ranged from 12.2 ± 2 g CO 2 -C m −2 day −1 in July 2017 to 14.2 ± 3 g CO 2 -C m −2 day −1 during July 2016 (Figure 7b). Total monthly GEP ( Figure S2b 439 g CO 2 -C m −2 month −1 (July 2014). Mean (±SD) daily TER (Figure 7c) was in the range 0.9 (February 2015) to 8.5 ± 0.6 g CO 2 -C m −2 day −1 (July 2014). Total monthly TER ( Figure S2c) was between 26 g CO 2 -C m −2 month −1 (February 2015) and 269 g CO 2 -C m −2 month −1 (July 2014). In spring and early summer of 2014, GEP and TER were both lower than during the corresponding periods of later years, with the magnitude of both fluxes increasing following full canopy closure in July 2014. GEP and TER both showed strong reductions during the P. vulgatissima outbreak in 2015. For example, total monthly GEP in August 2015 was only ca. 38% of values estimated for corresponding month of other years. The magnitude of TER was also lower during the outbreak compared with respective months of other years, but the relative reductions in TER were of a lesser magnitude when compared to relative reductions in GEP.

| Net ecosystem exchange
Estimates of daily NEE ranged from −12.5 ± 0.6 g CO 2 -C m −2 day −1 in June 2016 to 5.0 ± 0.2 g CO 2 -C m −2 day −1 (Figure 7a). The largest net daily losses were observed immediately after harvest in 2017 when TER remained high but GEP was terminated with the removal of the active aboveground biomass. The site switched from a net CO 2 sink to a net source on a number of days when GEP was reduced under low light conditions (Figure 7a). The number of days with net CO 2 uptake were 122 and 120 in 2014 and 2015, increasing to 196 and 172 for 2016 and 2017, correspondingly. Monthly NEE ranged from a net gain of −195 ± 14 g CO 2 -C m −2 month −1 in June 2016 to a net loss of 65 ± 4 g CO 2 -C m −2 month −1 during October 2014 ( Figure S2a).
Cumulative plots (Figure 8) of daily NEE show net CO 2 -C losses were highest during the early part of 2014 and the site did not become a cumulative net sink until the canopy closed in July of that year. Net CO 2 -C uptake started earlier during spring of all subsequent years, with the highest early season CO 2 -C uptake observed during spring in 2017. Net CO 2 -C gains were observed until September during 2014 and until harvest in 2017, and continued into October during the warm autumn period of 2016. By contrast, net CO 2 -C accumulation had ended by August in 2015, largely in response to low rates of GEP during the P. vulgatissima outbreak. In 2016, lower net CO 2 -C uptake during spring (compared to 2017) was compensated by higher net C gains during summer and autumn.

| Net ecosystem carbon balance
C balance terms for each year of the observation period are summarized in Table 3. Total accumulated NEE during the complete study period was −1,268 ± 303 g CO 2 -C m −2 . Annual GEP and TER were lowest during 2015 and highest during 2016. TER/GEP ratios showed that around 90% of CO 2 -C assimilated during photosynthesis was respired during 2014 and 2015, decreasing to 75% during 2016 and 2017. The least negative annual NEE was estimated for the first year of the rotation in 2014. In 2015, annual NEE was of similar magnitude to the previous year, with all of the observed net C uptake in that year occurring before the onset of insect defoliation. The magnitude of annual NEE for 2016, and for the January to November period of 2017 was more than three times more negative than for the two preceding years (Table 3). Dry biomass yield at harvest in September 2017 was 24 Mg/ha based on farm records, or 6 Mg/ha when annualized using four growing seasons. C export at harvest was 1,183 g C ha −1 for the whole production cycle, or 296 C m −2 year −1 when annualized. On this basis, the NECB indicated the plantation functioned as a small net sink of −85 ± 303 g CO 2 -C m −2 over the production cycle, or approximately −21 g CO 2 -C m −2 year −1 (e.g. approximately CO 2 -C neutral) when annualized.

| DISCUSSION
For the first time, the C balance of all years of a commercial rotation cycle at a mature SRC willow plantation in the UK has been observed using the EC technique. The plantation accumulated −12.7 ± 3 Mg CO 2 -C m −2 ha −1 as NEE during the production cycle, with 11.8 Mg C ha −1 removed as harvested biomass. The annualized dry biomass yield (6 Mg ha −1 yr −1 ) was at the lower end of the range (4 to 15 Mg ha −1 year −1 ) reported for SRC willow production in the temperate zone (Searle & Malins, 2014), most likely reflecting the combination of damage by willow beetle (Bjorkman et al., 2000) and the age class of the plantation.
The broadly neutral (annualized) NECB reported in this study (−21 g CO 2 -C m −2 year −1 ) supports the longer term perspectives provided by studies of soil C, which have shown limited (Rowe et al., 2016) or no change (Walter et al., 2015) in soil organic C following land use change to SRC willow systems. The evidence from this research therefore suggests the primary C benefit of bioenergy production at this location is one of displaced CO 2 emissions from fossil energy sources, rather than net removal of CO 2 from the atmosphere, at least for the time period and conditions encountered during this study. The current paucity of observational studies at SRC willow plantations currently limits direct comparisons with sites of similar age class in different climatic regions, and with different soil types, land use history and management regimes. However, the range of annual NEE reported in this study was less negative than that reported for a younger willow plantation in Sweden (−331 to 818 g CO 2 -C m −2 year −1 ) used for bioenergy production and waste water filtration (Grelle et al., 2007). The more negative annual NEE at the Swedish site is likely explained by the combination of fertilization by nutrients in waste water and irrigation, as well as broad geographical and climatic influences.
The results of this study contrast with the findings of the only other previous EC study of SRC willow in the UK. Harris et al. (2017) reported (annualized) NEE and NECB values of −620 ± 18 g CO 2 -C m −2 year −1 and −221 ± 66 g CO 2 -C m −2 year −1 , respectively, based on observations obtained during two years of a four-year rotation at a young plantation established on former grassland in southern England. Differences may be expected between the sites due to past land use history (Rowe et al., 2016), as well as more favorable growing conditions at their more southerly site and/or improved varietal productivity at the younger plantation. Additionally, the willow beetle mediated defoliation in this study may have also played a role, with peak season daily NEE values in 2016 (before any notable defoliation in that year) being of broadly similar magnitude to peak growth phase values observed by Harris et al. (2017). The approximately neutral NECB at the mature SRC willow plantation in this study compared to younger sites underlines Net ecosystem exchange (NEE) is net ecosystem carbon dioxide exchange (g CO2-C m −2 ); GEP is gross ecosystem production (g CO 2 -C m −2 year −1 ); TER is total ecosystem respiration (g CO 2 -C m −2 year −1 ). The values shown for 2017 represent the period 2017-01-01 to 2017-11-26.

T A B L E 3 Estimates of time integrated
carbon dioxide fluxes obtained between 2014 and 2017 at a mature short rotation willow plantation the importance of quantifying all life stages and management practices when assessing and projecting the C and GHG balance of perennial energy crops, and creates a strong case for extended observational studies designed to track C and GHG dynamics throughout the entire economic lifetime of bioenergy production systems. The observed defoliation by willow beetle in this study provided an unexpected and novel opportunity to explore such a biotic disturbance event, adding to the limited number of EC studies to have captured CO 2 flux dynamics in treed ecosystems during defoliating insect outbreaks (e.g. Clark, Skowronski, Gallagher, Renninger, & Schäfer, 2014;Clark, Skowronski, & Hom, 2010;Olsson, Heliasz, Jin, & Eklundh, 2017;Wilkinson, Eaton, Broadmeadow, & Morison, 2012). Although it is speculated that measures to control willow beetle and maximise biomass yield (e.g. Bjorkman & Eklund, 2006;Dalin, Kindvall, & Björkman, 2009) would likely translate into increased ecosystem C capture, the lack of comparative data from a non-affected production cycle and/or reference site limits a full quantitative analysis at this time. Despite this data gap, it was clear that photosynthesis, ecosystem respiration and net CO 2 gains were lower during the most severe defoliation event in 2015 compared with corresponding months of other years of the production cycle. Whereas lower photosynthesis and net C uptake were obviously related to reductions in green leaf area during defoliation, lower rates of ecosystem respiration are most likely explained by reductions in foliar respiration rates, as observed at other treed ecosystems experiencing disturbance by defoliating insects (Clark et al., 2010). It is possible that changes to the timing and quality of litter and frass inputs to the soil also served to enhance soil heterotrophic respiration rates (Clark et al., 2010;Grüning, Simon, Rennenberg, & L-M-Arnold, 2017;L-M-Arnold et al., 2016), although any such increase was clearly outweighed by reductions to foliar respiration. Future research should aim to quantify how biotic disturbance events alter ecosystem C fluxes, the partitioning between auto-and heterotrophic components of total ecosystem respiration, and subsequent influences on ecosystem C storage and biomass yield.
This study represents the first EC observations of net ecosystem CO 2 exchange to be conducted at a mature SRC willow plantation in the UK. Whilst the beetle outbreak clearly influenced the C balance of the plantation during this study, the ecosystem still functioned as a net in situ sink for atmospheric CO 2 in all years within the range −147 ± 70 to −502 ± 84 g CO 2 -C m −2 year −1 . However, the NECB of −21 g CO 2 -C m −2 year −1 shows that the SRC willow plantation was close to CO 2 -C neutral after harvested biomass was accounted for. Whilst this important finding may appear less favorable than a net increase in ecosystem C storage when considered within the context of decarbonization and climate change mitigation, the results support the assumption that these crops are approximately C neutral as is often adopted in bioenergy life cycle assessments (Rowe et al., 2011). This finding also supports claims that bioenergy with C capture and storage (BECCS) may prove viable in delivering negative GHG emissions required to meet domestic and international emissions reduction targets as and when such technology becomes viable (Energy Technologies Institute, 2015;Rosen, 2018). Understanding and modeling the C cycle implications of bioenergy cultivation over larger geographical areas clearly requires additional observational data from a broader range of locations and over longer time frames. Further work is also needed to address the impacts of insect defoliation on yield and the cycling of C and other nutrients. Despite these knowledge gaps, this study represents an important step forward in understanding the C balance of SRC willow production and provides a valuable data resource for future work on this topic.