## 2 Data and Methods
### Methods
#### 2.2.2 Parameter Optimization  
The derived GPP in 85 calibration site-years was used to optimize parameters in MOD17, TL-LUEn, and TL-LUE models. The optimization was implemented using the Markov chain Monte Carlo (MCMC) method [97, 98]. In the optimization, three models were driven using the same locally measured meteorological data and smoothed MODIS LAI. The MCMC simulation, as a stochastic simulation method, is based on Bayesian Theory, in which parameters are random variables instead of deterministic, but unknown constants in the classic thoughts. The fundamental formula of Bayesian Theory is:  
\[\pi\left(\theta\,|\,x\right)=\frac{p(x\,|\,\theta)\pi(\theta)}{\int\limits_{- \infty}^{\infty}p(x\,|\,\theta)\pi(\theta)d\theta} \tag{9}\]  
where \(\pi(\theta\,|\,x)\) is the posterior density of parameter \(\theta\) (a term distribution under the condition of given sample \(x\)); \(\pi(\theta)\) is the prior distribution of parameter \(\theta\) (the knowledge possessed before measurement); and \(p(x\,|\,\theta)\) is a likelihood function. To determine the posterior density \(\pi(\theta|x)\), the prior density and the likelihood function should be given in advance. We specified the prior density function as a uniform distribution over the following ranges:  
\[0<\varepsilon_{m}<20\] \[0<\beta<1500 \tag{10}\] \[0<\varepsilon_{m},\varepsilon_{m},\varepsilon_{m},\varepsilon_{ m}<10\]  
The lower and upper limits of \(\varepsilon_{m}\) (unit: g C MJ\({}^{-1}\)) and \(\beta\) (unit: \(\mu\)g C m\({}^{-2}\) s\({}^{-1}\)) were set according to the values compiled from 100 published datasets by [PERSON] _et al._[99]. The upper limits of \(\varepsilon_{\text{msu}}\), \(\varepsilon_{\text{msh}}\) and \(\varepsilon_{\text{max}}\) (unit: g C MJ\({}^{-1}\)) were assigned on the basis of previous findings [32, 100]. In the optimization, these parameters were assumed uniformly distributed in the given limits with equal probability for all possible values. The likelihood function was specified according to the distribution of simulation errors, which were assumed following a multivariate Gaussian distribution with a zero mean. This assumption is commonly made in many studies [101, 102, 103]. With this assumption, the likelihood function can be written as:  
\[p\left(x\,|\,\theta\right)=\prod_{i=1}^{n}\frac{1}{\sqrt{2\pi\sigma_{i}}}e^{- \left(P_{i}-Q_{i}\right)^{2}/2\sigma_{i}^{2}} \tag{11}\]  
where \(O_{i}\) and \(P_{i}\) are the tower-derived GPP and simulated GPP, respectively; \(\sigma_{i}\) are the standard error of tower-derived GPP. The sampling of parameters was implemented using the Metropolis-Hastings (M-H) algorithm. To find an effective proposal distribution, we first made a test run of the algorithm with 50,000 simulations. Based on the test run, a Gaussian distribution \(N(0,\,\text{cov}^{0}(\theta))\) was constructed (cov\({}^{0}(\theta)\) is a diagonal matrix with its diagonal elements equal to the estimated variances of parameters \(\theta\)). Then, the following proposal distribution was adopted to execute the consecutive MCMC simulations formally for 30,000 times:  
\[\theta^{k}=\theta^{k-1}+N(0,\text{cov}^{0}(\theta)) \tag{12}\]  
where \(\theta^{k}\) is the new parameters generated from its predecessor \(\theta^{k-1}\). The running means and standard deviations of parameter samples need time to approach stable. For bettering the statistical analysis of parameters, we discarded the initial 10,000 samples in the burn-in period and only used the remaining 20,000 samples for further analysis of each parameter. The histograms of the samples for each parameter indicate these parameters were well constrained in most situations because the posterior density functions were near the normal distribution (see Figure A1 in Appendix). Uncertainties of estimated parameters were quantified with a 95% highest-probability density interval.