**Changes in the Water Temperature of Rivers Impacted by the Urban Heat Island: Case Study of Suceava City**

#### **Andrei-Emil Briciu 1,, Dumitru Mihăilă 1, Adrian Graur 2, Dinu Iulian Oprea 1, Alin Prisăcariu <sup>1</sup> and Petru¸t Ionel Bistricean <sup>3</sup>**


Received: 26 March 2020; Accepted: 6 May 2020; Published: 9 May 2020

**Abstract:** Cities alter the thermal regime of urban rivers in very variable ways which are not yet deciphered for the territory of Romania. The urban heat island of Suceava city was measured in 2019 and its impact on Suceava River was assessed using hourly and daily values from a network of 12 water and air monitoring stations. In 2019, Suceava River water temperature was 11.54 ◦C upstream of Suceava city (Mihoveni) and 11.97 ◦C downstream (Ti¸său¸ti)—a 3.7% increase in the water temperature downstream. After the stream water passes through the city, the diurnal thermal profile of Suceava River water temperature shows steeper slopes and earlier moments of the maximum and minimum temperatures than upstream because of the urban heat island. In an average day, an increase of water temperature with a maximum of 0.99 ◦C occurred downstream, partly explained by the 2.46 ◦C corresponding difference between the urban floodplain and the surrounding area. The stream water diurnal cycle has been shifted towards a variation specific to that of the local air temperature. The heat exchange between Suceava River and Suceava city is bidirectional. The stream water diurnal thermal cycle is statistically more significant downstream due to the heat transfer from the city into the river. This transfer occurs partly through urban tributaries which are 1.94 ◦C warmer than Suceava River upstream of Suceava city. The wavelet coherence analyses and ANCOVA (analysis of covariance) prove that there are significant (0.95 confidence level) causal relationships between the changes in Suceava River water temperature downstream and the fluctuations of the urban air temperature. The complex bidirectional heat transfer and the changes in the diurnal thermal profiles are important to be analysed in other urban systems in order to decipher in more detail the observed causal relationships.

**Keywords:** diurnal thermal profile; urban tributaries; wavelet; covariance

#### **1. Introduction**

An increase in water temperature of rivers all over the world due to global warming is depicted in numerous studies [1,2]. At the same time, the increasing urban population worldwide leads to a greater and more territorially focused impact of cities on the environment. This impact, enhanced by the current climate changes, is exerted on urban rivers through multiple anthropogenic stressors such as the increasing imperviousness of catchments or the diminishing of areas with riparian vegetation [3,4]. The land use changes associated with the urbanization of an area (e.g., widespread paving) [4,5] lead to alterations in the thermal regime of the air in cities, and urban heat islands occur [6]. The warmer urban surfaces generate temperature surges in streams after localized rainstorms [3,7]. The impervious

urban areas reduce rainwater infiltration in the catchment and, as a consequence, the high waters are more intense and the baseflow is weaker [8,9]. The shallower waters during baseflow are more prone to being heated by the urban air [3]. The higher temperature of urban stream waters reduce their self-purification capacity by affecting the aquatic biota and the amount of dissolved oxygen in the water [10]. Changes in fish populations' structure and diversity are easily identifiable in rivers switched to an urban regime, and species that better tolerate warmer waters replace the previous ones [11,12].

Rivers passing through cities may have mean temperatures above or below the corresponding mean air temperature [13]. In both cases, discussions may arise concerning not only the influence of the city on the river, but also the inverse case. Colder rivers can have a cooling effect on a city, which may be enhanced by air circulation along the stream corridor [14–16]. At the opposite side, the warmer rivers have an effect of enhancing the urban heat island [17].

Due to the important role of water temperature in aquatic ecosystems, river water temperature is included in water quality indices [18]. Also, some software solutions were proposed to model the thermal impact of cities on rivers in order to find pollution causes and improve city management [19,20]. Overall, the increase in the urban stream water temperature is recognised as an essential aspect of the urban stream syndrome [21].

As an effect of climate change, river temperatures in north-eastern Romania are expected to increase by 1–2 ◦C at the end of this century, according to some models [13]. The predicted increases in air temperature in this part of Europe [22] will be higher in urban areas. The temporal evolution of stream water temperature of rivers in Romania was analysed only in a few studies [23–27]. Fewer are the studies about river water temperature in cities [26,27]. In this study, we aim to describe the thermal impact of Suceava city on its stream waters, especially Suceava River, by using high temporal resolution data and multiple analysis methods. The syntheses obtained from our data are relevant for a usual city–river relationship in Romania. Specific objectives of this study are (1) to find if and how the urban heat island alters the water temperature of Suceava River and (2) to identify the role played by urban tributaries and wastewaters in the observed variations of Suceava River temperature.

#### **2. Materials and Methods**

#### *2.1. Study Area*

Suceava city is the centre of a metropolitan area of approximately 150,000 people and is located in the north-eastern part of Romania, in Suceava county. It is a city divided by Suceava River, which flows from NW towards SE and crosses the industrial area in the middle of the city. Terrain elevation in the city ranges from 265 m above sea level (a.s.l.), in the floodplain/south-eastern part of the city, to 435 m a.s.l., in the hills of the northern limits. Climate is temperate continental. Winters record small amounts of precipitations and frequent negative temperatures, while summers are warm and rainy, with frequent torrential rainfalls. Mean air temperature is 8.16 ◦C, whilst the average annual sum of precipitation is 621 mm.

Suceava River has an average flow rate of 16.87 m3/s at I¸tcani station, located into the city [27]. It also has an annual flow regime imposed by the regional rainfalls and records frequent high waters (over 50 m3/s) in the late spring and early to mid-summertime interval. This river has well mixed waters due to its average maximum depth of under 1 m of water and due to numerous transversal hydraulic infrastructures, that generate hydraulic jumps (e.g., downstream of the I¸tcani and Burdujeni bridges). The main stream water monitoring stations used in this study are positioned on Suceava River upstream and downstream of Suceava city (Figure 1, Table 1). The other 3 monitoring stations are located on 3 small streamflows, typical for the tributaries received by Suceava River between Mihoveni (upstream) and Ti¸său¸ti (downstream). The cumulative streamflow of Suceava River tributaries in the study area is roughly 0.5 m3/s. The Collector streamflow is a semi-artificial streamflow which represents a mixture of urban rainwater and groundwater drainage and, in a small proportion, wastewater. It was created in the communist era as part of the measures to drain water from Suceava floodplain. Untreated wastewaters are also present in the flow of ¸Scheia and Cetă¸tii Creek, included in our monitoring network, and in the other streamflows, as a result of illicit domestic wastewater discharges. The wastewater treatment plant of Suceava city is located on the left bank of Suceava River and, in 2019, it discharged an average of 0.34 m3/s of water. This is also one of the 5 stations used to collect air temperature and relative humidity in Suceava floodplain (Figure 1, Table 1). These latter stations are part of a wider monitoring network (19 stations) that we used prior to our hydrological analysis in order to obtain maps of the urban heat island and its meteorological effects that help us to better understand our study area (Figure 2, Supplementary Figure S1).

Major and static sources of heat are a few and are represented by the thermoelectric plant (near the wastewater treatment plant), the cellulose factory and 5 plants that belongs to the food industry. Many (tens of thousands) minor static sources are represented by domestic water and air heaters (such as boilers). Tens of thousands of cars are a daily source of mobile thermal pollution. In the official build-up area of Suceava city, the industrial areas represent 11.16%, other buildings (mostly domestic) represent 24.37%, roads and pavements represent 6.14% and the remaining space is constituted of semi-natural spaces, such as private yards, parks or bare lands [27].

The urban heat island of Suceava is easily identifiable in summer (Figure 2b) and produces a distinct area of lower relative humidity in winter (Figure 2c).

**Figure 1.** Map of the study area (centred on Suceava city) and location of the monitoring stations.


**Table 1.** Details of the monitoring stations used in this study.

**Figure 2.** Urban heat island (UHI) above Suceava metropolitan area as revealed by hourly air temperature measurements in 2019 during winter—(**a**) (December, January, February)—and summer months—(**b**) (June–August)—UHI represents the area most impacted by the urban heat island.

#### *2.2. Instruments*

From June 28 until September 28, 2018, water temperature data used in this study for the analysis of Mihoveni (upstream) and Ti¸său¸ti (downstream) monitoring stations was measured with TruBlue 585 CTD instruments (accuracy: ±0.2 ◦C, resolution: 0.01 ◦C). From September 2018 until the end of 2019, water temperature and level measurements in our monitoring stations on Suceava River were obtained with 2 AquaTROLL500 instruments (accuracy: ±0.1 ◦C, resolution: 0.01 ◦C). The temperature measurements of the other streamflows were obtained with DS1922L-F5 iButton instruments (accuracy: ±0.5 ◦C, resolution: 0.0625 ◦C). Different stream waters, e.g., Suceava River or its tributaries, were measured with instruments with an accuracy identical for the entire streamflow type.

Air temperature and relative humidity (RH) in the stations used to create the maps in Figure 2 (with the exception of the official weather station of Suceava city) were measured with CEM DT-171 thermohygrometers (accuracy: ±1 ◦C and ±3.0% to ±5.0% RH, resolution: 0.1 ◦C and 0.1% RH) placed in miniature weather stations, made of wood, with good air circulation, preventing the direct solar radiation, with waterproof roof and fixed on wooden pillars at 2 m from the ground.

#### *2.3. Data*

We used hourly and daily time series with different lengths from 2018 and 2019 at UTC+2. The year 2019 was selected for in-depth analyses due to the broader data availability for all parameters. Data used in this study is divided into data obtained from our own instruments and data provided by various Romanian institutions. Data from our own instruments is represented by measurements of water temperature (5 sites) and, secondarily, water level (we used here, for a complementary analysis, only level data measured at Ti¸său¸ti with AT500 instruments); also, from the air environment, we obtained hourly values of temperature and relative humidity in 18 sites, of which 4 (located in Suceava floodplain—Table 1) were selected for further integration in the meteorological-hydrological analysis and comparison. Water temperature at Mihoveni (upstream) and Ti¸său¸ti (downstream) was recorded from 28 June 2018 until the end of 2019. Water level at Ti¸său¸ti used here was recorded for the entire year of 2019. Water temperature of the other streamflows was recorded from 20 September 2019 until the end of 2019. The sensors were placed in areas with persistent water flow. In the case of Suceava River, sensors were permanently under at least a 0.33 m column of water; for the tributaries, the minimum water column had a height of 0.1 m (maximum possible during low flow). Our meteorological data were recorded during the entire year of 2019 and above mostly grassy active surfaces. Water temperature data of Suceava River at Mihoveni (upstream) and Ti¸său¸ti (downstream) and air temperature and relative humidity data of Suceava floodplain are available at http://water.usv.ro/data.php.

Data from Romanian institutions is represented by data provided by ANAR (Romanian Waters National Administration—discharge and temperature data of Suceava River and atmospheric precipitation, both at I¸tcani—daily values for the entire year of 2019), ACET Suceava (the company responsible for water supply and sewerage in Suceava metropolitan area—discharge and temperature data of Suceava wastewater treatment plant effluent—daily data of 2019), ANM (National Meteorological Administration—air temperature, relative humidity and precipitation at Suceava Weather Station—hourly and daily data of 2019) and ANPM (National Environment Protection Agency—air temperature, relative humidity and solar radiation at SV2 station—hourly and daily data for 2019).

In 2019, the average yearly flow rate of Suceava River was 19.5 m3/s at I¸tcani station, which is only 15.6% higher than the multiannual discharge of 16.87 m3/s (6.9 m3/s standard deviation) [27] and, therefore, 2019 can be considered an average hydrological year. But, on the other side, the atmospheric environment above the city during the same year was hotter (+2.09 ◦C) and drier (−86 mm; −13.8%) than the last 60 years' average, the mark of a continuous climate change and increased urbanization of Suceava city metropolitan area. We suggest that 2019 can be taken into account as a case study year, which is valuable especially because studies using hourly data on thermal pollution of urban rivers are missing in Romania.

#### *2.4. Analysis Methods*

The maps that used interpolation of stations/points with measurements of air temperature and relative humidity were created in ArcGIS with a hybrid method that combines the multiple regression method with the surface interpolation of stations/points using an inverse distance weighted (IDW) method [28]. The spatial tendency of a parameter is calculated by applying the multiple regression. The difference between the real values and those estimated by regression (residues) are then interpolated using the IDW. Finally, the raster with the interpolated residues is added to the spatial trend in order to correct the spatial representation of the analysed parameter.

In this study, boxplots, seasonal adjustment, smoothing, normalising, surface plots and scalograms are done in MATLAB. The surfaces plot was created in MATLAB by using the "surf" function. The surfaces plot in this study used normalised matrices composed of vectors extracted from time series that had the seasonal oscillations removed by subtracting from the raw time series the smoothed variant of the selected time series. The subtraction acted as a high-pass filter. The smoothing used a moving average filter with a span of 49 values (in order to include 2 consecutive days). The vectors/columns of the matrices represent average hourly values and the normalisation process computed the z-score values of each vector (centre 0, standard deviation 1). The normalisation was necessary for representing on the same scale (for comparison purposes) the diurnal thermal amplitudes from both semesters of 2019 and from both environments, even if the real amplitude is higher in the warm semester and in the air.

Scalograms represent the graphic result of continuous wavelet analysis or wavelet coherence analysis. The wavelet transform involved the Morlet mother wavelet, and the continuous wavelet transform analysis used a rectified power spectrum [29]. We have chosen Morlet as the mother wavelet because it is frequently used in the analysis of time series in geosciences due to the good time–frequency localisation of events. The wavelet coherence analysis used the methodology of Grinsted et al. [30], which continued to implement the Morlet mother wavelet in the study of hydrological time series. For a detailed description of the wavelet analysis and equations, see our complementary studies [31,32].

XLSTAT (Statistical Software for Excel) was used for ANCOVA (analysis of covariance), which involved one quantitative dependent variable (parameter) and multiple quantitative explanatory variables (other parameters) for calculating the goodness-of-fit coefficients of the linear regression model, obtaining the Fisher's F test result (estimates the risk in assuming that the null hypothesis is wrong/that there is no effect of the explanatory variables) and producing the Type III SS Table (Sum of Squares analysis). This table estimates the contribution of each explanatory variable to the model by removing one variable at a time and evaluating its effect on the quality of the computed model (this removes the undesirable effect of the order in which the variables are selected in the process of model calculation). ANCOVA assumptions are: the dependent variable and covariate variables are measured on a continuous scale, there are no significant outliers, residuals are approximately normally distributed for the independent variable, there is a homogeneity of variances and there is a homogeneity of regression slopes.

An energy balance of the water system in this study was calculated according to the following formula:

$$\text{temperature}\_{\text{SVT}} = \left( \text{mass}\_{\text{SVM}} \cdot \text{temperature}\_{\text{SVM}} + \text{mass}\_{\text{TR}} \cdot \text{temperature}\_{\text{TR}} + \text{mass}\_{\text{WW}} \cdot \text{temperature}\_{\text{VW}} \right) \tag{1}$$

$$\left( \text{mass}\_{\text{SVM}} + \text{mass}\_{\text{TR}} + \text{mass}\_{\text{WW}} \right)$$

where mass is the mass of the water flows calculated from discharge, SVT is Suceava River downstream (Ti¸său¸ti), SVM is Suceava River upstream (Mihoveni), TR is the sum of Suceava River tributaries between the upstream and downstream monitoring stations on Suceava River and WW is the wastewater from the wastewater treatment plant of Suceava city.

#### **3. Results and Discussion**

In the monitored period of 2018 and 2019, water temperature of Suceava River ranged from 0 to 28 ◦C. It became supercooled during the 2018/2019 winter (minimum: −0.19 ◦C at Mihoveni (upstream) and −0.13 ◦C at Ti¸său¸ti (downstream)), when frazil ice was forming, and warm during the summers (maximum: 28.57 ◦C at Mihoveni and 28.65 ◦C at Ti¸său¸ti), due to heat waves in the atmosphere (Supplementary Figure S2). As the extreme values already revealed, Suceava River was generally warmer downstream of Suceava city (average: 11.97 ◦C upstream and 12.4 ◦C downstream). One can observe especially in the November–February months, that when water temperature upstream is near the freezing point, the same parameter downstream still exhibits diurnal cycles. These cycles in streamflows are created by similar cycles in air temperature and other mechanisms [33] and are still present downstream of Suceava city during cold days due to the heat input that the river received inside the city.

In order to find out if this increase in temperature is caused by the warmer air above Suceava city, we compared water and air temperature during an entire calendar year, 2019 (Figure 3). We selected two atmospheric datasets for this comparison. The first dataset is represented by the air temperature time series from Suceava Weather Station. Data from this monitoring station are less influenced by Suceava city due to its location outside the built-up area and are representative for the semi-natural environment that surrounds Suceava city and which represents the Suceava River catchment in Suceava Plateau.

**Figure 3.** Comparison of water temperature at Mihoveni (upstream) and Ti¸său¸ti (downstream) and air temperature in Suceava city area in 2019 (hourly data).

The second dataset is represented by time series obtained by averaging hourly values from the selected monitoring stations located in the urban floodplain of Suceava River (SV2, ITC, TRN, AMB, ACT—Table 1). This dataset is representative of the urban air of Suceava, as it was recorded in a complex residential, industrial and commercial area.

As expected, the thermal amplitude is higher in the air and the water temperature follows the air temperature (exceptions are the days with persistent negative air temperature, when water temperature stays around 0 ◦C due to the high energy loss necessary for freezing). In 2019, the average air temperature at Suceava Weather Station was 10.35 ◦C, while in the floodplain it was 11.05 ◦C. From the thermal difference of these stations (0.7 ◦C), only 0.5 ◦C might be explained by the vertical temperature gradient in the air due to the fact that the weather station is placed on a plateau at approximately 83 m above the mean elevation of the monitoring stations in the floodplain. Therefore, a surplus of 0.2 ◦C exists in the middle of the city as a result of the urban surfaces and activities. Water temperature was 11.54 ◦C upstream of Suceava city (Mihoveni) and 11.97 ◦C downstream (Ti¸său¸ti). Studies on urban river temperature increase often find increases above 1 ◦C downstream of urban areas (e.g., 4–5 ◦C increase at some stations in North Carolina [34]). The amplitude of the temperature increase of our studied river is more similar to those found by Zeiger and Hubbart [35] in Missouri (0.2–0.7 ◦C). The difference between the mean values of the monitoring stations upstream and downstream in 2019 was 0.43 ◦C. The accuracy of the temperature sensors of the AT500 instruments used for measuring

Suceava River water temperature is ±0.1 ◦C. The standard error of the mean is often used to assess the error of continuous measurements [36,37], and the standard error of the mean of Suceava River water temperature measurements is 0.085 ◦C for the upstream time series and 0.084 ◦C for the downstream time series. Therefore, we conclude that there is a relevant water temperature increase downstream.

We can observe that Suceava River is warmer than the air above the city at both moments, of entering and of exiting its metropolitan area. A calculus done for the warm semester (April–September, average value from hourly data), in order to exclude winter months with water persisting around 0 ◦C when the air temperature is frequently negative, has revealed a similar context: air at weather station = 16.67 ◦C, air in the floodplain = 17.7 ◦C, water upstream = 17.84 ◦C and water downstream = 18.19 ◦C. Suceava River is heated above the air temperature as a result of its warmer banks, which absorbed the solar radiation, and this reason is available for its tributaries too, whose shallower water is more prone to overheating. This heating is stronger in constructed areas (Suceava city in the floodplain is 0.7 ◦C warmer than at Suceava Weather Station) due to the numerous man-made structures (buildings, roads, channels) that transfer their excess heat to the water nearby. Small urban streams have a low inertia and are more susceptible to urban influence [38].

During the communist era, a network of pipelines was built in order to transport hot water from the thermal power station towards the city inhabitants. Those pipelines are still used today not only for heating (during winter), but also for hot showers (in some homes). The pipelines transport warm water from the thermal power station (located in the floodplain, near the wastewater treatment plant) to various urban users. These pipelines are on land surface or underground and transfer some of their heat to the nearby soil and groundwater (hot water leaves the plant with 70–80 ◦C and arrives at the consumers, sometimes after more than 10 km of pipe, with 48 ◦C or less). The warm water is then discharged into the sewerage system in order to be collected at the wastewater treatment plant. The sewerage network often contains wastewater that is warmer than the groundwater. Warmer groundwaters diffusely recharge Suceava River and its urban tributaries in the floodplain (this recharge direction is proven by the small springs that appear along Suceava Riverbanks). No studies exist about the temperature of the urban groundwaters in our study area. Groundwater tends to be cooler than air in the summertime and warmer in wintertime, but even during summer, some groundwaters may become warmer because hot urban surfaces may transfer their heat through soils into phreatic waters.

It appears that the average heat transfer between Suceava River and Suceava city is from water to air after the water has been heated by the urban soil and constructions, but the analysis of the average diurnal profiles of water and air in our study area indicates that the water and air relationship is bidirectional. One can observe in Supplementary Figure S3a,b, that the diurnal cycle of the air temperature in the floodplain has steeper slopes than at Suceava Weather Station and the moments of the maximum and minimum values occur earlier. During the late evening and the night, air is colder in the floodplain. Air temperature in the urban floodplain was warmest, compared to Suceava Weather Station air temperature, at 2 p.m. (2.46 ◦C difference), and air was colder by a maximum of 0.63 ◦C at 9 p.m. All these differences indicate that the air in the floodplain is strongly influenced by the urban land use (constructed surfaces heat up faster and cool more easily that a vegetated area). Similarly, after Suceava River flows through the city, the diurnal profile of Suceava River water temperature has steeper slopes and earlier moments of the maximum and minimum temperatures than upstream. Water downstream is warmer than upstream, especially in the afternoon when land and air become hotter. Water downstream was warmer than upstream with a maximum of 0.99 ◦C (4 p.m.) and a minimum of 0.05 ◦C (3 a.m.). Minima of the air temperature, at Suceava Weather Station and in the floodplain, and of water temperature, upstream and downstream, occur at 6 a.m., 5 a.m., 8 a.m. and 7 a.m., respectively. Moments of the average daily maxima for the same sites are 3 p.m., 2 p.m., 7 p.m. and 4 p.m., respectively. These values show that the peak position in the diurnal cycle of water has been strongly shifted towards the moment of the maximum air temperature because the heat gain during the day is an active process, in opposition to the passive process of heat loss during the night.

As previously observed, the heat transfer between air and water is bidirectional, both during an average day and during a year. This led to the question of whether this relationship alters the warming of Suceava River downstream. A simple statistical analysis revealed that 26.4% of the hourly values of Suceava River at the downstream station are actually lower than at the upstream station, and 11.2% of the days in 2019 had an average value that was lower downstream than upstream. We analysed the hourly differences in the diurnal regime of temperature during opposite seasons (summer and winter) and found that a constantly higher water temperature downstream is representative for winter, while, during summer night-time, Suceava River is colder downstream of the city (Figure 4a,b).

**Figure 4.** Hourly differences in the diurnal regime of temperature in 2019 during summer (left column) and winter (right column) for: Suceava River downstream—Suceava River upstream (**a**,**b**), floodplain air—Suceava Weather Station (**c**,**d**), floodplain air—Suceava River downstream (**e**,**f**).

During winter, the downstream–upstream thermal differences are lower than in summer. During the daytime, the heating of the active surface and of the air in the floodplain is much stronger than at Suceava Weather Station because of the built environment, but, during night-time, the floodplain becomes a space with frequent air temperature inversions of medium–low intensity (Figure 4c,d). These inversions are able, during most days in summer, to decrease the river water temperature downstream. During both seasons, the water of Suceava River is a heating source for the atmosphere of Suceava city, especially in winter (during 3/4 of an average day—Figure 4e,f).

Monthly average values in Figure 5a shows that, during March–June, air temperature in the urban floodplain of Suceava River is warmer than the stream water in both stations (or similar to it, in April). This might be attributed to the melting snows in the mountain area of Suceava River catchment and to the important rainfalls specific to this period of the year—all these create the colder water of Suceava River in the first half of the year. It is to be taken into consideration that 2019 was a dry year, with an annual sum of precipitation of 535.03 mm at Suceava station and 486.5 mm at I¸tcani station (in the floodplain). We can assume that, during years with average or excess precipitation, the impact of the urban heat island on Suceava River temperature is easier to detect.

**a b Figure 5.** Comparisons of water and air temperatures in the Suceava city area in 2019: (**a**) annual regime of temperature, (**b**) average diurnal profiles.

Studies on other city–river pairs indicate that a river can be warmer than the urban air even during winters with frozen surface waters [17]. Even small, linear water bodies such as a river may have a thermal effect on the surrounding environment, as proven by the 22 m wide UK river in the study of Hathway and Sharples [39]. Suceava River is usually 50–60 m wide inside the city, during baseflow. A warm river enhances the urban heat island [17].

On average, Suceava River receives heat from the urban air from ~9 a.m. until ~7 p.m. (Figure 5b). During the daylight, Suceava River water temperature inside the city increases quickly, but, during the night, it decreases slowly due to the increasingly higher heat loss necessary to lose 1 ◦C at lower and lower temperatures (in this case, temperatures under 30 ◦C; see isobaric specific heat [38]). In this way, during mid to late summer and during the autumn, the stream water gains energy from one day to another because the heat loss during the night-time is smaller than the energy gain during the daytime.

The most intense heat absorption from the atmosphere to the aquatic environment is done during the long summer daytimes by active evaporation processes from the water surface and the floodplain in an urban atmosphere drier than the surroundings. The atmosphere loses heat at the level of its basal layer—the heat is transferred to water through evaporation (Supplementary Figure S1). The most intense heat transfer from the aquatic environment and floodplain to the urban atmosphere takes place during the long winter nights, foggy mornings and evenings, when the condensation or sublimation of excess vapor in the urban atmosphere of Suceava creates a caloric surplus given to the atmosphere. This contributes to mitigating night-time and morning radiative losses (Supplementary Figure S4).

In both situations (both in the summer—through the remarkable intensity of the evaporation process—and in winter—through the long duration of the foggy/humid air interval in Suceava valley), the presence of water vapor in the atmosphere plays a very important and, at the same time, an ambivalent role in the thermo–caloric relationship between the urban atmosphere of Suceava city and the water of Suceava River. Similar strong thermal connections between water and air were also noted in previous studies [40].

Boxplots of hourly water and air time series show the higher, urban variability of the air temperature in the constructed floodplain, especially in August (Figure 6). In January, the interquartile interval of the water time series downstream of the city is bigger due to the less extreme negative temperatures in the floodplain air and the contribution of various heat sources. The observed distribution of air temperatures in 2019 is one of a warm year—the average temperature was 10.35 ◦C at Suceava Weather Station (in comparison, at the same station, the average temperature during 1950–2012 was 7.86 ◦C [27]).

**Figure 6.** Boxplots of water and air temperature time series (hourly data) for the following time intervals (the red line represents the median): (**a**) year 2019, (**b**) January 2019, (**c**) August 2019 (A1—Suceava Weather Station, A2—Suceava city in the floodplain, W1—Suceava River upstream, W2—Suceava River downstream).

The average diurnal profiles of water and air vary in their shape from month to month and from one monitoring station to another (Figure 7). The annual variation of the monthly average diurnal profiles is similar for the air time series (Figure 7a,b), but the variation in water time series is different (Figure 7c,d). Suceava River downstream recorded a high variability of the moment when the diurnal maximum value occurs. That moment varied from 3 p.m. to 8 p.m., and the urban air recorded only a fluctuation of ±1 h of the moment of the diurnal maximum. One can observe that the annual variation of Suceava River water has shifted from an original variation towards a variation that is specific to that of the local air temperature. In water time series upstream and downstream, the biggest differences in the position of maxima values occur during the warm semester. During winter, probably due to the atmospheric forcing that moves water temperatures to 0 ◦C for a long time, the differences in the hourly position of maxima values is minimal (the same is true for minima values).

A shape change of the stream water average diurnal profile from upstream to downstream of Suceava city similar to that observed in our study was also revealed by a previous study [27] that analysed data from a shorter time interval (October 2012–January 2013): the hour of the maximum temperature shifted from the evening (upstream) to the afternoon inside the city and also downstream of Suceava city.

**Figure 7.** Surface plots of seasonally adjusted and normalised average hourly values of air and water temperature in Suceava city area in 2019: (**a**) air at Suceava Weather Station, (**b**) air in the floodplain, (**c**) Suceava River upstream and (**d**) Suceava River downstream.

The scalograms of the continuous wavelet transform analysis indicate that the diurnal cycle occurs in more consecutive days in the monitoring station from Ti¸său¸ti—downstream (Figure 8). This persistence is sustained by heat transfer from the city into the river—this supplementary heat appears to attenuate the impact of external factors that tend to disrupt the repetition of a uniform diurnal cycle. The time interval with frequent significant (0.95 confidence level) diurnal cycles is longer downstream. During the cold semester, the diurnal cycles are less significant, especially in winter when, upstream, a diurnal periodicity may lack for many consecutive days. Also, the higher power of the diurnal periodicity downstream indicates that the diurnal cycle is easily shaped by the high thermal amplitude (day–night difference).

The wavelet coherence analysis of Suceava River water temperature at Mihoveni and Ti¸său¸ti (Figure 9a) reveals very frequent similar diurnal oscillations in both monitoring stations (that co-vary in a non-linear way). The diurnal periodicities are generally in phase, as indicated by the phase arrows pointing right. The scalograms of the wavelet coherence analysis of one station in air and the other station in water (Figure 9b,c) show that there is a better covariance of the diurnal periodicities in the urban floodplain air and Suceava River water downstream, at Ti¸său¸ti, than between the air relevant to the upstream catchment and Suceava River upstream, at Mihoveni. It is interesting to observe that the coherence between air and water downstream is stronger at the diurnal band of frequencies than the coherence between water upstream and water downstream. This indicates that the shape and timing of the diurnal water profile at Ti¸său¸ti is imposed by the diurnal distribution of the air temperature, rather than the characteristics of water at Mihoveni. However, at lower frequencies, the relationship between water time series in both monitoring stations is stronger, as indicated by the power spectrum and the fact that the continuous significant coherent area (dark red area topped by the thin black line of the statistical significance area limit) starts from ~5 days (128 h), instead of ~21 days (512 h). The wavelet analysis is useful for the analysis of the non-linear covariation of the diurnal cycle. The higher variability of the diurnal cycle does not have an impact only on a scale of a few consecutive days, but also, often, at a monthly scale. Thus, for example, the hourly position of minima values in the average diurnal profiles of the air temperature at Suceava Weather Station and in the floodplain and of the water temperature at Mihoveni (upstream) and Ti¸său¸ti (downstream) in July/August was 6/7 a.m., 5/6 a.m., 8/8 a.m. and 6/7 a.m., respectively.

**Figure 8.** Continuous wavelet transform analysis of hourly water temperature time series of Suceava River in 2019: (**a**) Mihoveni, (**b**) Ti¸său¸ti.

In order to find punctual sources of heat input for Suceava River, we analysed the temperature of some of its urban tributaries and the discharge characteristics of the wastewater treatment plant of Suceava city (the wastewater treatment plants are well-known sources of thermal pollution in rivers, worldwide [41] and in Suceava city [27]). With the exception of ¸Scheia River, the urban tributaries of Suceava River recorded water temperatures significantly greater than that of Suceava River (Figure 10) in the case study time interval.

During September 20–December 31, 2019, the average water temperature of (in upstream–downstream order; from hourly measurements) Suceava upstream, ¸Scheia, Cetă¸tii Creek, Collector and Suceava downstream was 8.54, 8.36, 11.64, 11.45 and 9.19 ◦C, respectively. Unlike the other urban tributaries, ¸Scheia flows mostly through a less urbanized landscape, before entering the city. This is why it has a temperature similar to that of Suceava River before entering the city. ¸Scheia lowers the temperature of Suceava River, as shown by measurements done ~0.4 km downstream, at I¸tcani station, where the annual water temperature is 12.4 ◦C and the average of the case study time interval is 7.9 ◦C (calculated from daily data, used with precautions because the daily values are based on measurements done only twice per day).

**Figure 9.** Wavelet coherence analysis of hourly water and air temperature time series of Suceava city in 2019: (**a**) Suceava River water temperature upstream and downstream, (**b**) air temperature at Suceava Weather Station and Suceava River water temperature upstream, (**c**) air temperature in the floodplain and Suceava River water temperature downstream.

**Figure 10.** Comparisons of water temperature of Suceava River and its urban tributaries during September 20–December 31, 2019: (**a**) temporal evolution of hourly data, (**b**) average diurnal profiles of water temperature.

The temperature of the other minor streamflows that discharge into Suceava River indicate how sensitive these waters are to the urban alteration of their streambed and catchment, to the urban heat island and to the disposal of wastewaters. Concrete embankments are frequent inside the city—these are an aquatic grey infrastructure that removes the natural processes from the rivers, including the interaction with the subsurface flow through the hyporheic exchange [6].

In 2019, the wastewater treatment plant of Suceava city discharged approximately 0.34 m3/s of wastewaters with an average of 16.31 ◦C. During September 20–December 31, 2019, it discharged 0.32 m3/s with an average of 15.6 ◦C. The effluent of the wastewater treatment plant is colder than Suceava River only during summer months, when the tap water that becomes wastewater is colder than Suceava River (the main water supply for Suceava city is from Berchi¸se¸sti groundwater, located at the plateau border with the mountain area). The average temperature of the discharged wastewater was 20.71 ◦C during summer (with ~1.16 ◦C lower than Suceava River at Ti¸său¸ti in the same time interval) and 12.24 ◦C during winter (with ~10.2 ◦C higher than Suceava River).

During September 20–December 31, 2019, we estimated the flow rate of Suceava River at I¸tcani as 3.94 m3/s. The average temperature of the tributaries that Suceava River receives between Mihoveni and Ti¸său¸ti is computed from the average of the three measured minor streamflows (10.48 ◦C). This average is combined in a weighted average with the wastewater temperature (cumulative discharge is 0.57 m3/s from 0.25 m3/s streamflow + 0.32 m3/s wastewater), resulting in a combined flow with a temperature of 13.35 ◦C (the discharge of tributaries is estimated as half of the 0.5 m3/s estimate by using the discharge of ¸Scheia during September 20–December 31, 2019, which is approximately half of the average specific to this time interval of the year). The combined flow temperature is added in a weighted average (discharge used for weighting) to Suceava River temperature upstream (8.54 ◦C) and results in a new temperature of Suceava River downstream of 9.14 ◦C.

The measured temperature in the field was 9.19 ◦C, resulting in this calculus/the surface water input into Suceava River explaining 0.6 ◦C (92.3%) from the total 0.65 ◦C difference between Ti¸său¸ti and Mihoveni. The remaining heat gain has to be attributed to the urban heat island (directly) and the groundwater input. It is to be remembered that the final heat gain is the result of complex processes of gaining or losing heat (e.g., through water input or exchange with atmosphere) and the estimated values are specific to the studied period of the end of 2019 only. Elements not taken into consideration (e.g., temperature of the raindrops) may have some role in the variability of Suceava River temperature. The indirect impact of the urban heat island on Suceava River temperature is exerted through urban tributaries. The thermal impact of the tributaries represents 31.75% (part of the above-mentioned 92.3%). It remains to be discovered how much of the temperature increase of tributaries over Suceava River temperature is caused only by the urban heat island, but a rough estimate is at least 50%. We can also estimate that the 0.05 ◦C left unexplained in the calculus of water temperature increase from upstream to downstream is also caused at least 50% by the urban heat island (direct influence). The measured impact of the wastewaters on the water temperature of Suceava River is not a surprise because the municipal treatment plants are considered as being probably the most important heat source for urban rivers [41]. Moreover, there was a constant increase in the temperature of wastewaters in cities reported, e.g., Tokyo [42].

The discharge of the wastewater treatment plant varies considerably, and the highest values are created by stormwaters collected during important rainfalls (Figure 11a). This is due to the fact that an exclusively pluvial drainage exists only in some small areas of the city. Most stormwaters are collected in the main (combined) sewerage and flow through the wastewater treatment plant. Therefore, the calculated impact of the wastewater is also representative of the impact of the stormwater runoff on the water quality of Suceava River. Temperature surges caused by stormwaters that washed warm city surfaces are difficult to detect, as most pluvial waters enter the combined sewers. Any detectable temperature surges caused by urban stormwater runoff must alter the average diurnal profile of Suceava River water temperature difference between the downstream station (Ti¸său¸ti) and the upstream station (Mihoveni) (Figure 11b). An analysis of the hourly temperature difference between these stations shows that the standard deviation is 0.76 ◦C and the maximum occurs at 2 p.m. In order to identify possible water temperature surges, we analysed the days with a temperature difference greater than 1 ◦C and with atmospheric precipitation greater than 10 mm (this amount generates consistent runoff, especially during high-intensity rains). From the 15 times when the rainfalls exceeded 10 mm/day, only 2 times generated a runoff that strongly impacted the stream water temperature downstream of Suceava city (Figure 11c,d). During the surge event in August 2019, the temperature difference was 1.64 ◦C and the temperature increase induced by urban stormwater merged the diurnal peaks of 2 days into a larger peak with a duration of 24 h. The surge in July 2019 lasted 13 h and recorded the maximum temperature difference of 2019 (4.14 ◦C, which is 1.1 ◦C higher than the second highest value). The temperature surges occurred during warm summer days, similar to many cases reported in other studies [3,35]. Our surge values are similar to those of Rice et al. [34] (1.9–3.27 ◦C), Zeigler and Hubbart [35] (4 ◦C) or Nelson and Palmer [3] (2–7 ◦C). The surge durations in our study are greater than those reported previously (up to 10–10.5 h) [34,35].

**Figure 11.** Links between rainfalls, wastewaters and the water temperature of Suceava River: (**a**) correlation between the rainfall amount and the wastewater treatment plant discharge, (**b**) average diurnal profile of Suceava River water temperature difference (downstream–upstream), (**c**) Suceava River water temperature difference (downstream–upstream) versus atmospheric precipitation at Suceava Weather Station during 3–4 August 2019 (48 h), (**d**) same as in (**c**), but during 14–15 July 2019 (24 h).

A correlation matrix (with Pearson coefficients) that took variables as daily averages during the entire year of 2019 into consideration (Supplementary Table S1) revealed very good correlations (>0.95) between Suceava River water temperature and air temperature or wastewater temperature. A good correlation in this linear estimate is between water temperature and solar radiation (>0.6). The good correlation of the solar radiation indicates the potentially high impact of other, unknown/unmeasured parameters (such as water turbidity) on the variability of water temperature and highlights the limits of any mathematical model.

ANCOVA tests, applied to some relevant parameters in order to find out which natural and anthropogenic factors influence the water temperature of Suceava River at Mihoveni (upstream) and Ti¸său¸ti (downstream), revealed the results displayed in Table 2. The goodness-of-fit (R2) of the linear model in explaining the upstream water temperature is 0.917, and it is 0.998 for the downstream case (in both cases, the risk in assuming that the selected explanatory variables explain the stream water temperature is <0.0001).

For the upstream case, the parameters that contribute more to the water temperature are air temperature, relative humidity and water level, but the contribution of the other parameters is not negligible. For the downstream case, the most important factors are water temperature from the upstream station and wastewater temperature, seconded by solar radiation and water level at the downstream monitoring station.


**Table 2.** Analysis of covariance (ANCOVA) of Suceava River water temperature upstream or downstream (dependent variable) versus other parameters (explanatory variables) - type III sum of squares analysis.

Abbreviations: WS—Suceava Weather Station, M—Mihoveni, I—I¸tcani, T—Ti¸său¸ti, fl.—floodplain, WTP—wastewater treatment plant.

In order to detect non-linear dependencies between multiple parameters, we applied a multiple wavelet coherence analysis. This uses the same characteristics of the previously discussed (simple) wavelet coherence and is a reliable method for the analysis of the non-linearly dependent processes [30]. The scalograms indicate the wavelet power of the first parameter explained by the other parameters (Figure 12). As in the case of ANCOVA, at the upstream station, the strongest links were detected between water temperature upstream and air temperature at Suceava Weather Station + water level upstream. However, at the downstream station, the most important factors linked to the evolution of the water temperature are air temperature in the floodplain and water temperature upstream (the latter has a stronger role than the temperature of the wastewaters, which was detected as relevant by ANCOVA). One can observe that the explanatory variables better explain the explained variable at the downstream station (the dark red areas are dominant and indicate stronger coherence, especially at the weekly and monthly scales).

An analysis of the scaled correlation between various time series was also conducted (Figure 13). The Pearson correlation coefficient was computed for every day (24 h/values generate 1 correlation coefficient and the result is a list of 365 correlation coefficients, 1 for each calendar day). We observed that, during 2019, the scaled correlation between water temperature at Ti¸său¸ti (downstream) and air temperature in the floodplain is better (the average of the 365 correlation coefficients is 0.6188) than the scaled correlation between water temperature at Mihoveni (upstream) and air temperature at Suceava Weather Station (0.3407). The classical correlation between the entire time series gave an irrelevant difference (0.8819 versus 0.8898); therefore, the scaled correlation is more useful in the process of finding differences between upstream and downstream monitoring stations.

Figure 13b shows that the scaled correlation between the downstream/floodplain pair and air temperature in the floodplain is better when the air temperature is higher, for almost every month. Figure 13c indicates that the same scaled correlation is lower when water level is high, especially during the warm season, when high waters occur frequently.

All data indicate that Suceava River temperature increases downstream of Suceava city because of numerous factors whose importance vary in time and in space. This increase is acquired effectively in the lower 2/3 of the distance between Mihoveni and Ti¸său¸ti (11.6 km straight line) and is caused by Suceava city, including its urban heat island.

**Figure 12.** Multiple wavelet coherence analysis scalograms between Suceava river water temperature and air parameters: (**a**) water temperature upstream versus air temperature at Suceava Weather Station and water level upstream, (**b**) water temperature downstream versus air temperature in the floodplain and water temperature upstream.

The increase in water temperature has an impact on other parameters of the river. For example, lower dissolved oxygen concentrations were reported downstream of Suceava city, at Ti¸său¸ti (compared to Mihoveni), partly caused by the increase in water temperature [31]. The lowered concentrations of dissolved oxygen are certainly diminishing the self-purification capacity of Suceava River that was attributed partly to this element [43,44].

**Figure 13.** Comparative temporal evolutions of various scaled correlations' time series (smoothed with a span of 1 week/7 values) and water/air parameters of Suceava city area in 2019: (**a**) scaled correlations of upstream versus downstream time series, (**b**) downstream scaled correlations versus air temperature in the floodplain, (**c**) downstream scaled correlation versus water level at Ti¸său¸ti.

**c** 

#### **4. Conclusions**

In 2019, the average air temperature was 10.35 ◦C at Suceava Weather Station and 11.05 ◦C in Suceava River floodplain/centre of the city. Suceava River water temperature was 11.54 ◦C upstream of Suceava city (Mihoveni) and 11.97 ◦C downstream (Ti¸său¸ti). The analysis of the average diurnal profiles of water and air in our study area indicates that the heat transfer between water and air is bidirectional. The peak position in the diurnal cycle of water has been strongly shifted towards the moment of the maximum air temperature of city. Monthly average values show that, in the first half of the year, air temperature in the urban floodplain of Suceava River was higher than the stream water temperature in both stations.

During a detailed spatial analysis (last quarter of 2019), the water temperature of some urban tributaries was ~2 ◦C higher than the water temperature of Suceava River. These tributaries are strongly impacted by the urban heat island and, together with the effluent of the wastewater treatment plant, increase the temperature of the main river.

The wavelet coherence analysis shows that there is a better covariance of the diurnal periodicities in the temperatures of urban floodplain air and Suceava River water downstream, at Ti¸său¸ti, than between the air temperature nearby the city and Suceava River upstream, at Mihoveni. Simple and multiple wavelet coherence analyses and ANCOVA indicate that the evolution of Suceava River water temperature downstream of Suceava city is controlled mainly by urban air temperature, water temperature upstream, wastewater temperature and solar radiation.

Last, but not least, 2019 was a dry and warm year (atmospheric environment) and Suceava River had warm waters. Further analyses using high temporal resolution data are needed during average or wetter/colder years in order to better understand the relationship between the urban heat island and Suceava River. For a better understanding of the distinct influence of the urban heat island and wastewaters on the changes of Suceava River, a monitoring station on Suceava River immediately upstream of the wastewater treatment plant is necessary for future studies. Data provided from this new station will help in assessing the gradual changes in the diurnal thermal profile of Suceava River from the uppermost monitoring station towards downstream of Suceava city. Data collected in this and future studies could be integrated into computer models that might help local authorities mitigate the impact of the urban heat island on urban stream waters. Best management practices should be promoted for Suceava floodplain. Restoration treatments may lower a river temperature and increase its dissolved oxygen concentration, as proven in other urban areas [6]. Because riparian and floodplain forests keep the river temperature low [45], we recommend a forestation with willows (*Salix* sp.) of Suceava River floodplain inside Suceava city in order to mitigate the urban stream syndrome.

**Supplementary Materials:** The following are available online at http://www.mdpi.com/2073-4441/12/5/1343/s1, Figure S1: Urban heat island above Suceava city as revealed by hourly air relative humidity measurements in 2019 during winter - a - and summer - b - (values used in maps are averages per season) – UHI represents the area most impacted by the urban heat island. Figure S2: Comparative evolution of Suceava River water temperature upstream (Mihoveni) and downstream (Ti¸său¸ti) of Suceava city (hourly data, June 28th, 2018 – December 31st, 2019). Figure S3: Comparisons of water temperature at Mihoveni (upstream) and Ti¸său¸ti (downstream) and air temperature in Suceava city area in 2019: a. average diurnal profiles of water temperature, b. average diurnal profiles of air temperature. Figure S4: Relative humidity diurnal regime of the air at Suceava Weather Station and in the floodplain in 2019 during summer (a) and winter (b) (ΔRH is the difference between the floodplain values and those of Suceava Weather Station). Table S1: Correlation matrix of selected water and air parameters in Suceava city in 2019.

**Author Contributions:** Conceptualization, A.-E.B.; Data curation, A.-E.B. and D.M.; Investigation, A.-E.B. and D.M.; Methodology, A.-E.B. and D.M.; Resources, A.-E.B., A.P., D.I.O. and P.I.B.; Supervision, A.G.; Visualization, D.I.O., A.P. and P.I.B.; Writing—original draft, A.-E.B. and D.M.; Writing—review and editing, A.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by CNCS–UEFISCDI, project number PN-III-P1-1.1PD-2016-2106.

**Acknowledgments:** Due to important field and laboratory work, Dumitru Mihăilă is to be considered, together with Andrei-Emil Briciu, as one of the principal authors. Wavelet coherence software was provided by A. Grinsted. Thanks for the data provided is addressed to three Romanian institutions, ANAR (Siret Water Basin Administration), ANM, ANPM, and to ACET Suceava. This study was supported by data obtained within the research project SQRTDA (Streamwater Quality Real-Time Data Analysis). This work was supported by a grant of Ministry of Research and Innovation, CNCS–UEFISCDI, project number PN-III-P1-1.1PD-2016-2106, within PNCDI III.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

## **Tracking Lake and Reservoir Changes in the Nenjiang Watershed, Northeast China: Patterns, Trends, and Drivers**

#### **Baojia Du 1,2, Zongming Wang 1,3, Dehua Mao 1,\*, Huiying Li <sup>4</sup> and Hengxing Xiang 1,2**


Received: 17 March 2020; Accepted: 31 March 2020; Published: 13 April 2020

**Abstract:** In terms of evident climate change and human activities, investigating changes in lakes and reservoirs is critical for sustainable protection of water resources and ecosystem management over the Nenjiang watershed (NJW), an eco-sensitive semi-arid region and the third-largest inland waterbody cluster in China. In this study, we established a multi-temporal dataset documenting lake and reservoir (area <sup>≥</sup> 1 km2) changes in this region using an object-oriented image classification method and Landsat series images from 1980 to 2015. Using the structural equation model (SEM), we analyzed the diverse impacts of climatic and anthropogenic variables on lake changes. Results indicated that lakes experienced significant changes with fluctuations over the past 35 years including obvious declines in the total area (by 42%) and number (by 51%) from 1980 to 2010 and a slight increase in the total lake area and number from 2010 to 2015. More than 235 lakes in the size class of 1–10 km2 decreased to small lakes (area < 1 km2), while 59 lakes covering 243.75 km2 disappeared. Total reservoir area and number had continuous increases during the investigated 35 years, with an areal expansion of 54.9% from 919 km2 to 1422 km2, and a number increase by 65.3% from 78 to 129. The SEM revealed that the lake area in the NJW had a significant correlation with the mean annual precipitation (MAP), suggesting that the MAP decline clarified most of the lake shrinkage in the NJW. Furthermore, agricultural consumption of water had potential impacts on lake changes, suggested by the significant relationship between cropland area and lake area.

**Keywords:** lakes; reservoirs; Nenjiang watershed; climate change; Landsat

#### **1. Introduction**

Inland waterbodies (lakes and reservoirs) are closely related to human life due to their diverse ecosystem services [1,2]. Lakes and reservoirs in the drylands, although covering only a small proportion of the landscape, play irreplaceable roles in fragile environments and for local residents [3]. Meanwhile, they are sensitive to climate change and human disturbances [4]. Monitoring changes in lakes and reservoirs and investigating their driving forces are thus of great significance to sustainable water resource management and regional economic development.

Due to diverse driving forces, lakes around the world have experienced changes in both area and number during the past decades. Thaw and "breaching" of permafrost have caused a widespread decline in the Arctic lake number and area [5]. Because of warmer air temperatures, which allow higher evaporation rates, as the world's largest surface freshwater system, the total area of the Great Lakes of North America experienced a continuous decrease from the 1980s to 2005 [6]. Under the background of global warming, glaciers in China's Tibetan Plateau (TP) showed a general retreat, while the glacial lakes expanded notably [7]. However, a drier climate and degraded permafrost have led to lake shrinkage in some basins of the TP [8]. Lakes over Mongolia experienced areal decline mostly related to a drier climate, while lakes in Inner Mongolia had notable shrinkage which was attributed mainly to human consumption of water resources, particularly coal mining [9]. Moreover, the lake area in the Jianghan Plain of China decreased dramatically from the 1950s to 1998, which was mainly caused by agricultural cultivation.

During the past 60 years, the number of reservoirs has increased notably, and there are about 50,000 large reservoirs nowadays around the world [10]. For example, China constructed nearly 45,000 reservoirs in the Yangtze River Basin to meet the large demand for water resources. By 2013, 98,000 reservoirs had been built in China [11]. Previous studies [12,13] indicated that most of large river systems and lakes in China were affected by reservoirs. Due to the construction of large reservoirs in upstream areas, water was impounded upstream, which seriously affected the water supply of lakes in the middle and lower reaches [14]. Therefore, it is necessary to track and analyze the changes in lakes and reservoirs at different scales and regions to form a scientific management response [15,16].

Optical images from different satellite sensors were widely used to monitor spatiotemporal variations of inland waterbodies at multiple geographic scales [17–22]. The Moderate Resolution Imaging Spectroradiometer (MODIS) data have been widely used to assess the water extent at daily to 16-day timescales despite a spatial resolution at 250 or 500 m [23]. Yet, areal changes of small lakes or reservoirs with irregular shapes were not accurately delineated due to the coarse resolution of source data. Other optical sensors, such as Quickbird (DigitalGlobe, Longmont, Colorado, USA) and IKONOS (DigitalGlobe, Longmont, Colorado, USA), provide finer images comparable to aerial photography for the extraction of lake or reservoir boundaries [24]. However, those data are limited in application at broad scale due to the high costs, narrow swath size, and so on [25,26]. Landsat series images have a fine spatial resolution (30 or 80 m) compared with previously mentioned data and have provided the longest temporal and spatial records for surface observations since their first launch in 1972 [27]. Consequently, Landsat imagery has been widely-used remote sensing data in examining changes in lakes or reservoirs.

As one of the third-largest waterbody clusters, an eco-environmentally fragile area, and an important base for grain production, the Nenjiang watershed (NJW) plays an important role in ecological conservation and national grain security in China [28]. However, the changes in lakes and reservoirs across the NJW during recent decades have rarely been examined and their correlations with climate change and human disturbances have rarely been quantitatively investigated. Doing so is critical to understanding the regional water cycle and sustainable water resource management. Therefore, this study aims to (1) employ Landsat 8 images to investigate the current status of lakes and reservoirs in the NJW, (2) evaluate changes in area and number of lakes and reservoirs using consistent Landsat images (i.e., 1980, 1990, 2000, 2010, and 2015), and (3) quantify the roles of climatic factors and artificial variables in driving lake changes.

#### **2. Materials and Methods**

#### *2.1. Study Area*

The NJW, covering an area of 297 <sup>×</sup> 103 km2, is located in the core region of Northeast China (Figure 1), with latitudes from 44◦1 48" N to 51◦42 1" N and longitudes from 119◦12 1" E to 127◦54 2" E [28]. Obvious terrain variances can be found in this area with elevations ranging from 120 to 1740 m above sea level. The highest elevation is in the northwest at the Greater and Lesser Khingan Mountains, and the lowest elevation is in the southeast at the Songnen Plain. This is why the lakes and reservoirs are mainly observed in the southeast. The NJW is dominated by a temperate

semi-arid continental climate with an mean annual precipitation of around 470 mm and a mean annual air temperature of 4 ◦C [29]. Most precipitation (82%) occurs in months from June to September which is a critical growth period for vegetation in this region. As an important grain base in China, the NJW plays an critical role in promoting regional economic development and ensuring national food security [30,31]. Moreover, this region provides important habitats for threatened species and migratory waterfowls.

**Figure 1.** Geographic location and general situation of the Nenjiang watershed.

#### *2.2. Data Source and Processing*

#### 2.2.1. Satellite Data

In this study, multi-temporal Landsat images from multispectral scanner (MSS), thematic mapper (TM), enhanced thematic mapper plus (ETM+), and operational land imager (OLI) sensors were acquired from the United States Geological Survey (USGS) to investigate the changes of lakes and reservoirs from 1980 to 2015. Specifically, we used 28, 27, 27, 29, and 26 scenes of images to extract the lakes and reservoirs for five dates of 1980, 1990, 2000, 2010, and 2015, respectively, (Figure 2). A total of 137 images acquired for months from June to September were used for the classification because the vegetation in this period is easier to be identified and the lakes and reservoirs have the largest amount. All of the images have a little cloud cover (less than 5%) and lakes/reservoirs on these images are clearly visible. In addition, in order to analyze the impact of floods on the area and the number of lakes and reservoirs in 1998 and 2013, we used 30 scenes of TM and 30 scenes of OLI images to extract the lakes and reservoirs in 1998 and 2013, respectively. Prior to image classification, data preprocessing including geometric, topographic, and radiometric corrections was performed for all the images using the ENVI 5.3 (Exelis Visual Information Solutions, Inc. Boulder, USA) software package. To ensure the data consistency, all images were re-projected to the 1984 WGS UTM zone 51N projection.

**Figure 2.** Landsat images selection in optimum months for different dates.

#### 2.2.2. Meteorological Data

In order to analyze the relationships of lake/reservoir area with climatic factors, the daily climatic data, including extreme air temperature and precipitation, mean wind speed, and sunshine hours during 1980–2015, were collected from the meteorological records of the China Meteorological Data Service Center. Figure 1 shows the locations of the meteorological stations in the NJW. Spatial patterns of the mean annual air temperature (MAAT) and mean annual precipitation (MAP) were interpolated from those meteorological records using an Anusplin software considering the elevation differences [32]. We used the daily meteorological data and the Food and Agriculture Organization of the United Nations (FAO) Penman–Monteith model [8,33] to calculate the actual evapotranspiration (ET). ET is calculated as follows:

$$ET\_0 = \frac{0.408\Delta (R\_n - G) + \gamma \frac{900}{T + 273} \mathcal{U}\_2 (\mathbf{e}\_s - \mathbf{e}\_a)}{\Delta + \gamma (1 + 0.34 \mathcal{U}\_2)} \tag{1}$$

$$ET = \; 9.78 + 0.0072 \times ET\_0 \times PPT + 0.051 \times PPT \times LAI \tag{2}$$

where *ET*<sup>0</sup> is the potential evapotranspiration, *Rn* is the net radiation, *G* is the soil heat flux, *es* is the saturation vapor pressure, *ea* is the actual vapor pressure, *es* − *ea* is the saturation vapor pressure deficit of the air, *T* is the air temperature at 2 m height, *U*<sup>2</sup> is the wind speed at 2 m height, Δ is the slope of the saturation vapor pressure temperature relationship, γ is the psychrometric constant, *PPT* is precipitation, and the calculation of *LAI* is referenced from the method of Lu et al. [34]. These parameters are directly calculated or derived from the average daily maximum and minimum temperature, the daily average temperature, the daily actual vapor pressure, the daily wind speed data, the actual duration of sunshine hours, the relative humidity data, and other empirical metrics.

In order to analyze the potential impact of climate change in the next 35 years (2015–2050) on regional lake area and number changes, we downloaded the Representative Concentration Pathways 5 (RCP 5) datasets from the China Agrosys Platform (http://stdown.agrivy.com). The future climate change data were generated based on the Fifth Generation Coupled Global Climate Model (CGCM 5) from the Canadian Centre for Climate Modeling and Analysis. The datasets constructed using original meteorological observations for each station included daily mean temperatures, daily precipitation, and daily potential evaporation. We converted daily mean temperature, daily precipitation, and daily potential evapotranspiration into annual mean temperature, annual mean precipitation, and annual potential evapotranspiration for our analysis.

#### 2.2.3. Other Data

The areal data of the cropland area and gross domestic product (GDP) index collected from local statistical yearbook (http://www.stats.gov.cn).

#### *2.3. Data Analysis*

#### 2.3.1. Extracting Lakes and Reservoirs

We focused on lakes and reservoirs greater than 1 km2 in terms of the 30 and 80 m resolutions of the satellite images [35]. An object-based image analysis (OBIA) was used rather than the traditional pixel-based classification method, because OBIA classifies objects instead of individual pixels [36–38]. In the process of OBIA classification, the spectrum, spatial information, texture, and geometric features characterized by remote sensing images were fully utilized. The lake and reservoir extents of the NJW in 1980, 1990, 1998, 2000, 2010, 2013, and 2015 were extracted in eCognition Developer 8.6 [14,39]. The normalized difference water index (*NDWI*) is the most popular used index for automated inland waterbodies delineation [35,40,41]. In addition, we tested different versions of *NDWI* (*mNDWI*, *NDRW*, and *NDWI*) [21,42] and found that the selected *NDWI* is more effective in our study area. Therefore, *NDWI* was applied to extract lakes and reservoirs, which was defined as:

$$NDWI = \left(Green - NIR\right) / \left(Green + NIR\right) \tag{3}$$

where *Green* and *NIR* represent the reflectance of the green and Near Infrared (*NIR*) bands, respectively [25]. The data processing steps for the lake and reservoir inventory are shown in Figure 3. The extraction of lakes and reservoirs consisted of three major steps: image multi-resolution segmentation, *NDWI* threshold testing, and classification rule designing. By defining the length/width index and the rectangular fit index built into the software, rivers and artificial ponds were removed. Then, we used visual interpretation to extract the reservoirs.

**Figure 3.** Flowchart for the lake and reservoir mapping process.

To validate the accuracy of extracted lakes and reservoirs, we calculated error matrices based on 1492 validation samples selected from Google Earth and Landsat sensors. The detailed number of validation samples for the five dates is given in Table 1. The overall accuracies of the lake and reservoir classification for the five dates were respectively evaluated and the standard errors bars (uncertainties) were estimated [43]. The overall accuracies of lakes and reservoirs were larger than 90%, and the kappa coefficients were larger than 0.87.


**Table 1.** Collected validation samples for the five dates.

MSS—multispectral scanner, TM—thematic mapper, ETM+—enhanced thematic mapper plus, OLI—operational land imager.

#### 2.3.2. Temporal Analysis of Lakes and Reservoirs

To fully understand the changes of lakes and reservoirs of different sizes, the lakes and reservoirs were categorized into four classes: 1–10 km2, 10–50 km2, 50–100 km2, and >100 km2. In order to examine the areal change rate of the lake or reservoir in different periods, an indicator of the lake or reservoir area dynamic degree, as shown Equation (4), was used to analyze their changes [44,45].

$$K = \frac{\mathcal{U}\_b - \mathcal{U}\_a}{\mathcal{U}\_a} \times \frac{1}{T} \times 100\% \tag{4}$$

where *K* is the dynamic indicator for the lake or reservoir area; *Ua* and *Ub* are the area of the lake and reservoir at the start date and end date, respectively; and *T* is the time scale under consideration.

#### 2.3.3. Assessing the Roles of Climatic Factors and Anthropogenic Causes in Lake Changes

The roles of climatic factors and artificial variables in changes of lakes during the study period were quantified using structural equation modeling (SEM) by the AMOS 22 software (IBM, Armonk, NY, USA) [39]. This paper analyzed the changes in MAAT, MAP, and ET to examine the influences of climate change on lake changes. The statistical data of cropland area and GDP were selected to investigate the impacts of artificial variables on lake changes [46,47]. Table 2 illustrates the optimum values for these indicators necessary for the SEM.

**Table 2.** Measures used to test the goodness of model.


#### **3. Results**

#### *3.1. Spatial Pattern of Lakes and Reservoirs in 2015*

Figure 4 shows the lake and reservoir distribution in 2015 across the NJW. Lakes and reservoirs were identified dominantly in the southeastern part of the NJW. A total of 233 lakes (area <sup>≥</sup> 1 km2) covering an area of 2110 <sup>±</sup> 53 km2 in 2015 were extracted from satellite images. Most of the lakes (89.4%) had an area of smaller than 10 km2. There were three lakes with area greater than 100 km2, and their accumulated area accounted for approximately 33.8% of the total lake area in the NJW. The Chagan Lake was the largest lake with an estimated area of 292 km2.

There were 129 reservoirs (area <sup>≥</sup> 1 km2) in the NJW with a total area of 1422 <sup>±</sup> 56 km2 in 2015. Of these, 108 reservoirs had an area smaller than 10 km2. There were two reservoirs with an area larger than 100 km2, with the accumulated area accounting for approximately 41.5% of the total reservoir area in the NJW. The Nierji Reservoir was the largest reservoir with an estimated area of 360 km2.

**Figure 4.** Spatial distribution of lakes and reservoirs in the Nenjiang watershed (NJW) in 2015.

#### *3.2. Temporal Changes of Lakes and Reservoirs from 1980 to 2015*

Figure 5 suggests that dramatic changes in the total area and number of lakes and reservoirs over the NJW occurred from 1980 to 2015. During the observed 35 years, the total lake area decreased by 38.7% from 3440 <sup>±</sup> 72 km2 to 2110 <sup>±</sup> 53 km2, whilst the total number of lakes decreased by 233 from the original 468 in 1980. Specifically, lake changes had significant fluctuations over the past 35 years, including obvious declines in total area (42%) and number (51%) from 1980 to 2010 and slight increases in the total lake area and number from 2010 to 2015. Reservoirs in the NJW experienced continuous expansion during 1980–2015. The total number of reservoirs increased from 78 to 129 with the total area expansion being 55%, from 919 <sup>±</sup> 53 km2 to 1422 <sup>±</sup> 56 km2.

**Figure 5.** Decadal variations in the total area and number of lakes and reservoirs (area <sup>≥</sup> 1 km2) in the NJW from 1980 to 2015. The error bars represent the 95% confidence level.

#### 3.2.1. Lake Changes

For a further understanding of the temporal changes of lakes, the spatial heterogeneity of lake area changes across the NJW was investigated (Figure 6). Changes in the lake area presented clear variations. During the first period, 1980–1990, most of the expanded lakes were mainly distributed in the southeastern part of the NJW, while the lakes that shrunk were mainly distributed in the central part (Figure 6a). A rapid shrinkage of the total lake area and a decline in lake number occurred during 1990–2000. The larger lake shrinkage occurred in the southern NJW (Figure 6b), while lakes in the Chagan Lake, Yangsha Lake, and Tumuji Lake zones exhibited expanding trends during this period. During 2000–2010, lakes in the NJW with large areal loss were distributed mainly in the Dalong Lake zones, while expanded lakes were identified mainly at the intersection point between the Nenjiang River and Taoerhe River (Figure 6c). During 2010–2015, it is noteworthy that the characteristic of lake expansion was relatively obvious in the NJW compared to the other periods. Lakes with shrinkage were mainly distributed in the southeastern part, whereas lake expansions mostly occurred in the eastern part (Figure 6d).

**Figure 6.** Spatial variations of lake changes in the Nenjiang watershed during different periods. Red and blue dots present lake shrinkage and expansion, respectively, while the dot extent denotes change proportion in the corresponding period. (**a**–**d**) represent the periods 1980–1990, 1990–2000, 2000–2010, and 2010–2015, respectively.

The detailed lake changes in each class are shown in Table 3. From 1980 to 1990, the largest areal decline of lakes was observed for the 50–100 km2 size class. Interestingly, the Dalong Lake, which was in the size class of 50–100 km2, separated into five lakes in the size class of 10–50 km2 during this period. Specifically, 10 lakes with a size of 1–10 km<sup>2</sup> covering a total area of 34 km2 in the central region disappeared in this period. During 1990–2000, lake shrinkage occurred with the total lake area declining by 18%, and 53 lakes (area <sup>≤</sup> 50 km2) vanished during this decade. During 2000–2010, both the total number and area of lakes in the size class of 50–100 km<sup>2</sup> increased, while both the total number and area of lakes in other size classes decreased. From 2010 to 2015, both the total area and number of lakes in the size classes 1–10, 10–50, and 50–100 km2, increased.


**Table 3.** Lake number, area, and changes in different classes between 1980 and 2015.

Note: \* denotes change at a significant level of 0.05, \*\* denotes significant at 0.01.

#### 3.2.2. Reservoir Changes

Figure 7 shows the spatiotemporal patterns of reservoir changes in the NJW during different periods. During 1980–1990, the reservoirs that shrunk were distributed mainly in the central and northeast parts, whereas the expanded reservoirs were mainly distributed in the southern and eastern parts of the NJW (Figure 7a). During 1990–2000, the expanded reservoirs were mainly distributed in the northeast and southern parts, whereas the reservoirs that shrunk were mainly distributed in the central and southern parts of the NJW (Figure 7b). During 2000–2010, significant reservoir expansions were identified across the NJW, while the larger reservoirs that shrunk were mainly distributed in the southern part of the NJW (Figure 7c). After 2010, the reservoirs that shrunk were mainly distributed in the southern and eastern parts. Larger expanded reservoirs were mainly identified in the northern part and the Taoer River basin (Figure 7d).

**Figure 7.** Spatial distribution of reservoir changes in the Nenjiang watershed during different periods. Purple and blue dots present reservoir shrinkage and expansion, respectively, while the dot extent denotes change proportion in the corresponding period. (**a**–**d**) represent the periods 1980–1990, 1990–2000, 2000–2010, and 2010–2015, respectively.

During the period of 1980–1990, the total area of reservoirs in size classes of 10–50 km2 and 50–100 km<sup>2</sup> increased by 121%, with an areal increase from 149 <sup>±</sup> 10 to 328 <sup>±</sup> 16 km2, and by 133% from 52 <sup>±</sup> 14 to 121 <sup>±</sup> 17 km2, respectively. The total areas of reservoirs in the size classes of 1–10 km2 and <sup>&</sup>gt;100 km<sup>2</sup> decreased from 205 <sup>±</sup> 9 to 200 <sup>±</sup> 6 km2 (−2%) and from 513 <sup>±</sup> 20 to 151 <sup>±</sup> 14 km2 (−71%), respectively (Table 4). During 1990–2000, the total areas of reservoirs in the size classes of 1–10 km2 and 50–100 km<sup>2</sup> increased by 23%, with an areal increase from 200 <sup>±</sup> 6 to 247 <sup>±</sup> 5 km2, and by 123% from 121 <sup>±</sup> 17 to 270 <sup>±</sup> 15 km2, respectively. The total area of reservoirs in the size classes of 10–50 km2 decreased from 328 <sup>±</sup> 16 to 263 <sup>±</sup> 8 km2 (−20%). During the period of 2000–2010, the total area of reservoirs increased from 780 <sup>±</sup> 28 to 1086 <sup>±</sup> 58 km2. The total areas of reservoirs in the size classes of 1–10 km2 and 10–50 km2 increased by 30% and 24%, respectively. The total area of reservoirs in the size class of 50–100 km<sup>2</sup> decreased by 50%. During the period of 2010–2015, the total reservoir area increased by 31% with an increase in area from 1086 <sup>±</sup> 58 to 1422 <sup>±</sup> 56 km2. Detailed reservoir changes for these different classes are shown in Table 4.

**Table 4.** Reservoir number, area, and changes in different classes between 1980 and 2015.


Note: \* denotes change at a significant level of 0.05, \*\* denotes significant at 0.01.

#### *3.3. Roles of Climatic Factors and Artificial Variables in Driving Lake Changes*

As shown in Figure 8, the model passed the reliability test, convergent validity test, and discriminant validity test. The dominant climatic factor that influenced lake changes was the MAP (β = 0.66, *p* < 0.001), followed by the ET (β = 0.39, *p* < 0.05) and MAAT (β = 0.15, *p* < 0.1). Agricultural consumption of water had a significant effect on lake changes, suggested by the significant relationship between cropland area and lake area (β = 0.17, *p* < 0.1).

**Figure 8.** Structural model and path coefficient. Arrows show the effect path and direction. Numbers adjacent to arrows are path coefficients, (β). The path coefficients (β) characterize the extent of the effect of the examined variable on the lake changes, while the larger value indicates a strong positive or negative effect. P represents a significant level. \* represents *p*-value < 0.1, \*\* represents *p*-value < 0.05, and \*\*\* represents *p*-value < 0.01. MAP denotes the mean annual precipitation, MAAT represents the mean annual air temperature, ET represents the evapotranspiration, and GDP represents the gross domestic product. The "e" represents the error (residual term) of path analysis of observed variables in the structural model. The numbers on the arrows are the values of the standardized regression weights of the model.

#### **4. Discussion**

Lake and reservoir mapping is affected by the resolution of the used satellite data [26]. We focused only on lakes and reservoirs with area greater than 1 km2 to reduce as much as possible the errors induced by a coarse resolution of 30 and 80 m. Although some images out of the optimal season were used, these data mainly covered the Greater and Lesser Khingan Mountains where few lakes and reservoirs were identified. This did not yield large uncertainties in our analysis. This study integrated OBIA and visual interpretation instead of automatic classification to extract lakes and reservoirs, which ensured data accuracy and effective analysis.

Both climate change and human activities contributed to the changes of lakes in the NJW. On the one hand, the climate over the study area is a semi-arid continental climate. Water supply and output for these lakes dominated with precipitation and evapotranspiration, respectively. SEM analysis revealed that the lake shrinkage in the NJW had a significant correlation with the MAP (β = 0.66, *p* < 0.001) (Figure 8), followed by the ET (β = 0.39, *p* < 0.05), and the MAAT (β = 0.15, *p* < 0.1). This suggests that precipitation had a significant statistical relationship with lake area and potentially had the largest impacts on the lake area. During the investigated 35 years, a warmer climate was identified for the NJW with a significant increase of MAAT (*p* < 0.05) (Figure 9a). A reduced water supply from precipitation and output by increased ET characterized most of the lake shrinkages in the NJW. Specifically, lakes in the NJW showed areal changes with a decrease from 1980 to 2010 and then an increase from 2010 to 2015 (Figure 5), which is consistent with the changed trend of MAP during the 35 years (Figure 9a). The MAP could be regarded as the main climatic factor to explain the lake shrinkage in the NJW from 1980 to 2015 in terms of its decline.

**Figure 9.** Past changes and future projections of climatic factors in the NJW. (**a**) Changes of mean annual air temperature (MAAT), mean annual precipitation (MAP), and evapotranspiration (ET) from 1980 to 2015 and (**b**) changes of future MAAT, MAP, and PET from 2015 to 2050.

Future projections of climate change using the Intergovernmental Panel on Climate Change (IPCC) models (RCP 5) indicate that the warming and drying trend will continue in the NJW (Figure 9b). Both the MAAT and PET show increasing trends with a rate of 0.03 °C yr−<sup>1</sup> and 1.64 mm yr−<sup>1</sup> (*p* < 0.05), while the MAP exhibits a significant decline with a rate of <sup>−</sup>1.63 mm yr−<sup>1</sup> (*<sup>p</sup>* <sup>&</sup>lt; 0.01). If this trend continues, the total area of lakes in the NJW may continue to decrease, and some small lakes will disappear.

This study found that flood events can markedly affect the lake area and number, especially for the small lakes (1–10 km2). According to the hydrological records, the most serious floods during the recent century occurred in 1998 and 2013 [51]. Compared with the total lake area in 2010 (normal flow year), the total lake area in the NJW respectively increased by 2946 km2 in 1998 and 895 km2 in 2013 during the flooding events. It is clear that extreme precipitation events (floods) have accelerated the expansion of the total lake area (Figure 10). In arid and semi-arid regions, a multidimensional view on the prevention and exploitation of floods is required. Lakes and reservoirs have huge water storage capacity, and thus they could serve as hydrological buffers to prevent floods and provide irrigation water for sustainable agricultural development [10].

**Figure 10.** Increased lake area due to floods in 1998 and 2013.

On the other hand, the study area is an important grain production base in China. Agricultural demand for water had potential impacts on lakes. Besides climate change, human activities have imposed marked impacts on lakes and reservoirs [52]. While our study revealed that climatic factors drive striking lake and reservoir changes, human-induced changes should be responded to as quickly as possible. In this study, a clear correlation (β = 0.17, *p* < 0.1) was observed between lake area and crop area (Figure 8). Lake shrinkage caused by artificial variables is mainly attributed to the agricultural water consumption in the NJW [50,53]. Due to the construction of large reservoirs upstream water was therefore impounded upstream, which seriously affected the water supply of lakes in the middle and lower reaches. In particular, the largest reservoir, Nierji Reservoir, constructed over the upstream area of the Nenjiang River exerted a marked influence on downstream lakes. Due to reclamation from lakes and agricultural development during the past 35 years, some lakes in the size class of 1–10 km2 decreased to small lakes with areas smaller than 1 km2. For example, a large area of shallow waters was reclaimed for planting rice (Figure 11). Since 1990, the area of paddy fields in the NJW has increased significantly by 90.26 km<sup>2</sup> from 3.46 to 93.73 km2. The total cropland area increased from 459.42 to 948.18 km2. With the population increase and agricultural land expansion in this region, more and more open water resources have been applied to agricultural irrigation. It is evident that the human impacts on lake and reservoir changes could not be ignored.

**Figure 11.** Image examples for lakes and reservoirs affected by the expanded paddy fields.

Lake shrinkage and disappearance lead to ecological and environmental degradation, such as aggravating the degree of sandstorms and desertification and reducing the number of wild animals [4,54]. To optimize water distribution programs through scientific water conservancy projects is very important. Therefore, appropriate measures need to be implemented by managers to further reduce the decline of lake areas [55], such as enhancing the drainage capacity of reservoirs in drought years but the storage capacity of reservoirs in flood years. In addition, we should control the expansion of paddy fields to relieve water stress and guarantee regional sustainable development in the NJW [56].

#### **5. Conclusions**

In this paper, we established a multi-temporal dataset of lakes and reservoirs in the NJW using long time-series Landsat images from 1980 to 2015 to document their changes on a decadal scale and quantified the contribution degree of MAAT, MAP, ET, cropland, and GDP to change in lake area. A notable decline in the total lake area by 1330 km<sup>2</sup> in the period of 1980 to 2015 was identified, while the lake number decreased contemporaneously. In contrast to lake shrinkage, the total area and number of reservoirs in the NJW experienced continuous increases. We identified 51 newborn reservoirs with

total area of 504 km2 in 2015 compared to 1980. SEM analysis revealed that decrease of the MAP is the dominant factor driving the changes of lakes, followed by ET and MAAT, especially for those of small size. Furthermore, the human impacts on lake and reservoir changes could not be ignored. Timely and appropriate policies and measures are required to reduce lake shrinkage and respond to the degraded environment. The results and analysis in this study are expected to provide guidance for sustainable management of water resource in the NJW.

**Author Contributions:** B.D. and D.M. designed the analytical framework of this study. B.D. performed the data analysis and drafted the manuscript. H.L. and H.X. provided methodological advice. D.M. and Z.W. made major revisions of the manuscript. All authors have read and agreed to the published version of the manuscript.

**Funding:** This study was jointly funded by the National Key Research and Development Program of China (2019YFA0607100, 2016YFC0500408, 2016YFC0500201, and 2016YFA0602301), the National Natural Science Foundation of China (41771383), and the funding from Youth Innovation Promotion Association, Chinese Academy of Sciences (2017277, 2012178) and National Earth System Science Data Center of China (http://northeast.geodata.cn).

**Acknowledgments:** We appreciate the facility that make Landsat 8 OLI images accessible through the USGS. We thank the three anonymous reviewers for the constructive comments and suggestions, which help improve the quality of this manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
