Towards an Improved Environmental Understanding of Land Surface Dynamics in Ukraine Based on Multi-Source Remote Sensing Time-Series Datasets from 1982 to 2013

Ukraine has experienced immense environmental and institutional changes during the last three decades. We have conducted this study to analyze important land surface dynamics and to assess processes underlying the changes. This research was conducted in two consecutive steps. To analyze monotonic changes we first applied a Mann–Kendall trend analysis of the Normalized Difference Vegetation Index (NDVI3g) time series. Gradual and abrupt changes were studied by fitting a seasonal trend model and detecting the breakpoints. Secondly, essential environmental factors were used to quantify their possible relationships with land surface changes. These factors included soil moisture as well as gridded air temperature and precipitation data. This was done using partial rank correlation analysis based on annually aggregated time-series. Our results demonstrate that positive NDVI trends characterize approximately one-third of Ukraine’s land surface, located in the northern and western areas of the country. Negative trends occurred less frequently, covering less than 2% of the area and are distributed irregularly across the country. Monotonic trends were rarely found; shifting trends were identified with a greater frequency. Trend shifts were seen to occur with an increased frequency following the period of the 2000s. We determined that land surface dynamics and climate variability are functionally interdependent; however, the relative influence of the drivers varies in different locations. Among the factors analyzed, the air temperature variable explains the largest portion of NDVI variability. High air temperature/NDVI correlation coefficients (r = 0.36  ́ 0.77) are observed over the entire country. The soil moisture content is of significant influence in the eastern portion of Ukraine (r = 0.68); precipitation (r = 0.65) was most influential in the central regions of the country. These results increase our understanding of ecosystem responses to climatic changes and anthropogenic activities.


Introduction
Ukraine has experienced immense environmental and institutional changes during the last three decades due to both human-induced and environmental processes such as socio-economic transformation, ongoing urbanization, increased climate variability and frequency of hazardous events [1,2].The drastic social, institutional and economic alterations at the end of 1991 as a result of the Soviet Union collapse induced significant changes in distribution and extent of land cover, land use intensity, economic productivity in the agricultural sector and shifts in the land surface phenology [2,3].These socio-economic changes had radical results and one of the consequences was the widespread farmland abandonment especially in the northern and western regions [1,[4][5][6][7].Besides the ongoing socio-economic transformations, several studies discussed the changes in climate conditions.In the past 30-50 years, droughts have become more frequent and intense, covering up to half of the area of Ukraine every 10-12 years, and up to 20% every 2-3 years [8].All these described alterations had an impact on the ecosystem that can be reflected on the land surface.To understand the effect of these phenomena, long-term monitoring of land surface changes is required.
There is a limited number of studies on long-term land surface change analysis for all of Ukraine, as previous studies focused on specific applications of land cover change (e.g., forest loss, farmland abandonment) in particular areas of the country (Carpathian ecoregion, the river basin of the Dnieper river) [1,2,9].Moreover, little attention has been paid to the understanding of the main drivers of these changes as most of the studies focused on the factors influencing land abandonment [1,10].
Satellite image time series observations of land surface reflectance and vegetation indices, such as Normalized Difference Vegetation Index (NDVI) calculated from them, are often used for land surface change analysis as NDVI uses the spectral differences between red and near infrared reflectance and as a result exposes the properties of ecosystem [11].NDVI-based trend analysis can assist the investigation of changes of land surface over a period of time and help to understand the mechanisms behind the observed changes.This is because NDVI is not only related to a variety of vegetation properties (structural properties of plants, vegetation productivity) [12,13], but it can also indirectly represent the ecosystem conditions [14].For instance, several studies used NDVI as a proxy of land cover change [15] and land dynamics [16], land degradation [17], farmland abandonment [18], and phenological changes [19].Based on these complex characteristics of NDVI, we assume that it can be used to assess land surface dynamics over time.In the context of this study, we define land surface dynamics as the sudden or gradual transformation or modification of earth surface including vegetation, soil, water and its ecosystem functioning in the unit of time (e.g., day, month, year) including fast (e.g., land use change) and slow processes (e.g., land degradation), which can indicate the response to both anthropogenic and environmental factors [12,13,20].
The detection of significantly increasing or decreasing trends is often not enough to assess the environmental impact or to be able to attribute the change to factors behind it.Abrupt changes or non-linearities may occur within the long-term time series, as the trends in complex systems usually are not monotonously increasing or decreasing [20][21][22].Furthermore, the interpretation of these abrupt changes without using additional datasets or expert knowledge is often challenging [23].
In this regard, the overall aim of this research was not only to analyze the land surface dynamics over all of Ukraine for the time period from 1982 to 2013 but also to assess underlying processes and drivers of the changes.The specific objectives were: (1) to assess long-term land surface variability; (2) to integrate data from different sensors and investigate the suitability of different methodologies; and (3) to estimate the effects of the essential environmental variables namely soil moisture, air temperature and precipitation on land surfaces.Our study is, therefore, the first study that analyzed spatiotemporal changes of land surface over all of Ukraine for the long time period of 30 years to provide integrated and spatially continuous knowledge of the changes and its triggers using remote sensing approaches.

Study Area
Ukraine is located in Eastern Europe in the southeastern part of the East European plain, between 44 ˝20 1 N and 52 ˝20 1 N latitudes and 22 ˝51 E and 41 ˝15 1 E longitude (Figure 1).It covers 603,550 km 2 and consists of 24 administrative units (oblast) [24].Most of Ukraine is composed of vast plains stretching north from the Black and Azov Seas.There are also two mountainous areas in the west and south of the country, the Carpathians and Crimean mountains, respectively.The climate in the East European Plain, the medium Ukrainian Carpathians, and the Crimean Mountains is mostly temperate continental, with hot summers (t > 18 ˝C) and cold winters (t < ´5 ˝C); the southern coast of the Crimea has a subtropical Mediterranean climate [25].The area of Ukraine is divided into three main agro-ecological zones: steppe (southern and eastern regions, 40% of the area), semi-steppe (central regions) and mixed forest (north and north-west) [7,26] (Figure 1d).European Plain, the medium Ukrainian Carpathians, and the Crimean Mountains is mostly temperate continental, with hot summers (t > 18 °C) and cold winters (t < −5 °C); the southern coast of the Crimea has a subtropical Mediterranean climate [25].The area of Ukraine is divided into three main agro-ecological zones: steppe (southern and eastern regions, 40% of the area), semi-steppe (central regions) and mixed forest (north and north-west) [7,26] (Figure 1d).Despite the prevailing moderate continental climate, the study area is characterized by fairly significant differences in the humidity, temperature and length of the growing season [27].Hot and relatively dry summers, together with relatively severe winters, create good conditions for the formation of chernozems (black soils) that are highly fertile for most agricultural crops [28].The average winter temperature ranges from −8 °C to −12 °C.In the southern regions, the mean winter temperature is 0 °C.The mean summer temperature ranges from 18 °C to 25 °C, although it can reach above 35 °C.Despite the prevailing moderate continental climate, the study area is characterized by fairly significant differences in the humidity, temperature and length of the growing season [27].Hot and relatively dry summers, together with relatively severe winters, create good conditions for the formation of chernozems (black soils) that are highly fertile for most agricultural crops [28].The average winter temperature ranges from ´8 ˝C to ´12 ˝C.In the southern regions, the mean winter temperature is 0 ˝C.The mean summer temperature ranges from 18 ˝C to 25 ˝C, although it can reach above 35 ˝C.
The area of Ukraine is characterized by significant regional differences in rainfall and its distribution throughout the year with a summer maximum.The largest amplitude of the volume is characteristic to the south.Most of the precipitation falls in the Ukrainian Carpathians (more than 1500 mm per year).The precipitation ranges from 650 mm (in the west) to 400 mm (in the southeast).In dry years, precipitation is significantly reduced: in the coastal areas of the Azov and Black Seas up to 100 mm, in the steppe up to 150-200 mm, and in forest-steppe up to 250-350 mm.In winter, permanent snow cover occurs almost everywhere in the country.One of the key features of the climate of Ukraine are the droughts which in recent years have increased frequency and affect large areas [8,27,29].

Normalized Difference Vegetation Index
As a primary data source for the regional study the Global Inventory Monitoring and Modeling System (GIMMS) NDVI3g dataset was chosen as it is the only currently available dataset that spans since 1982 and includes a dense time series needed for the analysis of land surface dynamics.The GIMMS NDVI3g dataset is an improved 8 km NDVI dataset produced by the National Oceanic and Atmospheric Administration (NOAA) Advanced Very High Resolution Radiometer (AVHRR) instruments that extend from July 1981 to December 2013 [12,30].
The third generation GIMMS NDVI is a 15-day Maximum Value Composite that was acquired from seven different NOAA satellites (7,9,11,14,16,17,18), which have been processed using an adaptive Empirical Mode Decomposition.NDVI3g is appropriate for long-term studies of land surface trends in vegetation, seasonality and coupling between climate variability and vegetation over the last three decades [31,32].
The NDVI3g has different quality flags [17].Based on these flags, prior to analysis, all values with reduced quality were excluded.We used pixels with minimum 75% valid data excluding interpolated observations, possible snow covered pixels.In addition, the images for the year of 1981 were also excluded to have full-year observations.The maximum value compositing approach was applied in the NDVI3g dataset.

Environmental Variables
Several studies discussed the effects of climate variability on land surface change [33][34][35].The monthly mean temperature (TMP) and precipitation total (PRE) datasets were obtained from the Climate Research Unit (CRU).The datasets cover the period 1901 to 2014 with a spatial resolution of 0.5 ˝ [36].
Recent studies also show the impacts of soil moisture on vegetation at various temporal scales [34,37,38].They discuss the co-variation of NDVI and soil moisture, stating coherent trend changes [39].Furthermore, other authors showed that a soil moisture NDVI pixel-wise residual trend indicates land degradation [17].Soil moisture influences hydrological and agricultural processes, runoff generation, drought development and other processes.A soil moisture dataset [38,40] was chosen for the study, which was a merged product of active and passive microwave soil moisture products [40].The dataset covering the period 1978-2013 provides daily surface soil moisture with a spatial resolution of 0.25 ˝and depicts the soil moisture in around 2 cm layer depth [40].Therefore, we chose the following environmental variables precipitation, temperature and soil moisture for the analysis.
In addition, we used a land cover map [41] of 1990 as a reference (the main institutional changes happened after 1990) to characterize the possible associations of the environmental variables and NDVI in main vegetation formations: cropland, grassland, and forest, which have different seasonal characteristics and response to the disturbances [41].

Methods
Prior to processing, the raster time series datasets were reprojected to World Geodetic System (WGS) 84 coordinate reference system and resampled to a common 8 km grid, using a bilinear interpolation algorithm.The datasets were further smoothed with Savitzky-Golay filter as the datasets were not extremely noisy, and we intended preserving the local variations [42].We extracted the temporal window (1982-2013) of the environmental time series to match the same period with NDVI dataset.The water bodies and build up areas were masked out.All the analyses were conducted using R statistical software [43].

Trend Analysis
After the pre-processing of the image time series (1982-2013), we explored the trends of annually aggregated NDVI series.There are several methods described in the literature for time series analysis of satellite imagery [12,13,20].The simplest and most common method is the ordinary least-squares (OLS) regression, where the main assumption is that the land surface is changing linearly and gradually over time [44].The drawbacks of this method have been depicted by several authors [23,43].In particular linear regression assumes normality and independence of data (neglecting the temporal autocorrelation).To overcome this, we used the Mann-Kendall nonparametric trend test, which does not depend on the distribution of data and can be used with data with serial dependence.Tüshaus et al. [45] demonstrated the performance of trend analysis with the aforementioned method for the land surface dynamics monitoring.This non-parametric test evaluates the occurrence of a monotonic single direction trend in the time series [45].Furthermore, non-stationarity (time dependency) was excluded due to a significant Augmented Dickey-Fuller test and inspection of the autocorrelation graphs and the related tests of annually aggregated time series.
However, the land surface dynamics are not always a gradual process and often the changes are abrupt due to factors such as extreme weather conditions or changes in land management.In general, the changes that can be monitored using remotely sensed data can be classified as follows: (1) a seasonal or cyclic change; (2) a gradual change over time consistent in direction (monotonic); and (3) an abrupt shift at a specific point in time [20,44].
To study sudden and gradual changes in land surface, the Break For Additive Season and Trend (BFAST) algorithm was used, [43] which decomposes the series into trend, seasonal and reminder components to evaluate gradual and seasonal dynamics occurring within indices derived from satellite image time series such as NDVI [46][47][48].The season-trend model [47] (i.e., linear trend and harmonic component) was fitted to pre-processed monthly NDVI time series, and the stability of the model was assessed using a test that determines the significance of structural trend breaks.As we were focusing on trajectories of the changes, rather than the number of the breakpoints, we chose the model with either 0 or 1 breakpoints.These breakpoints were assumed to catch the most influential trend shifts in the time series, which represent drastic changes in ecosystem functioning [47,48].We classified the trends based on the modified methodology described by de Jong et al. [47].The authors used the term greening/browning for categorizing the changes based on the signal before and after the breakpoint.Instead of the above-mentioned terms, we used the concept of increase or decrease of the NDVI trends.

Correlation Analysis
We used correlation analysis to investigate the effects of essential environmental variables, namely soil moisture, air temperature and precipitation on variability of NDVI.Correlation analysis is one of the commonly used methods for assessing the associations between vegetation indices and independent variables.In the case of simple correlation, the linear relationship between two variables is determined.The disadvantage of such an approach is that the effect of other variables is ignored.For example, when the correlation between soil moisture and NDVI is studied, the impact of other factors, such as precipitation and temperature, on NDVI is disregarded.To reflect the internal relation between two variables, a partial correlation coefficient should be calculated.Partial correlation coefficients reveal the relation between two variables when controlling the impact of other variables.
As a result, when two variables are both linked to the third, only the relation between the former two variables is estimated while the impact of the third variable is removed [49,50].Partial Spearman correlation is nonparametric in the form of relationship, and it is robust against outliers and invariant against monotone transformation.

Trends of NDVI Series
Figure 2 shows the results of the non-parametric Mann-Kendall trend analysis of annually aggregated NDVI time series from 1982 to 2013 for Ukraine.Areal statistics (Figure 2b) depict the occurrence of positive and negative trends in each administrative unit proportional to its area.As changes of environmental conditions are most likely to affect the land-cover classes in a different way [33], we interpreted the relative influence of each factor based on the spatial context of the observed associations i.e., land cover class and the agro-ecological zone.coefficients reveal the relation between two variables when controlling the impact of other variables.As a result, when two variables are both linked to the third, only the relation between the former two variables is estimated while the impact of the third variable is removed [49,50].Partial Spearman correlation is nonparametric in the form of relationship, and it is robust against outliers and invariant against monotone transformation.

Trends of NDVI Series
Figure 2 shows the results of the non-parametric Mann-Kendall trend analysis of annually aggregated NDVI time series from 1982 to 2013 for Ukraine.Areal statistics (Figure 2b) depict the occurrence of positive and negative trends in each administrative unit proportional to its area.As changes of environmental conditions are most likely to affect the land-cover classes in a different way [33], we interpreted the relative influence of each factor based on the spatial context of the observed associations i.e., land cover class and the agro-ecological zone.Overall, around one-third of the area of Ukraine is characterized by significant positive trends, clustered mainly in northern and western parts, which are mostly in the relatively humid and temperature limited area with forests and forest-steppe zones.These areas are also characterized by less frequent crop cultivation and observed land abandonment processes [18].The secondary succession on the former croplands is captured by a positive trend signal in our analysis.Around 30% of cropland is characterized by positive trends, which also indicates changes in agricultural practices.A large cluster of the positive trend was identified in the northwest from Kiev, where the Chernobyl zone is located.This area has been closed since 1986 after the occurrence of the Chernobyl disaster, which is one of the major anthropogenic disasters of the last century.Within the Chernobyl zone, active vegetation growth and environmental recovery are observed and documented [4], which is also reflected in the results of the trend analysis.The observed positive trends in forest cover 32% of its area.This is evident in the areas such as mountainous regions in the Carpathian Mountains in the west of Ukraine and in the Crimean Mountains [9].
Statistically significant negative trends were less frequent (<2% of the total area) and are mostly scattered across the country.This could probably be explained by the fact that the area was more dynamic, and, instead of a monotonic trend, land surface was varying over the entire study period.Most of the steppe zone does not reveal significant trends.The steppe zone of Ukraine, sometimes also referred to as the 'bread basket' of Ukraine, is the area where most intensively cultivated lands are located.The absence of the big clusters of significant trends within this part of Ukraine may indicate that, despite the overall socio-economic and environmental processes that have occurred in the country since 1982, these areas were continuously managed in a similar way, which has not led to noticeable land surface trends.
It should be emphasized that the analysis of long-term earth observation time series is affected by several factors, such as the spatial resolution of the imagery, the level of temporal aggregation, number of observations, pre-processing steps, and the analysis method itself.The temporal resolution of the time series has an effect on the trend significance estimation.In comparison to the full-resolution time series, the use of annually aggregated data reduces the number of observations, which results in the underestimation of the trend significance.Nevertheless, annual aggregation decreases the risk of detecting false positive trends [12].Furthermore, the aggregation is an effective method to account for autocorrelation: dependence between successive observations in a time series [50].
Although monotonic trend analysis can reveal the general tendency of land surface changes, the trends are not always continuously increasing or decreasing but can fluctuate over time [21].To account for these changes, a breakpoint analysis was conducted.Figure 3 illustrates the spatial distribution of change categories.Pixels with insignificant trends (p ě 0.05) were masked out.For each category of land surface change, representative points were selected, and the temporal profiles of NDVI response and the observed trend were studied in detail (Figure 3a-d).
In general, NDVI varied and pixels with trend shifts prevailed (Figure 3).Increasing trends with the shifts is the dominant type covering 28% of Ukraine (Table 1).Although the hotspots of the monotonic increase can be noticed (3.67% of the area), the results of the analysis of the change type reveal a higher number of interrupted trends (i.e., increase with setback) and reversal trends (increase to decrease).The interrupted increasing trends could probably be conditioned by human interventions and by the effect of extremely dry periods (observed in 10.24% of the area).The reversal trend indicating trend change from increase to decrease stated the decline of ecosystem functioning.We suggest that these types of changes can be associated with the large socioeconomic disturbances followed by the collapse of the Soviet Union and associated processes of land abandonment after 1991 and land re-cultivation after the 2000s [6,7].
We also evaluated the percentage of change types in each land cover class (Table 1).The estimates reflect the percentage of the trend type within the dominant class, using the land cover information from the 1990 reference year.Generally, the interpretation of these trends is challenging as the trend shifts can be explained by several factors [47].The reversal trend can be attributed to long-term increasing human pressures [48], particularly the cultivation of abandoned croplands, which was also observed by Smaliychuk et al. [7].Another abundant category of change is the increasing trend with an abrupt negative change, which can be related to the extreme events such as drought.In Ukraine, frequent drought events have been observed during the last 16 years [29].The drought in 2003 in most of the pixels was not detected as a major break and this can be explained by the fact, that in long-term series from 1982-2013, the drought did not have the biggest break in the NDVI trajectory.Furthermore, during the consecutive drought in 2007-2010, the break was not detected in the first year of the drought.The effects of droughts should be further investigated.Another challenge for the interpretation is the spatial resolution of the dataset, as one pixel of NDVI3g dataset can contain diverse land cover types, due to the fragmented landscapes in some regions and the size of the land management units [51,52].Despite of this sub-pixel heterogeneity, the dataset is still representing the large-scale land use/land cover patterns and their spatiotemporal variance over this vast area.The validation of these trends is not straightforward because the field observations are not temporally dense and mostly were restricted to a few-time steps.In addition, the scale compatibility was another issue, as the datasets might be not representative for an 8 km grid [53].To overcome this, we compared our results with the results of the previous studies [1,[4][5][6][7][8][9] and investigated the validity of the results by local experts.

Characteristics of the Main Time Periods of the Trend Shifts in Ukraine
To understand the land surface dynamics and the main drivers of its changes, it is essential not only to look at the spatial distribution of the trend shifts but also the timing of these changes.We classified the time of the NDVI changes in three time periods of 1982-1992, 1993-2002 and 2003-2013 as the socio-economic pressures and environmental conditions were different during these periods (Figure 4).According to our results, the decrease with the abrupt positive break predominated during the first period of 1982-1992 within the steppe zone of Ukraine, which is the main agricultural zone of the country.This period was the last period before the collapse of the Soviet Union and was still characterized by extensive cultivation and agricultural land expansion in the country.
drought.In Ukraine, frequent drought events have been observed during the last 16 years [29].The drought in 2003 in most of the pixels was not detected as a major break and this can be explained by the fact, that in long-term series from 1982-2013, the drought did not have the biggest break in the NDVI trajectory.Furthermore, during the consecutive drought in 2007-2010, the break was not detected in the first year of the drought.The effects of droughts should be further investigated.Another challenge for the interpretation is the spatial resolution of the dataset, as one pixel of NDVI3g dataset can contain diverse land cover types, due to the fragmented landscapes in some regions and the size of the land management units [51,52].Despite of this sub-pixel heterogeneity, the dataset is still representing the large-scale land use/land cover patterns and their spatiotemporal variance over this vast area.The validation of these trends is not straightforward because the field observations are not temporally dense and mostly were restricted to a few-time steps.In addition, the scale compatibility was another issue, as the datasets might be not representative for an 8 km grid [53].To overcome this, we compared our results with the results of the previous studies [1,[4][5][6][7][8][9] and investigated the validity of the results by local experts.

Characteristics of the Main Time Periods of the Trend Shifts in Ukraine
To understand the land surface dynamics and the main drivers of its changes, it is essential not only to look at the spatial distribution of the trend shifts but also the timing of these changes.We classified the time of the NDVI changes in three time periods of 1982-1992, 1993-2002 and 2003-2013 as the socio-economic pressures and environmental conditions were different during these periods (Figure 4).According to our results, the decrease with the abrupt positive break predominated during the first period of 1982-1992 within the steppe zone of Ukraine, which is the main agricultural zone of the country.This period was the last period before the collapse of the Soviet Union and was still characterized by extensive cultivation and agricultural land expansion in the country.During the second period (1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002), most of the changes occurred in the northwest of Ukraine (Figure 4).The first post-socialistic years were described with high abandonment rates as well as with prevailing forest disturbances.The agricultural activities were mostly affected after the collapse of the Soviet Union due to the ongoing transformation processes and structural changes in the agricultural sector at that time.Almost 70% of cropland abandonment occurred within the first 10 years of the transition from a state-command to a market-driven economy [5].Since the late 1990s and early 2000s, due to the economic recovery, the re-cultivation of formerly abandoned land emerged [5].The last period (2003-2013) was also associated with strong environmental impacts (i.e., droughts and extreme weather conditions) [8,27].
Although the effect of sensor change that affected the quality of the AVHRR data is significantly reduced in NDVI3g [17,54], some research found that GIMMS3g NDVI data are still showing an inconsistency between sensors in humid, dry-sub humid, semi-arid and hyper-arid regions [54].As the AVHRR platform changes (1985, 1988, 1994, 1995, 2000 and 2004) could potentially cause the breakpoints in the trends [47], we have conducted the corresponding tests.We compared these changes with the frequency of the detected changes and the detected breakpoint did not show the same result, suggesting that our results are unlikely to be affected by the platform changes.

Factors of Land Surface Dynamics in Ukraine
We explored the correlations between NDVI time series and environmental variables (mean annual temperature, total precipitation and mean soil moisture) using the partial rank correlation method.Prior to correlation analysis, we evaluated multicollinearity by obtaining VIF (variance inflation factor) for each variable.Non-stationarity (time dependency) could be excluded due to a significant the Augmented Dickey-Fuller test and inspection of the autocorrelation graphs and the related tests of annually aggregated time series.As the VIF was smaller than two, we continued the partial ranked correlation analysis of the driving factors and NDVI.
Figure 5 illustrates the spatial distribution of the extent of the effect of each factor.We used a subtractive model for visualization, where the primary factors are illustrated with a pure color and the overlapping areas of influence are depicted by the intersection of colors.
the Soviet Union due to the ongoing transformation processes and structural changes in the agricultural sector at that time.Almost 70% of cropland abandonment occurred within the first 10 years of the transition from a state-command to a market-driven economy [5].Since the late 1990s and early 2000s, due to the economic recovery, the re-cultivation of formerly abandoned land emerged [5].The last period (2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013) was also associated with strong environmental impacts (i.e., droughts and extreme weather conditions) [8,27].
Although the effect of sensor change that affected the quality of the AVHRR data is significantly reduced in NDVI3g [17,54], some research found that GIMMS3g NDVI data are still showing an inconsistency between sensors in humid, dry-sub humid, semi-arid and hyper-arid regions [54].As the AVHRR platform changes (1985, 1988, 1994, 1995, 2000 and 2004) could potentially cause the breakpoints in the trends [47], we have conducted the corresponding tests.We compared these changes with the frequency of the detected changes and the detected breakpoint did not show the same result, suggesting that our results are unlikely to be affected by the platform changes.

Factors of Land Surface Dynamics in Ukraine
We explored the correlations between NDVI time series and environmental variables (mean annual temperature, total precipitation and mean soil moisture) using the partial rank correlation method.Prior to correlation analysis, we evaluated multicollinearity by obtaining VIF (variance inflation factor) for each variable.Non-stationarity (time dependency) could be excluded due to a significant the Augmented Dickey-Fuller test and inspection of the autocorrelation graphs and the related tests of annually aggregated time series.As the VIF was smaller than two, we continued the partial ranked correlation analysis of the driving factors and NDVI.
Figure 5 illustrates the spatial distribution of the extent of the effect of each factor.We used a subtractive model for visualization, where the primary factors are illustrated with a pure color and the overlapping areas of influence are depicted by the intersection of colors.Although the studies have discussed the impact of climatic variation on land surface changes [33,55,56], very few of them have focused on the impact of environmental factors on land surface trends in Ukraine as opposed to management and land tenure factors [1][2][3][4][5][6][7][8].Our results show that the changes in land surface can be partially explained by the combination of thermal and water conditions.
It is important to mention, that we do not raise any causal argument for the statistical results of presented correlations.The only aim we had was to assess the effect of an estimator and describe possible associations between environmental variables and land surface changes.Among all of the analyzed factors, air temperature explained most of the NDVI variability (Table 2).The impacts of air temperature were observed in 31.5% of the area, with a positive correlation coefficient (0.36-0.77).This can be explained by the fact, that in western and northern regions, the temperature is the main limiting factor for the vegetation growth, and with the increase of temperature over the study period, the correlation becomes higher.The percentage of the area in each land cover class with NDVI/Air temperature correlation depicts that grasslands were the most sensitive to temperature.Soil moisture content is influential in eastern (r = 0.68), semi-steppe, and steppe zones, covering around 12% of the area.High correlation coefficients between NDVI and soil moisture is observed in grasslands (Table 3).Furthermore, we observed small clusters of negative correlations between NDVI and soil moisture that might be referred to anthropogenic influences in this area and unfavorable socio-economic conditions.Finally, precipitation is a dominant factor (r = 0.65) in central regions of the country (Table 4).Relatively low correlation coefficients were estimated between NDVI and precipitation in forests, which were conditioned by lower water limitation compared to other classes.Again, there is a small cluster of a negative correlation between NDVI and precipitation in the northwest of Ukraine.Since it is not driven by climate variability, this could be an indicator of human-induced changes such as land transformations [57] (not represented on a map).Moreover, the results of NDVI/soil moisture partial correlation analysis show a higher number of pixels with a significant association within the study area than NDVI/rainfall.Although theoretically the precipitation is considered as an essential forcing factor for soil moisture variability, the latter is also affected by conditions such as soil hydraulic properties, land cover, evapotranspiration, and topography.This is in line with several studies [17,58] where soil moisture had better performance for the explanation of land surface dynamics.In addition, when evaluating the partial rank correlations between variables, the effect of the other confounders is eliminated by adjusting the target correlation for any other monotonic association with the variables of this study.Although the additional datasets give insight about the processes that can affect the land surface changes, further analysis should consider the use of complementary factors and include datasets representing anthropogenic activity.Several studies discussed the factors that affect the land abandonment in Ukraine, and they depicted that factors such as soil type, yield decrease, unfavourable socio-economic conditions and population density demonstrate a significant association with land abandonment and cultivation rate [1,10].The observed increasing and decreasing trends and their interruptions are only partially explained by temperature, precipitation and soil moisture.Other factors-such as CO 2 fertilization, agricultural management practices and the intensity of land use and fertilizer use may affect the land surface changes.To improve our analysis, this should be considered in the context of overall land surface dynamics.

Conclusions
The current study addressed the assessment of long-term land surface dynamics over a 32 -year period between 1982 and 2013 in Ukraine.The interpretation of the influences of environmental variables on the observed land surface dynamics were supplemented and improved by employing additional datasets.Fundamentally, we concluded that long-term increasing NDVI trends are widespread across the country.Decreasing trends with abrupt positive breaks predominated during the initial 1982-1992 study period within the steppe zone of Ukraine.This period is characterized by extensive cultivation and agricultural land expansion in the country.During 1993 and 2002, most of the changes occurred in northwest Ukraine.The first post-socialistic years in the country were characterized by high rates of land abandonment as well as widespread forest disturbances.Following the collapse of the Soviet Union, agricultural activities and practices were principally affected by transformational processes and structural changes in this sector.A general economic recovery has been occurring in Ukraine since the late 1990s and early 2000s.These conditions have stimulated re-cultivation of formerly abandoned land.This process has emerged mainly in the central portion of the country and is reflected by the monotonic increase and increase with a setback that is seen to predominate during the recent years of 2003-2013.Critical, high impact environmental events (including droughts and other extreme weather conditions) also occurred with an increased frequency during these years.Based on our results, we conclude that the combined trend and breakpoints analysis captures and reflects the spatial patterns of these trends as well as describes the timing of those changes.This capability will allow increased understanding of dynamic land surface processes over extensive areas.We found that the GIMMS NDVI3g dataset was extremely useful for regional scale trend analysis.In particular, a data record spanning more than three decades enabled a markedly improved trend characterization.
With respect to the effects of environmental factors on NDVI change, we conclude that partial correlation analysis is useful for determining the individual contribution of each specific environmental factor on NDVI dynamics.Correlation analysis shows that temperature accounts for approximately 30% of overall NDVI variance.Grassland and cropland environments show similar responses to the drivers, demonstrating high sensitivity to soil moisture and precipitation as well as temperature.
Additional analyses should include complementary causal factors along with datasets that accurately characterize anthropogenic activity.There is also a requirement for analyses at higher spatial resolutions.This will enhance interpretation of land surface trends and provide for assessments at scales complimentary to land management practices.
We are confident that this study demonstrates that the effective combination of different trend analysis techniques, integrated multiple datasets, and effective statistical modelling allows derivation of valuable, previously unavailable spatially explicit information on land surface dynamics and their causal factors.Finally, such information will enable an improved understanding of land surface processes and environmental changes ongoing over space and through time that are specific to Ukraine.

Figure 1 .
Figure 1.Study area, the Walter-Lieth climate graphs for different locations in the country (a-c) and the agro-ecological zones (d) [26].

Figure 1 .
Figure 1.Study area, the Walter-Lieth climate graphs for different locations in the country (a-c) and the agro-ecological zones (d) [26].

Figure 2 .
Figure 2. Spatial distribution of annually aggregated Normalized Difference Vegetation Index (NDVI) trends for the period 1982 to 2013 (a); and the areal statistics for each region (oblast) (b).

Figure 2 .
Figure 2. Spatial distribution of annually aggregated Normalized Difference Vegetation Index (NDVI) trends for the period 1982 to 2013 (a); and the areal statistics for each region (oblast) (b).

Figure 3 .
Figure 3. Spatial distribution of the change categories.The pixels with significant trends in both segments and significant monotonic trends are represented (p < 0.05).Temporal profiles show the most frequent change types.(a) Increase to decrease (b) Decrease with burst (c) Increase with setback (d) Increase to decrease.

Figure 3 .
Figure 3. Spatial distribution of the change categories.The pixels with significant trends in both segments and significant monotonic trends are represented (p < 0.05).Temporal profiles show the most frequent change types.(a) Increase to decrease (b) Decrease with burst (c) Increase with setback (d) Increase to decrease.

Figure 4 .
Figure 4. Spatial distribution of the time of the shift classified in three major time intervals.Figure 4. Spatial distribution of the time of the shift classified in three major time intervals.

Figure 4 .
Figure 4. Spatial distribution of the time of the shift classified in three major time intervals.Figure 4. Spatial distribution of the time of the shift classified in three major time intervals.

Figure 5 .
Figure 5. Dominant factors of Normalized Difference Vegetation Index (NDVI) change: pixel-wise correlation coefficients between NDVI and environmental variables.The graycolor denotes pixels without statistically significant correlations (p < 0.05).

Figure 5 .
Figure 5. Dominant factors of Normalized Difference Vegetation Index (NDVI) change: pixel-wise correlation coefficients between NDVI and environmental variables.The graycolor denotes pixels without statistically significant correlations (p < 0.05).

Table 1 .
Percentage of change types in each land cover class.

Table 1 .
Percentage of change types in each land cover class.

Table 2 .
Percentage of the area in each land cover class with Normalized Difference Vegetation Index (NDVI)/Air temperature correlation.

Table 3 .
Percentage of the area in each land cover class with Normalized Difference Vegetation Index (NDVI)/Soil Moisture correlation.

Table 4 .
Percentage of land cover class with significant correlation coefficients Normalized Difference Vegetation Index (NDVI)/ precipitation.