Suitability of the boreal ecosystem simulator (BEPS) model for estimating gross primary productivity in hemi-boreal upland pine forest

The Soontaga flux tower (58°01′24 ″N, 26°04′15 ″E) is located in a dry hemi-boreal forest dominated by Scots pine (Pinus sylvestris L.). Soontaga area is a coniferous upland forest site of Vaccinium type (Lõhmus, 2004) where the stand age, where surrounding stands fall into this range, was approximately 60–210 years with the maximum canopy height of 30 m and the existing second layer of Picea abies (L.) H. Karst. with an average height of 15 m. Mean annual temperature was 7.2 °C, with mean annual sum of precipitation of about 728 mm in 2019 (Estonian Environment Agency, 2020). The location of the study area and tower is shown in Figure 1.

Figure 1

Following Román et al. (2009) and Wang et al. (2012, 2014), surface heterogeneity was evaluated using the surface albedo retrievals from 30 m spatial resolution Landsat data as an input to the model. The surface albedo was obtained with the narrow-band-to-broadband conversion (Liang et al., 2003; Smith, 2010): (1) ∝=0.356p1+0.130p3+0.373p4+0.085p5+0.072p7−0.00180.356+0.130+0.373+0.085+0.072, \propto = {{0.356{p_1} + 0.130{p_3} + 0.373{p_4} + 0.085{p_5} + 0.072{p_7} - 0.0018} \over {0.356 + 0.130 + 0.373 + 0.085 + 0.072}}, where are respective Landsat-7 ETM+ sensor bands. The spatial representativeness at different scales was assessed by calculating the gamma variance with 500 m, 1 km, 1.5 km, and 2 km footprint areas centered at the flux tower (Figure 2A).

Figure 2

The gamma variance estimator, was used to attain half of the average square difference between albedo values, which are within a specific distance, classes, or bins, and these are defined by the multiplication of 30 m (Román et al., 2009): (2) yE(h)=0.5∑i=1N(h)(Zxi−Zxi)2N(h). {y_E}\left( h \right) = 0.5{{\sum\nolimits_{i = 1}^{N\left( h \right)} {{{\left( {{Z_{xi}} - {Z_{xi}}} \right)}^2}} } \over {N\left( h \right)}}. Here, refers to the surface albedo at the pixel location of ; and corresponds to the surface albedo at another pixel which should be within a lag distance h. N(h) is the number of paired data at a distance of h. The maximum lag distance which is used in each gamma variance is constrained by the half maximum distance of the given subset and as a rule of thumb it has to be a factor of the minimum lag (i.e., 30 m) (Román et al., 2009). Thus, for a 1.0 km² subset = 690 m, for a 1.5 km² subset = 1050 m, and for a 2.0 km² subset = 1410 m.

The BEPS model consists of an advanced canopy radiation sub-model for quantifying the effects of canopy architecture on the radiation distribution and photosynthesis in the canopy (Feng et al., 2007). BEPS computes the total photosynthetically active radiation absorbed by the canopy (Liu et al., 1999). The sunlit and shaded leaf separation approach strategy (Chen et al., 1999) has been used in the BEPS model (Liu et al., 2002). The stratification strategy is preferred in the BEPS model over other strategies (Bonan, 1995) because it captures the key radiation variation inside the canopy which allows effective temporal integration (Chen et al., 1999) during a given time.

The total GPP can be calculated with the separation of sunlit and shaded leaf groups (Liu et al., 2002): (3) GPP=(AsunLAIsun+AshadeLAIshade). GPP = \left( {{A_{sun}}LA{I_{sun}} + {A_{shade}}LA{I_{shade}}} \right). In Equation (3), the subscripts ‘sun’ and ‘shade’ denote sunlit and shaded components, A is the photosynthesis rate, and LAI is leaf area index.

Photosynthesis rates are calculated as: (4) Ac,i=VmCi−ΓCi+Kc(1+O/K0), {A_{c,i}} = {V_m}{{{C_i} - \Gamma } \over {{C_i} + {K_c}\left( {1 + O/{K_0}} \right)}}, (5) Aj,i=JCi−Γ4(Ci+2Γ), {A_{j,i}} = J{{{C_i} - \Gamma } \over {4\left( {{C_i} + 2\Gamma } \right)}}, where and represent the Rubisco-limited and light-limited gross photosynthesis rates, represents the maximum carboxylation rate, is the radiation-dependent electron transport rate, and represent the intercellular CO₂ and the oxygen concentration in the atmosphere, respectively, stands for the CO₂ compensation point without dark respiration, and represent the Michaelis-Menten constants for CO₂ and O₂, respectively, and calculations are done separately for sunlit leaves (i = 1) and shaded leaves (i = 2).

The photosynthesis rate is calculated as the minimum of Rubisco-limited and light-limited photosynthesis rates: (6) Ai=min(Ac,i,Aj,i). {A_i} = \min \left( {{A_{c,i}},{A_{j,i}}} \right). The leaf area of sunlit (LAI_sun) and shaded (LAI_shade) leaves is estimated as follows: (7) LAIsun=2 cos θ(1−exp(−0.5ΩLAI/cosθ), LA{I_{sun}} = 2\,\cos \,\theta \left( {1 - \exp \left( { - 0.5\Omega LAI/\cos \theta } \right),} \right. (8) LAIshade=LAI−LAIsun. LA{I_{shade}} = LAI - LA{I_{sun}}. θ represents the solar zenith angle, represents the clumping index. The clumping index describes the aggregation of the foliage elements within a canopy. The random distribution would equal 1; a clumping index with values less than one indicates more clumped vegetation (Nilson, 1971). BEPS, its daily parameterization, and structure are described in detail elsewhere (Ju et al., 2006; Luo et al., 2018).

The input data required by BEPS are LAI, land cover, clumping index, and daily meteorological data (Feng et al., 2007).

MODIS product MCD15A3H.006 with a 500 m resolution and 4-day step was used to provide LAI input data. The original LAI values were interpolated linearly to daily steps. The MODIS LAI values occasionally reached the maximum value of 6, which is not realistic for the Soontaga tower area. Since LAI value of 3 was previously found to be a maximum value at a comparable Scots pine stand in Järvselja, Estonia (Pisek et al., 2011), the input LAI values were subsequently capped at a value of 3 to provide a more realistic scenario encountered at the Soontaga site. Evergreen needleleaf land cover type class and the value of 0.69 as the clumping index provided for the location in the global map by, He et al. (2012) were used in this study.

The necessary meteorological data for the BEPS model include temperature, humidity, radiation, precipitation, and wind speed information (Liu et al., 1997). All the meteorological data used in this study were measured directly at the site. The BEPS model was run with an hourly step. Any missing values in meteorological data were interpolated and converted to the hourly values from half-hourly values.

CO₂ concentration and wind parameters were continuously measured in the tower, at the height of 39 m (5 m above the canopy) in 10 Hz using an enclosed-path infrared CO₂ gas analyzer (Li-7200, LI-COR Biosciences, Lincoln, USA) and 3D-sonic anemometer (uSonic Class A, METEK GmbH, Elmshorn, Germany). Also, environmental parameters, such as air temperature and relative humidity (HC2A-S3, Rotronic AG, Bassersdorf, Switzerland) and solar radiation (LI-190SL, LI-COR Biosciences, Lincoln, USA) were measured at the site at a sampling frequency of 1 Hz. All the readings were stored on a data logger (CR3000, Campbell Scientific, Inc., Logan, UT, USA). Tower GPP values were obtained using the nighttime-data-based flux partitioning method in REddyProcWeb online tool (Wutzler et al., 2018; Reichstein et al., 2005). Nighttime (solar radiation less than 10 W) net ecosystem exchange (NEE) values are considered to be fully represented by nighttime ecosystem respiration (ER) and GPP is set to zero. Daytime GPP is then obtained as a difference between measured daytime NEE and modeled daytime ecosystem respiration (ER). Lloyd & Taylor (1994) regression model was fitted to the nighttime ER data using air temperature as the driving factor to fill the missing values. (9) ER=ERref×eE0×(1Tref−T0−1T−T0), ER = {\rm{E}}{{\rm{R}}_{{\rm{ref}}}} \times {e^{{E_0} \times \left( {{1 \over {{T_{ref}} - {T_0}}} - {1 \over {T - {T_0}}}} \right)}}, where R_ref (μmol C m⁻² s⁻¹) is the respiration at the reference temperature; E₀ (kJ mol⁻¹) is the activation energy; T (°C) is the measured air temperature. T_ref was set to 15°C, and T₀ was kept constant at −46.02 °C as in Lloyd & Taylor (1994).

Model parameters estimated from the nighttime data with a sliding 7-day window and daytime air temperature were then used to model daytime ER (Lloyd & Taylor, 1994). The tower footprint area was calculated using Kljun et al. (2015) model.

Figure 2B presents the geostatistical assessment of the site heterogeneity assessed at four different scales (0.5 km, 1 km, 1.5 km, 2 km). The gamma variance for the 0.5 km footprint area (the nominal spatial resolution of the satellite input data) around the tower reaches a plateau at 0.0003. According to Wang et al. (2017), if the gamma variance levels off at a value less than 0.0005, the site can be deemed spatially representative. The Soontaga area within the 500 m distance from the tower can be considered spatially homogeneous at the tower footprint (20.88 ha) (estimated using the Kljun model (Kljun et al., 2015)) as well as the nominal resolution of the used remote sensing data. However, the gamma variance values for the distances of 1 km, 1.5 km and 2 km level off at values > 0.0005 (Figure 3). Based on the gamma variance values, the Soontaga site might not be representative for a greater area (≥ 1 km) as the gamma variance is over the suggested homogeneity threshold by Wang et al. (2017). It shall also be noted that the MODIS 500 m resolution products are generated from multi-angular observations that may come from a larger area (Wang et al., 2017) and could be thus influenced by the heterogeneity at the greater area around Soontaga tower.

Figure 3

BEPS model validation and comparison The flux tower GPP_tower measurements were used to validate the BEPS result for each year. Table 2 provides the general overview of the GPP_BEPS values obtained with the BEPS model and the eddy covariance data measured at the site. A very close match was observed in 2017 and 2019. There was an excellent agreement (only 0.33% difference) in the total GPP_BEPS obtained with BEPS (1293) and tower measurements (1288) in 2019 with the root mean square error (RMSE) of 1.47. Similarly, only a 3.96% difference was observed in 2017, with GPP estimated by BEPS at 1214 and tower GPP_tower at 1166 and RMSE of 1.51. Compared to the tower measurements, the BEPS model overestimated GPP by 14.96% in 2018 with RMSE of 2.25 and underestimated it by 23.12% with RMSE of 2.20 in 2016.

First, as mentioned in the Methods section, there may be a possible footprint mismatch between the tower eddy covariance measurements and scale resolution of BEPS input data from remote sensing. The footprint area of the Soontaga flux tower is 20.88 ha. The remote sensing data that are used to drive the BEPS model (LAI, clumping index) are provided at nominal 500 m resolution. However, the actual signal may come from a bigger area, and may also not always align with the footprint of the flux tower. The GPP_tower measurements also have uncertainty introduced during calculation steps (Liu et al., 2016). The nighttime-data-based GPP_tower measurement method employed in the eddy covariance technique relies heavily on the fluxes estimated during the periods with a limited amount of quality data due to the low turbulence conditions (Reichstein et al., 2005).

The meteorological weather condition of the specific year can impact the BEPS model’s performance. Table 3 provides an overview of warm/cold and dry/wet conditions at the site during the studied period.

The expected GPP_BEPS rate was higher for 2016 as it had cold, and wet weather conditions based on the meteorological data presented in Table 3. The results were less convincing compared to other years with the 23.12% underestimation of the annual sum of GPP_BEPS estimates by BEPS compared to tower-based estimates (Figure 3A). While there was a close relationship with a very little difference between BEPS and tower estimates in the second half of the growing season (DOY > 210; Figure 3A), there was an unrealistic drop in BEPS GPP_BEPS retrievals from the middle of June till the third week of July (DOY 180–205). Here the used meteorological data was from the field data and that had some missing values in humidity calculation and those missing humidity values caused this drop in GPP_BEPS values during that period. A more precise and accurate dataset can improve the result. Additionally, clearer evidence can be seen if the daily GPP_BEPS values are plotted against the tower GPP_tower values with the dropped humidity input data and without the dropped humidity input data. All days of the year included in Figure 4A show a diverging relationship from the 1:1 line and the value is 0.7, while the omission of days with missing humidity inputs improves the value to 0.81 (Figure 4A). On the other hand, the GPP_BEPS was over-estimated with very high values from the last week of May to the second week of June. It was caused by the apparently over-estimated MODIS LAI values at the beginning of the growing season when the LAI compared to the probable actual LAI values at the site. From the beginning of the year until mid-February, the LAI values were zero (Figure 5A), resulting in zero GPP_BEPS estimates from BEPS during that period. Except for the LAI impact and the humidity input issue, the second half of the year showed BEPS GPP_BEPS values to closely match the expected tower-based observations. GPP_BEPS course also mirrors low temperature and high precipitation (Figure 5A) episodes well during this period.

Figure 4

Figure 5

Based on the meteorological data, 2017 was relatively cold and wet (Table 3), with a much higher precipitation rate in autumn compared to other years used in this study (Estonian Weather Service, 2021). Again, there was a close relationship with very little difference between GPP_BEPS and GPP_tower values throughout the whole year in 2017 (Figure 3B). Compared to tower-based estimates, the BEPS model overestimated the yearly GPP_BEPS sum by 4%. There were occasional small differences from the last week of May to the third week of August (Figure 3B). The GPP_BEPS trajectory by BEPS closely followed the seasonal temperature profile (Figure 5C). Also, higher precipitation clearly affected higher GPP_BEPS values from both approaches (Figure 5C). Similarly to the situations in 2018 and 2019, the occasional small differences in BEPS GPP_BEPS from the last week of May to the third week of August were caused by the apparent LAI overestimation by the MODIS LAI product during this part of the growing season. The differences in GPP_BEPS (Figure 3B) match the LAI saturation timing (Figure 5B).

The year 2018 experienced drought during July in Soontaga area (Estonian Weather Service, 2021) and the year was facing higher temperatures than usual (Table 3). According to the tower-based measurements, this caused relatively lower GPP_tower production as there is less water to support the process because dry weather and high temperature caused high evaporation. The GPP_BEPS matched GPP_tower values very well at the beginning and end of the growing season in 2018 (Figure 3C). The BEPS model provided much higher GPP_BEPS estimates than the tower GPP_tower estimates for most of the growing season. The GPP_BEPS course matched the radiation and temperature profiles very well (Figure 5C). Figure 5C indicates that low precipitation or dryness is behind the low GPP_BEPS values. Similarly, in 2019 later on, the over-estimated LAI from the MODIS LAI product caused the over-estimation in the result of GPP_BEPS from May to the third week of August (DOY 150–230). The first peak of BEPS overestimation started from 7 May (DOY 127), which matched the sudden increase in MODIS LAI values beyond the threshold LAI value of 3, which lasted with few exceptions till 17 August (DOY 229). Results from the 2018 season confirm clear, very high sensitivity BEPS to LAI input.

Based on meteorological data in Table 3, 2019 was a relatively warm and wet year at this site compared to the other years included in this analysis. The site had enough water to support the photosynthesis, but at the same time, high temperature caused a high evaporation rate. The two effects combined had a neutral effect and the GPP production was as expected. There was a good match between BEPS and tower estimates throughout most of the seasonal course in 2019 (Figure 3D). Compared to GPP_tower values obtained with tower eddy covariance measurements, BEPS underestimated GPP_BEPS from January to April and November to December. Figure 5A shows that the LAI input values for BEPS are mostly zero during the first few months (till the middle of April) as well as the last two months of the year. Since this is an evergreen needle-leaf site, the actual LAI values were higher and allowed the photosynthesis process to start immediately under suitable conditions, captured by the eddy covariance measurements. In contrast, close to zero input LAI values provided to BEPS did not allow to match the observed tower GPP_tower values during these periods.

The GPP_BEPS values were higher than the tower estimates from the middle of May to the middle of June 2019 (DOY 130–151) (Figure 3A). Compared to the beginning and towards the end of the year discussed above, input LAI values for BEPS did not suffer from underestimation during this period. The sharp increase in GPP_BEPS values around 10 May (DOY 130) coincided with the moment when the maximum LAI value, as provided by the MODIS LAI product, was reached. The period from the middle of May till the end of June may point to the opposite effect compared to the start of the season. While the MODIS LAI product underestimated LAI over the site at the beginning of the season and consequently underestimated the GPP, the overestimated GPP_BEPS from the middle of May till the beginning of July was caused by the apparent LAI overestimation in the MODIS LAI product during this section of the growing season. Previously, the maximum LAI for Scots pine stands was observed to be reached later in the season at the beginning of July (Heiskanen et al., 2012), compared to the earlier dates observed by the MODIS LAI product (Figure 5A). This is also supported by an excellent agreement between the BEPS and tower based GPP estimates later in the season despite the similarly high MODIS LAI values being observed until the end of August, which also agrees with the observations for similar Scots pine stands by Heiskanen et al. (2012).

Finally, Figure 6 is used to crosscheck the LAI impact. The model was tested with constant LAI values of 1 and 2 as input throughout 2019. A closer agreement between GPP_BEPS and GPP_tower estimates using tower observations is observed from the beginning of the season till the beginning of July (DOY 182), confirming the LAI overestimation during the period by the MODIS LAI product was the cause of the disagreement between the GPP_BEPS and GPP_tower estimates. The sudden, concurrent drops in GPP values predicted by both BEPS and tower observations are caused by cloudy conditions, limiting the incoming irradiation below optimal levels for GPP production.

Figure 6