## 2 Materials and Methods
### Datasets  
Site-specific CO\({}_{2}\) flux and climate data measurements were conducted by the EC flux tower, which was installed at a flat site in the alpine wetland in October 2007 to continuously measure the net ecosystem CO\({}_{2}\) exchange (NEE) between the vegetation and atmosphere. The footprint of this flux tower is approximately 200-250 m. A three-axis sonic anemometer with a fast response open path infrared CO\({}_{2}\)/H\({}_{2}\)O gas analyzer (IRGA, LI 7500, LI-COR Inc., Lincoln, NE, USA) was  
Figure 1: Location of the study site in the Zoige alpine wetland, China. installed at a height of 2.2 m above ground level to measure the fluctuations in three wind components (\(w\), \(u\), \(v\)) and CO\({}_{2}\) fluxes. Data quality control was also implemented to reduce measurement-induced uncertainties. Daily GPP and NEE data during 2008-2009 were then used to analyze the dynamics of wetland phenology. The averaged eight-day GPP and climate data observed in 2008 and 2009 were utilized to support model simulations and validate the results of the VPM model. The daily weather dataset for the study site location was obtained from the climate database of the Zoige meteorological station (2000-2007, 2010-2011). The meteorological station is located about 15 km from the EC tower. More details on the instruments and the data quality control process are described in [PERSON] _et al._ and [PERSON] _et al._ (2019); [PERSON] _et al._ (2019). Site-specific vegetation indices were derived from MODIS datasets. We downloaded the eight-day composite MOD09A1 datasets for the 2000-2011 period from the Oak Ridge National Laboratory's Distributed Active Archive Center (DAAC) website ([PERSON] _et al._, 2001). Then, reflectance values of four spectral bands (blue (459-479 nm), red (620-670 nm), near infrared (NIR 841-876), and shortwave infrared (SWIR1628-1652) over 2000-2011 were used to calculate three vegetation indices: normalized difference vegetation index (NDVI), enhanced vegetation index (EVI), and land surface water index (LSWI) ([PERSON] _et al._, 2017). \[NDVI=\frac{\rho_{nir}-\rho_{red}}{\rho_{nir}+\rho_{red}} \tag{1}\]  
\[EVI=G\times\frac{\rho_{nir}-\rho_{red}}{\rho_{nir}+(C_{1}\times\rho_{red}-C_{2 }\times\rho_{blue})+L} \tag{2}\]  
\[LSWI=\frac{\rho_{nir}-\rho_{swir}}{\rho_{nir}+\rho_{swir}} \tag{3}\]  
where \(G=2.5\), \(C_{1}=6\), \(C_{2}=\)7.5, \(L=1\), and \(\rho_{nir}\), \(\rho_{red}\), \(\rho_{blue}\), and \(\rho_{swir}\) are the reflectance of NIR, red, blue, and SWIR bands, respectively ([PERSON] _et al._, 2017; [PERSON] _et al._, 2019; [PERSON] _et al._, 2019; [PERSON] _et al._, 2019; [PERSON] _et al._, 2019). The three vegetation indices calculated from surface reflectance were affected by noise from clouds and cloud shadows ([PERSON] _et al._, 2019). Therefore, we used a simple gap filling method to gap-fill vegetation index values for which cloudy pixels existed in the time series ([PERSON] _et al._, 2019). Based on the two-step gap-filling procedure reported in earlier studies ([PERSON] _et al._, 2019), we first selected a three-point time-series filter (X(t \(-\) 1), X(t), and X(t + 1)) and used values of non-cloudy pixels in this window to replace cloudy pixels. If both the X(t \(-\) 1) and X(t + 1) pixels were reliable and X(t) was unreliable, the average value of X(t \(-\) 1) and X(t + 1) was used to replace X(t). If only one pixel (either X(t \(-\) 1) or X(t + 1)) was cloud-free and reliable and X(t) was unreliable, we used that pixel value to replace X(t). If the algorithm did not result in a three-point time series filter, we then extended it to a five-point time series filter, X(t \(-\) 2), X(t \(-\) 1), X(t \(+\) 1), X(t \(+\) 2), using the same procedure as the above three-point time series filter ([PERSON] _et al._, 2019).