Understanding the effect of disturbance from selective felling on the carbon dynamics of a managed woodland by combining observations with model predictions

The response of forests and terrestrial ecosystems to disturbance is an important process in the global carbon cycle in the context of a changing climate. This study focuses on the effect of selective felling (thinning) at a managed forest site. Previous statistical analyses of eddy covariance data at the study site had found that disturbance from thinning resulted in no significant change to net ecosystem carbon uptake. In order to better understand the effect of thinning on carbon fluxes, we use the mathematical technique of four‐dimensional variational data assimilation. Data assimilation provides a compelling alternative to more common statistical analyses of flux data as it allows for the combination of many different sources of data, with the physical constraints of a dynamical model, to find an improved estimate of the state of a system. We develop new observation operators to assimilate daytime and nighttime net ecosystem exchange observations with a daily time step model, increasing observations available by a factor of 4.25. Our results support previous analyses, with a predicted net ecosystem carbon uptake for the year 2015 of 426 ± 116 g C m−2 for the unthinned forest and 420 ± 78 g C m−2 for the thinned forest despite a model‐predicted reduction in gross primary productivity of 337 g C m−2. We show that this is likely due to reduced ecosystem respiration postdisturbance compensating for a reduction in gross primary productivity. This supports the theory of an upper limit of forest net carbon uptake due to the magnitude of ecosystem respiration scaling with gross primary productivity.


Introduction
The response of forests and terrestrial ecosystems to disturbance (e.g., felling, fire, or insect outbreaks) is one of the least understood components in the global carbon cycle [Ciais et al., 2014]. Current land surface models fail to represent the effect of disturbances on long-term carbon dynamics [Running, 2008], although these RESEARCH ARTICLE 10.1002/2017JG003760 Key Points: • No change was found in forest NEE after selective felling • Data assimilation was used to investigate the effect of this disturbance • Concurrent reductions in respiration and GPP allow for unchanged NEE Supporting Information: • Supporting Information S1 • Figure S1 Correspondence to: E. M. Pinnington, e.m.pinnington@pgr.reading.ac.uk Thinning is a silvicultural practice used to improve ecosystem services or the quality of a final tree crop and is globally widespread. The effect of thinning on carbon budgets has largely been ignored [Liu et al., 2011]. Thinning has been shown to increase the basal growth increment of remaining trees due to increased light and water availability which may indicate increased net primary productivity in subsequent years [Bréda et al., 1995;Martín-Benito et al., 2010]. However, Misson et al. [2005] found that the immediate effect of thinning can change an ecosystem from a sink to a source of CO 2 , due to reduced gross primary productivity (GPP) following a reduction in total leaf area and unchanged or heightened ecosystem respiration. Other studies, analyzing flux tower eddy covariance records, find no significant change in the observed net ecosystem exchange (NEE) of CO 2 after thinning [Vesala et al., 2005;Moreaux et al., 2011;Dore et al., 2012;Saunders et al., 2012;Wilkinson et al., 2015]. These studies suggest that this is due to increased light availability and reduced competition allowing ground vegetation to display increased GPP and compensate for an increase in heterotrophic respiration postdisturbance.
Other studies have shown a significant reduction in the carbon content of rhizosphere soils following tree felling [Hernesmaa et al., 2005]. It has been shown that tree roots provide a rhizosphere priming effect, greatly increasing the rate of soil organic carbon decomposition [Dijkstra and Cheng, 2007]. This is consistent with previous work carried out at the study site in this paper, where it has been shown that the magnitude of ecosystem respiration is strongly coupled to the magnitude of GPP [Heinemeyer et al., 2012]. Predictions made by Kurz et al. [2008] about the impacts of mountain pine beetle outbreaks in Northern American forests suggested a switch from sink to source of carbon following this disturbance. However, the analysis of a diverse set of observations for an area with approximately 70% infested trees by Moore et al. [2013] revealed little change in net CO 2 flux, due to concurrent reductions in GPP and ecosystem respiration. Similar results were also found from large-scale tree girdling experiments [Högberg et al., 2001], where 1-2 months after girdling a 54% decrease in soil respiration was observed.
Here we used data assimilation, which is a mathematical technique for combining observations with prior model predictions in order to find the best estimate of a dynamical system. Functional ecology models have been combined with many different observations relevant to the carbon balance of forests [Quaife et al., 2008;Fox et al., 2009;Zobitz et al., 2011;Richardson et al., 2010;Zobitz et al., 2014;Niu et al., 2014;Pinnington et al., 2016], leading to improved estimates of model parameter and state variables and reduced uncertainty in model predictions. There have been many efforts to model the effect of disturbance on forest ecosystems [Thornton et al., 2002;Seidl et al., 2011], with a growing number of dynamic global vegetation models [Sitch et al., 2008], some of which explicitly model the impact of disturbance, e.g., Moorcroft et al. [2001]. However, the use of data assimilation has been limited to a few examples, all of which used satellite data [Hilker et al., 2009;Kantzas et al., 2015]. The authors are not aware of any studies assimilating site level data to quantify disturbance effects. By assimilating observations relevant to postdisturbance ecosystem carbon dynamics with prior model predictions of ecosystem behavior, we can analyze the retrieved parameters after data assimilation to find the model-predicted effects of disturbance.
In this paper we investigate the effect of thinning on the carbon dynamics of the Alice Holt flux site , a deciduous managed woodland, following an event in 2014, when one side of the site was thinned and the other side left unmanaged. We present new observation operators for the assimilation of daytime and nighttime NEE observations with a daily time step model, in this case the Data Assimilation Linked Ecosystem Carbon (DALEC2) model [Bloom and Williams, 2015]. These methods require no model modification. We combine all available observations for 2015 with prior model predictions to find two sets of optimized model parameter and initial state values, corresponding to thinned and unthinned sides of the forest. We then use these two versions of the model to seek to explain why the net uptake of carbon remains unchanged even after removing a large proportion of the trees from one side. We find a net ecosystem carbon uptake for the year 2015 of 426 ± 116 g C m −2 for the unthinned forest and 420 ± 78 g C m −2 for the thinned forest, despite a reduction in GPP of 337 g C m −2 for the thinned forest when compared to the unthinned forest. We find that reduced ecosystem respiration for the thinned forest allows for this unchanged net carbon uptake.

10.1002/2017JG003760
The data assimilation techniques presented in this paper could be applied for similar analyses at other sites and provide a novel method to help elucidate the reasons behind ecosystem responses.

Alice Holt Research Forest
Alice Holt Forest is a research forest area managed by the UK Forestry Commission located in Hampshire, SE England. Forest Research has been operating a CO 2 flux measurement tower in a portion of the forest, the Straits Inclosure, continuously since 1998. The Straits Inclosure is a 90 ha area of deciduous broadleaved plantation woodland located on a surface water gley soil and has been managed for the past 80 years. The majority of the canopy trees are oak (Quercus robur L.), with an understory of hazel (Corylus avellana L.) and hawthorn (Crataegus monogyna Jacq.), but there is a small area of conifers (Pinus nigra ssp. laricio (Maire) and P. sylvestris L.) within the tower measurement footprint area depending on wind direction. Further details of the Straits Inclosure site and the measurement procedures are given in Wilkinson et al. [2012], together with analysis of stand-scale 30 min average net CO 2 fluxes (NEE) from 1998 to 2011.
As part of the management regime, the Straits Inclosure is subject to thinning, whereby a proportion of trees are removed from the canopy in order to reduce competition and improve the quality of the final tree crop. At the Straits, an intermediate thinning method is used with a portion of both subdominant and dominant trees being removed from the stand to stimulate the growth of remaining dominant trees [Kerr and Haufe, 2011]. The whole of the stand was thinned in 1995. Subsequently, the eastern side of the Straits was thinned in 2007 and then the western side in 2014. The flux tower at the site is situated on the boundary between these two sides. This allows for the use of a footprint model to split the flux record and thus analyze the effect of this disturbance on carbon fluxes at the site. In Wilkinson et al. [2015] a statistical analysis of the eddy covariance flux record found that there was no significant effect on the net carbon uptake of the eastern side after thinning in 2007. In this paper we focus on the effect of disturbance on the western side after thinning in 2014. We therefore refer to the western side as "thinned" forest and the eastern side as "unthinned" forest, although the "unthinned" forest was previously thinned in 2007 and so will have a different structure to a completely unmanaged forest.

Observations
In order to assess the effect the 2014 thinning had on the Straits Inclosure, an intensive field campaign was undertaken by the lead author in 2015 to measure leaf area index and also to estimate standing woody biomass. From the site we also have a long record of flux data, as discussed in section 2.1. These observations span both the thinned and unthinned sides of the forest.

Leaf Area Index
To assess the impact of the 2014 thinning, three transects were established in the Straits Inclosure for intensive sampling during 2015. A total of 435 sampling points were marked at 10 m apart, using a GPS and fluorescent tree spray paint. Measurements of peak leaf area index (LAI) (July 2015 to September 2015) were made using both a ceptometer and hemispherical photography. The transects were walked twice with the ceptometer taking readings at every sampling point, giving 870 readings in total. Hemispherical photographs were taken every 50 m as shown in Figure 1, giving 89 photographs in total.
We measured below-canopy photosynthetically active radiation (PAR) using the ceptometer while logging above-canopy PAR using a data logger and PAR sensor positioned outside the canopy. We then estimated LAI using the above-canopy and below-canopy PAR readings [Fassnacht et al., 1994]. For the hemispherical photographs, we used the HemiView software [Rich et al., 1999] which calculates the proportion of visible sky as a function of sky direction (gap fraction) which it then uses to calculate LAI [Jonckheere et al., 2004].
Six litter traps were also established at points along the transects (positions shown in Figure 1) allowing for comparison with the other methods. These were sampled throughout the leaf fall season in 2015. We found that the LAI derived from the litter traps was always greater than LAI estimated from optical methods, as expected [Bréda, 2003]. From the sampling of the litter traps we also have estimates of leaf mass per leaf area for use in data assimilation. As the six litter traps are not enough to describe the LAI for the research site [Kimmins, 1973], we used estimates from the ceptometer and hemispherical photographs for data assimilation. We took the weighted average (dependent on number of observations taken of each type) of the hemispherical photograph and ceptometer estimated LAI and derived an LAI of 4.42 with a standard error of 0.07 for the eastern unthinned forest and an LAI of 3.06 with a standard error of 0.07 for the western thinned forest. We assimilated the mean of 299 LAI observations in the unthinned and 225 in the thinned section of forest. Consequently, the appropriate representation of error for data assimilation is the standard error of the mean. From our litter trap observations we find a leaf mass per area of 29 g C m −2 free soluble carbohydrates for both sides of the forest.

Woody Biomass
The method of point-centered quarter (PCQ) was used to conduct a biomass survey as specified in Dahdouh-Guebas and Koedam [2006]. Along the three transects 114 points were sampled measuring the diameter at breast height (DBH) and the density of trees. We then used allometric relationships between DBH and total aboveground biomass and coarse root biomass, found in work carried out by Forest Research and in McKay et al. [2003], to find an estimate of total woody and coarse root carbon (referred to as C woo in the DALEC2 model). These observations are shown in Table 1.
Forest Research has carried out its own mensuration studies at the site. One such study of the western thinned forest (at a similar time to our own PCQ measurements) found a tree density of 225 ha −1 and an average DBH of 32 cm, which are in close agreement to the estimates in Table 1. This gives us confidence that earlier measurements taken by Forest Research before the thinning are representative of the methods we have used. Measurements of the same section of forest from 2009 found a tree density of 418 ha −1 and an average DBH of 28 cm. This suggests that approximately 46% of trees have been removed during the 2014 thinning. From these estimates we can also see the effect that thinning has on the type of trees found at the site. The trees per hectare has dropped dramatically after thinning, but the mean DBH has increased, because the smaller subdominant trees have been removed. The greater mean DBH of the eastern unthinned section, 34 cm, indicates that the thinning that took place in 2007 has allowed the dominant trees to grow as a result of reduced competition.

Flux Tower Eddy Covariance
The Straits Inclosure flux tower provides half-hourly observations from January 1999 to December 2015. These consist of the NEE fluxes and meteorological driving data of temperature, irradiance, and atmospheric CO 2 concentration for use in the DALEC2 model. The NEE data were subject to u * filtering (with a value of 0.2m s −1 ), where data corresponding to low friction velocity values were removed from the data set, and quality control procedures as described by Papale et al. [2006] but were not gap filled. The resultant half-hourly NEE data set was then split between observations corresponding to the western thinned and eastern unthinned sides of the site using a flux footprint model. This model is dependent on wind speed and direction to calculate the location that the flux tower is sampling. Partitioning the NEE data set in this way reduces the total number of available NEE observations (see Wilkinson et al. [2016] for more details).
To match the time step of our model, we computed daily NEE observations by taking the mean over the 48 measurements made each day, selecting only days where there are no missing data. As we have been strict on the quality control of the flux record and not used any gap filling, this presented a problem in terms of the number of daily NEE observations available. By further splitting the flux record between two sides, we retrieved very few total daily observations of NEE for either side. In order to address this, we computed daytime and nighttime NEE fluxes (NEE day and NEE night , respectively) for use in data assimilation. We used a solar model to define whether NEE measurements were made at daytime or nighttime. We then took the mean over the half-hourly day or nighttime measurements, again only taking periods where there were no gaps in the data so that we were only considering true observations. This provided us with a factor of 4.25 more observations of NEE for assimilation, as seen in Table 2. Because we are averaging over shorter time periods, we have a smaller probability of gaps and erroneous data. We see that we have more daytime NEE observations than nighttime, as we tend to have much more turbulent air mixing in daylight hours. Times of low turbulent mixing lead to an underestimation in flux values. In section 2.3.2 we give details of how we relate these twice daily observations of NEE to a daily time step model.
The errors in observations of daily NEE were assumed to be constant and set at 0.5 g C m −2 day −1 by Williams et al. [2005], whereas Braswell et al. [2005] found these errors to be closer to 1 g C m −2 day −1 . A more recent study finds a mean value of 0.8 g C m −2 day −1 for NEE flux errors [Post et al., 2015]. However, Richardson et al. [2008] show that flux errors are heteroscedastic. To account for the heteroscedastic nature of NEE errors, we define an error function that scales between 0.5 and 1 g C m −2 day −1 based on the magnitude of the observation. This function is defined as 0.5 + 0.04|NEE i day | g C m −2 day −1 , where |NEE i day | is the magnitude of the daytime NEE observation on day i; this function is derived by considering the maximum and minimum |NEE i day | values. Raupach et al. [2005] comment that nighttime measurements of NEE are much more uncertain than daytime measurements. This is difficult to quantify, but here we assume that nighttime flux errors are 3 times the magnitude of daytime errors. We therefore have the error function of 1.
night | is the magnitude of the nighttime NEE observation. We also include correlations in time between the errors in our observations of NEE, as discussed in Pinnington et al. [2016] (see Figures S6 and S7 in the supporting information).

Model and Data Assimilation 2.3.1. DALEC2 Ecosystem Carbon Model
The DALEC2 model is a simple process-based model describing the carbon dynamics of a forest ecosystem [Bloom and Williams, 2015]. The model is constructed of six carbon pools (labile (C lab ), foliage (C fol ), fine roots (C roo ), woody stems and coarse roots (C woo ), fresh leaf and fine root litter (C lit ), and soil organic matter and coarse woody debris (C som )) linked via fluxes. The aggregated canopy model (ACM) [Williams et al., 1997] is used to calculate daily gross primary production (GPP) of the forest, taking meteorological driving data and the modeled leaf area index (a function of C fol ) as arguments. Figure 2 shows a schematic of how the carbon pools are linked in DALEC2; full model equations can be found in section A1.

Data Assimilation
We implement four-dimensional variational data assimilation (4D-Var) with the DALEC2 model for joint parameter and state estimation [Navon, 1998]. In 4D-Var we aim to find the parameter and initial state values such that the model trajectory best fits the data over some time window, given some prior information about the system. This prior information takes the form of an initial estimate of the parameter and state variables of the model, x b , valid at the initial time. This prior is assumed to have unbiased, Gaussian errors with known covariance matrix B. Adding the prior term ensures that our problem is well posed and that we can find a locally unique solution [Tremolet, 2006]. We aim to find the parameter and initial state values that minimize the weighted squared distance to the prior and the weighted squared distance of the model trajectory to the observations, over a time window of length N, with individual time points t 0 , … , t N . We do this by finding the state x a at time t 0 that minimizes the cost function where x 0 is the vector of parameter and initial state values to be optimized, x i is the vector of model variables at time t i , h i is the observation operator mapping the parameter and state values to the observations, y i is the vector of observations at time t i , and R i is the observation error covariance matrix. The time step, i, is 1 day in this case. Further details of the implemented data assimilation scheme and specification of prior and observational errors can be found in Pinnington et al. [2016].
In this paper we assimilate daytime and nighttime NEE in order to increase the number of observations available to us and also better partition our modeled estimate of GPP and total ecosystem respiration. As the DALEC2 model runs at a daily time step, this requires us to relate the daily parameter and state values from the model to the twice daily observations of NEE. We do this by writing two new observation operators, one relating the model state and parameters to daytime NEE and the other to nighttime NEE. The NEE of CO 2 at any given time is the difference between GPP and ecosystem respiration. For an observation of total daily NEE on day i we have where Ψ represents meteorological driving data used in the calculation of GPP, f auto is the fraction of autotrophic respiration, lit is the litter carbon turnover rate, som is the soil and organic carbon turnover rate, Θ is the temperature dependence exponent factor, and T i is the mean temperature over 24 h. Further description can be found in section A1. The first term in equation (2) represents gross primary productivity, the second autotrophic respiration, and the third and fourth terms heterotrophic respiration.
For total daytime NEE we have where day is number of daylight hours 24 , and T i day is the mean temperature over daylight hours. Here we still have the same term for GPP as in equation (2) as all photosynthesis occurs during daylight hours. We have made the assumption that respiration is spread uniformly in time; therefore, the respiration terms are scaled by the fraction of daylight hours. For nighttime NEE we have where night is number of night hours 24 , and T i night is the mean nighttime temperature. In equation (2) we do not have a term for GPP as no GPP will occur during the night. The respiration here is scaled by the fraction of nighttime hours. The length of day and night are calculated using a solar model.
These new observation operators allow for assimilation of daytime/nighttime NEE without the need for altering the model and can be applied to other ecosystem models to allow for the assimilation of eddy covariance data at a finer temporal resolution.

Experimental Setup
In order to assess the information content of the three available data streams (described in section 2.2) and their impact on the effect of disturbance as predicted by the model, we conducted a data denial procedure. This involved assimilating different combinations of observations, in three experiments, as shown in Table 3. In our first experiment, Experiment A, we used only the eddy covariance data, as this is the data type most commonly used in data assimilation studies. In the second, B, we added the observations relating to leaf mass and area, and finally in the third experiment, C, we added the observations of woody biomass, as NEE observations have been shown to be unable to constrain this [Fox et al., 2009]. In each experiment we used the prior model as specified in the appendix in Table A1 and Figure S3 (supporting information). This prior model was found by assimilating daytime and nighttime NEE, leaf mass per area, and LAI observations from 2012 and 2013 before the thinning occurred. More information on the data assimilation methods used to find this prior model can be found in Pinnington et al. [2016].
In each experiment we ran the assimilation for both the thinned forest and the unthinned forest, using the distinct data for each side. This allowed us to retrieve a unique set of parameter and initial state values for each section of forest. We analyzed the optimized parameter and initial state values for the thinned and unthinned forest and also the model predictions of different variables for each side postdisturbance. This allowed us to judge the effect the thinning in 2014 had on the carbon dynamics of the forest in 2015.
It would be expected that we will retrieve different estimates for each of the experiments outlined in Table 3, with our most confident estimate being when all observations types are assimilated together in experiment C. This would allow us to see how much information each data stream provides and assess whether NEE data alone are enough to understand the effect of disturbance.

Results
In Figures 3 and 4 we show the observations and model trajectories after assimilation for the thinned and unthinned forest for experiments A and C, respectively. We can see that performing the data assimilation has allowed the model to fit all the assimilated observations well for both experiments (in experiment A only NEE observations are assimilated). This can also be seen by the reductions in root-mean-square error compared to the prior in Table 4. From Figures 3a and 3b we see that the modified observation operators presented in section 2.3.2 have allowed our model to represent both daytime and nighttime NEE well.

Comparison of Experiments
In experiment A we have only assimilated NEE observations. From Table 4 we can see that we improve the fit to the assimilated observations for both the unthinned and thinned forest when compared to the prior model. The root-mean-square error (RMSE) is within the specified observation error for both daytime and nighttime NEE after assimilation. By only assimilating observations of NEE, we have not been able to accurately predict LAI. Although we have improved the fit of the model to LAI after assimilation for the thinned forest (see Table 4), we have significantly degraded the fit of the model to LAI for the unthinned forest. Partitioning the NEE data set between the thinned and unthinned forest (as described in section 2.2.3) has resulted in a gap in the observations for the unthinned forest during the period of greatest carbon uptake (June 2015 to August 2015); see Figure 3a. This is due to the prevailing wind in this period being from the southwest. This gap is potentially causing an underestimation of NEE for the unthinned forest. From Figure 3d and Table 4 we can see that NEE observations alone do not give us enough information to recover a value of C woo with the DALEC2 model. This is also found in previous work [Fox et al., 2009].
In experiment B we have assimilated observations of NEE, LAI, and leaf mass per area. From Table 4 we see that including the extra observations has allowed the model to fit LAI well for both the unthinned and thinned forest, and although the fit of the model to the NEE observations is slightly degraded compared to experiment A, it is still well within the specified observation error from section 2.2.3. Table 4 also shows that including these extra observations in experiment B still does not allow us to recover an accurate value of C woo . Further results from experiment B can be found in the supporting information ( Figures S1 and S2).
In experiment C we assimilate all available observations. This gives us very similar results as in experiment B, except that including the observations of C woo in the assimilation allows the model to fit this observation well, as seen in Table 4. We see from Figures 4a and 4c that including observations of LAI in the assimilation removes the issue of underprediction of NEE for the unthinned forest, as discussed for experiment A. The distinct difference in stand structure is now clear in Figure 4, with reduced LAI and woody carbon for the thinned forest. Table 5 shows the cumulative annual fluxes for the year 2015 for the three experiments. All three experiments show only small differences in the net ecosystem carbon uptake between the thinned and unthinned forest. The area common to the distributions of the thinned and unthinned forest NEE in Table 5 was found to be 89%, 79%, and 81% in experiments A to C respectively, calculated using the Weitzman overlap measure [see Inman and Bradley, 1989]. However, the partitioning of this carbon uptake between GPP and total ecosystem respiration (TER) is markedly different (distribution overlap ≪ 1%), with experiment A predicting increased TER and GPP after thinning and experiment C (and B) predicting reduced TER and GPP after thinning. This can be seen more clearly in Figure 5. The difference between the results of experiments A and C highlights the issue that NEE is the difference between two large fluxes (NEE = −GPP + TER) and we can therefore find an accurate prediction of NEE despite underestimating/overestimating both GPP and TER. Therefore, care should be taken when interpreting model results based solely on NEE data, especially in this case, as we have seen that the partitioning of the NEE data between the thinned and unthinned forest has introduced a bias into our data set. If we were to base our analysis on experiment A, we would assume that the thinning had caused an increase in ecosystem respiration and that this had been compensated for by an increase in GPP. This is the   opposite conclusion to the one we find in experiment C when we include observations relating to the structure of the forest. Table 5 also shows that in some cases adding more observations reduces the uncertainty in model-predicted annual fluxes. However, in other cases we see the opposite effect. We believe that this is due to the assimilation overfitting to the observations (see section 4.1) and therefore underpredicting the modeled uncertainty in NEE.

Partitioning of Carbon Fluxes
In Figure 6 we show the partitioning of modeled cumulative ecosystem respiration for the year 2015 between total autotrophic respiration and heterotrophic respiration from litter and soil for both the unthinned and thinned forest in experiment C. The DALEC2 model represents autotrophic respiration as a constant fraction of GPP. From Figure 6 we can see the strong dependance of autotrophic respiration on GPP with the growth rate being much greater between June 2015 and September 2015 (when GPP will be of greater magnitude). For heterotrophic respiration the growth rate is more constant throughout the whole year. Total ecosystem respiration is reduced by 331 g C m −2 for the thinned forest when compared to the unthinned forest, with reductions in both heterotrophic and autotrophic respiration of 169 g C m −2 and 162 g C m −2 , respectively.   Figure 7 shows the change in parameter and initial state values for the thinned and unthinned forest after assimilating all observations in experiment C. It is important to note that this is the difference when compared to our prior model estimate, which was found by assimilating only eddy covariance, LAI, and leaf mass per area observations from 2012 and 2013. We therefore expect changes in parameter and state values for both the thinned and unthinned forest, as we are assimilating new data streams. This is particularly noticeable in the initial carbon pool state variables in Figure 7. Constraints on the carbon pool state variables are provided by the assimilated observations of woody biomass and coarse roots (C woo ), LAI, and leaf mass per area (c lma ). LAI and c lma give us a constraint on foliar carbon (C fol ) as LAI = C fol c lma . We can see that the values for the model-predicted carbon pools are as we might expect with the thinned forest having less carbon in all pools when compared to the unthinned forest. For litter carbon (C lit ) we expect a reduction in input of leaf litter for the thinned forest, and, although there might be increased woody debris after thinning, this is much less readily decomposed and so possibly has little impact in the year after thinning . The difference in predicted soil carbon content (C som ) between the thinned and unthinned forest is consistent with studies analyzing soil carbon contents after felling [Hernesmaa et al., 2005]. For the parameters the biggest changes appear to be in the litter carbon turnover rate parameter ( lit ), with the retrieved parameter being significantly reduced for the unthinned forest when compared to the thinned forest. However, we still see reduced total litter respiration in Figure 6 for the thinned forest compared to the unthinned forest. This is due to the significant difference in litter carbon content (C lit ) for both sides, with the unthinned forest having a much higher litter carbon content than the thinned forest. The large change in the lit parameter between the two sides is therefore compensating for an overestimated difference in litter carbon content between the two sides.

Twin Experiments
Experiment C should give us the best possible results as we have assimilated all available data. To ensure that this is the case, we have run three "twin" experiments (described in the supporting information), where we aim to estimate a set of known model parameters and initial state variables (referred to as the model "truth") by assimilating synthetic observations generated from a normally distributed random sample around the known model mean. We use the same combination of observations as in experiments A to C but generate the observations from the model truth. The results from these experiments are shown in Figure S5 and Table S1 for the unthinned and thinned forest. Explanation of parameter and state variable symbols in "" Table A1. in the supporting information. From these experiments we find the smallest error in parameter and initial state values in twin experiment C where all synthetic observations are assimilated. The error in parameter and initial state values is reduced by 28% in twin experiment C, compared to the results when only NEE observations are assimilated in twin experiment A. Although this does not prove that experiment C will give us the correct results (as twin experiments are based on modeled observations which will not directly reflect physical observations and so have limitations; see Errico et al. [2013]), it lends confidence that the best results will be achieved when all observations are used.

Discussion
In this paper we have investigated the possible explanations for the observation that a thinning event, where approximately 46% of trees were removed from the study site, had no impact on net ecosystem carbon uptake. We used data assimilation to combine observations and prior model predictions of ecosystem carbon balance in order to understand how the state of an ecosystem might be altered after a disturbance event. We have confidence in the optimized model prediction of NEE as we have demonstrated previously that assimilating a single year of data can accurately forecast the carbon uptake of the site for a long time period (15 years) [Pinnington et al., 2016].

Impacts of Increased Data Steams
We conducted three experiments assimilating different combinations of available data streams. For all experiments we find no significant change in net carbon uptake for the studied ecosystem following stand thinning. This is consistent with other studies of ecosystem carbon dynamics following thinning [Vesala et al., 2005;Moreaux et al., 2011;Dore et al., 2012;Wilkinson et al., 2016].
We find different reasons for this unchanged carbon uptake dependent on which data streams are assimilated. When only assimilating NEE observations in experiment A, we find increased ecosystem respiration and increased GPP postdisturbance. These results are unreliable due to bias introduced into the NEE data set from partitioning between the thinned and unthinned forest. It is likely that the gap in the NEE data set for the unthinned forest when observations would have been of highest magnitude has caused our model to underpredict the carbon uptake for the unthinned forest. This is a potential explanation for the large underprediction of LAI for the unthinned forest in experiment A, as seen in Figure 3c. It is also likely that in experiment A we are overfitting to the NEE data as there are no additional independent observational constraints. This means that unrealistic parameter values can be retrieved in order to find the best possible fit to the observations. This overfitting is potentially giving us an unrealistically low posterior uncertainty for our optimized parameters and model-estimated carbon fluxes (shown in Table 5). Introducing other independent observations in experiment B has reduced the problem of overfitting. MacBean et al. [2016] have shown that for nonlinear data assimilation, assimilating all available data concurrently is superior to assimilating individual data streams in sequence.
We can see from Table 5 that both experiments B and C predict very similar cumulative fluxes, suggesting that the assimilated observations of C woo have not had much impact on the model carbon dynamics for this time period. Because the rate parameters controlling this pool are relatively slow, it is likely that observations of C woo will become much more important over longer time scales. Here we have only assimilated a single observation of C woo for either side of the forest; if multiple observations of C woo were available throughout time, this would give us an estimate of the rate of woody biomass accumulation, providing an important constraint on the carbon assimilation of the forest. For experiment C the time of senescence in LAI predicted by the model is consistent with PhenoCam observations made by Forest Research at the site, as shown in the supporting information ( Figure S4). However, the time of green-up in LAI predicted by the model is later than the PhenoCam observations. We hypothesize that this is due to the model predicting the photosynthetically effective LAI implicitly, rather than the LAI related to canopy green index measured by the PhenoCam. The PhenoCam will show a peak in green index LAI before new leaves become competent at photosynthesis [Reich et al., 1991;Morecroft et al., 2003;Klosterman et al., 2014].
From our most confident estimate, where all available observations are assimilated, the model shows that reductions in GPP, following a decrease in total leaf area postthinning, are being offset by simultaneous reductions in ecosystem respiration. This is in contrast to current suggestions that reduced canopy photosynthesis is compensated for by increased GPP by ground vegetation postthinning [Vesala et al., 2005;Moreaux et al., 2011;Dore et al., 2012;Wilkinson et al., 2016]. It is important that more independent data relating to both productive and respiration fluxes are sought to further verify the results of this study.

GPP and Respiration Are Closely Linked
Our results show a decrease in both autotrophic and heterotrophic respirations following thinning. We follow the definition of Heinemeyer et al. [2012] and characterize belowground autotrophic respiration as respiration from roots, mycorrhizal fungi, and other microorganisms dependent on the priming of soils with labile carbon compounds from roots. Heterotrophic respiration is respiration by microbes not directly dependent on autotrophic substrate; however, the largest fraction of heterotrophic respiration is based on the decomposition of young organic matter (e.g., leaves and fine roots) whose availability also depends on the GPP of an ecosystem [Janssens et al., 2001]. We find similar decreases in both heterotrophic and autotrophic respiration for the thinned forest when compared with the unthinned forest. While it has been shown that heterotrophic respiration can decrease after disturbance events [Bhupinderpal et al., 2003], it is possible that we overestimate the reduction in heterotrophic respiration and underestimate the reduction in autotrophic respiration. This is understandable as we have assimilated no data on this partitioning. Also, our model description of autotrophic respiration is simple (described as a constant fraction of GPP), and therefore the heterotrophic respiration component of the model might compensate and in this instance describe the behavior of mycorrhizal fungi and other microbes commonly categorized in the autotrophic component of respiration.
In a study measuring soil CO 2 fluxes over 4 years at the Straits Inclosure (the study site in this paper) Heinemeyer et al. [2012] showed a large contribution (56%) of autotrophic respiration (characterized as root and mycorrhizal respiration) to total measured soil respiration. Heinemeyer et al. [2012] also suggested that mycorrhizal fungi play a role in priming the turnover of soil organic carbon by other microbes, with evidence from Talbot et al. [2008]. Högberg and Read [2006] find similar figures for the autotrophic contribution to total soil respiration, with around half or more of all soil respiration being driven by recent photosynthesis. Heinemeyer et al. [2012] discuss the possibility of this tight coupling between GPP and ecosystem respiration leading to an upper limit for forest CO 2 uptake due to increased GPP leading to increased respiration, which is also discussed by Heath et al. [2005]. Our results support this hypothesis, as ecosystem respiration scales with GPP after approximately 46% of trees are removed from the study site, meaning that we find no significant change in net ecosystem carbon uptake after thinning.
Studies analyzing eddy covariance flux records also find no significant change in the net ecosystem exchange of CO 2 after thinning [Vesala et al., 2005;Moreaux et al., 2011;Dore et al., 2012;Wilkinson et al., 2016]. Vesala et al. [2005] used a model of light interception and ground vegetation photosynthesis to show that the unchanged NEE was due to increased GPP by ground vegetation (following increased light availability and reduced competition). This compensated for increases in heterotrophic respiration and reduced canopy photosynthesis postthinning. Similar conclusions were drawn by Moreaux et al. [2011] who destructively sampled ground vegetation and showed an increase in biomass postthinning. We do not find evidence to support such conclusions and instead suggest that reduced ecosystem respiration is the most important component for the unchanged NEE of the forest following thinning. However, it is important to note that our measurements of LAI are made at approximately 1 m above the forest floor, which means that our measurements of LAI do not account for ground vegetation. Therefore, any effect of this ground vegetation is not simulated by our model. Despite this, observations made during multiple walks of the three established transects find no evidence of increased ground vegetation in the year after thinning. In fact, much of the ground vegetation and subcanopy was removed during thinning and did not appear to have recovered in the following year. At longer time scales regrowth of the subcanopy and ground vegetation will play an important role in increased productivity. Our results suggest that this increased productivity would also be met with subsequent increases in ecosystem respiration.

The Use of Data Assimilation to Predict Disturbance Effects
Data assimilation provides a valuable alternative to more common statistical analyses of flux data records to calculate constituent ecosystem carbon fluxes, as it allows for the combination of many distinct data streams and the dynamics of a physically meaningful model to construct our solution. We have shown that basing our results on a single data stream (NEE) can give us much different conclusions than when all data streams are assimilated. Through the use of data assimilation, we find different results for the effect of thinning on ecosystem carbon fluxes than in previous statistical analyses of flux tower data at the same study site by Wilkinson et al. [2016]. Wilkinson et al. [2016] found increased ecosystem respiration after a thinning event at Alice Holt in 2007 and suggested that this was counteracted by an increase in GPP by ground vegetation. After assimilating all available data streams, we instead find that a thinning event in 2014 led to reductions in both TER and GPP.
In this work we have been strict on the quality control procedures for NEE data, and have allowed no gap filling, to ensure that we base our results on only the best quality true observations. This and the partitioning of the NEE data set between the thinned and unthinned areas of forest has resulted in a distinct data gap during the growing season for the thinned site NEE, as discussed in section 3. From the data assimilation experiments conducted we have shown that the combination of multiple distinct data steams has helped to reduce the impact of this data gap on assimilation results.
The effect of disturbance is poorly characterized in current land surface and global climate models [Running, 2008]; it is important to better understand how parameters and carbon pools might change following disturbance. DALEC2 and many other ecosystem models assume that respiration rates are proportional to carbon pool size. It has been suggested that although this assumption works well in equilibrium conditions, it may not allow such models to predict ecosystem carbon dynamics following disturbance [Schimel and Weintraub, 2003]. The data assimilation techniques in this paper present a way for these simple models to cope with step changes in ecosystem behavior, by allowing parameters and carbon pools to be updated following disturbance events.

Conclusion
In this work we have investigated the response of a managed forest ecosystem to the disturbance of selective felling by using data assimilation. Assimilating all available data streams after an event of disturbance with a prior model prediction allows us to assess changes to model parameter and state variables due to this disturbance. We have also created modified observation operators to allow for the assimilation of daytime and nighttime NEE observations with a daily time step model. This negated the need for model modification and increased the number of observations by a factor of 4.25.
Our modeled estimates show no significant change in net ecosystem carbon uptake after a thinning event in 2014 where approximately 46% of trees were removed from the studied area. Similar results were also found following a thinning activity in 2007 . From our optimized model we find that reduced ecosystem respiration is the main reason for this unchanged net ecosystem carbon uptake. Therefore, even for a decrease in GPP following thinning, there is no significant change in NEE. We hypothesize that this reduction in ecosystem respiration is due to reduced input of autotrophic substrate following thinning, meaning that both autotrophic and heterotrophic respiration are reduced. These results support work suggesting that GPP is the dominant driver for ecosystem respiration [Janssens et al., 2001;Bhupinderpal et al., 2003;Högberg and Read, 2006;Heinemeyer et al., 2012;Moore et al., 2013]. This has implications for future predictions of land surface carbon uptake and whether forests will continue to sequester atmospheric CO 2 at similar rates or if they will be limited by increased GPP leading to increased respiration.