Forest Carbon Flux Simulation Using Multi-Source Data and Incorporation of Remotely Sensed Model with Process-Based Model

Abstract

Forest carbon flux is critical to climate change, and the accurate modeling of forest carbon flux is an extremely challenging task. The remote sensing model (the MODIS MOD_17 gross primary productivity (GPP) model (MOD_17)) has strong practicability and is widely used around the world. The ecological process (the Biome-BioGeochemical Cycles Multilayer Soil Module model (Biome-BGCMuSo)) model can describe most of the vegetation’s environmental and physiological processes on fine time scales. Nevertheless, complex parameters and calibrations pose challenges to the application and development of models. In this study, we optimized all the input parameters of the MOD_17 model for the calibration of the Biome-BGCMuSo model to obtain GPP with continuous spatiality. To determine the contribution of input parameters to the GPP of different forest types, an Extended Fourier Amplitude Sensitivity Test (EFAST) was performed on the Biome-BGCMuSo model firstly. Then, we selected the sample points of each forest type and its different ecological gradients (30 for each type), using the GPP simulation value of the optimized MOD_17 model corresponding to the time and space scale to calibrate the Biome-BGCMuSo model, to drive the calibrated Biome-BGCMuSo, and we simulated the different forest types’ net primary productivity (NPP). According to dendrochronological measurements, the NPP simulation results were verified on the whole regional scale. The results showed that the GPP values of different forest types were highly sensitive to C:N_leaf (C:N of leaf), SLA1 (canopy average specific leaf area in phenological phase 1), and FLNR (fraction of leaf N in Rubisco). The coefficient of determination (R²) between the simulated forest NPP and the measured NPP was 0.64, and the root-mean-square (RMSE) was 26.55 g/C/m²/year. Our study aims to reduce uncertainty in forest carbon fluxes simulated by the Biome-BGCMuSo model, providing feedback for understanding forest ecosystem carbon cycling, vegetation productivity, and climate change.

1. Introduction

Forests are an important part of terrestrial ecosystems and play an important role in the global carbon cycle [1,2]. The forest ecosystem is the main carbon pool on land, and its carbon content accounts for more than half of the entire terrestrial ecosystem [3]. However, due to the influence of climate change, human activities, and the natural growth and degradation of forest vegetation, it shows a dynamic change law [4,5] and leads to uncertainties in the simulation of carbon fluxes in different regions and forest types, so it is always a challenge to model forest carbon fluxes accurately [6,7]. Vegetation productivity is an important indicator for evaluating carbon flux including gross primary productivity (GPP), net primary productivity (NPP), and net ecosystem exchange (NEE) [8,9].

The current methods for quantitative estimation of forest ecosystem productivity can be summarized as: inventory method, eddy correlation method, and model method. The forest inventory method usually estimates the carbon storage in the forest by measuring the change in biomass, and then obtains the net primary productivity through the conversion factor [10]. The advantage of the inventory method is to directly measure the carbon storage of vegetation and soil at the sample scale, but it has a long inventory cycle and suffers great uncertainty from sample to area. The eddy covariance method has been proven to be accurate in detecting changes in ecosystem carbon fluxes on fine time scales [11,12], but since the sites are scarce and costly that makes it difficult to extend to regional-scale applications. With the development of remote sensing technology, people can take advantage of the macroscopic and dynamic characterizations of remote sensing and its unique application potential in large-scale surveys that are spatially continuous. However, remote sensing can only obtain instantaneous observation data and cannot provide continuous temporal characterization; moreover, pure remote sensing information is also difficult for reflecting the interaction between the forest ecosystem and the outside world [13,14].

For vegetation productivity simulation, current widely used models include two types: light use efficiency (LUE) models and ecological process models. The LUE model based on remote sensing can obtain regional-scale and time-continuous information. Representative models include CASA [15], GLO-PEM [16], MODIS MOD17 model [17], C-Fix [18], etc. However, such models cannot reflect the internal mechanism of the forest ecosystem or the interaction between the forest ecosystem and the surrounding environment [13]. The ecological process model can describe most of the physiological and ecological processes of vegetation on fine time scales [19], such as photosynthesis, respiration and transpiration. Representative models include TEM [20], BEPS [21], Biome-BGC [22], etc. However, the ecological process models require a large amount of detailed ground information as input data, and a series of related initialization states and vegetation physiological ecology; these will result in great uncertainties especially in the regional application.

Many scholars in previous studies have used the incorporation of LUE model and ecological process model to reduce the uncertainty of simulations of carbon fluxes. Chiesi et al. [23] proposed a method for parameter calibration of the process model based on the incorporation of the remote sensing model (C-Fix)-process model (Biome-BGC); the mothed suppressed the uncertainty of the simulation results, but the mothed was only applied to a single ecological gradient in the Mediterranean region and its generalizability remains unknown. The optimized MOD_17 model was used to couple the process model Biome-BGC in the region, and the simulation results were demonstrated more accurately and reliably [24,25,26], but this method was only suitable for single dominant tree species (Picea crassifolia) in this region. Sánchez-Ruiz et al. [27] simulated forest-year NPP in the Spanish peninsula based on the Monteith method and the Biome-BGC model incorporation, but it is limited to a temperate oceanic climate area.

According to the research gaps mentioned above, our research aims were: (1) to explore the sensitivity parameters for different forest types in the study area using the extended Fourier analysis method; (2) to strengthen the robustness and applicability of the incorporation of the remote sensing model and Biome-BGCMuSo model in obtaining forest carbon fluxes; and (3) to reveal the temporal and spatial trends of forest NPP in typical forest-covered areas such as the cold temperate zone of China, and the impact of climate change on forest carbon fluxes as well.

2. Materials and Methods

2.1. Study Area

The Greater Khingan Mountains in Inner Mongolia (119°36′26″–125°24′10″E, 47°03′26″–53°20′00″N) are located in the northeast of Inner Mongolia Autonomous Region (Figure 1) and are one of the four major state-owned forest areas in China. The main ecological function area is 1.067 × 105 km² (about 46% of the entire Greater Khingan Mountains), and the main forest coverage rate is 79.56%. The region plays an irreplaceable and important role in water conservation, oxygen production, carbon sequestration, soil conservation, biological protection, and genetic diversity. This area belongs to the northern forest belt of Eurasia. It is an important barrier for China’s northern homeland security system and a key forest carbon pool in China. The soil types in the study area are mainly brown coniferous forest and swamp soil, with a thickness of about 40–60 cm [28]. It has an extremely important strategic position of ecological construction. This area is located in the high latitude and alpine zone, with an altitude of 400–1500 m. It belongs to the cold temperate continental monsoon climate. The annual average temperature is −3.5 °C, extreme temperatures reach −50.2 °C, and the annual precipitation levels from 300 to 450 mm are usually concentrated in the range from July to September; it is one of the most sensitive areas in the world to respond to climate change and the soil layer is poor.

2.2. Data and Preprocessing

2.2.1. Meteorological Data and Dendrochronological Measurements

Model input meteorological files of Biome-BGCMuSo include daily minimum and maximum temperature, daily average temperature, daily precipitation, daily shortwave download radiation (Sard), daily average water vapor pressure difference (VPD) and sunshine hours (Daylen). Therefore, we collected the ground meteorological element-driven data (CMFD) with high spatial and temporal resolution from 1992 to 2018 in China [29]. This dataset is produced by combining the existing international Princeton reanalysis data, GLDAS data, GEWEX-SRB radiation data, and TRMM precipitation data as the background field, and combining the conventional meteorological observation data of the China Meteorological Administration, with a spatial resolution of 0.1° [30,31]. To keep the same size as the Biome-BGCMuSo model simulation unit, we first used python (version 3.10) to convert 75,920 (8 × 365 × 27) scene CFMD data from network Common Data Form (NetCDF) to Tag Image File Format (TIFF) and then use the nearest neighbor interpolation method to resample them into a resolution of 1 km (Albers projection, WGS-84 coordinate system). Then, we extracted daily minimum temperature, daily maximum temperature and daily precipitation data from the CMFD data. Finally, VPD, Sard and Daylen were calculated using the MT-CLIM 4.3 algorithm for the Biome-BGCMuSo model-driven.

In order to verify the feasibility of parameters of the Biome-BGCMuSo model and the accuracy of the simulation results, 46 sample plots were selected from the Greater Khingan Mountains in Inner Mongolia during August 2013 and August 2016. The main forest types included Larix gmelinii and Betula platyphylla pure forest, Larix and Betula mixed forests, etc., and forest structural parameters and dendrochronological measurements were carried out. Collected samples were packaged back to the laboratory after the dendrochronological measurements, and the annual ring width was measured by a WinDENDRO annual ring analyzer with a measurement accuracy within 0.001 mm. The COFECHA algorithm was used for manual cross-dating, and the DBH value of each standard sample tree was obtained year by year. According to the biomass growth equations of each forest type [32,33], the annual forest biomass increment was calculated, and also the average ratio of the biomass to the carbon of Larix gmelini (0.5211) and Betula platyphylla (0.491) [34]; the forest NPP for each plot was obtained, then they were used to verify the model simulation results.

2.2.2. Remote Sensing Data

In this study, the time series (2012–2015) of the GLASS LAI/fPAR product [35] was obtained to optimizeMOD_17 model. GLASS LAI is a high-quality remote sensing data product with a long time series, no missing data, and global coverage. During the production procedure of this product, cloud removal, snow removal, missing value filling, filtering and other processes had been carried out to reduce some errors and increase the practical applications and the possibility of obtaining data, and they were provided free of charge to users.

The land-use data were produced by using the fine resolution observation and global land cover monitoring (FROM-GLC-seg) product (http://0-data-ess-tsinghua-edu-cn.brum.beds.ac.uk/index.html (accessed on 25 June 2021)) from Landsat TM and ETM+ data, using the Random Forest (RF) classifier with high overall accuracy. The CO₂ concentration data are downloaded from the observational dataset from the Mauna Loa Observatory in Hawaii (https://gml.noaa.gov/obop/mlo/programs/esrl/co2/co2.html (accessed on 21 July 2021)). The soil texture dataset is provided by the Resource and Environmental Science Data Center (RESDC) of the Chinese Academy of Sciences (https://www.resdc.cn (accessed on 29 July 2021)); it is a 1:1 million soil-type map which was mapped on the second soil survey in China profile data. The soil types were divided into three categories: sand, silt, and clay; their values were in percentages, and the spatial resolution is 1 km. The DEM data were obtained from ASTER GDEM (https://www.gscloud.cn (accessed on 5 September 2021)) and it was used to calculate the slope and aspect datasets. In order to maintain consistency with the simulated spatial scale of the Biome-BGCMuSo model, we reprojected and resampled the above data into Albers projection, WGS-84 coordinate system with the resolution of 1 km.

2.3. Methods

To assess the spatiotemporal dynamic changes in forest NPP, the method we used included four steps: first, the MOD_17 model was optimized using CMFD data, GLASS LAI/fPAR products, and the value of ${\epsilon }_{\mathit{max}}$ for the main forest types in the study area; second, Biome-BGCMuSo was subjected to sensitivity analysis using the EFAST method, and the most sensitive parameters were retained for the next calibration step; third, sample plots representing different forest types were selected in the region, and the process model’s sensitive parameters for each forest type were calibrated based on the model–model coupling method; and fourth, the simulated forest NPP was validated using dendrochronological measurements, together with the spatiotemporal dynamic changes in forest NPP and its relationship with climatic factors were analyzed. Figure 2 shows the overall flowchart of the approach used in this study. Details of each step are given in subsequent sections.

2.3.1. The MOD_17 Model

The MOD_17 model [36] was developed according to the radiation conversion efficiency theory of Monteith [37], a remote sensing-based light energy utilization model. The MODIS GPP algorithm can be expressed as:

$$\mathit{GPP}{=\epsilon }_{\mathit{max}}\text{ × }\mathrm{f}\left({\mathit{T}}_{\mathit{min}}\right)\text{ × }\mathrm{f}\left(\mathit{VPD}\right)\text{ × }\mathrm{APAR}$$

(1)

$$\mathrm{APAR}=\mathrm{PAR}\text{×}\mathit{f}\mathrm{PAR}$$

(2)

where ${\epsilon }_{\mathit{max}}$ is the maximum LUE; $\mathrm{PAR}$ is the incident photosynthetically active radiation; $\mathit{f}\mathrm{PAR}$ is the fraction of PAR absorbed by vegetation; $\mathrm{f}\left({\mathit{T}}_{\mathit{min}}\right)$ and $\mathrm{f}\left(\mathit{VPD}\right)$ are the scalars of the minimum air temperature and VPD, respectively. The values of $\mathrm{f}\left({\mathit{T}}_{\mathit{min}}\right)$ and $\mathrm{f}\left(\mathit{VPD}\right)$ vary from in the range of [0, 1]. The expressions of $\mathrm{f}\left({\mathit{T}}_{\mathit{min}}\right)$ and $\mathrm{f}\left(\mathit{VPD}\right)$ are as follows:

$$\mathrm{f}\left({\mathit{T}}_{\mathit{min}}\right)=\{\begin{array}{c}{0,\mathit{T}}_{\mathit{min}}{\text{<}\mathit{T}}_{\mathit{minmin}}\\ \frac{{\mathit{T}}_{\mathit{min}}{-\mathit{T}}_{\mathit{minmin}}}{{\mathit{T}}_{\mathit{minmax}}{-\mathit{T}}_{\mathit{minmin}}}{,\mathit{T}}_{\mathit{minmin}}{\text{<}\mathit{T}}_{\mathit{min}}{\text{<}\mathit{T}}_{\mathit{minmax}}\\ {1,\mathit{T}}_{\mathit{min}}{\text{>}\mathit{T}}_{\mathit{minmax}}\end{array}$$

(3)

$$\mathrm{f}\left(\mathit{VPD}\right)=\{\begin{array}{c}{0,\mathit{VPD}\text{>}\mathit{VPD}}_{\mathit{max}}\\ \frac{{\mathit{VPD}}_{\mathit{max}}-\mathit{VPD}}{{\mathit{VPD}}_{\mathit{max}}{-\mathit{VPD}}_{\mathit{min}}}{,\mathit{VPD}}_{\mathit{min}}{\text{<}\mathit{VPD}\text{<}\mathit{VPD}}_{\mathit{max}}\\ {1,\mathit{VPD}\text{<}\mathit{VPD}}_{\mathit{min}}\end{array}$$

(4)

where $\mathit{VPD}$ and ${\mathit{T}}_{\mathit{min}}$ are 8-day averaged values, and ${\mathit{T}}_{\mathit{minmin}}$, ${\mathit{T}}_{\mathit{minmax}}$, ${\mathit{VPD}}_{\mathit{max}}$, and ${\mathit{VPD}}_{\mathit{min}}$ are parameters derived from Biome Properties Look-Up Table, and their values vary with vegetation types.

In this study, we first collected ${\epsilon }_{\mathit{max}}$ for the test site according to the literature studied at similar study areas with similar latitudes consistent with the cold-temperate vegetation of the Greater Khingan Range (Table S1). Then, we used the processed meteorological data to optimize $\mathrm{f}\left({\mathit{T}}_{\mathit{min}}\right)$ and $\mathrm{f}\left(\mathit{VPD}\right)$. Finally, the GLASS LAI/fPAR products were used to simulate the different forest-type GPP values from 2012 to 2015 and the values were used to calibrate the Biome-BGCMuSo model.

2.3.2. Biome-BGCMuSo Model

The Biome-BGCMuSo model was developed from the Biome-BGC model series [38], and we use the latest version of the model, Biome-BGCMuSo 6.2. Compared with Biome-BGC, Biome BGCMuSo improves the model’s ability to simulate carbon and water cycles in managed ecosystems and provides options for managing cropland, grassland, and forest, among which the most important modules of the model are the carbon flux module, phenology module and soil flux module. The GPP of vegetation is calculated by Farquhar’s [39] photosynthetic reaction mechanism model in the carbon flux module, and maintenance respiration is a function of reactive nitrogen content. Leaf development affects the accumulation of carbon and nitrogen in leaves, stems, roots, and litter, and is calculated by the phenology module. The soil flux module describes the decomposition of dead plant material and soil carbon pools [40].

The time scale of the Biome-BGCMuSo model is daily, and at least four input files are used for each execution: they are initialization file (INI file), meteorological data file (MET file), soil property file (SOIL file), and physiological ecological parameter file (EPC file). The simulation process of the model is divided into two stages: the first phase is the spinup simulation, where the model runs from a lower soil carbon and nitrogen level and it stops until a steady state is reached; the second phase is the normal simulation stage, which uses the initial values obtained in the initialization stage to simulate NPP of GPP [38]. To bring the model to an equilibrium state, daily meteorological data from 1992 to 2018 were selected and used in this study; in the procedure, the first 20 years were used for the spinup simulation phase and the last 7 years for the normal simulation phase.

2.3.3. Sensitivity Analysis

The EFAST method is a global sensitivity analysis method of model variance analysis proposed by Sailteli et al. [41], which combines the advantages of the Sobel method and the Fourier amplitude sensitivity test method. It has the advantages both of stability and fast calculation. Considering the complexity of the physiological and ecological parameters of the Biome-BGCMuSo model and the close interaction between the parameters, the influence of each physiological and ecological parameter on the output of the model is mainly represented by the first-order sensitivity index and the total sensitivity index. The first-order sensitivity index mainly reflects how much the variance of the model output can be reduced when a certain parameter value is fixed as the true value. The total sensitivity index refers to the sum of the first-order sensitivity index and the interaction sensitivity index between this parameter and other parameters. The total sensitivity index is suitable for global, quantitative, and model-independent sensitivity analysis.

The EFAST method often uses a search function to search in the multidimensional space of parameters, and obtains the sensitivity index by estimating the contribution rate of each data input Xi to the variance of the result Y, which can be expressed as:

$$\mathit{Y}=f\left(X\right)=f\left(X_1,X_2,\cdots ,X_n\right)$$

(5)

where X_i (n = i) indicates that each input parameter has its own value range to indicate its uncertainty. The total variance of the model output can be expressed as:

$${\mathit{V}}_{\mathit{Y}}={\displaystyle \sum }_{\mathit{i}}{\mathit{V}}_{\mathit{i}}+{\displaystyle \sum }_{\mathit{i}}{\displaystyle \sum }_{\mathit{j}\text{>}\mathit{i}}{\displaystyle \sum }_{\mathit{k}\text{>}\mathit{j}}{\mathit{V}}_{\mathit{ijk}}{+\text{…}+\mathit{V}}_{1,2,\text{…}.\mathit{n}}$$

(6)

where ${\mathit{V}}_{\mathit{i}}=\mathit{V}\left(\mathit{E}\left({\mathit{Y}/\mathit{X}}_{\mathit{i}}{=\mathit{x}}_{\mathit{i}}^{*}\right)\right)$, ${\mathit{V}}_{\mathit{ij}}=\mathit{V}\left(\mathit{E}\left(\frac{\mathit{Y}}{{\mathit{X}}_{\mathit{i}}}{=\mathit{x}}_{\mathit{i}}^{*}{,\mathit{X}}_{\mathit{i}}{=\mathit{x}}_{\mathit{j}}^{*}\right)\right){-\mathit{V}}_{\mathit{i}}{-\mathit{V}}_{\mathit{j}}$, and other variance calculation methods can be deduced in turn. $\mathit{E}\left({\mathit{Y}/\mathit{X}}_{\mathit{I}}{=\mathit{x}}_{\mathit{i}}^{*}\right)$ represents the expected value of Y when X_i takes a fixed value xi. The first-order sensitivity index of the parameter ${\mathit{S}}_{\mathit{i}}{=\mathit{V}}_{\mathit{i}}/\mathit{V}$ has a value range between 0 and 1, and the sensitivity index quantifies the influence of the parameter X_i on the results. The total sensitivity index STi is the sum of the first-order sensitivity index of this parameter and the interaction sensitivity index of this parameter with other parameters: ${\mathit{S}}_{\mathit{Ti}}{=\mathit{S}}_{\mathit{i}}{+\mathit{S}}_{\mathit{ij}}{+\mathit{S}}_{\mathit{ijk}}{+\mathit{S}}_{\mathit{ijk}}{+\mathit{S}}_{\mathit{ijk}\text{…}\mathit{m}}$ [42].

Sensitivity analysis of the Biome-BGCMuSo model was carried out in SimLab2.2. A total of 51 parameters were selected in the EPC input file (

Table 1. Selected input parameters of Biome-BGCMuSo used in the sensitivity analysis.

Number	Parameter	Description	Unit
1	TGP	transfer growth period as fraction of growing season	prop.
2	LGS	litterfall as fraction of growing season	prop.
3	LFRT	annual leaf and fine root turnover fraction	1/yr
4	LWT	annual live wood turnover fraction	1/yr
5	FM	annual fire mortality fraction	1/yr
6	WPM	whole-plant mortality fraction in vegetation period	1/yr
7	C:Nleaf	C:N of leaves	kgC/kg N
8	C:Nlit	C:N of leaf litter, after retranslocation	kgC/kg N
9	C:Nfr	C:N of fine roots	kgC/kg N
10	C:Nlw	C:N of live wood	kgC/kg N
11	C:Ndw	C:N of dead wood	kgC/kg N
12	DMCleaf	dry matter carbon content of leaves	(kgC/kgDM)
13	DMClit	dry matter carbon content of leaf litter	(kgC/kgDM)
14	DMCfr	dry matter carbon content of fine roots	(kgC/kgDM)
15	DMCf	dry matter carbon content of fruit	(kgC/kgDM)
16	DMCs	dry matter carbon content of soft stem	(kgC/kgDM)
17	DMClw	dry matter carbon content of live wood	(kgC/kgDM)
18	DMCdw	dry matter carbon content of dead wood	(kgC/kgDM)
19	Llab	leaf litter labile proportion	DIM
20	Lcel	leaf litter cellulose proportion	DIM
21	FRlab	P fine root labile proportion	DIM
22	FRcel	fine root cellulose proportion	DIM
23	Flab	fruit litter labile proportion	DIM
24	Fcel	fruit litter cellulose proportion	DIM
25	DWcel	dead wood cellulose proportion	DIM
26	Wint	canopy water interception coefficient	1/LAI/d
27	k	canopy light extinction coefficient	DIM
28	SPLR	all-sided to projected leaf area ratio	DIM
29	LAIall:pro	ratio of shaded SLA:sunlit SLA	DIM
30	FLNR	fraction of leaf N in Rubisco	DIM
31	gsmax	maximum stomatal conductance (projected area basis)	m/s
32	gcl	conductance (projected area basis)	m/s
33	gbl	boundary layer conductance (projected area basis)	m/s
34	SW	stem weight corresponding to maximum height	(kgC)
35	Rdmax	maximum depth of rooting zone	(m)
36	GR	growth resp per unit of C grown	(prop.)
37	MRpern	maintenance respiration in kg C/day per kg of tissue N	(kgC/kgN/d)
38	NSC:Scmax	theoretical maximum prop. of non-structural and structural carbohydrates	(DIM)
39	NSCMR	of non-structural carbohydrates available for maintenance respiration	(DIM)
40	SWClim2	minimum of soil moisture limit2 multiplicator (full anoxic stress value)	prop
41	VPDs	vapor pressure deficit: start of conductance reduction	Pa
42	VPDc	vapor pressure deficit: complete conductance reduction	Pa
43	TRwsl	turnover rate of wilted standing biomass to litter	prop
44	TRcwl	turnover rate of non-woody cut-down biomass to litter	prop
45	SLA1	canopy average specific leaf area in phenological phase 1	m2/kg
46	SLA2	canopy average specific leaf area in phenological phase 2	m2/kg
47	SLA3	canopy average specific leaf area in phenological phase 3	m2/kg
48	SLA4	canopy average specific leaf area in phenological phase 4	m2/kg
49	SLA5	canopy average specific leaf area in phenological phase 5	m2/kg
50	SLA6	canopy average specific leaf area in phenological phase 6	m2/kg
51	SLA7	canopy average specific leaf area in phenological phase 7	m2/kg

), and the variation range of each parameter was fluctuated around 20% [43]. To evaluate the selected 51 input parameters, we first sampled each parameter uniformly by the Monte Carlo method according to its range (Tables S2, S3 and S4), with a sampling frequency of 3315 (51 × 65) for each forest type. The total sampling frequency was 9945 (3315 × 3). Next, the Biome-BGCMuSo model was run in batches with the generated multiple sets of input parameters, and the statistical annual GPP averages of different forest types from 2012 to 2018 were input into SimLab2.2 for sensitivity analysis. Finally, the sensitivity indices (Si and STi) were divided into 3 groups: the indices with values greater than 0.2 were set as high sensitivity indices; the values between 0.1 and 0.2 were set as medium sensitivity indices; and the indices with values less than 0.1 were set as low sensitivity indices.

2.3.4. Incorporation of the MOD_17 and Biome-BGCMuSo

After sensitivity analysis by the EFAST method, the sensitive parameters of the Biome-BGCMuSo model (sensitivity index greater than 0.1) were retained, the sensitivity parameters were sorted according to the level of the parameter sensitivity, and the sensitivity parameters were calibrated for each forest type with respect to specific environmental factors. First, we selected sample points (30 for each type) of each forest type and its different ecological gradients (site conditions, meteorological conditions, stand structure, productivity, etc.) to simulate GPP from 2012 to 2015 (the number of samples is: 30 × 4 years × 23 layers × 3 = 8280). Then, the optimized MOD_17 model was used to fit the GPP simulation values of the Biome-BGCMuSo model corresponding to the spatiotemporal scale. During the fitting process, the parameters were adjusted according to the order of the sensitivity index. The globally representative calibration parameters in the process model are obtained when the fitting results are stable.

3. Results

3.1. Optimization of the MOD_17 Model

The optimized MOD_17 model was used to estimate the forest GPP in the Greater Khingan Mountains forest area of Inner Mongolia. Figure 3 shows the distribution of GPP in different forest types in 2012, 2013, 2014, and 2015. The highest GPP values was acquired by deciduous coniferous forest (DNF), next was deciduous broad-leaved forest (DBF), then the third was mixed coniferous (MF).

The highest value for DBF GPP was in 2013 and the highest value for DNF and MF GPP was in 2014.

Table 2. GPP statistical results of the MOD_17 model simulation after optimization.

Year	Forest Type	Count	Mean	STD	Min	Max	CV (%)
2012	DBF	41,583	635.48	130.39	176.43	1269.95	20.52
	DNF	7651	977.62	160.46	378.92	1585.82	16.41
	MF	22,631	641.17	114.07	160.94	1067.52	17.79
2013	DBF	41,580	737.82	118.83	108.4	1242.4	16.1
	DNF	7651	1008.21	144.6	357.08	1571.69	14.34
	MF	22,630	657.27	102.53	167.75	989.16	15.6
2014	DBF	41,575	646.63	139.75	138.71	1312.03	21.61
	DNF	7650	1000.89	158.98	388.92	1643.15	15.88
	MF	22,628	665.82	114.89	175.5	1055.16	17.25
2015	DBF	41,578	561.34	125.7	139.76	1160.08	22.39
	DNF	7649	927.53	157.95	353.94	1509.93	17.03
	MF	22,628	611.4	110.91	163.67	1012.6	18.14

shows the descriptive statistical results of the annual mean value of GPP for different forest types, including year, forest type, number of pixels (Count), mean (Mean), maximum value (Max), minimum value (Min), standard deviation (STD), and coefficient of variation (CV).

The uncertainty (CV) fluctuation of DBF had the largest range, spanning from 16.1% to 22.39%, and the range of GPP value was 108.4–1312.03 g/C/m²/year. The smallest uncertainty fluctuation range of DNF ranged from 14.34% to 17.03%, and the GPP value range was from 353.94 to 1643.15 g/C/m²/year. In general, the uncertainty fluctuation range of GPP was [14.34%, 22.39%], the annual average value of GPP ranged from 561.34 to 1008.21 g/C/m²/year, and the range of GPP value was [108.4, 1643.15] g/C/m²/year.

3.2. Sensitivity Analysis of Biome-BGCMuSo

The mean annual GPP of different forest types in the study area and the first-order sensitivity index (Si) and total sensitivity index (STi) of the selected 51 parameters are shown in Figure 4, and the sensitive parameters are on the right side of the red solid line. The total sensitivity index was greater than the first-order sensitivity index, indicating that the interaction among the 51 parameters was much larger than the influence of a single parameter. The most sensitive first-order sensitivity and total sensitivity index of the annual GPP mean of DNF, DBF and MF were C:N_leaf and SLA1. The first-order sensitivity index was 0.3, 0.43 and 0.25, respectively, and the second-order sensitivity was 0.81, 0.68 and 0.69, respectively.

The mean annual GPPs of the different forest types are all sensitive to five parameters, namely C:N_leaf, SLA1, k, FLNR and LWT. It is worth noting that the sensitivity parameters were also different between the forest types, where MF and DBF were both sensitive to TR_wsl and VPD_c, while DNF was not sensitive to these two parameters.

3.3. Calibration of Biome-BGCMuSo Mode

After the EFAST method was used to analyze the sensitivity of the physiological and ecological parameters of the Biome-BGCMuSo model, the results of calibrating the sensitive parameters of the model for different forest types were shown in Figure 5. The GPP values simulated by the optimized MOD_17 model are in agreement with the GPP values simulated by the calibrated Biome-BGCMuSo model. The best effect is DNF, R² = 0.84, RMSE = 24.16 g/C/m²/16d; followed by MF, with R² = 0.81 and RMSE = 17.25 g/C/m²/16d; the last one is DBF, with R² = 0.76 and RMSE = 15.39 g/C/m²/16d.

Through the dendrochronological measurements of 46 sample plots collected in the Greater Khingan Mountains area of Inner Mongolia in 2013 and 2016, some samples with serious interference and large measurement errors were screened. After conversion to forest NPP, they were verified with the calibrated Biome-BGCMuSo model. The results are as follows (Figure 6). The verification results showed that the calibrated Biome-BGCMuSo model simulates the forest NPP better, with R² = 0.64 and RMSE = 26.55 g/C/m²/year.

3.4. The Spatial and Temporal Dynamics Analysis of Forest NPP

To analyze the temporal and spatial trends of forest NPP in the study area, the calibrated Biome-BGCMuSo simulated forest NPP in the study area and calculated the average forest NPP from 2012 to 2018. The results are shown in Figure 7a. The annual mean value of forest NPP ranged from 26.83 to 718.27 g/C/m²/year, which was low in the middle and high in the east and west in terms of spatial distribution. Figure 7b shows the forest NPP variation trend from 2012 to 2018. It revealed that the temporal variation trend of forest NPP was not obvious. The variation range of most areas was between −42.25 and 2.5 g/C/m²/year, while the regional change trend was obvious and showed a negative growth trend. Figure 8 shows the time distribution trend of NPP in different forest types in the growing season from 2012 to 2018. The highest NPP value in the growing season is coming from deciduous coniferous forest (DNF), next was coniferous and broad-leaved mixed forest (MF), the third was deciduous broad-leaved forest (DBF). During the seven-year period, the NPP in the growing season did not change significantly and showed a slow growth trend.

Figure 9 shows the correlation analysis results of the annual mean value of NPP and climatic factors in the forest growing season in the Greater Khingan Mountains forest area of Inner Mongolia. The correlation coefficient matrix (Figure 9) shows that the annual mean NPP of the forest growing season from 2012 to 2018 has an upward trend with a correlation coefficient of 0.54. The annual mean value of NPP in the forest growing season was positively and strongly correlated with VPD and T_mean, which were 0.82 and 0.81, respectively. The weak correlation with Sard was 0.11 and the negative correlation with PRE was −0.66.

Figure 10 shows the statistical results of the annual average forest NPP dynamics and climatic factors in the Greater Khingan Mountains forest area of Inner Mongolia. The statistical data were only calculated for the forest growing season. From 2012 to 2018, the average temperature of the growing season in the study area was between 10.1~12.3 (°C), the annual precipitation was between 230~770 (mm), the average saturated vapor pressure was between 660~830 (Pa), and the annual shortwave radiation was 1.14~1.95 (106 W/m²). Among them, the change trend of forest NPP was basically consistent with the change trend of average vapor pressure and average temperature, which was opposite to the change trend of annual precipitation and the change trend of annual shortwave radiation was not obvious.

4. Discussion

This study showed that the incorporation strategy of the MOD_17 model and the Biome-BGCMuSo model can effectively reduce the uncertainty during the simulation of carbon fluxes in different forest types at the Greater Khingan Mountains, Inner Mongolia, in the cold temperate zone. By calibrating the MOD_17 model with the optimized Biome-BGCMuSo model, the correlation coefficient between the simulated forest NPP and the measured forest NPP was 0.64. Although the model was not validated with observational data such as flux sites in the study area, the simulation results were validated with two sample plot data points and tree dendrochronological measurements. Compared with the observation data such as flux sites, the sample plot data and dendrochronological measurement data had the characteristics of various disturbance scenarios and multihabitat conditions, strong spatial continuity, long-term scales, and many samples; the phenomenon revealed the effectiveness of the simulation through the incorporation model strategy and confirmed the reliability and global representation of forest carbon fluxes simulation.

The MOD_17 model has strong practicability and has been widely used at regional, national and even global scales, but the parameters of the Biome Properties Look-Up Table have significant uncertainties, which depend on the quality of remote sensing data, model parameterization schemes and empirical knowledge [44,45]. Therefore, we optimized each parameter. Among them, ${\epsilon }_{\mathit{max}}$ had not been recalibrated but the ${\epsilon }_{\mathit{max}}$ consistent with the cold temperate forest type had been collected through the literature data. The $\mathit{f}({\mathit{T}}_{\mathit{min}})$, $\mathit{f}\left(\mathit{VPD}\right)$, and PAR were calculated from CMFD data with a time-resolve of 3 h [31], providing a basis for obtaining carbon fluxes with temporal and spatial continuity. The fPAR was derived from the GLASS LAI/fPAR product, which had been proven to be reliable [35].

Biome-BGCMuSo is a newly developed process model of Biome-BGC that is widely used to simulate forest NPP, but the sensitivity analysis of many parameters and model calibration are the key steps in forest carbon flux simulation. The results of our sensitivity analysis of the Biome-BGCMuSo model showed that there were five sensitive parameters that had important effects on different forest types. However, many parameters in the Biome-BGCMuSo model were extremely insensitive to GPP, such as TGP (transfer growth period as a fraction of growing season), LGS (litterfall as a fraction of growing season), DMC_leaf (dry matter carbon content of leaves), DMC_f (dry matter carbon content of fine roots), etc. The parameters that controlled the carboxylation rate of photosynthesis, such as C:N_leaf, SLA1 and FLNR, had the greatest impact on the GPP of different forest types, and the results were consistent with White et al. [43] and Ren et al. [46]. These three parameters determined the size of the maximum carboxylation rate, and the Rubisco enzyme directly affected the carboxylation of vegetation to fix atmospheric CO₂. If the nitrogen content of Rubisco enzyme increased its activity would be promoted, thereby enhancing photosynthesis and promoting organic substance accumulation [47,48].

A strategy based on the incorporation of parametric models and ecological process models to simulate carbon fluxes with temporal and spatial continuity had been proven to be reliable and effective for Mediterranean forests and single dominant tree species in cold temperate zones [23,24,25,26], but feasibility experiments had not been carried out on different forest types in the cold temperate zone. Due to the heterogeneity of the growth environment, different forest types and the same forest type under different habitats had different physiological and ecological adaptation mechanisms, so it was difficult to calibrate the parameters of the Biome-BGCMuSo model on a regional basis, and it would result in the model distortion simulation [23]. Therefore, by incorporating the optimized MOD_17 model and the calibrated Biome-BGCMuSo model, it was effective to suppress the simulation of different forest types and it was possible to simulate carbon flux with spatial and temporal continuity.

Although we have reduced the uncertainty of the carbon flux simulation of different forest types in the Greater Khingan Mountains forest area of Inner Mongolia through the model–model coupling strategy, the simulation scale of the Biome-BGCMuSo model is 1 km, so the problem of pixel mixing cannot be avoided. The key to solving this problem is to use high-resolution remote sensing products. Although the classification products used in this study are fine classification products with a resolution of 30 m, the problem of pixel mixing cannot be avoided after resampling to 1 km. In addition, due to differences in model structure, parameter values, and driving data among different models, even if the coupling between models can improve the global representation of sensitivity parameter calibration, there are still errors in the simulation of forest ecological processes in large-area time series. In the model simulation, state parameters (such as LAI, soil moisture, freeze-thaw state, etc.) have large uncertainties, which will cause error transmission during the simulation process. Data assimilation is an objective and quantitative analysis method that integrates observations with heterogeneous, irregular distributions and different precisions on various types of spatiotemporal scales into the model, and obtains the best guesses for model parameters and state variables as much as possible [49].

The response of forest NPP to climate change is complex, and different forest types vary greatly. Climate controls the formation of forest NPP by changing environmental conditions to affect the physiological structure and process of forests, so changes in NPP can reflect the response of forest ecosystems to the environmental climate [50]. Generally, forest NPP has a positive response to precipitation changes, and has both positive and negative responses to temperature changes, but the response to temperature is stronger than that to precipitation. This study found that from 2012 to 2018, the forest NPP in the Greater Khingan Mountains, Inner Mongolia, had a positive response to the average temperature and average saturated vapor pressure, which was consistent with the results of He et al. [51] and Li et al. [52]. However, the difference is that we found a negative response between forest NPP and precipitation, which may have the following two reasons: First, the response of temperature to climate in the study area is much stronger than that of precipitation, especially in cold places. The annual average temperature in our selected study area is −3.5 °C and the extreme temperature reaches −50.2 °C, so the forest NPP may be more sensitive to temperature. Second, it may be due to the mutual interference of different climatic factors and the joint action of multiple factors. As shown in Figure 10, the minimum value of forest NPP in the study area was in 2013, and the average temperature, average saturated vapor pressure and annual shortwave radiation were the minimum values in 7 years, while the annual precipitation was the highest value in 7 years. Through the joint action of multiple factors, the stomatal conductance of plants is reduced or even closed, and both plant transpiration and photosynthesis are significantly reduced, thereby reducing the accumulation of dry matter.

From 2012 to 2018, the annual mean value of forest NPP ranged from 26.83 g to 718.27 g/C/m²/year, which was positively and strongly correlated with VPD and T_mean: 0.82 and 0.81, respectively. The weak correlation with Sard was 0.18, and the correlation with Pre was negatively correlated at −0.66. In terms of spatial distribution, it is low in the middle and high on the four sides; in terms of temporal trend, the change trend is not obvious, although it shows slow growth. The results showed that the GPP values of different forest types were highly sensitive to C:N_leaf (C:N of leaf), SLA1 (canopy average specific leaf area in phenological phase 1), and FLNR (fraction of leaf N in Rubisco).

In summary, this study also has some shortcomings, which can reduce the uncertainty sources of model simulation from the following aspects: (1) For the Biome-BGCMuSo model, site-specific input parameter values can achieve the best simulation of forest carbon flux, but some input parameters cannot be observed (such as Lcel, Llab, etc.), and it is difficult to obtain the variation range of parameters. How to collect and determine the range of these parameters needs to be further considered and explored; (2) The Biome BGCMuSo model improves the ability to simulate the carbon and water cycles in managed ecosystems. The simulation process is based on grids and does not consider the flow of matter and energy between grids; (3) We only used the sample plot data and dendrochronological measurement data to verify the results of forest carbon flux simulation in the study area, with low time resolution; (4) Although we simulated forest carbon flux with satisfactory spatial and temporal resolution through the strategy of model–model coupling, we did not consider disturbances such as forest diseases and insect pests, forest land changes, and forest fires.

5. Conclusions

Based on the strategy of incorporating the optimized remote sensing model and the ecological process model, this study simulated the carbon flux of different forest types in the Greater Khingan Mountains forest region of Inner Mongolia from 2012 to 2018, and then analyzed the spatial and temporal dynamic changes in forest carbon flux over the past 7 years. The results showed that C:N_leaf, SLA1 and FLNR were the most sensitive parameters of DNF, DBF and MF. Using the incorporation of optimized MOD_17 model and Biome-BGCMuSo, the simulated forest carbon flux can effectively suppress the distortion simulation of different forest types, which is verified by dendrochronological measurements data, and the effect is good.