Evaluating gridded BIOME-BGC for simulating LAI at Kasilian watershed-Iran

Abstract Leaf Area Index is a vegetation characteristic which affects water cycle in watersheds by intercepting water. Hence, we decided to simulate LAI with high spatial resolution for different land covers in Kasilian watershed located in Iran from 2004 to 2013. For this purpose, the gridded BIOME-BGC was developed with 500 m spatial resolution which includes 319 grids within Kasilian watershed domain. In this study, The particle swarm optimization (PSO) algorithm was applied to find a combination of 3 parameters including maximum stomatal conductance (Gs_max), new fine root C: new leaf C (FRC: LC), and canopy average specific leaf area which contribute to the best fit between BIOME-BGC LAI and MODIS LAI in each pixels. According to the PSO results, this algorithm indicated an appropriate performance during the calibration and the root mean square error (RMSE) between BIOME-BGC LAI and MODIS LAI converged to an optimum point during the primary iterations in all pixels. Moreover, the computed per cent error (PE) values in all pixels were less than approximately 30% demonstrating that LAI was simulated by BIOME-BGC with a reasonable accuracy for three-dominant land covers Deciduous Broadleaf Forest, Shrub, and C3 grasses in the case study.


Introduction
LAI is a dimensionless ratio of the leaf area per unit of the ground area (m 2 /m 2 ). It controls many ecological and hydrological processes such as photosynthesis, evapotranspiration, and run-off in watersheds. Among several land cover indices, the daily or monthly LAI is used to compute the accurate amount of water interception in the watershed hydrological modelling (Gigante, Iacobellis, Manfreda, Milella, & Portoghese, 2009;Tesemma, Wei, Peel, & Western, 2015;Wu, Gu, Luo, Shi, & Yu, 2014).
Despite its importance, investigating the relationship between LAI and other factors such as atmospheric and ecological conditions is so complicated. LAI is dependent on water, carbon, and energy fluxes in terrestrial ecosystems (White, Thornton, Running, & Nemani, 2000). In order to simulate LAI, an ecosystem model simulator like BIOME-BGC has been used in many studies (Wang et al., 2016;Zhang, Huang, & Lian, 2015). For example, Biome-BGC was used to simulate carbon fluxes in crop farms in china (Wang, Watanabe, & Ouyang, 2005). They demonstrated that BIOME-BGC was capable to simulate net primary production (NPP) and LAI. We chose BIOME-BGC model, because of its appropriate ability to simulate LAI for the several types of vegetation such as the different kinds of forests, C3 and C4 grasses, and Shrubs (Kang et al., 2014).
Obviously, it is necessary to have observed LAI to calibrate the BIOME-BGC model but it is rarely measured in some study areas, especially at the watershed scale and long time periods. In the absence of field LAI measurements, remotely sensed LAI is the best way to provide LAI in long term period and large area such as a catchment (Dong et al., 2016). Here, we decided to use 8-day MODIS LAI as a proper satellite source which its capability to derived LAI have been indicated by an array of studies in the past few years (Fang, Wei, & Liang, 2012;Yan et al., 2016). MODIS LAI is used as an input data for WRF/CMAQ modelling system (Ran et al., 2016) and as observed values for calibrating the Aqua-Crop model (Trombetta, Iacobellis, Tarantino, & Gentile, 2016). In addition, using high-resolution satellite land cover maps could be a proper alternative way to provide ecosystem modelling input types of vegetation. Among several satellite land cover maps such as MODIS, SPOTVEGETATION, and MERIS, the MODIS annually land cover products which define the type of land cover at different pixel size have been used in numerous studies (Pouliot, Latifovic, Zabcic, Guindon, & Olthof, 2014;Sun, Liang, Xu, Fang, & Dickinson, 2008). The overall objective of this work was to examine the ability of BIOME-BGC for simulating LAI in a watershed located in the mountainous land. Specifically, our main focuses are as follows: (1) developing gridded BIOME-BGC with 500 m spatial resolution, (2) using MODIS LAI to evaluate the BIOME-BGC, and (3) using the PSO to find the best fit between simulated and observed LAI.

Case study
Kasilian watershed covers an area of approximately 68 km 2 located in the mountainous lands of Mazandaran province in the north of Iran. It has a cold climate in southern parts and humid climate in central and northern parts. The highest and the lowest elevations are 3349 and 1120 m, and the annual average rainfall is about 700 mm, and the annual average of temperature is about 12 °C. The location of watershed and climatological station are shown in Figure 1.

BIOME-BGC model
BIOME-BGC is a biogeochemical model which simulates main physiological processes such as photosynthesis, evapotranspiration, respiration, and decomposition within the terrestrial ecosystems (Running & Hunt, 1993;Thornton et al., 2002). One of the model outputs is LAI produced at daily, monthly, and annually time steps. BIOME-BGC excel version 5 was used in this study that has the ability to simulate LAI for seven types of land covers including Deciduous Broadleaf Forest (DBF), Evergreen Broadleaf Forest, Deciduous Needleleaf Forest, Evergreen Needleleaf Forest, Shrubs, C3 grasses, and C4 grasses. BIOME-BGC is performed for each land cover separately. Hence, due to this fact that land cover is spatially distributed in large areas (e.g., watersheds), we developed gridded BIOME-BGC in this study. Due to lack of high-resolution land cover map and LAI field data in the case study, we used MODIS land cover map and MODIS LAI to address this problem.

MODIS land cover
This annually land cover map has provided the types of vegetation for all parts of the world at 500 m spatial resolution since 2001 (Friedl et al., 2010). This product has five land cover classification scheme including the International Geosphere-Biosphere Programme (IGBP) classification scheme (Friedl et al., 2002), the University of Maryland classification scheme (UMD) (Hansen, DeFries, Townshend, & Sohlberg, 2000), the LAI/FPAR Biome scheme (Myneni, Ramakrishna, Nemani, & Running, 1997), the Biome Classification Scheme (Running, 1994), and the Plant Functional Type Scheme (Bonan et al., 2002). Land cover types derived from The LAI/FPAR Biome scheme is more consistent with BIOME-BGC land cover types; therefore, the LAI/FPAR Biome scheme is selected as appropriate land cover classification in this study. This land cover map which involves 319 pixels was obtained in 2004 ( Figure  2). According to this map, 60% of Kasilian watershed is covered by DBF. Also, 15% is covered by savanna and Shrubs, and 25% is covered by crops like broadleaf and cereal crops.

MODIS LAI
This product has produced average LAI values at the 500 m spatial resolution (same as MODIS land cover) and at 8-day temporal resolution since February 2000 (Myneni, Knyazikhin, & Park, 2015). The algorithms used in the MODIS LAI products contain the main algorithm (LUT) and the backup algorithm. The main algorithm uses MODIS spectral data of red (648 nm) and infrared (858 nm) which reflect from the Earth's surface and simulate the LAI values by 3D radiative transfer model (Knyazikhin, Martonchik, Myneni, Diner, & Running, 1998). The backup algorithm estimates LAI values using empirical relations between Normalized Differences Vegetation Index and LAI and Fraction of Photo-synthetically Active Radiation (FPAR).

Data collection for BIOME-BGC
In order to run the BIOME-BGC model, it needs three major information including meteorological data, site characteristics, and ecophysiological parameters. In this study, 319 pixels derived from MODIS land cover with 500 m spatial resolution were considered as BIOME-BGC simulation sites. MODIS land cover reclassified into BIOME-BGC database land cover which is shown in Table 1. The daily maximum temperature, minimum temperature, and precipitation from January 2004 to December 2013 were provided from the climatological station located in the centre of Kasilian watershed ( Figure 1). The total required meteorological data for each pixel were simulated by Mountain Climate Simulator (MT-CLIM) model (Thornton, Hasenauer, & White, 2000). Information of the site characteristics such as effective soil depth and soil texture (the per cent of clay, silt, and sand) was obtained from soil map and soil profiles provided by the agricultural ministry of Iran in 1995. Furthermore, because there were no measurements for ecophysiological parameters in the case study, the default ecophysiological parameters values of BIOME-BGC model were used (White et al., 2000).

Calibration phase
Output variables of BIOME-BGC model can be adjusted by an array of the ecophysiological parameters. The aim of this study is simulating LAI; therefore, three effective parameters on LAI were selected for tuning the BIOME-BGC LAI output. Here, we selected three parameters including Maximum stomatal conductance (Gs_max), new fine root C: new leaf C (FRC: LC), and canopy average Specific Leaf Area (SLA) for calibrating the BIOME-BGC model. Then, in order to find the best combination of these parameters that contributes to the minimum Root Mean Square Error (RMSE), the PSO algorithm was chosen (Equation (1)). Also, the per cent error index (PE) was used for demonstrating optimization performance results (Equation (2)).

PSO algorithm
One of the well-known optimization methods developed by Kennedy and Eberhart in 1995 is the PSO algorithm (Eberhart & Kennedy, 1995). This algorithm consists of two main loops.

First loop: the production of initial population
Firstly, random values were generated in the range of each parameter. These parameter values introduce to BIOME-BGC and RMSE between MODIS LAI and BIOME-BGC LAI were obtained for each random combination of parameter values. In this loop, the current positions of parameter values are the personal best position (P BEST) and the position of parameter values that lead to minimum RMSE among all results considered as global best position (G BEST). The range of calibration parameters is shown in Table 2.

Second loop: iteration loop (updating the position of parameter values)
In this loop, regarding to the previous parameter value position (X(j,i)) and new parameter value velocity (V(j,i + 1)), the new parameter value position (X(j,i + 1))  R2 coefficients are random values between 0 and 1. In the first loop, the iteration number (i) and the particles velocity value (V(i,j)) are 0. The inertia weight (W(i)) is the coefficient of previous parameter value velocity which defines how much the previous velocity is effective on current parameter value velocity. One of the proper ways to adjust this coefficient is a linear reduction. In this equation, the total number of iterations is i max, Wmax and Wmin are 0.9 and 0.4, respectively (Equation (5)).
These new parameter values are introduced to BIOME-BGC model. Then, the position of the personal best (PBEST) and the global best (GBEST) are updated, and this process continues until the algorithm achieves the total number of iterations. In this study, we have set the initial population equal to 100, number of iterations equal to 12, and the value of both convergence coefficients (C1 and C2) equal to 0.5.
was computed with Equation (3). The iteration number and value parameter are i and j, respectively.
With regard to the inertia weight factor (W(i)), parameter value distance to personal best position (Pbest (j,i) − X(j,i)), and parameter value distance to global best position (Gbest (j,i) − X(j,i)), the new velocity of parameter value (V(j,i + 1)) is updated (Equation (4)).
In Equation (4), C1 and C2 are coefficients that control converging speed to the final result, and R1 and (1) (2)  values in iteration loops. Based on the results, it is clear that convergence speed is noticeably high and RMSE values converge to the lowest RMSE in approximately 5th iteration. Also, PSO results for other pixels are the same with these results, and it is demonstrated that LAI can be well adjusted by the three calibration parameters (Gs_max, FRC: LC, and SLA) used in this study for all pixels. Generally, this study indicates that PSO algorithm has an appropriate ability to optimize BIOME-BGC.

PSO algorithm results
BIOME-BGC model was calibrated by PSO algorithm for all pixels from 2004 to 2013. As example, PSO algorithm results for 3 of 319 pixels are shown in Figure  3-5. According to these graphs, the iteration number 0 shows the 100 RMSE values in initial population loop, and the iteration number 1-12 shows the 100 RMSE  Figure 6. comparing simulated and observed laI (site 1 -dBF land cover). to simulate LAI at a watershed located in the north of Iran. Initial BIOME-BGC land cover and LAI data were extracted from MODIS land cover and MODIS LAI for 319 pixels within the case study domain. Subsequently, the BIOME-BGC model was calibrated well by the optimized parameter values resulted from PSO algorithm for different kind of land covers in all pixels. We demonstrated that the PSO algorithm can be used as a robust optimization algorithm for BIOME-BGC calibration in the future studies. Also, evaluation results of BIOME-BGC showed the good capability of this model to simulate LAI for three-dominant land covers (DBF, C3 grasses, and Shrubs). Moreover, our results demonstrate that BIOME-BGC has a good ability to capture the phenology of plants.
Nevertheless, it should be noted that MODIS LAI products uncertainty and the inability of BIOME-BGC model to simulate the rapid growth of LAI in the early growing season increase the total error in LAI simulation. This problem may be improved by changing codes of BIOME-BGC model that needs further studies.
Therefore, PSO algorithm can be used to calibrate the BIOME-BGC model in subsequent studies. Moreover, it should be noted that other ecophysiological parameters can be used for calibrating the BIOME-BGC model with LAI, GPP, and NPP in further studies to examine the ability of the PSO algorithm.

BIOME-BGC evaluation results
In this section, BIOME-BGC performed with calibrated parameters computed by PSO algorithm. Consequently, per cent error values in all pixels become less than 30%. This indicates that LAI was simulated by BIOME-BGC with reasonable accuracy. As example, the calibrated value of parameters and per cent error values for 14 of 319 pixels are shown in Table 3.
Moreover, as an example, the graphs which compare 8-day LAI resulted from calibrated BIOME-BGC and 8-day MODIS LAI (2004-2013 for site number 1, 167, and 43 are shown in Figure 6-8, respectively. Based on these graphs, it is evident that simulated and observed LAI trends are generally similar and the seasonal timing of growing season is quite captured. However, these trends show discrepancies in early growing season in all the years. This is due to the fact that Biome-BGC phenology does not allow rapid growth of LAI in the early growing season. To address this problem, further studies should be considered to modify the internal BIOME-BGC codes such as the phenology codes and the logic for C-allocation. In addition, by comparing BIOM-BGC and MODIS annual watershed average LAI from 2004 to 2013, it is indicated that the general magnitude of LAI was estimated with good accuracy in all pixels from 2004 to 2013 ( Figure 9). However, it is obvious that the BIOME-BGC LAI was under estimated. Finally, these results show that LAI was estimated with reasonable accuracy in all pixels with three-dominant land covers (DBF, Shrubs, and C3 grasses) within Kasilian watershed domain.
Simultaneous estimation of daily solar radiation and humidity from observed temperature and precipitation: An application over complex terrain in Austria.