## 2 Materials and Methods  
### Data Sources  
Annual ET data used in this study were obtained from [PERSON] et al. (2018) [32], which estimated through machine learning approach (model tree ensemble, MTE) integrating eddy covariance ET and remote-sensing data from 1982 to 2015 at spatial resolution (hereafter MTE-ET, spatial resolution: 0.1\({}^{\circ}\)) [33; 34]. Previous studies have demonstrated the capacity of MTE-ET to predict ET, even in regions where training data are not enough [32; 35]. We used annual GLASS GPP product from 1982 to 2015 (hereafter GLASS-GPP, spatial resolution: 0.05\({}^{\circ}\)) produced by National Earth System Science Data Sharing Infrastructure, National Science and Technology Infrastructure of China ([[http://www.geodata.cn](http://www.geodata.cn)]([http://www.geodata.cn](http://www.geodata.cn))). It was estimated based on the Bayesian algorithm ensemble of eight light use efficiency models integrating the MODIS (Moderate Resolution Imaging Spectroradiometer) dataset, AVHRR (the Advanced Very High Resolution Radiometer) dataset and meteorological data [36]. The GLASS GPP was widely validated through 155 global eddy covariance measurements which contained nine terrestrial ecosystem types and applied in previous ecosystem carbon cycles studies at global and regional scales [36; 37].  
The widely used Standardized Precipitation-Evapotranspiration Dirac (SPEI, spatial resolution: 0.5\({}^{\circ}\)) developed by [PERSON] et al. (2010) was applied which combined both flexible time scales and the sensitivity to changes in evaporative demand compared with other drought indexes (e.g., SPI--Standardized Precipitation Index and PDSI--Palmer Drought Severity Index) [38; 39; 40]. The 12-month SPEIs were selected to define annual drought severity during the period of 1982-2015 which can be directly acquired from the SPEI base version 2.5 products ([[http://digital.csic.es/handle/10261/153475](http://digital.csic.es/handle/10261/153475)]([http://digital.csic.es/handle/10261/153475](http://digital.csic.es/handle/10261/153475))). Annual gridded precipitation data (spatial resolution: 0.1\({}^{\circ}\)) was downloaded from the National Meteorological Information Center of the China Meteorological Administration ([[http://data.cma.cn/](http://data.cma.cn/)]([http://data.cma.cn/](http://data.cma.cn/))) at 613 meteorological stations in China, which was interpolated to raster data in order to be consistent with the resolution of ET and GPP through the Inverse Distance Weight (IDW) method which was weighted average interpolation based on the distance between points. Land cover map of 2015 (spatial resolution: 300 m ) from the European Space Agency Climate Change Initiative ([[http://maps.elie.ucl.ac.be/CCI/](http://maps.elie.ucl.ac.be/CCI/)]([http://maps.elie.ucl.ac.be/CCI/](http://maps.elie.ucl.ac.be/CCI/))) was also aggregated to 0.1\({}^{\circ}\) by using the nearest neighbor method and also used to classify vegetation types (we reclassified the 38 land cover types into nine categories in this study, namely croplands--CRO, grasslands--GRA, evergreen broadleaf forests--EBF, deciduous broadleaf forests--DBF, evergreen needle-leaf forests--ENF, deciduous needle-leaf forests--DNF, mixed forests--MF, shrub lands--SHR and non-vegetation regions--NOV, Figure 1a) over China [41].  
In addition, the study was carried out at nine sub-regions according to the Chinese climatic and administrative boundaries [42; 43] and four eco-climatic zones according to multiyear average precipitation in China, as shown in Figure 1b. Moreover, GLASS-GPP and SPEI datasets were also uniformly aggregated to 0.1\({}^{\circ}\) using the nearest neighbor resampling method.