*2.3. Description of MOD17A2H Algorithm*

The MOD17A2H algorithm is based on light-use efficiency (LUE) approach [49,50], which provides global GPP estimates of 8 day temporal and 500 m spatial resolution [15]. The MODIS GPP product is calculated from the following equation:

$$\text{GPP} = \varepsilon\_{\text{max}} \times 0.45 \times \text{SW}\_{\text{rad}} \times \text{FPAR} \times \text{f}(\text{T}\_{\text{min}}) \times \text{f}(\text{VPD}) \tag{1}$$

where εmax is the maximum LUE obtained from the Biome-specified Parameters Look Up Table (BPLUT) on the basis of vegetation type. The BPLUT contains values specifying minimum temperature and VPD limits, specific leaf area and respiration coefficients for the standard land cover classes [48]. SWrad is shortwave solar radiation of which 45% is photosynthetically active radiation (PAR), FPAR is the fraction of PAR absorbed by vegetation and the scale factors f(Tmin) and f(VPD) reduce εmax under unfavorable conditions of low temperature and high VPD. The forcing data such as SWrad, Tmin and VPD in the MOD17A2H GPP product were implemented by the Global Modeling and Assimilation Office (GMAO) Reanalysis data. The MODIS GPP algorithm is described in detail in previous literature [12,15,18].

#### *2.4. Parameter Optimization and Uncertainty Analysis*

The current MOD17 BPLUT is too general for local regional application [20]. The same set of parameters was applied indiscriminately to diverse types of the same ecosystems, introducing large uncertainties for the simulation of GPP in the arid region. To improve the accuracy of the GPP estimation in desert–oasis–alpine ecosystems in the arid region, we calibrated the parameters of the MOD17 model based on in situ flux tower observations using Bayesian model-data fusion approach. The model parameters were calibrated against GPP time series from the flux tower measurement network through a Bayesian data model synthesis [33,38]. According to Bayesian theory, posterior probability density functions (PDFs) of model parameters (θ) given the existing data (D), denoted P(θ|D), can be obtained from prior knowledge of the parameters and information generated by comparison of simulated and observed variables, and can be described as:

$$\mathbf{P}(\boldsymbol{\theta}|\mathcal{D}) = \frac{\mathbf{P}(\boldsymbol{\theta})\mathbf{P}(\mathcal{D}|\boldsymbol{\theta})}{\mathbf{P}(\mathcal{D})} \tag{2}$$

where P(D) is the probability of observed GPP and P(D|θ) is the conditional probability density of observed GPP with prior knowledge, also called the likelihood function for parameter θ.

Given a collection of *N* measurements, the likelihood function (L) can be expressed as:

$$\mathcal{L} = \prod\_{i=1}^{N} \frac{1}{\sqrt{2\pi}\sigma} e^{\frac{(\overline{\chi\_i} - \mu\_i)^2}{2\sigma}} \tag{3}$$

where *σ* represents the standard deviation of the data-model error, *Xi* represents the *i*th of *N* measurements, and *μ<sup>i</sup>* is the model-derived estimates of a measurement.

In our study, we assumed the parameter priors are uniform, and the posterior PDFs for the model parameters were generated from prior PDFs P(θ) with observation data by a Markov chain Monte Carlo (MCMC) sampling technique [33]. Herein, the Metropolis–Hasting algorithm [51,52] was adopted to generate a representative sample of parameter vectors from the posterior distribution. We ran the MCMC chains with 50,000 iterations each, and regarded the first 15,000 iterations as the burn-in period for each MCMC run. All accepted samples from the runs after burn-in periods were used to compute the posterior parameter statistics of the models.

In this study, the MOD17A2H GPP algorithm contains 5 parameters: Maximum light-use efficiency (εmax or LUEmax), temperature-constrained factors (Tmin\_min, Tmin\_max,) and the water-constrained factors (VPDmin, VPDmax). The lower and upper bounds of εmax (0.3–3.0 gC/m2/day/MJ APAR) were determined from the range of εmax used in PEMs [9,15,53]. Following the related References [8,15,50,54], we specified the initial bounds of these parameters: Tmin\_min ( ◦C), Tmin\_max ( ◦C), VPDmin (Pa) and VPDmax (Pa) as [−35,−2], [6,30], [60,1000], and [1500,6500], respectively.
