*2.1. Study Area*

The SNP is located in the central part of Northeast China, in a range of 121◦38 to 128◦33E, 42◦49 to 49◦12N, with a total area of 22.35 × 10<sup>4</sup> km<sup>2</sup> (Figure 1). It is an alluvial plain situated in the central Songliao Basin between the Xiaoxing'an and Changbai Mountains, through which the Songhua River and the Nenjiang River flow [10]. The SNP belongs to a temperate continental semi-humid and semi-arid monsoon climate zone, characterized by four seasons with a hot, rainy summer and a cold, dry winter and with significant windy days. The annual precipitation is between 350 and 800 mm [11]. Soils are fertile with chernozem, meadow soil, and black soil widely distributed. From west to east, the natural ecosystem type is typical grassland, meadow steppe and forest steppe, respectively. It is an important ecological protection barrier in the Northeast.

**Figure 1.** Study Area. FL: farmland; DBF: deciduous broad-leaved forest; MF: mixed forest; GRA: grassland and WET: wetland.

#### *2.2. Data and Processing*

In this study, CUE at the monthly and seasonal scale derived from MODIS data products and the ancillary data were used to explore the spatial-temporal variations of CUE of natural ecosystems and their responses to climate and LSP changes. The main steps are as follows: (1) estimating NPP at monthly scale by the CASA (Carnegie-Ames-Stanford approach) model; (2) calculating monthly CUE of different ecosystems and performing the trend analysis; (3) extracting LSP metrics and analyzing the effects of phenology and climatic factors on the variations of ecosystem CUE by the correlation and partial correlation analysis methods. Figure 2 illustrates the technical approach of this study.

**Figure 2.** The flow chart of the research approach.

#### MODIS and Meteorological Data

The MOD17A2H version 6 GPP product is a cumulative 8-day composite of values based on the radiation use efficiency concept, which can be potentially used as an input to data models to compute energy, carbon, water cycle processes and biogeochemistry of vegetation [12]. Monthly Normalized Difference Vegetation Index (NDVI) spatial distribution data set was obtained from the MOD13A2. For the calculation of NPP, the Maximum Value Composite (MVC) was used to synthesize the monthly NDVI data. The phenological parameter extraction was based on 16-day data. Each time period of the MODIS data included 6 images.

The MODIS land cover type dataset (MCD12Q1) was downloaded from the Land Processes Distributed Active Archive Center (LP DAAC) (https://lpdaac.usgs.gov/) [13]. In this study, land cover data from 2001 to 2015 were adopted. According to the International Geosphere-Biosphere Program

(IGBP) class scheme and the regional characteristics of the SNP, natural ecosystem included the following four types, deciduous broad-leaved forest (DBF), mixed forest (MF), grassland (GRA) and wetland (WET).

The total monthly radiation was obtained from 7 radiation stations (Figure 1). The remaining meteorological data of 29 meteorological stations within and near the SNP from 2001 to 2015 were acquired from the National Meteorological Information Center (http://data.cma.cn/). Most studies sugges<sup>t</sup> that precipitation and temperature be the main meteorological factors affecting CUE [3,13]. In this study, monthly precipitation (mm) and temperature (◦C) were computed from May to November in order to explore the relationship between CUE and climate factors in growing seasons. MODIS, meteorological and other data were all resampled at a resolution of 1 km. Table 1 summarizes and describes the characteristics of the data and sources.



#### *2.3. Estimating NPP with the CASA Model*

We estimated the monthly vegetation NPP of the SNP based on the CASA model, which considers the physiological and ecological characteristics of vegetation and the environmental conditions related to growth [14]. Vegetation NPP estimation was derived by using vegetation cover type, NDVI, monthly average temperature, total precipitation and solar radiation [15]. The basic calculation formula of the CASA model is as follows [21]:

$$\text{NPP}(\mathbf{x}, \mathbf{t}) = \text{SOL}(\mathbf{x}, \mathbf{t}) \times \text{FPAR}(\mathbf{x}, \mathbf{t}) \times 0.5 \times \text{T}\_{\varepsilon 1}(\mathbf{x}, \mathbf{t}) \times \text{T}\_{\varepsilon 2}(\mathbf{x}, \mathbf{t}) \times \mathcal{W}\_{\varepsilon}(\mathbf{x}, \mathbf{t}) \times \varepsilon\_{\text{max}} \tag{1}$$

where SOL(x,t) is the total solar radiation at pixel x for month t. FPAR(x,t) is the fraction of photosynthetically active radiation absorbed by vegetation. 0.5 indicates the proportion of solar active radiation (0.4–0.7 μm) that can be utilized by vegetation to the total solar radiation. Tε1(x,t) and Tε2(x,t) represent temperature stress coefficients, Wε(x,t) is the coefficient of water stress, and εmax is the maximum light use efficiency under ideal conditions [22]. Comparison between the calculated NPP and the reported study conducted in Northeast China and the western part of Jilin Province [23,24] as the cross checking and validation of the analysis. Mao et al. [23] verified the NPP by comparing the simulated value with the flux observation data. The simulated value was close to the measured value, and the error was within 25%.

#### *2.4. Calculation of CUE*

The CUE of ecosystem describes the relationship between photosynthesis and respiration, which is an important indicator of the ability of plants to transfer carbon [16]. As one of the key controlling factors of ecosystem carbon storage, CUE is defined as follows [19]:

$$\text{CUE} = \frac{\text{NPP}}{\text{GPP}} \tag{2}$$

where GPP represents the ability to capture energy and carbon through photosynthesis plants and the total amount of carbon assimilation. NPP reveals the energy of plants stored after losing carbon from GPP through autotrophic respiration [17]. The higher CUE means the greater the proportion of GPP kept by ecosystems after self-consumption. However, uncertainty issues have been recognized by studies using public-domain data, e.g., with respect to water use efficiency (WUE) [25].

#### *2.5. Extraction of Land Surface Phenology Metrics*

We used the dynamic threshold method to extract metrics of LSP. The polynomial method was used to fit and reconstruct the NDVI time series data from 2000 to 2016. The software TIMESAT with a seasonal parameter of 0.5, an adaptation strength of 2.0, a Savitzky–Golay window size of 2, and an amplitude of 20% was run in MATLAB R2015b (The Mathworks, Inc., Natick, MA, USA). The parameters were set according to Qi et al. [18]. The start of growing season (SOS) and the end of season (EOS) for each year were calculated, and the length of season (LOS) was obtained as the difference between SOS and EOS in each grid.
