Impact of Nitrogen Fertilization on Forest Carbon Sequestration and Water Loss in a Chronosequence of Three Douglas-Fir Stands in the Pacific Northwest

To examine the effect of nitrogen (N) fertilization on forest carbon (C) sequestration and water loss, we used an artificial neural network model to estimate C fluxes and evapotranspiration (ET) in response to N fertilization during four post-fertilization years in a Pacific Northwest chronosequence of three Douglas-fir stands aged 61, 22 and 10 years old in 2010 (DF49, HDF88 and HDF00, respectively). Results showed that N fertilization increased gross primary productivity (GPP) for all three sites in all four years with the largest absolute increase at HDF00 followed by HDF88. Ecosystem respiration increased in all four years at HDF00, but decreased over the last three years at HDF88 and over all four years at DF49. As a result, fertilization increased the net ecosystem productivity of all three stands with the largest increase at HDF88, followed by DF49. Fertilization had no discernible effect on ET in any of the stands. Consequently, fertilization increased water use efficiency (WUE) OPEN ACCESS Forests 2015, 6 1898 in all four post-fertilization years at all three sites and also increased light use efficiency (LUE) of all the stands, especially HDF00. Our results suggest that the effects of fertilization on forest C sequestration and water loss may be associated with stand age and fertilization; the two younger stands appeared to be more efficient than the older stand with respect to GPP, WUE and LUE.


Introduction
The factors affecting the physiological processes controlling the amount of terrestrial carbon (C) sequestered in terrestrial ecosystems mainly include atmospheric CO2 concentration, climatic variability, land use change and nitrogen (N) fixation [1][2][3].In terrestrial ecosystems, the important processes of determining whether an ecosystem is a C sink or source are C fixation and release through photosynthesis and respiration, respectively.For N-limited terrestrial ecosystems, additional N supply can affect these processes, and in turn affect the strength of C sinks or sources [4][5][6][7][8][9].
Research has shown that the response of C exchange to N fertilization in terrestrial ecosystems is positive or negative depending on the ecosystem N status [10][11][12][13].Many researchers have studied the effects of N enrichment on net ecosystem productivity (NEP) in the forest ecosystem, but it remains uncertain as to how much N addition contributes to the magnitude of terrestrial NEP [1,14,15].
Many studies have found that N addition increases NEP [16][17][18][19][20][21].Magnani et al. [22] found NEP was highly correlated with N deposition in forest ecosystems and the response of NEP to N in wet deposition was approximately 725 kg C kg −1 N.However, Sutton et al. [23] found that the response of NEP to N was about 50-75 kg C kg −1 N based on a model re-analysis across 22 European forest sites and considering the impacts of climatological differences among stands.de Vries et al. [24] also found that the response of NEP to N would be 30-70 kg C kg −1 N based on a multi-factor analysis of European forest measurements at nearly 400 intensively monitored forest plots in Europe.The results of both Sutton et al. [23] and de Vries et al. [24] are much smaller than the estimates of Magnani et al. [22].Therefore, there is still some controversy over the magnitude and sustainability of NEP resulting from N enrichment and its potential mechanisms [25][26][27].
However, other studies have reported that N addition may cause a relatively small increase in NEP [12,28], have no effect [29], or have a negative effect by altering plant and microbial communities, including threatened and endangered species [30].Nadelhoffer et al. [31] pointed out that the influence of N enrichment on ecosystem C storage should be minor in forest ecosystems in the northern temperate zone.Currie et al. [32] stated that small increases in C storage occurring primarily in living and dead wood may result from elevated N deposition over the next few decades.Using meta-analysis from global N addition experiments, Liu and Greaver [33] concluded that the increases in short-term C sequestration below ground caused by N enrichment are due to increased C storage in the surface organic layer, but it appeared to be difficult to predict the response of C storage to N enrichment through a long-term experiment.Harpole et al. [34] reported that the ability of grassland ecosystems to sequester C late in the growing season was decreased by N enrichment because of increased growing season water use and earlier leaf senescence.
Furthermore, water and C cycles in terrestrial ecosystems are closely coupled.A global sensitivity analysis for an integrated forest ecosystem model [35] also showed the intrinsic coupling between water and the C cycle through stomatal conductance, which affect both water and C fluxes between the biosphere and atmosphere.Nutrients may not only affect productivity and foliar biomass but are also associated with evapotranspiration (ET) in forests and other ecosystems [36].Increased productivity and foliar biomass induced by increased nutrient availability are probably responsible for the increase in ET resulting from increased rainfall interception and transpiration.As a consequence, more attention should be paid to water use efficiency (WUE), defined as the ratio of C gain (usually gross primary productivity, GPP) to water loss, i.e., ET [37,38].Studies have reported that by increasing WUE, N addition could enhance plant productivity [39][40][41].On the other hand, other researchers have found that N enrichment has no impact on WUE [42][43][44] or a negative effect on WUE [45,46].Consequently, it is difficult to account for the differences in N-induced effects on WUE [47].Furthermore, fertilization may increase NEP through an increase in GPP, resulting from an increase in resource capture and resource use efficiency, e.g., increase in absorbed photosynthetically active radiation (APAR), light use efficiency (LUE, the ratio of GPP to APAR), and/or ecosystem carbon-use efficiency (CUEe, the ratio of NEP to GPP) [5].Although LUE models have been extensively used to estimate regional or global GPP through remote sensing techniques [48][49][50], the biophysical controls on LUE in terrestrial ecosystems, especially soil N, remain poorly understood, especially in Douglas-fir stands on the west coast of Canada.
The three different-aged Douglas-fir stands in this study (DF49, established in 1949, HDF88, established in 1988, and HDF00, established in 2000) are in the Pacific Northwest coastal forest region of Canada.Forests in this region between Oregon and Alaska cover approximately 10 5 km 2 and play an important role in the global C cycle.In this region, due to their considerable distances from heavy industrial activity, there is very little atmospheric N deposition, and N deficiency generally occurs in soils [51].Fertilization using a single dose of urea in early 2007 in all three of these Douglas-fir stands and on-going long-term measurements of C fluxes using the eddy-covariance (EC) technique offer an opportunity to examine the effects of N fertilization on C and water fluxes.Previous studies on the effects of N fertilization have focused on C fluxes at DF49 in the first post-fertilization year using a process-based model [52], and C fluxes, ET and WUE in all three stands in the two post-fertilization years using an empirical model (multiple linear regression, MLR) [53].
In this study, we have used the measured C fluxes (viz., GPP, ecosystem respiration (R), NEP (= GPP − R)), ET and climatic variables in the three stands from 1998 to 2010.The objectives of this study are threefold: (1) to estimate the effects of N fertilization in the three stands on the C fluxes, ET, WUE and LUE using two common data-driven models, MLR and an artificial neural network (ANN); (2) to examine the variation of the N fertilization effects during the 4-year post-fertilization period; and (3) to gain insight into the mechanisms of N fertilization controlling the C fluxes, ET, WUE and LUE in the three stands.

Site Descriptions
Three stands were located on the east coast of Vancouver Island, BC, between Campbell River and Denman Island.In 2010, the ages of the stands were 61, 22 and 10 years old, which spanned the typical ages of stands in this area.The three stands were dominated by Douglas-fir (Pseudotsuga menziesii (Mirb.)Franco) with relatively similar stand densities, soil, topography, elevation, and biogeoclimatic classification [54].On the whole, stand age and the corresponding stand structural characteristics were considered to be largely responsible for the differences among these stands in C uptake and ET.The oldest stand (49°52′7.8″N, 125°20′6.3″W, flux tower location), DF49, was planted with Douglas-fir seedlings in 1949 and occupied an area of 130 ha.The stand comprised 80% Douglas-fir, 17% western red cedar (Thuja plicata), and 3% western hemlock (Tsuga heterophylla).The mean annual temperature and precipitation measured during 1998-2010 in the stand are 8.1 °C and 1350 mm, respectively.The 22-year-old (pole-sapling) stand (49°32′10.49″N, 124°54′7.18″W), HDF88, 110 ha in size, was situated about 30 km southwest of DF49.The stand comprised 75% Douglas-fir, 21% western red cedar and 4% grand fir (Abies grandis).The mean annual temperature and rainfall in the stand, since measurements began in 2002, are 9.1 °C and 1579 mm, respectively.The youngest stand (49°52′1.08″N, 125°16′43.80″W), HDF00, was located about 2 km southeast of DF49.This 32-ha stand was harvested in the winter of 1999/2000 and in the following spring was planted with one-year-old seedlings (93% Douglas-fir, 7% western red cedar).Much of the leaf area at this site was due to the growth of pioneer and understory species from the previous stand.The mean annual temperature and rainfall at the site are 8.5 °C and 1321 mm, respectively.Additional details on stand history, vegetation and soil can be found in Humphreys et al. [54], Chen et al. [55], Krishnan et al. [56] and Morgenstern et al. [57].

Stand Fertilization
Urea fertilizer was spread aerially at 200 kg•N•ha −1 over about 110 ha at DF49 and 20 ha at HDF88, which included most of the respective tower flux footprints, on 13 January and 17 February 2007, respectively.Because of the young age of the planted trees and competing understory at HDF00, 80 g urea per tree (~60 kg•N•ha −1 ) was manually applied along the tree drip line of trees located in the 5-ha tower flux footprint on 13-14 February 2007.The fertilized areas around the respective towers accounted for more than 80% of the EC fluxes measured during unstable atmospheric conditions, which occurred during the daytime and nighttime [55].

EC and Climate Measurements
The EC technique has been used to measure continuous half-hourly CO2 and water vapor fluxes since 2000, 2001 and 1997 at DF49, HDF88 and HDF00, respectively.Climate variables that were also measured included precipitation, net radiation (Rn), photosynthetically active radiation (PAR), wind speed, air temperature (Ta) and humidity, soil temperature (Ts) and volumetric soil moisture content (θ).Details on EC instrument characteristics, measurements and calculation procedures for these stands are described in Humphreys et al. [54], Chen et al. [55] and Morgenstern et al. [57].Half-hourly EC-measured CO2 fluxes above the stand were corrected by adding the estimated rate of change in CO2 storage in the air column below the EC sensor height to obtain the net ecosystem exchange (NEE).NEP was obtained as NEP = −NEE [54].
The methodology of gap-filling for EC-measured NEP and its separation into GPP and R was adopted from the Fluxnet-Canada Research Network protocol [55,58].Missing nighttime values of R were obtained using an exponential relationship between NEE and Ts (Q10 model) using measured NEE values corresponding to friction velocity values greater than a threshold value.This model was also used to obtain daytime values of R. GPP during the growing season was obtained as the sum of daytime-measured NEP and daytime R and was set to zero during the cold season.An empirical light-response relationship was used to fill gaps in GPP.The details on NEP gap-filling and its separation into GPP and R as well as the uncertainties in the annual sums of NEP, GPP and R can be found in Chen et al. [59].

Artificial Neural Networks
An ANN model is characterized by flexible mathematical structures, which can be used to investigate complicated non-linear relationships between inputs and outputs [60].For training of time series prediction, the ANN model is conducted by propagating the input data and then back-propagating the error by means of self-adjustment of the weights using least squares residuals so that the simulated outputs best approximate the target outputs (observed data).The performance of an ANN model as a purely empirical non-linear regression model is commonly affected by the quality of the training dataset, network architecture and network training [61].Additional mathematical background on ANN models can be found in White [62], Jain et al. [63], Zhang et al. [64] and Basheer and Hajmeer [65].The back-propagation (BP) algorithm used in this study for network training has been successfully applied in ecosystem C flux estimation [66][67][68] because it is regarded as the most widely applied algorithm in the ANN literature [69,70].A three-layer BP neural network (BPNN) was employed with one hidden layer because Funahashi [71] and Cybenko [72] demonstrated that even with one hidden layer, a BPNN can approximate any continuous multivariate function with reasonable precision.In the calibration phase, the optimization approach applied was based on Levenberg-Marquardt [73], which is known to outperform the simple gradient descent and other conjugate gradient methods.In this three-layer BPNN, a hyperbolic tangent sigmoid transfer function (tansig) at hidden layer and a linear transfer function (purelin) at output layer was used.Aiming at obtaining the hidden node number and avoiding over-fitting during the training period, we used a trial and error method to select the optimal solution by altering the number of hidden nodes.

Comparison of ANN and MLR Models
Six environmental variables that included Ts at the 5-cm depth, Ta, downwelling PAR ( Q ), vapor pressure deficit (D), Rn, and θ in the 0-30 cm layer, and stand age were selected to train the ANN model using pre-fertilization flux measurements.The particular environmental variables that worked best for GPP, R and ET were chosen as inputs after trying several ANN networks.GPP, R and ET were simulated independently using different input variables for the ANN model.Specifically, three different variable groups, (Ta, θ, Q , D, stand age), (Ta, θ, Ts, stand age) and (Ta, θ, Rn, D, stand age), were then used to train the ANN model and predict monthly GPP, R and ET, respectively, for the three sites assuming the stands had not been fertilized.To ensure high precision in the period of model prediction, we used multi-year monthly values of climate variables, stand age and EC-measured C fluxes and ET before 2005 to train the ANN model and then verified the trained model with measurements in 2005 and 2006.When we were convinced that the optimized ANN model successfully simulated the multiyear seasonal variations in C fluxes and ET, the input values for the post-fertilization period were transferred into the trained model to predict the GPP, R and ET for 2007 to 2010.The resulting differences between the measurements and predictions were used to discern the impact of fertilization.
At the beginning of the training, inputs and outputs were normalized between 0 and 1, which is a common preprocessing method for variables with different values.Moreover, we also used an MLR model to simulate C fluxes and ET with the same variables to assess the advantages of the ANN modeling approach over traditional regression modeling techniques.Similar to the ANN model, we used a variable group (e.g., Ta, θ, Q, D, stand age) as the independent variables to establish an MLR equation and predict the dependent variable (e.g., GPP).Because of the differences in age among these three stands, it is important to note that both ANN and MLR models were run separately for each stand.The performances of the models were evaluated by comparing predicted with measured values in the calibration and verification periods.In this study, computer programming and data analysis were performed using MATLAB (version 7.13, R2011b).

Environmental Variables
Since all three stands were exposed to similar weather conditions (warm dry summers and wet cool winters), we describe only the seasonal and annual variations in environmental variables at DF49.The seasonal patterns of Ta, Ts, θ, D, Q and Rn during 2007-2010 were similar to those measured before fertilization from 1998 to 2006 (Figure 1).Compared with the 9-year means during the pre-fertilization period, the annual average values of Ta during 2007-2010 were lower by about 0.7, 1.1, 0.4 and 0.8 °C, respectively, and the annual average values of Ts were higher by about 0.3, 0.1, 0.1 and 0.6 °C, respectively (Figure 1a,b).Annual values of total Q in 2007, 2008 and 2010 were lower by about 9%, 3% and 6%, respectively, suggesting possible growth limitation in those years, but were higher by about 8% in 2009 (Figure 1f).Annual average values of θ in the 0-30 cm soil layer in the growing season (May-September) in the post-fertilization period during 2007-2010 were 0.22, 0.20, 0.19 and 0.20 m 3 •m −3 , respectively, which were higher, especially in the noticeably wet year of 2007, compared with the pre-fertilization 9-year mean (0.16 m 3 m −3 ) (Figure 1c).

Seasonal and Interannual Variations of Observed C Fluxes and ET before and after Fertilization
The seasonal variations in the monthly EC-measured C fluxes and ET in the three stands before and after fertilization are presented in Figure 2. Monthly GPP values during the 4 post-fertilization years were similar to the means during the 9 non-fertilized years (Figure 2a).The monthly values of R during the 4 post-fertilization years were close to the 9-year pre-fertilization means for Jan-May but lower for the later months except for June and July in 2009 (Figure 2b).Therefore, the monthly NEP values during the 4 post-fertilization years were greater than the 9-year pre-fertilization average values for March-December, especially in June, August and September (Figure 2c).Similarly, for HDF88, monthly NEP values during the 4 post-fertilization years at HDF88 were generally greater than the pre-fertilization average values for 2002 to 2006 (Figure 2g).In the case of HDF00, compared with the 2001-2006 means, the largest increases in monthly values of NEP during 2007-2010 occurred in June, July and August (Figure 2k).The monthly GPP values during 2007-2010 at HDF88 and HDF00 were higher than the mean values for the pre-fertilization periods (Figure 2e,i).The monthly R values for HDF00 were also generally higher than the monthly means in the pre-fertilization period (Figure 2j).Furthermore, at DF49, the maximum NEP occurred in April and May (Figure 2c), while maximum NEP at HDF88 and HDF00 occurred in June and July (Figure 2g,k).Most of the monthly ET values during 2007-2010 at DF49 and HDF88 were close to the monthly means in the pre-fertilization period (Figure 2d,h).However, at HDF00, monthly ET values for May-July during 2009 to 2010 were significantly higher than the means for the pre-fertilization period (Figure 2l).Moreover, for DF49, the monthly values of ET during January to April in 2010 were higher than the previous 9-year average values (Figure 2d).
To better examine the interannual variations of EC-measured C fluxes, ET, WUE and LUE, the annual measured values are compared in Figure 3.The results indicate that EC-measured NEP during the post-fertilization four years increased in all three stands (Figure 3c).Annual NEP at DF49 remained stable over post-fertilization years with values about 200 g C m −2 year −1 higher than in the pre-fertilization period, despite the post-fertilization years being cooler and wetter than the pre-fertilization years (Figure 1a,c).However, at HDF88, annual NEP consistently increased with increasing stand age (Figure 3c).
The interannual variations in GPP and R at DF49 since 1998 show that R decreased while GPP showed little change during 2007-2010 (Figure 3a,b).However, for the other two stands, annual R increased in the first post-fertilization year and decreased in the next three years (Figure 3b).In addition, compared with the pre-fertilization (9-years) mean CUEe (17%) in DF49, annual average values of CUEe during 2007-2010 were higher by about 9%, 8%, 9% and 12%, respectively.
Annual ET during the post-fertilization period did not show an apparent difference from the previous 9-year means in all three stands (Figure 3d).Annual values of WUE varied greatly among the three stands (Figure 3e).Annual WUE of DF49 was the highest, ranging from 4.5 to 5.

Verification of the ANN Model and Its Comparison to the MLR Model
The differences between EC-measured, ANN-and MLR-modeled monthly GPP, R, and NEP in the three stands are compared in Figure 4   MLR simulations for DF49 for the model calibration period, however, explained about 97%, 92% and 76% of the variance of monthly GPP, R, and NEP, respectively (Figure 4a and Table 1), which are somewhat lower than for the ANN simulations.The linear regression analysis showed RMSE = 21, 29, and 19 g•C•m −2 •month −1 for GPP, R and NEP, respectively, which are significantly higher than those for the comparison of ANN-simulated and EC-measured values.In addition, for the model verification period, linear regression between MLR-simulated and observed values showed RMSE = 30, 34 and 24 g•C•m −2 month −1 for GPP, R, and NEP, respectively, which are all significantly higher than those for the calibration period.
For the other two stands, the performance of the ANN model was also superior to that of the MLR model in the estimates of GPP, R and NEP (Table 1).Regarding ET, similar results were observed in the calibration and verification periods in all three stands (Table 1).As a consequence, we focused on describing and discussing the results of ANN modeling in the prediction period to ascertain the responses of C fluxes and ET to fertilization.

Effects of N addition on C exchange Fluxes, ET, WUE and LUE in a Chronosequence
In order to provide insights into the influence of fertilization on the C fluxes, ET as well as WUE and LUE during the period of prediction, monthly measured and ANN-modeled values, and their differences are shown for the four post-fertilization years (2007-2010) in Figures 5 and 6.The annual effects of N fertilization during this period are summarized in Tables 2 and 3. Figure 7 compares the annual measured and ANN modeled values of the C fluxes and ET, and shows ± 1 standard deviation estimated using the optimized model parameters.
Over the four post-fertilization years at DF49, most of the ANN modeled monthly values of GPP were lower than the measured values (Table 1 and Figure 5a), but were higher for R (Table 1 and Figure 5d) and consequently, modeled monthly values of NEP were lower than the observed values (Table 1 and Figure 5g).
At the annual time scale over the four years, N fertilization at DF49 caused a 6% to 11% increase in annual GPP and a 6% to 10% decrease in R, resulting in an 84% to 164% increase in NEP (Table 2).For HDF88, N fertilization caused an 18% to 29% increase in annual GPP, an 11% increase in R in 2007, and a 6% to 14% decrease during 2008-2010, resulting in a 208% to 321% increase in NEP (Table 2).For HDF00, N fertilization caused a 43% to 70% increase in annual GPP and an 8% to 21% increase in R, resulting in a 27% to 50% increase in NEP (Table 2).Overall, N fertilization had a positive impact on GPP for all three stands with the response being the highest at HDF00, followed by HDF88 (Table 2).
No consistent effects on annual ET were evident in the four post-fertilization years in all three stands (Table 3).For annual WUE, N fertilization of DF49 caused an increase of 0.30 to 0.  3).For annual LUE, N fertilization had positive impacts for all three stands, especially at HDF00 with the greatest response ranging from 0.0039 to 0.0057 mol•C•mol −1 photons (Table 3).

Table 1. Comparisons of artificial neural network (ANN) and multiple linear regression (MLR) model performances among the calibration, verification and prediction periods for monthly gross primary productivity (GPP
and evapotranspiration (ET, mm•year −1 ) in the three Douglas-fir stands.Table 2. Effects of nitrogen fertilization on annual gross primary productivity (GPP, g and net ecosystem productivity (NEP, g•C•m −2 •year −1 ) for the four post-fertilization years 2007-2010 in the three Douglas-fir stands.

Effects of N Fertilization on Gross Primary Productivity
N fertilization generally results in GPP and net primary productivity (NPP) increases in terrestrial ecosystems by promoting growth and biomass accumulation in terrestrial plants [7,21,25,74,75].Our results showed that N fertilization generally increased GPP in all three stands, especially in the younger stands, HDF00 and HDF88.N addition resulted in a GPP increase of 8% at DF49 in the first post-fertilization year, which is in excellent agreement with an 8% increase obtained by Chen et al. [52] using a model-data synthesis approach and a 10% increase calculated by Jassal et al. [53] using a simple empirical model (Table 4).This is lower than the 14% increase estimated by Grant et al. [76] using a process-based model.In addition, our results, which showed that N fertilization led to an increase in GPP in the three stands during the first two years after fertilization, are consistent in sign with the findings by Jassal et al. [53], although the magnitude of the increases in GPP are somewhat different.This difference may partly result from the reported uncertainty due to the use of the empirical relationships used for gap-filling, resulting in an uncertainty of ±50 g•C•m −2 •year −1 in EC-measured GPP, and ±75 and ±25 g•C•m −2 •year −1 for measured R and NEP, respectively [54].Overall, N fertilization had a positive effect on GPP and the effectiveness of fertilization varied with stand age.

Effects of N Fertilization on Ecosystem Respiration
Our results indicated that N fertilization of DF49 reduced R by approximately 101 g•C•m −2 •year −1 in the first post-fertilization year, which is in good agreement with the 93 g•C•m −2 •year −1 decrease calculated by Chen, et al. [52] (Table 4).During the next three years (2008-2010), modeled R remained greater than measured R by 114 to 170 g•C•m −2 •year −1 , implying that N fertilization suppressed R.However, the increase in R induced by N fertilization at HDF00 during the first two years after fertilization agrees with the findings reported in Jassal et al. [53].Our ANN modeling indicated that this positive effect on R at HDF00 continued for the third and fourth post-fertilization years (2009 to 2010) (Table 2).It should be noted that Jassal et al [53] partitioned NEE by estimating daytime R using the relationship between daytime NEP and PAR rather than by using the relationship between nighttime NEE and Ts.Such variable effects of N addition on R in different stands suggest the need to investigate the mechanisms.

Effects of N fertilization on Net Ecosystem Production
N fertilization significantly increased NEP during the four post-fertilization years in all the three stands, especially HDF88, followed by DF49, largely due to N-induced GPP increases and R decreases.These results are consistent with those reported by Chen, et al. [52] in the first post-fertilization year, and are also consistent in sign with the results obtained by Jassal et al. [53] both at DF49 and HDF88 in the first two post-fertilization years (Table 4), although the magnitude of the increase in NEP was less than that in our study.However, at HDF00, a substantial increase in NEP in the four post-fertilization years was primarily due to GPP increasing more than R (Table 2).
The C:N response [kg C (sequestered) kg −1 N added] has been widely used to estimate the response of C sequestration to N addition [22][23][24]77].Our results suggest that C:N response varied with stand age.In the first post-fertilization year, DF49 and HDF88 had approximately the same C:N response (13 kg C kg −1 N, followed by HDF00 (8 kg C kg −1 N).Our results are consistent with those of Jassal et al. [53], de Vries et al. [24], Högberg [26] and Sutton et al. [23], which are significantly smaller than the assessment made by Magnani et al. [22], with an extremely high C:N response (approximately 725 kg C kg −1 N) to wet deposition.Furthermore, it seems that the 60 kg•N•ha −1 application to individual trees at HDF00 was less efficient with respect to NEP than the larger application (200 kg N ha −1 ) at DF49 and HDF88.

Effects of N fertilization on ET, WUE and LUE
N application showed no consistent effect on ET among the three stands.N fertilization slightly decreased ET at HDF88 in the last three years and had no clear effect on ET at HDF00 over the four post-fertilization years (Table 3).On the other hand, increased nutrient availability (e.g., N fertilization) had significantly positive effects on WUE at HDF88 and HDF00 due to its strong positive effect on GPP compared with its small effect on ET, but had no significant effect at DF49.Ripullone et al. [47] found that N addition increased both WUE and biomass production.They concluded that this was mainly due to the positive response of photosynthesis to N availability with it having no impact on either stomatal conductance or transpiration.Consequently, improved tree nutrition (e.g., N fertilization) could cause an increase in forest productivity [78] without leading to an increase in demand for water.
Furthermore, N fertilization caused an increase in the 4-year mean WUE of 0.34 g•C (kg•water) −1 for DF49, implying a relatively conservative response to N fertilization, and increases of 0.78 and 1.19 g•C (kg•water) −1 for HDF88 and HDF00, respectively, likely due to their younger stand ages and less developed root system [56].Wu et al. [79] suggested that low or appropriate N application might improve the survivability of Sophora davidii seedlings under water deficit conditions, resulting in an increase in WUE and stimulating growth and biomass production.Our results indicated that the supply of 60 kg N ha −1 applied to individual trees at HDF00 may be more efficient in terms of WUE compared with the higher application of 200 kg N ha −1 at HDF88 and DF49, especially in the latter stand.
N fertilization in our study had positive impacts on LUE in all three stands in the four post-fertilization years, especially at HDF00 followed by HDF88 with the smallest increase occurring in 2009 in all three stands accompanied by the highest annual PAR.This suggests that there was a reduction in the increase in LUE induced by N fertilization under high PAR, which has also been reported by Ibrom et al. [80] and Schwalm et al. [81].Furthermore, N fertilization caused a 0.002•mol C•mol −1 photon increase in the 4-year mean LUE over the post-fertilization period for DF49, implying a relatively conservative plant productivity response to N fertilization, while higher increases of 0.0032 and 0.0047 mol•C•mol −1 photons for HDF88 and HDF00, respectively, were likely due to their younger stand age and rapid stand development.Overall, stand development had a significant impact on plant productivity, LUE and WUE as well as nutrient demand (e.g., N fertilization) [82].

Modeling Uncertainty and Limitations
Several studies have found that there is observable interannual variability of C component fluxes and ET in these stands, which to a large extent, results from interannual variability in climate [42,54,56,59].We used both ANN and MLR approaches for removing the impacts of climatic variability on C fluxes and ET in the three stands during the post-fertilization period.In addition, a significant effect of stand development on C fluxes in the two younger stands (HDF88 and HDF00) was demonstrated in Figure 3.Our results show that pre-fertilization monthly C fluxes and ET affected by climate and stand age were better described using the ANN approach than the MLR approach.
The determination of BPNN model parameters including the hidden layer node number, momentum coefficient and learning rate and its shortcomings, such as the slow rate of convergence and local minimum, have considerably restricted its application and have led to uncertainties in the estimate of the N effect in this study.Nevertheless, we believe this study contributes significantly to quantifying the effects of N fertilization on C fluxes and ET and reducing the uncertainty in model simulation.Considering the limitations of the BPNN model, in order to further improve the estimate of the N effect, more importance should be attached to a hybrid BPNN model with optimized algorithms such as the genetic algorithm and particle swarm optimization.

Conclusions
1. Pre-fertilization monthly C fluxes and ET were better described using the ANN approach than the MLR approach.
2. N fertilization increased GPP in all three Douglas-fir stands over the four post-fertilization years with the greatest increase in the youngest stand (HDF00) followed by HDF88, suggesting that the effectiveness of fertilization on GPP may be associated with stand age.In addition, fertilization decreased R over the four years after fertilization at DF49, and over the last three years at HDF88, but increased R over the four years at HDF00.As a result, N fertilization increased NEP in all three stands, with the greatest increase in HDF88 and followed by DF49.
3. N fertilization had no discernible effects on ET among the three stands.This result led to a significant increase in WUE in the two younger stands (HDF88 and HDF00) but a small increase in only the first three years for the oldest stand (DF49).
4. N fertilization increased LUE in all three stands, especially the youngest stand (HDF00).

Figure 2 .
Figure 2. Comparisons of mean monthly values of eddy-covariance tower measured gross primary productivity (GPP), ecosystem respiration (R), net ecosystem productivity (NEP) and evapotranspiration (ET) in all three stands for the pre-fertilization and post-fertilization years (2007-2010).The means of measured carbon component fluxes and ET at DF49, HDF88 and HDF00 are shown for 1998 to 2006, 2002 to 2006 and 2001 to 2006, respectively.(a-d) for DF49; (e-h) for HDF88; and (i-l) for HDF00.The error bars are ± 1 SD for all the pre-fertilization data.
and Tables 1-3.ANN simulations for the model calibration period (1998-2004) at DF49 showed strong agreement with the observed values and explained about 99%, 99% and 96% of the variance in monthly GPP, R, and NEP, respectively (Root mean square error (RMSE) = 2, 8, and 8 g•C•m −2 •month −1 , respectively).These results indicate that by optimizing the ANN parameters, the trained ANN model achieved reasonable simulation of monthly GPP, R and NEP at DF49.The linear regression analysis between ANN-modeled and measured values for the model verification period (2005-2006) for DF49 showed R 2 = 0.99, 0.99 and 0.98 and RMSE = 3, 6 and 6 g•C•m −2 •month −1 for GPP, R, and NEP, respectively (Table1).

Figure 5 .
Figure 5. Artificial neural network (ANN) simulated and eddy-covariance tower measured monthly gross primary productivity (GPP), ecosystem respiration (R) and net ecosystem productivity (NEP) in all three stands for the four post-fertilization years 2007-2010.(a-c) for GPP; (d-f) for R; and (g-i) for NEP.The effects of N fertilization on the carbon (C) component fluxes were estimated as the differences between the measured C fluxes and their corresponding modeled values.

Figure 6 .
Figure 6.Artificial neural network (ANN) simulated and eddy-covariance tower measured monthly water-use efficiency (WUE) and light-use efficiency (LUE) in all three stands for the four post-fertilization years 2007-2010.(a-c) for WUE and (d-f) for LUE.The effects of N fertilization on WUE and LUE were estimated as the differences between the measured WUE and LUE values and their corresponding modeled values.

Figure 7 .
Figure 7. Comparisons of measured and artificial neural network (ANN) simulated annual gross primary productivity (GPP), ecosystem respiration (R), net ecosystem productivity (NEP) and evapotranspiration (ET) at DF49 from 1998 to 2010.(a) for the ANN calibration years 1998-2004; (b) for the model validation years 2005-2006; and (c) for the four post-fertilization years 2007-2010.The error bars with ± 1 SD are the standard errors for ANN-modeled annual values estimated using the optimized model parameters.