Remote Sensing Greenness and Urbanization in Ecohydrological Model Analysis: Asia and Australasia (1982–2015)

Linking remote sensing information and ecohydrological models to improve understanding of terrestrial biosphere responses to climate and land use change has become the subject of increased interest due to the impacts of current global changes and the effect on the sustainability of human lifestyles. An application to Asia and Australasia (1982–2015) is presented, revealing the following results: (i) The broad distribution of regions with the enhanced vegetation greenness only follows the general pattern as for the whole, without obvious dependence on regional or climate fluxes ratios. That indicates a prevailing increasing greenness over land due to both the impacts of current global changes and the sustainability of human lifestyles; (ii) regions with vegetation greenness reduction reveal a unique distribution, concentrating in the water-limited domain due to the impacts of external (climatically “dry gets drier and wet gets wetter”) and internal (anthropogenically increased evaporation) changes; (iii) the external changes of dryness diverge at the boundary separating energy from water-limited regimes, and the internal changes indicate large-scale afforestation and deforestation) that occur mainly in China and Russia due to a conservation program and illegal logging, respectively, and a massive conversion of tropical forest to industrial tree plantations in Southeast Asia, leading to an increased evaporation.


Introduction
Changes in vegetation greenness have been reported at regional and continental scales on the basis of forest inventory and satellite measurements [1][2][3][4]. Long-term changes in vegetation greenness are driven by multiple interacting direct factors (human land-use management) and indirect factors (such as climate change, CO 2 fertilization, nitrogen deposition, and recovery from natural disturbances) [5][6][7][8][9][10].
Recent satellite data reveal a greening pattern induced by direct factors (human land use), which enhances the vegetation greenness as it is strikingly prominent in China and India, notwithstanding the impact of indirect factors [11]. An enhanced vegetation growth has broad implications for surface energy, water and carbon budgets, and ecosystem services across multiple scales [7,[12][13][14]. Linking remote sensing and climate and ecosystem models to provide improved understanding of terrestrial biosphere responses to climate and land use change has become the subject of increased interest due to the impacts of current global changes and sustainability of human lifestyles [2,11,15].
The reliable detection and attribution of changes in vegetation growth are prerequisites for the development of strategies for the sustainable management of ecosystems [3]. However, human exploitation of land will remain a complex dynamic endeavour, and the underlying mechanisms are not yet fully understood [5,11]. To quantify impacts of climate change and anthropogenic activities on land surface dynamics, an ecohydrological modelling approach has been introduced [16,17], which is based on two nondimensional flux ratios as ecohydrological variables. They separate energy and water supply to distinguish climate from land use-induced change effects.
The application of the ecohydrological conceptual model has been proved and expanded to regional scales by combining remote sensing vegetation information and reanalysis data in an ecohydrological state space, spanned by relative excess energy U and excess water W [18,19]. Here, the diagnostic is used, combined with nighttime lighted city distribution (urban) to answer the following questions: Do the increased/decreased regions of greenness show a general pattern? Are there any differences between the whole study regions, significant ecohydrologically changed regions (significant U,W changes), and urban areas? How do climate changes and anthropogenic activities contribute to those regions? The aim of the following analysis is to attribute and quantify these changes to climatic or human-induced causes across Asia and Australasia. Our evaluation comprises ecohydrological surface information of the changing rainfall-runoff chain jointly with remote sensing observations of vegetation and urbanization, which, to our knowledge, has not been performed.

Data and Method
Data on climate, vegetation greenness, and urbanization were provided by ERA-Interim, GIMMS NDVI (1982-2015 averages, see [20,21]) and DMSP/OLS nighttime light (2012, see [22,23]). The spatial resolutions of the ERA-Interim and GIMMS vegetation greenness were 0.75 • by 0.75 • and 8 km by 8 km, respectively. The global DMSP/OLS nighttime stable light from NOAA National Geophysical Data Center provided the information for urbanization, given the brightness range 0 < DN < 63 and using the spatial resolution of about 1 km or 30 arc-seconds. The spatial resolutions were resampled according to the relatively coarse resolution of ERA-Interim into 0.75 • by 0.75 • by bilinear resampling.

Data Preprocessing
Unlike traditional sociodemographic research using administrative boundaries, population, size or density, and economic indicators to define a city, remote sensing techniques are employed in this study to analyze continuous spatial variation in nighttime light intensity of development or degree of modification. To reduce the effects of overglow [22,24], only spatially contiguous lighted pixel units (subscript "i") with DN i ≥ 12 are identified as cities (for detailed descriptions of data quality control and threshold selection, see Small et al. [24]). The area of each pixel within spatially contiguous lighted areas, a(j), is defined as a i , where j corresponds to the sequence of spatially contiguous areas exceeding the brightness threshold, and i corresponds to the number of pixels in each of these areas. Thus, a city (numbered by "j") size is characterized by size (j) = i a i ( j).

Ecohydrological State Space Analysis
Ecohydrological analysis is a physical approach to diagnose the controlling factors of the rainfall-runoff chain on watershed scale [16,17]. Water and energy are supplied by precipitation P and net radiation N or potential evapotranspiration. In a climatological mean, these fluxes are balanced by the partitioning of evapotranspiration E and runoff Ro and by evapotranspiration E and sensible heat flux H, respectively. The energy flux units are in water flux equivalents, m yr −1 .
Energy and water balance equations are supplemented by an empirical equation of state, which, in the framework of Budyko (1974 [25]), characterizes the rainfall-runoff chain. Here, Schreiber's equation ( [26]; for a stochastic interpretation, see [27]) was used.
For ecohydrological analysis, water and energy fluxes are separated by introducing relative excess energy U and relative excess water W. These flux ratios characterize energy and water fluxes, which are unused by the ecosystem (and, therefore, in excess); that is, excess energy is available for atmospheric heating, and excess water is available for geomorphological formation [28].
The flux ratio of the net radiation and precipitation yields Budyko [25] dryness index: Surface climates are conveniently projected into the (U,W) space as an ecohydrological state space spanned by flux ratios to separate energy from water flux-related excesses ( Figure 1). The dryness index D is displayed in the (U,W) space, separating energy (0 < D < 1) from water (D > 1) limited regimes. In addition, biomes can be suitably identified ranging from tundra (D < 0.3) via forests (0.3 < D < 1) and steppe/savanna (1 < D < 2) to semidesert (2 < D < 3) and desert (D > 3) [15,25]. Energy and water balance equations are supplemented by an empirical equation of state, which, in the framework of Budyko (1974 [25]), characterizes the rainfall-runoff chain. Here, Schreiber's equation ( [26]; for a stochastic interpretation, see [27]) was used.
For ecohydrological analysis, water and energy fluxes are separated by introducing relative excess energy U and relative excess water W. These flux ratios characterize energy and water fluxes, which are unused by the ecosystem (and, therefore, in excess); that is, excess energy is available for atmospheric heating, and excess water is available for geomorphological formation [28].
The flux ratio of the net radiation and precipitation yields Budyko [25] dryness index: Surface climates are conveniently projected into the (U,W) space as an ecohydrological state space spanned by flux ratios to separate energy from water flux-related excesses ( Figure 1). The dryness index D is displayed in the (U,W) space, separating energy (0 < D < 1) from water (D > 1) limited regimes. In addition, biomes can be suitably identified ranging from tundra (D < 0.3) via forests (0.3 < D < 1) and steppe/savanna (1 < D < 2) to semidesert (2 < D < 3) and desert (D > 3) [15,25].  [29,30]). Squares and circles represent the ecohydrological reference, states (squares) representing the first period followed by subsequent second period (circles). Directions and lengths of the vectors (arrows) connecting first and second period determine quality and magnitudes of the causes of change (see text), as described by squares of different colors. Lines of constant dryness or net radiation vs. precipitation and the graph of an ideal rainfall-runoff chain, U = 1 + (1 − W)/log(W) (see [29,30]). Squares and circles represent the ecohydrological reference, states (squares) representing the first period followed by subsequent second period (circles). Directions and lengths of the vectors (arrows) connecting first and second period determine quality and magnitudes of the causes of change (see text), as described by squares of different colors. In the (U,W) space, ecohydrological changes are visualized as pieces of trajectories (or vectors, see Figure 1) whose origin and end characterize the (U,W) space representing the first and the subsequent second averaging period. The causes of change can be attributed to a forcing induced by contributions from external or climate and from internal or anthropogenic origin [16][17][18]: (1) External (climate-induced) processes change the climate forcing, which can be visualized as the vector components in the direction of the negative diagonal towards the second yellow or the fourth dark blue quadrant. For example, given E = constant, for increasing (decreasing) aridity dD > 0 (< 0), one obtains an increasing (decreasing) energy excess U = 1 − E/N and decreasing (increasing) water excess, W = 1 − E/P; with magnitudes given by the inverse aridity ratio (Equation (7)). The (inverse) initial slope at (U 1 ,W 1 ) point (corresponding to the D = 1 − line) defines the attribution coordinates. Vector components of ecohydrological states changing from a first to a second period, which are aligned along D-lines (perpendicular to the D-lines), are attributed to internal anthropogenic control (external climate forcing). Related orthogonalization provides quantitative measures of the external and internal causes of change. In this sense, the attribution of the causes of change is quantified by the lengths of the vector components projected onto (perpendicular to) the dryness or D = 1 − line ending at (U 1 ,W 1 ) (for detailed calculation, see [18]).

Application of Ecohydrological Analysis: Asia and Australasia
Changes of the surface energy and water fluxes, which represent the changing rainfall-runoff chain, are evaluated in the following by employing the attribution analysis to classify and quantify the causes of change observed in Asia and Australasia. Depending on vegetation greenness increase or decrease, two domains were selected for comparative analysis: (i) The continental domain covering the administrative region of Asia and Australasia with significant ecohydrological (U,W) change, and (ii) the embedded smaller scale regions affected by prominent anthropogenic activity (that is, cities or urbanized areas). The periods of 1982-1998 and 1999-2015 are suggested by the climate warming tendency showing high and low values during the first and second periods, respectively [30,31].

Land Surface States in (U,W) Space: Means and Changes
The climate mean state is presented in (U,W) space comprising the geographical distribution of relative excess energy U and excess water W in terms of frequencies of pixels, each of which is associated with a (U,W) pair. Thus, the physical states in the (U,W) space were characterized by a number (density) of area units, associated with long-term (U,W) climate means (1982-2015) obtained from the ERA-Interim dataset: (1) Administrative Asia and Australasia count 9782 pixels on the ERA-interim scale (Figure 2a  Vegetation increase (Figure 2, middle row): The frequency distributions of the three land surface states associated with regions where vegetation greenness increased in the second period followed a general pattern as the whole regions ( Figure 2, upper row) showed no obvious dependence on regional or climate fluxes. That is in agreement with previous research, which had shown a prevailing vegetation growth over the Earth's lands since the early 1980s when satellite data became available [8,11,29,32,33]. This indicates that regions of increasing vegetation greenness were affected by similar positive climate or human impact on greenness. states associated with regions where vegetation greenness increased in the second period followed a general pattern as the whole regions ( Figure 2, upper row) showed no obvious dependence on regional or climate fluxes. That is in agreement with previous research, which had shown a prevailing vegetation growth over the Earth's lands since the early 1980s when satellite data became available [8,11,29,32,33]. This indicates that regions of increasing vegetation greenness were affected by similar positive climate or human impact on greenness.
(iii) Vegetation decrease (Figure 2, bottom row): Unlike regions of increased vegetation greenness (Figure 2, middle row), climate change and anthropogenic activity contributed negatively to the three land surface modes associated with vegetation greenness reduction, and they revealed unique distribution patterns: Asia and Australasia were characterized by a trimodal distribution with one mode aligned along D = 1 and the other two modes separated by dryness differentiation (1 < D < 2 and D > 3). The dominant mode of Asian and Australasian greenness decrease was concentrated in the water-limited domain, and a similar distribution could be found in the nighttime lighted cities. That is, cities located in the water-limited domain were more likely to show decreasing vegetation greenness. The mass of significant (U,W) change was concentrated in the water-limited domain, but tended to move closer to the Schreiber (Equation (6)) curve.
In general, the pixel distributions of the three climate mean states in (U,W) space ( Figure 2, upper row) showed structurally similar patterns without considering vegetation greenness change and for increased greenness regions (Figure 2, middle row). However, decreasing vegetation greenness regions (Figure 2, bottom row) showed different energy-water dependence and urbanization dependence as noted in the above. Thus, it is reasonable to attribute state change trajectories in (U,W) space to external (or climate) and internal (or anthropogenic) causes over the decreasing vegetation greenness regions and diagnose the underlying negative climate vegetation and urbanization relations.

Attribution of Change: Regions of Decreased Greenness
Regions of greenness change associated with three land surface states, (i) administrative climate mean, (ii) ecohydrological significant change, and (iii) nighttime lighted urbanization, being attributed to external and internal causes are partitioned into four quadrants and first displayed geographically and statistically (Figure 3). Then, changes in (U,W) space are represented by pieces of trajectories to highlight and visualize vegetation greenness decreasing from the first period to the second period ( Figure 4).  (Figure 2, middle row), climate change and anthropogenic activity contributed negatively to the three land surface modes associated with vegetation greenness reduction, and they revealed unique distribution patterns: Asia and Australasia were characterized by a trimodal distribution with one mode aligned along D = 1 and the other two modes separated by dryness differentiation (1 < D < 2 and D > 3). The dominant mode of Asian and Australasian greenness decrease was concentrated in the water-limited domain, and a similar distribution could be found in the nighttime lighted cities. That is, cities located in the water-limited domain were more likely to show decreasing vegetation greenness. The mass of significant (U,W) change was concentrated in the water-limited domain, but tended to move closer to the Schreiber (Equation (6)) curve.
In general, the pixel distributions of the three climate mean states in (U,W) space (Figure 2, upper row) showed structurally similar patterns without considering vegetation greenness change and for increased greenness regions (Figure 2, middle row). However, decreasing vegetation greenness regions (Figure 2, bottom row) showed different energy-water dependence and urbanization dependence as noted in the above. Thus, it is reasonable to attribute state change trajectories in (U,W) space to external (or climate) and internal (or anthropogenic) causes over the decreasing vegetation greenness regions and diagnose the underlying negative climate vegetation and urbanization relations.

Attribution of Change: Regions of Decreased Greenness
Regions of greenness change associated with three land surface states, (i) administrative climate mean, (ii) ecohydrological significant change, and (iii) nighttime lighted urbanization, being attributed to external and internal causes are partitioned into four quadrants and first displayed geographically and statistically ( Figure 3). Then, changes in (U,W) space are represented by pieces of trajectories to highlight and visualize vegetation greenness decreasing from the first period to the second period (Figure 4).    Attribution of external and internal causes ( Figure 3): Changes in the three land surface states show both a similar percentage of area cover affected by increased versus decreased vegetation greenness (around 70% versus 30%) and regions predominantly affected by external processes (between 70% and 80%). Of the external controlling regions, the fourth dark blue quadrant and second orange quadrant represent decreasing and increasing aridity associated with increasing and decreasing W and decreasing and increasing U, respectively. Thus, geographically, most regions of Russia, from Southwest China to India, the islands in Southeast Asia, and Western Australia indicate a decreasing aridity.
Of the internally affected regions, large-scale afforestation (third light blue quadrant) occurs mainly in China, which is in agreement with previous research about the largest net gain of forest [11,34]. Deforestation (first pink quadrant) could be detected by the ecohydrological model in Russia where illegal logging is a serious issue, especially in remote areas. The extent of illegal logging is estimated to affect between 5 and 30 percent of the boreal forest [34,35]. For Southeast Asia, with a massive conversion of tropical forest to industrial tree plantations [36], this attribution analysis identifies few pixels in pink as deforestation (but more in light blue as afforestation). The possible reason is that industrial tree plantations evaporate more than tropical forests (see also Röll et al. [37]). Attribution of external and internal causes ( Figure 3): Changes in the three land surface states show both a similar percentage of area cover affected by increased versus decreased vegetation greenness (around 70% versus 30%) and regions predominantly affected by external processes (between 70% and 80%). Of the external controlling regions, the fourth dark blue quadrant and second orange quadrant represent decreasing and increasing aridity associated with increasing and decreasing W and decreasing and increasing U, respectively. Thus, geographically, most regions of Russia, from Southwest China to India, the islands in Southeast Asia, and Western Australia indicate a decreasing aridity.
Of the internally affected regions, large-scale afforestation (third light blue quadrant) occurs mainly in China, which is in agreement with previous research about the largest net gain of forest [11,34]. Deforestation (first pink quadrant) could be detected by the ecohydrological model in Russia where illegal logging is a serious issue, especially in remote areas. The extent of illegal logging is estimated to affect between 5 and 30 percent of the boreal forest [34,35]. For Southeast Asia, with a massive conversion of tropical forest to industrial tree plantations [36], this attribution analysis identifies few pixels in pink as deforestation (but more in light blue as afforestation). The possible reason is that industrial tree plantations evaporate more than tropical forests (see also Röll et al. [37]).   (Figure 4) and the light blue quadrants (Figure 3), which both characterize the conversion of forest to industrial tree plantations and the greening of cities; (ii) Climate (external) induced changes: The regions of significant (U,W) change revealed a general pattern of dry remaining dry or getting dryer and wet remaining wet or getting wetter, which is accompanied with a dryness change by crossing the (D = 1) threshold of energy to water-limited regimes. Associated with the changing density distribution in Figure 2f, it appears that the peak in the water-limited region is enhanced, and the less obvious peak aligned along D = 1 will move towards the energy-limited region. Unlike regions with significant (U,W) changes, the spatially contiguous lighted cities did not show obvious general patterns.
Attribution analysis employed for the 30% of the area of Asia and Australasia where vegetation greenness has changed provided the following results for decreasing vegetation: (i) Significant anthropogenic or internally induced land surface changes (light blue arrows, Figure 4a) were caused by large-scale afforestation and deforestation mainly in China and Russia supported by a conservation program and illegal logging, respectively, while the tropical Southeast Asia was subjected to a massive conversion of forestation of the industrial tree plantation. All of these changes were associated with increasing evaporation; (ii) climate-controlled or externally changed (Figure 4b) areas show that the yellow and dark blue arrows were tending to regimes of enhanced dryness or wetness, respectively. Furthermore, an excess water ration of W = 1 − E/P − 0.4 may be identified as a critical threshold from which the direction of climate change diverges, tending towards drier/wetter regimes (characterized by larger/smaller dryness D or less/more excess water, respectively). In this sense, climate change forces regions from a dry to a drier and from a wet to a wetter regime.

Conclusions
An attribution analysis of the causes of ecohydrological change was employed to determine (1) the climate-induced external impact in terms of water and energy supply as external forcing, which is distinguished from (2) the internal processes due to anthropogenic influence modifying the partitioning of the surface fluxes compensating the climate-induced forcing [18].
This study applied the ecohydrological model analysis in relation to vegetation greenness change and urbanization characterized by nighttime lighted cities. This diagnostic reveals how regions suffering vegetation greenness reduction show different patterns of energy-water dependence. The results of the attribution analysis of change, which are conditional to vegetation greenness change and urbanization, support a shift from traditional management of resources and risks to an integrated monitoring and holistic response across wide ranges of space-time scales by linking remote sensing, climate, and ecosystem model analysis. The following results are presented as samples for possible management: (1) Regions where vegetation greenness increased in the second period indicate similar distribution in (U,W) space (as for the whole region). That is, a prevailing increase of vegetation growth could be found over the whole of Asia and Australasia since the early 1980s and did not show regional or climate flux dependence. However, regions of decreasing vegetation greenness were concentrated in the water-limited regime, both in the administrative Asia and Australasia regions of significant (U,W) change and spatially contiguous nighttime lighted cities. (2) The attributions of change to external/climate and internal/human-induced effects indicated large-scale afforestation and deforestation occurring mainly in China and Russia, respectively, which is in agreement with previous research about the largest net gain and loss by illegal logging of forest [11,35,36]. Southeast Asia, where a massive conversion of tropical forest to industrial tree plantations [37] occurred, showed few pixels in pink as deforestation (but in light blue as afforestation). The possible cause is that industrial tree plantations there evaporate more than tropical forests, just like forests evaporate more than bare land as observed after afforestation in China (light blue). (3) Significant changes in (U,W) space, which were associated with decreasing vegetation greenness, showed that dry areas remained dry or got dryer, and wet regions remained wet or got wetter, and these changes were separated at the (D = 1) threshold line.
In general, complementing geographical presentation of remote sensing information and climate and ecosystem models provides improved understanding of terrestrial biosphere responses to climate and land use change. It will help shift traditional management of resources and risks to integrated monitoring and holistic responses across wide ranges of space-time scales.
Author Contributions: D.C. and K.F. designed the experiments; R.S. obtained the funding of RRI-KLOF201904; Y.G. and S.G. obtained the funding of XDA23100504; D.C. and C.Z. processed the data; D.C. analyzed the data and wrote the original paper, and all authors contributed to the revising of the manuscript.