Mapping the reduction in gross primary productivity in subarctic birch forests due to insect outbreaks

It is projected that forest disturbances, such as insect outbreaks, will have an increasingly negative impact on forests with a warmer climate. These disturbance events can have a substantial impact on forests’ ability to absorb atmospheric CO2, and may even turn forests from carbon sinks into carbon sources; hence, it is important to develop methods both to monitor forest disturbances and to quantify the impact of these disturbance events on the carbon balance. In this study we present a method to monitor insect-induced defoliation in a subarctic birch forest in northern Sweden, and to quantify the impact of these outbreaks on gross primary productivity (GPP). Since frequent cloud cover in the study area requires data with high temporal resolution and limits the use of finer spatial resolution sensors such as Landsat, defoliation was mapped with remote sensing data from the MODIS sensor with 250 m× 250 m spatial resolution. The impact on GPP was estimated with a light use efficiency (LUE) model that was calibrated with GPP data obtained from eddy covariance (EC) measurements from 5 years with undisturbed birch forest and 1 year with insect-induced defoliation. Two methods were applied to estimate the impact on GPP: (1) applying a GPP reduction factor derived from EC measured GPP to estimate GPP loss, and (2) running a LUE model for both undisturbed and defoliated forest and deriving the differences in modelled GPP. In the study area of 100 km2 the results suggested a substantial setback to the carbon uptake: an average decrease in regional GPP over the three outbreak years (2004, 2012, and 2013) was estimated to 15± 5 Gg C yr−1, compared to the mean regional GPP of 40± 12 Gg C yr−1 for the 5 years without defoliation, i.e. 38 %. In the most severe outbreak year (2012), 76 % of the birch forests were defoliated, and annual regional GPP was merely 50 % of GPP for years without disturbances. The study has generated valuable data on GPP reduction, and demonstrates a potential for mapping insect disturbance impact over extended areas.


Introduction
It is estimated that forests account for half of the global terrestrial net primary productivity and act as important sinks of atmospheric CO 2 (Bonan, 2008).Forests in the Northern Hemisphere contribute significantly to this sink, with the mid-and high-latitude ecosystems as major contributors (Goodale et al., 2002;Kurz et al., 2008b).The highlatitude forests are predicted to be among the ecosystems that are most strongly influenced by climate change (Kurz et al., 2008b); a warmer climate is likely to increase forest productivity (e.g.Nemani et al., 2003;Boisvenue and Running, 2006), and result in higher uptake of CO 2 from the atmosphere.On the other hand, it is projected that the impact of forest disturbances will increase with a warmer climate (Seidl et al., 2014), and there are indications that disturbances such as wind, fires, and insect outbreaks have led to saturation of the carbon sink in European forests (Nabuurs et al., 2013).One important forest disturbance agent is insects; it is projected that the temporal and spatial dynamics, as well as the intensities and ranges of insect herbivore outbreaks, will be influenced by global warming (Vanhanen et al., 2007;Battisti 2008;Jepsen et al., 2008;Netherer and Schopf, 2010).These insect outbreaks can severely disturb forest ecosystems, and have a strong impact on carbon dynamics (Kurz et al., 2008a;Jepsen et al., 2009;Heliasz et al., 2011).Quantitative effects of insect outbreaks on the carbon balance are, however, not well known (Clark et al., 2010;Schäfer et al., 2010;Hicke et al., 2012), and insect outbreaks are generally excluded in large scale carbon modelling, which may result in overestimation of forests' ability to act as carbon sinks (Kurz et al., 2008b;Hicke et al., 2012).Consequently, it is important to develop methods both to monitor the spatial extent of insect outbreaks and to quantify the impact of these outbreaks on the carbon balance.
One alternative to estimate the impact on forest productivity is modelling: the impact of a large-scale outbreak of the mountain pine beetle (Dendroctonus ponderosae Hopkins) in British Columbia, Canada, was studied with a forest ecosystem model by Kurz et al. (2008a).The impact on the carbon balance of gypsy moth (Lymantria dispar L.) defoliation in New Jersey, USA, was also modelled with both a canopy assimilation model (Schäfer et al., 2010) and a terrestrial biosphere model (Medvigy et al., 2012).The impact of spruce budworm (Choristoneura fumiferana Clem.) outbreaks in eastern Canada were modelled by Dymond et al. (2010), and Landry et al. (2016) developed a Marauding Insect Module (MIM) in the Integrated Biosphere Simulator (IBIS) that enables simulation of insect outbreak for three insect functional types.Another alternative to quantify the influence of an insect outbreak on the carbon balance is to apply eddy covariance (EC) measurements: Brown et al. (2010Brown et al. ( , 2012) ) studied how a mountain pine beetle outbreak influenced net ecosystem productivity (NEP) in British Columbia, Canada; Clark et al. (2010Clark et al. ( , 2014) ) studied differences in net ecosystem exchange (NEE) between undisturbed years and years with severe defoliation by the gypsy moth in New Jersey, USA; and Heliasz et al. (2011) estimated the reduction in NEE during the growing season due to outbreaks of autumnal moth (Epirrita autumnata Borkhausen) and winter moth (Operophtera brumata L.) in northern Sweden in 2004.Even though not explicitly studied, there was gypsy moth defoliation of holm oak (Quercus ilex L.) present in a time series of EC measurements in southern France (Allard et al., 2008).These methods generate valuable data on the impact of insect defoliation on the carbon balance; however, to quantify the total regional impact, data on the extent of defoliation events are required.
To generate wall-to-wall estimates of the disturbance effect on the carbon balance, remotely sensed data from satellites can be used.Several studies have demonstrated that satellite based remote sensing techniques can be applied to detect insect disturbances with high accuracy; see, for example, Wulder et al. (2006), Adelabu et al. (2012), andRullan-Silva et al. (2013) for reviews.In this paper we study outbreaks of autumnal moth and winter moth in subarctic mountain birch (Betula pubescens ssp.czerepanovii N.I.Orlova) forests in northern Sweden.These outbreaks often cover large areas, but are often followed by within-season recovery of the foliage in parts of the outbreak areas, which in combination with cloudy conditions can limit the possibility to map the outbreaks with remote sensing methods.Nevertheless, outbreaks of autumnal and winter moth have been mapped in northern Fennoscandia with high accuracy with Landsat data (Tømmervik et al., 2001;Babst et al., 2010).The low temporal resolution of Landsat (16 days revisit time) can, however, be a limitation; as an example, only fractions of the area included in this study were visible in Landsat data during the peak of a severe outbreak in 2013.An alternative to Landsat data is coarse spatial resolution data from, for example, the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor, which provides data with high (daily) temporal resolution and a spatial resolution of 250 m × 250 m or coarser.MODIS-derived normalized difference vegetation index (NDVI) has been used to map autumnal and winter moth outbreaks with high accuracy in northern Fennoscandia (Jepsen et al., 2009), and Olsson et al. (2016b) developed a method for near-real-time monitoring of insect-induced defoliation that facilitates monitoring of refoliation later in the growing season.
Furthermore, there is a large body of research demonstrating that vegetation primary productivity can be estimated with remotely sensed data and a light use efficiency (LUE) approach (e.g.Prince, 1991;Ruimy et al., 1994;Running et al., 2004;Xiao et al., 2004;Wu et al., 2010;McCallum et al., 2013;Gamon 2015).The LUE concept was introduced by Monteith (1972) and Monteith and Moss (1977), suggesting that the primary productivity of plants has a strong linear relationship to the absorbed amount of photosynthetically active radiation (APAR), i.e. solar radiation in the spectral range 400-700 nm that is absorbed by the plant canopy.Since near-linear relationships between satellite-derived vegetation indices and the fraction absorbed PAR (fAPAR) have been established (e.g.Asrar et al., 1984;Sellers, 1987;Goward and Huemmrich, 1992;Myneni and Williams, 1994;Olofsson and Eklundh, 2007), it is possible to create a LUE model driven by remote sensing data.Such a LUE model could be applied for large-area estimates of the impact of forest disturbance on the uptake component of the carbon balance.Bright et al. (2013) utilized Landsat data to map bark beetle damage in northern Colorado, USA, and MODIS GPP data, which are based on a LUE model, to quantify the impact of the damage on GPP.However, to the knowledge of the authors, no previous study has utilized remote sensing data and developed a LUE model to monitor and quantify the impact of defoliating insects' outbreak on GPP.
In this study we utilized EC measured GPP to develop a LUE model, driven by MODIS-derived NDVI, to quantify the regional impact on GPP of insect-induced defoliation, and to map the spatial extent of the defoliation.Our main study objective was to compare GPP for years with insect damage (2004, 2012 and 2013) with GPP for years without insect damage (2007, 2009, 2010, 2011 and 2014) in the birch forest of a subarctic valley of northern Sweden.The analysis was achieved with two methods: (1) finding GPP for undisturbed forest and estimate the impact of an insect outbreak with a common reduction factor derived from EC Biogeosciences, 14, 1703-1719, 2017 data, and (2) by applying a LUE model for both undisturbed and defoliated pixels and computing the differences.

Study area
The study area was the mountain birch (Betula pubescens ssp.Czerepanovii N.I.Orlova) forests in a valley south-west of the village of Abisko (68.35 • N, 18.82 • E), and along Lake Torneträsk, as illustrated in Fig. 1 (green).The area is located in the subarctic zone in northern Sweden with Lake Torneträsk at an altitude of 345 m a.s.l. and with the highest mountains reaching 1700 m a.s.l.(Interact, 2016).These birch forests are infested by the autumnal moth (Epirrita autumnata Borkhausen) and the winter moth (Operophtera brumata L.) in time intervals of 9-10 years (Bylund, 1995;Tenow et al., 2007).The first reported outbreaks by the autumnal moth in northern Fennoscandia are from the mid-1800s, and the winter moth has been observed in the northern parts of Fennoscandia since late 1800 (Tenow, 1972).These insect outbreaks strongly influence the birch forests (Ammunét et al., 2015): severe defoliation events may result in stem mortality, requiring decades of recovery (e.g.Tenow, 1996;Tenow and Bylund, 2000;Jepsen et al., 2013), and understorey vegetation can shift into more grass-dominated communities (Karlsen et al., 2013;Jepsen et al., 2013).Root-associated fungal communities can change (Saravesi et al., 2015), as can chemical and physical properties of the soil (Kaukonen et al., 2013).A warmer climate, especially a lower frequency of years with extremely cold winters, as reported by Callaghan et al. (2010), strongly influences birch moth populations (Babst et al., 2010).The autumnal moth outbreaks have expanded into colder, more continental regions, and the winter moth has reached further to the north-east into areas where the autumnal moth previously dominated (Jepsen et al., 2008).The latest outbreaks in the study area occurred in 2004, with a documented reduction in carbon sink strength of 89 % at an EC tower located in birch forest (Heliasz et al., 2011), and in 2012 and 2013 (Bengt Landström, county administrative board of Norrbotten, personal communication, 31 October 2013).These outbreak events were included in this study.

Remote sensing data and smoothing of time series
We used two Terra/MODIS satellite data products with 8day temporal resolution: (1) MOD09Q1 version 5, surface reflectance in the red and near-infrared (NIR) bands, including quality assurance (QA) information, with 250 m × 250 m spatial resolution, used mainly to derive NDVI (LPDAAC 2016a); and (2) MOD09A1 version 5, surface reflectance, as well as QA data, with 500 m × 500 m spatial resolution (LP-DAAC, 2016b), utilized due to the product's more comprehensive QA data.
NDVI was computed from the MODIS data as (Rouse et al., 1973;Tucker, 1979) where red is reflectance in the red wavelength band, and NIR is reflectance in the near-infrared wavelength band.We created time series for the period 2000-2014 for all pixels in the study area and processed in TIMESAT ver.3.2.TIMESAT is a software package used to reduce the influence of noise by fitting smoothed functions to time series of data (Jönsson andEklundh, 2002, 2004).In this study we applied the same fittings and weights as in Olsson et al. (2016b): Double logistic functions were used to smooth the raw NDVI data and QA data from both MOD09Q1 and the more comprehensive QA flags in MOD09A1 were utilized to estimate the quality of the NDVI observations.In this study we use the term NDVI DL to refer to the smoothed time series of NDVI.

Fraction of canopy-absorbed PAR and relationships with NDVI
The fraction of PAR absorbed by the canopy (fAPAR canopy ) was measured at a spectral tower located in birch forest north-west from Abisko (Fig. 1, black star).fAPAR canopy was obtained using the four-component method, i.e. measurements of incoming PAR above canopy, the total reflected PAR above the canopy, the transmitted PAR below the canopy, and the reflected PAR by the understorey vegetation and ground below the canopy.See Eklundh et al. (2011) for detailed information about the estimation of fAPAR canopy .All PAR sensors were calibrated at the field site following the procedure by Jin and Eklundh (2015), and fAPAR canopy at solar noon time was calculated and used in the final analysis.fAPAR canopy data were available for the years 2010 and 2011.
Average fAPAR canopy over 8-day periods, coinciding with the MODIS 8-day periods, were computed, and an ordinary least squares (OLS) regression was performed to find the relationship between fAPAR canopy and NDVI DL Myneni and Williams (1994).The linear equation derived was used in the LUE model to obtain fAPAR from the double logistic fitted NDVI.

Eddy covariance and meteorological data
The EC tower is situated in the eastern part of the study area (Fig. 1, black triangle) and located near the crossing point of four nominal MODIS pixels with 250 m × 250 m spatial resolution (Fig. 2).Vegetation in the four pixels is similar, with some open mires in the north-east pixel and a paved road crossing the two southernmost pixels.The tower's footprint is estimated to be about 200 m long, which is slightly smaller than a MODIS pixel.The prevailing wind directions are from the west and from the east; hence, the main footprint of the EC tower is to the west and east from the tower, where vegetation is most homogeneous.Time series of NDVI were extracted and mean values and standard deviations were computed for the four MODIS pixels to study whether there were any larger deviations in the pixels' NDVI signals.In Fig. 3, mean NDVI and standard deviation for the four pixels in the period 2010-2014 are shown.The low standard deviations indicate that there are minor differences in the NDVI signal between the pixels during the main growing season for both raw NDVI and NDVI DL both for years without disturbance and for outbreak years.Hence, we assume that a varying footprint of the EC tower due to varying wind directions and stability will have a limited influence on the EC measurements.
The EC measurements were made 8 m above ground, 3.3 m above canopy, using a three-dimensional sonic anemometer (Metek USA-1; METEK GmbH, Germany) and an open-path infrared gas analyzer (Licor 7500; LI-COR Inc., USA).The system was operated with a frequency of 20 Hz, and data were recorded by a data logger (CR1000; Campbell Scientific, Inc., USA).Additional measurements of air temperature (Vaisala WXT510; Vaisala, Finland) and incoming photosynthetic flux density (PPFD; JYP 1000, SDEC, France), used for flux partitioning and gap filling, were made at the tower.Data were obtained each year during the period 1 May to 30 September, which is from before the start of the growing season (Karlsson et al., 2003) until the late growing season; during the years included in this study GPP was approaching zero by the last week of September.For the years 2004 and 2013, temperature and PAR were obtained from Abisko Scientific Research Station (ANS, 2015); comparisons between data from ANS and the EC tower showed small differences for the years when data were available from both sources.
EC flux calculations were done with the EddyPro software ver.5.2.1 (LI-COR Inc., USA).Gaps caused by bad weather conditions, bad EC measuring conditions, or short breaks in instrument functioning were filled with the "Eddy covariance gap-filling & flux-partitioning tool" (http://www.bgc-jena.mpg.de/~MDIwork/eddyproc/).The main reasons for removing data were precipitation, as we used an openpath gas analyser, and atmospheric conditions that did not fulfil the turbulence conditions required.We considered the gap-filling approach suitable also for defoliated years since the gap-filling function is created based on data from short time windows, usually 7 days, and hence adjusts the fitting parameters for changing ecosystem conditions.A model from the same website was used to partition NEE into GPP and ecosystem respiration (R eco ).It was assumed that nighttime NEE is equal to night-time R eco .Accordingly, the accepted night-time data were fitted to the Lloyd and Taylor (1994) model based on air temperature.This model was also used to estimate R eco during daytime conditions.GPP was estimated as the residual after subtracting R eco from the measured NEE.Details about gap filling and flux partitioning are described in Reichstein et al. (2005).

Land cover and elevation data
Land cover data were obtained from the Swedish mapping, cadastral, and land registration authority (Lantmäteriet; Dnr: I2014/00579).These land cover data are based on a classification of Landsat TM data, were updated in the year 2000 as a part of the CORINE land cover project, and have a spatial resolution of 25 m × 25 m (Lantmäteriet, 2010).Birch forests in the study area were identified by extracting all pixels with broadleaved forest.Since birch is the dominating tree species with only a few sporadic individuals of other species (Sonesson and Lundberg, 1974), all forests were considered to be birch.These data were used to calculate the fraction forest cover per MODIS pixel.
Elevation data were obtained from Lantmäteriet (National survey of Sweden) as a digital elevation model (DEM) with 50 m × 50 m spatial resolution (Lantmäteriet; Dnr: I2014/00579).Mean elevation for each MODIS pixel was computed as the average altitude of all DEM pixels covered by a MODIS pixel.To adjust for altitudinal differences in temperatures across the study area, a mean summer temperature gradient of 0.5 • C per 100 m (Josefsson, 1990;Holmgren and Tjus, 1996) was applied to the temperature data from the EC tower.

Light use efficiency model
A LUE model with mean values of daily GPP in 8-day intervals (GPP lue ) (g C m −2 day −1 ), corresponding to the time interval of the MODIS data, was developed as where ε (g C MJ −1 ) is the light use efficiency, fAPAR 8day is fAPAR for a MODIS 8-day period derived from NDVI DL , and PAR 8day (MJ m −2 day −1 ) is mean daily PAR measured at the EC tower over the 8-day period.The light use efficiency varies between vegetation types, and variability in meteorological conditions is accounted for through reductions factors for temperature and vapour pressure deficit (e.g.Field et al., 1995;Prince and Goward, 1995;Potter et al., 1999;Turner et al., 2003).In this study the light use efficiency was computed as where ε max (g C MJ −1 ) (see Sect. 2.3.3) is the maximum efficiency applied in the model and f 8day is a reduction factor.We assumed that accounting for temperature only is sufficient in our study region, which is supported by Bergh et al. (1998) and Lagergren et al. (2005).Two models were created to describe f 8day , as in Lagergren et al. (2005) where T mean8 ( • C) is the mean temperature for a MODIS 8day period.The reduction factor was computed as where GDD thres (see Sect. 2.3.3) is a threshold applied to decide when temperature and frost events no longer influence ε, in a similar fashion to Bergh et al. (1998) and Lagergren et al. (2005).S GGD is a reduction factor influenced by GDD and frost events and computed as where P frost is a reduction factor controlled by frost events and computed as where T min8 ( • C) is the lowest temperature during a MODIS 8-day period.

Second part of the growing season
In the second part of the growing season, covering late June to September, f 8day is controlled by mean temperature only as where T thres ( • C) (see Sect. 2.3.3) is a temperature factor for controlling the influence of the 8-day mean temperature during the second part of the growing season.

LUE model optimization
The LUE model was optimized to find three factors: (1) the GDD threshold (GDD thres ), (2) the temperature factor (T thres ), and (3) the period to change from the first to the second seasonal model.These were found by minimizing the root mean square error (RMSE) and maximizing R 2 , based on GPP lue and daily mean values of GPP from the EC tower over MODIS 8-day periods (GPP EC ).To compute ε max , the mean value of the light use efficiency for all MODIS periods with maximum efficiency, i.e. f 8day = 1, was calculated, where the efficiency was computed as where GPP EC was derived from the EC tower.Two ε max values were computed: one including data from the 5 years (2007, 2009, 2010, 2011 and 2014) with undisturbed birch forest, and one (ε max, def ) for the year 2012 with insect defoliation.No data were available from 2008 due to equipment failure, and in 2013 the measurements were disturbed by larvae climbing the equipment.

LUE model uncertainty
A Monte Carlo approach was applied to evaluate the uncertainty of the LUE model by creating sets with 100 parameter values each for ε max and slope and intercept derived from the OLS regression between fAPAR canopy and NDVI DL .The standard deviation of ε max was estimated from all MODIS periods with maximum efficiency, as described in Sect.2.3.3, and a 95 % confidence interval for the regression line was estimated.The different sets of parameters were created randomly from a uniform distribution, and the Monte Carlo simulation was run for all possible combinations of parameter values for the 5 years with undisturbed forests and over 15 sets of 100 MODIS pixels with birch forest.Mean and standard deviation of LUE modelled GPP were estimated from these simulations.

Identifying MODIS pixels with defoliated birch forest
Defoliated MODIS pixels were identified for the 3 years with insect outbreaks with a near-real-time monitoring method based on Kalman filtering and cumulative sums (Olsson et al., 2016b).The method identifies a seasonal trajectory of NDVI representing birch forest during a year without disturbances, called stable season.A Kalman filter (Kalman, 1960) is applied to the raw NDVI observations from the year of study and deviations from the stable season are computed.A cumulative sum (CUSUM) filter (Page, 1954) is applied to these deviations, and a pixel is classified as defoliated when the cumulative sum of deviations reaches a given threshold.
In a near-real-time application the stable season can only be derived from years prior to the year of study.In this study we modified the method so that the stable season was derived from all years with available data.For high detection accuracy, the method requires that a MODIS pixel is covered by at least 50 % forest.Hence, based on the land cover data from Lantmäteriet, forest in pixels with lower forest cover was excluded, resulting in 100 km 2 of the in total 125 km 2 birch forest in the study area being included; the mean forest cover was 80 % per MODIS pixel.The method detected 74 % of the defoliated sampling areas in the study area with a misclassification of undisturbed areas of 39 % (Olsson et al., 2016b).

Annual GPP loss due to insect defoliation
GPP for years without insect defoliation was estimated for all pixels by applying the LUE model and computing the mean value for the 5 years without insect outbreak, and with data available from the EC tower.The 8-day average of incoming PAR (PAR 8day ), measured at the EC tower, was assumed to be valid for all pixels in the study area, which was also suggested by comparisons between PAR measured at the EC tower and ANS.
Two methods were applied to study the reduction in annual GPP due to the insect outbreaks: (1) a method based on a reduction factor derived from the EC data from 2012, when the birch forest in the footprint of the tower was severely defoliated and no refoliation occurred.This reduction factor was applied to all pixels in the study area and (2) a method where the LUE model was applied to all defoliated pixels with ε max, def computed for defoliated growing seasons, and where the loss in GPP was computed as the difference between undisturbed and defoliated years.

Method 1 -GPP reduction factor
The fraction of the measured annual GPP at the EC tower that was lost due to the insect outbreak in 2012 was computed as where GPP redfact is the reduction factor and GPP defoliated is annual GPP from the EC tower in 2012.GPP undisturbed is GPP from the tower representing a year without disturbances and computed as the mean of annual GPP for the 5 years without disturbances.
The reduction in annual GPP was computed for each pixel by applying the reduction factor to GPP for undisturbed years and multiplying with the area forest cover in the pixel.The same reduction factor was applied to all years with insect defoliation.The total impact of the defoliation was computed as the sum of GPP loss for all defoliated pixels in the study area, and for each year with insect outbreak.

Method 2 -LUE model for defoliated pixels
The LUE model, modified to model growing season with defoliation, was applied to all defoliated pixels in the study area to estimate annual GPP for each year with defoliation.Derivation of ε max, def was done with the same method as ε max , but only data from 1 year with insect outbreak (2012) were available to estimate ε max, def and to evaluate the performance of the defoliation LUE model.For each year with insect outbreak, the regional reduction in GPP was computed by summing, over all pixels identified as defoliated, the difference between GPP for years without outbreak and GPP for this specific outbreak year.

Influence of refoliation
We also studied how recovering foliage later in the growing season influenced the two methods.The assumption was that recovering foliage would result in slightly higher NDVI DL values, which would enable Method 2 to capture the refoliation and, hence, estimate GPP losses more accurately.All pixels that were detected as defoliated were classified as refoliated or non-refoliated with the defoliation monitoring method.The differences between GPP loss derived with methods 1 and 2 were computed as GPP loss method 1 minus GPP loss method 2. Finally, the mean differences for refoliated and non-refoliated pixels were derived.

Correlation between fAPAR and NDVI
There was a strong linear relationship between 8-day mean values of fAPAR canopy and NDVI DL for NDVI DL values ≥ 0.4 (Fig. 4).The influence of observations with NDVI DL values < 0.4 and with f 8day > 0 was small.For the www.biogeosciences.net/14/1703/2017/Biogeosciences, 14, 1703-1719, 2017 The 95 % confidence intervals for slope and intercept applied in the Monte Carlo simulation to estimate the LUE model's uncertainty were −0.05 ± 0.18 (intercept) and 0.60 ± 0.11 (slope).

Light use efficiency
Optimization resulted in a GDD thres of 32 growing degree days (Fig. 5, left) and a T thres of 8 • C (Fig. 5, middle).The optimal period to change the model for f 8day was after MODIS period 23, i.e. the last week of June (Fig. 5, right).
Light use efficiency for years with no disturbance and with f 8day = 1 (black line with error bars in Fig. 6) gave an ε max The correlation between GPP EC and GPP lue was strong with R 2 = 0.90 (Fig. 7).The intercept is −0.11 and the slope is 1.01, indicating that the LUE model performs well for years without outbreaks.The low GPP observations with several zero values for LUE modelled GPP are from May, before budburst for the birch forest.These low GPP values have little influence on annual GPP.The Monte Carlo simulation resulted in an estimated standard deviation of 30 % of the mean annual GPP.Hence, all annual GPP values derived from the LUE model are given with a standard deviation of 30 % of annual GPP.GPP measured from the EC tower and the 5 years with available data (Table 1) resulted in a mean annual GPP of 440 g C m −2 yr −1 .During the outbreak in 2012 annual GPP was 180 g C m −2 yr −1 , which resulted in a reduction in GPP compared to undisturbed conditions of 59 %.Hence, a reduction factor of 0.59 was applied to quantify the impact of the insect outbreak on GPP.

Method 2 -LUE model for defoliated pixels
The correlation between GPP EC and GPP lue for the year with defoliation ( 2012) and data available from the EC tower was weaker than for years without disturbances, with an R 2 of 0.83 (Fig. 8).The figure, with an intercept of −0.54 and a slope of 1.25, indicates that the LUE model underestimates GPP for lower values.The light use efficiency for the MODIS 8-day periods with f 8day = 1 (black circles in Fig. 6) gave an ε max, def of 0.98 ± 0.25 g C MJ −1 (±1 standard deviation), resulting in the following LUE model for defoliated pixels: GPP lue, defoliated = 0.98 × f 8day • (−0.05 + 0.60 × NDVI DL ) In Fig. 6, NDVI DL has higher values in the year with defoliation compared to undisturbed years in May (periods 16-18).These high NDVI DL values are due to poor fitting of the double logistic function during winter and early spring in 2012 (see Fig. 3, where NDVI DL increases earlier in 2012 compared to the other years).The impact on the result is, however, small, since these eight periods are in the early part of the growing season when f 8day is zero.The Monte Carlo simulation resulted in an estimated standard deviation of 35 % of the mean annual GPP for years with defoliation.Hence, all annual GPP losses estimated with Method 2 are given with a standard deviation of 35 %.

Defoliated areas and quantifying the insect outbreaks impact on annual GPP
In the year 2012, with the most widespread defoliation in this study, 76 % of the 100 km 2 forests were defoliated (Table 2 and Fig. 9).In 2004 and 2013, 53 and 55 % of the forests were defoliated, respectively.The mean annual reduction in regional GPP due to the insect outbreaks for the three outbreaks studied was 15 ± 5 Gg C yr −1 according to Method 1, with the largest outbreak in 2012 with a negative impact on regional GPP of 18 ± 5 Gg C yr −1 (Table 2).The average annual regional GPP in the study area, derived with the LUE model (Eq.12) and the 5 years without insect outbreak, was 41 ± 12 Gg C yr −1 , which gives a reduction in 2012 of 44 %.
The impacts of the outbreaks in 2004 and 2013 were reduction in regional GPP of 32 and 34 %, respectively.There were no differences in the GPP reduction per square metre between the outbreak years.When a LUE model was also applied to model GPP during defoliation events (Method 2) the mean annual decrease in regional GPP was 15 ± 5 Gg C yr −1 , which is the same estimate as with Method 1.The regional GPP loss in 2012 was 20 ± 7 Gg C yr −1 , which is slightly higher compared to Method 1.In the year 2004 the two methods resulted in similar decreases in GPP, while the GPP decrease was larger with Method 1 in 2013.Differences in GPP loss per square metre between the years were larger with Method 2: 190 ± 67 g C m −2 yr −1 in 2013 was the lowest GPP loss, and 270 ± 95 g C m −2 yr −1 in 2012 was the largest GPP loss.
Table 3. Differences in GPP loss (g C m −2 yr −1 ) between methods 1 and 2 for MODIS pixels with recovering foliage later in the season, and pixels with no refoliation according to the defoliation monitoring method.Higher GPP loss with Method 2 gives negative values.growing season (Table 3).For all years the mean differences in GPP loss (g C m −2 yr −1 ) between the methods were lower for pixels that recovered later in the growing season.These results suggest that Method 2 captured some of the refoli-ation, though the differences are small and within the error margin.

Discussion
This study has shown a substantial setback in GPP caused by insect defoliation in a subarctic deciduous forest in northern Fennoscandia.At the EC tower, GPP decreased by 260 g C m −2 yr −1 (59 %) during the outbreak in 2012 compared to the mean GPP of undisturbed years.In the entire study area annual mean values of decrease in GPP ranged from 190 ± 67 to 270 ± 95 ± g C m −2 yr −1 .The total decrease in regional GPP due to the three insect defoliation events studied here was estimated to be 45 ± 14 Gg C, which is of the same magnitude as the average annual regional GPP of 41 ± 12 Gg C yr −1 for single years with no disturbances.During the most severe outbreak year (2012), the annual regional GPP loss was nearly 50 % (20 Gg C yr −1 ), with 76 % of the 100 km 2 birch forests in the study area defoliated.In this study we have estimated the impact on GPP only but we noted that during the outbreak in 2012 the decrease in R eco was larger than the decrease in GPP during the growing season around the EC tower.Respiration is affected by insect outbreaks in two ways: (1) autotrophic respiration is reduced as defoliated trees cannot photosynthesize and (2) heterotrophic respiration increases when dead larvae decompose.The amount of carbon respired by larvae is likely to be the same as the amount of carbon in the eaten leaves, so we should only observe a shift of respiration in time.In addition, larvae transport nutrients from trees to fungi and bacteria living in soil, which further increase respiration.The increase in heterotrophic respiration did not offset decrease in autotrophic respiration, and R eco for the outbreak year decreased in comparison to non-disturbed years.This study also highlights the advantage of combining EC data and remote sensing data by applying data from an EC tower to calibrate a LUE model and applying satellite data to estimate the impact on GPP over larger areas.EC measurement alone cannot be extrapolated with high accuracy if the spatial and temporal extent of an outbreak is unknown, and the LUE model could not be developed without EC data.The combination facilitates wall-to-wall mapping of forest disturbances and quantitative estimates of the impacts on primary productivity.
There are, however, limitations in the study that must be considered.One major challenge is to establish baseline conditions for GPP in areas with reoccurring insect outbreaks, as in Abisko.As a comparison, Olsson et al. (2016a) tested a defoliation detection method on the outbreak in Abisko in 2013 and achieved the highest detection accuracies when the baseline conditions were based on the 6 years with highest NDVI values in the period 2000-2012.In this study the annual GPP for years without disturbances was estimated as the mean of the 5 years without insect outbreak and with available EC data.It is likely that some of these years were influenced by the insect outbreak in 2012.The 2 years prior to when the insect populations reached outbreak levels (2010 and 2011) had lower annual GPP than the years 2007 and 2009 (Table 1), and it is likely that GPP in 2014 was influenced by the insect defoliation in 2012 and 2013.Heliasz (2012) suggested that GPP reaches pre-outbreak levels 2-3 years after an outbreak, and Hoogesteger and Karlsson (1992) showed that leaf area index (LAI) returned to pre-defoliation levels 2 years after 100 % artificial defoliation even though tree ring width was lower than normal at least 3 years after the experiment.For the birch forests it may take decades to fully recover from severe outbreaks (Tenow and Bylund, 2000).
To get an indication of the potential influence on GPP by insect defoliation for the non-outbreak years we modelled GPP based on PAR for the years with data available from the EC tower and compared with EC-derived GPP (see Supplement).The result showed that measured GPP at the EC tower, and GPP modelled with PAR data, were similar in 2007 and 2009.In the 2 years prior to the outbreak (2010 and 2011), measured GPP was lower than PAR modelled GPP, indicating that there were signs of defoliation by growing larval population.Also in 2014 when the birch forests were recovering, measured GPP was lower than PAR modelled GPP.During the insect outbreak in 2012 measured annual GPP was 290 g C m −2 yr −1 lower than PAR modelled GPP, which is larger than the decrease of 260 g C m −2 yr −1 applied in this study.In addition, we ran the LUE model with meteorological data from ANS for the year 2008 to fill the gap in the time series with measured GPP and to study how well it agreed with the years 2007 and 2009.According to the LUE model, annual GPP at the EC tower was 440 g C m −2 yr −1 in 2008, which agrees with the GPP value for undisturbed years of 440 g C m −2 yr −1 applied in the study.However, since years that are influenced by pre-outbreak defoliation as well as a recovery year are included as undisturbed years, it is likely that the baseline GPP applied in this study is lower than GPP for undisturbed conditions.This is also indicated by the larger difference between PAR modelled and measured GPP in 2012 and suggests that the estimated decreases in GPP due to insect outbreaks in this study are on the lower side.
Another limitation is the assumption that no other factors than insect outbreaks influence annual GPP, even though it is likely that also meteorological conditions influence GPP.The comparison between EC-derived GPP and PAR modelled GPP suggests that only 2 years with EC data represent undisturbed forest; hence, the amount of data from the EC tower is too small to study correlations between EC-derived GPP and meteorological variables.Instead we studied correlations between NDVI and meteorological data from ANS, where we used the mean of the highest NDVI DL value of each year derived from 200 MODIS pixels with birch forest.To minimize the influence of insect-induced defoliation we excluded the outbreak years and years prior to and after outbreaks.No linear correlations between PAR and GPP were found.There were, however, negative correlations between temperature and seasonal maximums of NDVI DL , with the strongest correlation between NDVI and the mean temperature in May-June.The influence of temperature on NDVI was weak, and due to the estimated uncertainties of the LUE model of 30 % we did not include these correlations in the analysis.However, with data from the EC tower available for more years it would be a potential improvement of the method to include meteorological data when estimating the decrease in annual GPP.
There are also uncertainties in the LUE model.The relationship between fAPAR 8day and NDVI DL (Eq.11) was estimated from two growing seasons without disturbances.Due to larvae disrupting the PAR sensors there were no fA-PAR data available from the outbreak years; hence, Eq. ( 11) was used also for defoliation events.Furthermore, the relationship was derived from fAPAR obtained from the upper canopy, which may not be representative of the entire forest, since the relationship between fAPAR 8day and NDVI DL is likely to vary with understorey and forest densities in the study area.The relationship is also likely to vary with varying understorey responses due to defoliation, which may influence the estimated decreases in annual GPP.Accounting for these uncertainties would require more data on the fAPAR and NDVI relationship as well as more detailed land cover data, which would make the model more complex.Hence, we assume this limitation to be acceptable, and since the aim of the study was to estimate the influence of defoliation of the birch trees, we considered fAPAR canopy to be the most suitable variable.Another potential limitation is that the LUE model developed for years with defoliation seems to under-P.-O.Olsson et al.: Mapping the reduction in gross primary productivity estimate GPP for values lower than about 1.5 g C m −2 day −1 (Fig. 8).However, for the outbreak year with available EC data (2012) the underestimated values from the LUE model are mainly due to a cold spring that resulted in a large reduction factor (f 8day ).During the main growing season LUE modelled and EC-derived GPP agrees well, which increases confidence in the modelling.
It may also seem surprising that the difference in NDVI DL was comparably low in relation to the difference in light use efficiency.It is, however, known that NDVI saturates for high LAI and that small changes in NDVI can be associated with large changes in LAI (e.g.Myneni et al., 2002).The light use efficiency, on the other hand, can decrease substantially with lower LAI since more leaves will operate in the lightsaturated portion of the photosynthesis (e.g.Medlyn, 1998).There are also uncertainties in how well the EC tower footprint represents the entire study area.Heliasz (2012) utilized a permanent EC tower as reference and a mobile EC tower to study variability in carbon exchange in the birch forests around Abisko and concluded that there were only minor differences in GPP at seven sites during the peak growing season in 2008 and 2009.Hence, we consider the EC tower footprint to be representative for the study area.
The accuracy of the defoliation detection method also influences the results of the study.The method missed 26 % of the defoliated MODIS pixels and misclassified 39 % of the undisturbed pixels as defoliated in the evaluation data used by Olsson et al. (2016b).This implies that the defoliated areas in 2004 and 2013 were slightly overestimated, while the defoliated area in 2012 is likely underestimated, though the impact on the total numbers is likely small.It should also be considered that 20 % of the forests in the study area were excluded since they are located in MODIS pixels with < 50 % forests cover.Thus, the current estimate of a total reduction in GPP may be conservative.
A limitation with the developed LUE model for largearea estimates is that it includes observed meteorological data (temperature and PAR).An alternative for running the model over larger areas would be to use modelled meteorological data (Olofsson et al., 2007;Schubert et al., 2010).There are also uncertainties related to the temperature data utilized.The gradient applied to model mean temperatures depending on altitude is likely to give accurate estimates in the study area.However, minimum temperatures are more uncertain since cold air can drain downhill and accumulate in valleys and low areas, rather than decrease with altitude.Altogether, since the EC tower is located on a small ridge in the lower, flat parts of the study area, we anticipate that the temperatures there are not substantially lower than the area in general.We compared with lowest daily temperature from Abisko research station, which is located near the spectral tower 10 km to the west (Fig. 1), and at a slightly higher altitude than the EC tower.For all periods with frost events during the early season, i.e. when the lowest temperature influences f 8day , the mean value of absolute differences, with the coldest temperatures at the research station, was only 0.4 • C. With these small temperature differences and since frost events only influence GPP in the early growing season, the impact on annual GPP was considered minor.
The defoliation detection methods used in this study gives a time series of smoothed NDVI that captures the timing of the defoliation event as well as the potential refoliation.The LUE model, on the other hand, utilizes NDVI smoothed with double logistic functions.These functions do not capture the typical seasonal trajectory for years with refoliation.This is illustrated in Fig. 3, where raw NDVI stays around 0.6 during the entire growing season in 2012, when there was no refoliation around the EC tower.In 2013, when there was substantial refoliation around the EC tower, raw NDVI stays around 0.6 during June, but it increases to pre-outbreak levels in early July, when refoliation occurs.In 2004 the raw NDVI values has a pattern similar to that in 2013 with low values (around 0.6) until early August, when refoliation results in a later season peak in NDVI.This seasonal development of raw NDVI agrees well with GPP for the limited period with available EC data in the outbreak year 2004.NDVI DL does not capture this trajectory, with sharply increasing NDVI values that level off and start increasing again later in the season.However, even though the actual timing of the defoliation is not captured during years with refoliation the total growing season GPP is well modelled.A new version of TIME-SAT, currently developed and tested, will also capture more detailed seasonal trajectories with smooth fitting of curves.These new curve-fitting methods have a potential to improve the performance of the LUE model.
We applied two methods to quantify the impacts on GPP to study which methods performed better for refoliating birch forests.The assumption was that Method 2 would be more adaptive and adjust for differences in defoliation intensities between MODIS pixels.Since the level of defoliation, as well as understorey responses to the defoliation, is likely to influence NDVI DL , which in turn will influence fAPAR, it was anticipated that a method based on a LUE model to derive GPP during defoliation events would capture variability in defoliation levels and understorey responses between MODIS pixels.Method 1, on the other hand, with a common reduction factor, does not account for local differences between pixels and is similar to upscaling the local conditions at the EC tower, even though the method has the advantage that annual GPP for each pixel is derived with a LUE model, and hence should be more accurate than assuming that GPP for all MODIS pixels is identical to GPP at the EC tower.For the years 2004 and 2012, the two methods resulted in similar estimates of the GPP loss with slightly larger decrease in GPP for Method 2. In 2013, the difference between the methods was larger with the highest decrease in annual GPP for Method 1.One possible explanation for the smaller decrease in annual GPP according to Method 2 for the year 2013 is that the growing season seems to have been shorter and that refoliation started earlier and was stronger in 2013 compared to 2004; this is indicated by the seasonal developments of NDVI.It should also be noted that higher NDVI might be due to increasing growth of understorey grasses favoured by the changed light conditions due to defoliation (Karlsen et al., 2013) rather than recovering birch.
The impact of insect outbreaks on the carbon balance has been quantified in earlier studies: Heliasz et al. (2011) studied the impact on NEE of the autumnal moth and winter moth outbreak in Abisko in 2004, but these measurements started on 2 July, which was around 10 days after the larvae reached peak densities, which most likely resulted in an underestimated reduction in NEE.To facilitate a comparison between the outbreak years 2004 and 2012, we computed GPP for the period 2 July to 30 September for all years with EC data.This indicated that the two outbreak years had a similar impact on the carbon balance during the period studied with a GPP loss of 210 g C m −2 yr −1 in 2004 and 200 g C m −2 yr −1 in 2012 compared to years without disturbance.Furthermore, the loss of 200 g C m −2 yr −1 in the year 2012 and for the same time period as studied in the year 2004, compared to the GPP loss of 260 g C m −2 yr −1 for the entire growing season in 2012, suggests that the impact on NEE was underestimated by Heliasz et al. (2011).Clark et al. (2010) found the highest difference in NEE between undisturbed years and years with severe defoliation by the gypsy moth in New Jersey, USA, to be 266-480 g C m −2 yr −1 , and Clark et al. (2014) found that midday NEE during complete defoliation was 14 % of pre-defoliation rates.Allard et al. (2008) noted that cumulative NEE was lower during a year with insect defoliation compared to years without disturbances; however, the low NEE value might to a large extent have been caused by a dry spring.Brown et al. (2010) found that a mountain pine beetle outbreak turned a forest into a carbon source; no preoutbreak EC data were available to quantify the impact on NEP, but recovery after the outbreak was faster than anticipated (Brown et al., 2012).It should be noted that the mountain pine beetle feeds within the phloem and directly kills trees, while the moth species discussed above are defoliators that usually only kill trees in cases of severe and repeated outbreaks (Hicke et al., 2012).Modelling studies have also found that forests have changed from sinks into sources of carbon, in some cases for extended periods (Kurz et al., 2008a;Dymond et al., 2010;Schäfer et al., 2010;Medvigy et al., 2012).However, to our knowledge, this is the first study that has utilized remote sensing data and developed a LUE model calibrated with EC data to both quantify and map the spatial extent of the impact of defoliating insects' outbreaks on GPP.
The results of this study could help to reduce uncertainties in the impact of insect outbreaks on primary productivity as well as to improve carbon budgets by including insectinduced defoliation.For the mountain birch forests in this study the estimated reduction in annual GPP, compared to years without disturbances, was 50 % when there was limited refoliation in the study area.For years with widespread refoliation, the annual GPP losses were about one-third of GPP for years without disturbances.In addition, the spatial and temporal mapping of insect defoliation provided by remote sensing is important for accurate simulation of the carbon dynamics.Furthermore, the outbreak area included in this study is only a fraction of the 10 000 km 2 estimated to have been severely defoliated in northern Fennoscandia during the period 2000-2008 (Jepsen et al., 2009).Assuming that the conditions were similar over northern Fennoscandia, the insect defoliation over these vast areas would result in a potential total regional GPP loss for the time period of the magnitude 2-3 Tg C. Models not accounting for such recurring disturbance events would seriously overestimate the ability of these forests to absorb atmospheric CO 2 .

Conclusions
This study showed, with the aid of MODIS NDVI and eddy covariance data, a substantial loss in regional GPP due to insect-induced defoliation in subarctic deciduous forests in northern Fennoscandia.The estimated mean annual decrease in regional GPP for a year with insect outbreak was 15 ± 5 Gg C yr −1 ± in the study area of 100 km 2 .This should be compared with the average annual GPP of 41 ± 12 Gg C yr −1 for years without disturbances.In the most severe outbreak year (2012), 76 % of the birch forests were defoliated and annual GPP was merely 50 % of GPP for years without disturbances.
The study also demonstrated the use of remote sensing data both to monitor the spatial extent of the defoliation and to estimate the impact on the primary productivity of these defoliation events.The insect disturbance is shown to have major impacts on the primary production of the subarctic forest; consequently, the derived methods, based on combining remote sensing and eddy covariance measurements, are of major importance to support carbon balance estimates over large areas.The Supplement related to this article is available online at doi:10.5194/bg-14-1703-2017-supplement.

Figure 1 .
Figure 1.The studied birch forest (green) along the south-west part of Lake Torneträsk and in the valley to the south-west of the village of Abisko.The locations of the eddy covariance (EC) tower used to obtain GPP and the spectral tower used to obtain fAPAR data are also shown.Reference system is SWEREF99 TM and latitude and longitude are in WGS84.Source of background map: Lantmäteriet (Dnr: I2014/00579).

Figure 3 .
Figure 3. Mean (black) and standard deviation (grey) of NDVI 2010-2014 for the four pixels around the eddy covariance (EC) tower.Panel (a) is raw NDVI and (b) is NDVI fitted with double logistic functions in TIMESAT (NDVI DL ).Years 2012 and 2013 are those with insect outbreak.In 2013 the birch forest was refoliating later in the growing season.There is small secondary peak in raw NDVI (a) appearing each year.This peak appears during the winter, when there is no vegetation in the study area, and is hence removed from the smoothed data (b).

Figure 4 .
Figure 4. Correlation between ground-measured fraction of absorbed PAR by the canopy (fAPAR canopy ) and MODIS-derived NDVI smoothed with a double logistic function in TIMESAT (NDVI DL ) in 8-day intervals.Only NDVI DL values ≥ 0.4 were included in the OLS regression resulting in the black line.R 2 = 0.81 and N = 29.

Figure 5 .
Figure 5. Influence on RMSE and R 2 of GDD thres (a), T thres (b), and the optimal period to change from the first to the second f 8day model (c).RMSE is computed from mean of daily GPP over 8-day periods.

Figure 6 .
Figure 6.Light use efficiency (ε), NDVI fitted with double logistic functions (NDVI DL ) scaled ×2 (green), and PAR (orange) for the 6 years with data from the EC tower.Black lines with error bars and black circles are the light use efficiency values included when ε max and ε max, def were computed for undisturbed and defoliated years, respectively.The error bars are symmetric and 1 standard deviation higher or lower than the mean values.

Table 1 .
Annual GPP derived from the EC tower for the 5 years without insect outbreak and the year 2012 with insect outbreak.

Figure 9 .
Figure 9. Reduction in annual GPP (g C m −2 yr −1 ) due to the outbreak of autumnal moth and winter moth in 2012 computed with a LUE model also for defoliation (Method 2).One standard deviation of the GPP losses is estimated to 35 % of the given values.Areas with only the background map have a canopy cover less than 50 % or are outside the study area shown in Fig. 1.The reference system is SWEREF99 TM and latitude and longitude are in WGS84.Source of background map: Lantmäteriet (Dnr: I2014/00579).
Data availability.The fAPAR and EC data used in this study are available from the corresponding author upon request.The meteorological data from Abisko Scientific Research Station are available from the research station: http://polar.se/abisko-naturvetenskapliga-station/.The MODIS data used are available from the references LPDAAC (2016a, b).
Senn et al.'s (1992)pping the reduction in gross primary productivity2.3.1 First part of the growing seasonDuring the first part of the growing season, covering May to late June, f 8day depended on growing degree days (GDDs) and frost events, where GDD was computed with a base temperature of 5 • C, followingSenn et al.'s (1992)method applied to mountain birch development in northern Finland: : one model for the first part of the growing season and one model for the second part of the growing season.P

Table 2 .
Defoliated area (km 2 ) and annual reduction in GPP * (Gg C yr −1 ) for the 3 years with insect defoliation since the year 2000.The total area with forest cover is 100 km 2 .