## 2 Materials and Methods
### 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 (\(\theta\)) given the existing data (D), denoted \(\text{P}(\theta\,\text{I}\,\text{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:  
\[\text{P}(\theta|\text{D})=\frac{\text{P}(\theta)\text{P}(\text{D}|\theta)}{ \text{P}(\text{D})} \tag{2}\]where P(D) is the probability of observed GPP and P(D | \(\theta\)) is the conditional probability density of observed GPP with prior knowledge, also called the likelihood function for parameter \(\theta\).  
Given a collection of \(N\) measurements, the likelihood function (L) can be expressed as:  
\[\mathrm{L}=\prod_{i=1}^{N}\frac{1}{\sqrt{2\pi}\sigma}e^{\frac{(X_{i}-p_{i})^{2} }{2\sigma}} \tag{3}\]  
where \(\sigma\) represents the standard deviation of the data-model error, \(X_{i}\) represents the \(i\)th of \(N\) measurements, and \(\mu_{i}\) 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(\(\theta\)) 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 (\(\epsilon_{\mathrm{max}}\) or LUE\({}_{\mathrm{max}}\)), temperature-constrained factors (T\({}_{\mathrm{min},\mathrm{min}}\), T\({}_{\mathrm{min},\mathrm{max}}\)) and the water-constrained factors (VPD\({}_{\mathrm{min}}\), VPD\({}_{\mathrm{max}}\)). The lower and upper bounds of \(\epsilon_{\mathrm{max}}\) (0.3-3.0 gC/m\({}^{2}\)/day/MJ APAR) were determined from the range of \(\epsilon_{\mathrm{max}}\) used in PEMs [9; 15; 53]. Following the related References [8; 15; 54], we specified the initial bounds of these parameters: T\({}_{\mathrm{min},\mathrm{min}}\) (\({}^{\circ}\)C), T\({}_{\mathrm{min},\mathrm{max}}\) (\({}^{\circ}\)C), VPD\({}_{\mathrm{min}}\) (Pa) and VPD\({}_{\mathrm{max}}\) (Pa) as [\(-\)35\(,-\)2], [6; 30], [60; 1000], and [1500,6500], respectively.