Interannual and Seasonal Variations of Hydrological Connectivity in a Large Shallow Wetland of North China Estimated from Landsat 8 Images

Abstract

Hydrological connectivity is an important characteristic of wetlands that maintains the stability and functions of an ecosystem. This study investigates the temporal variations of hydrological connectivity and their driving mechanism in Baiyangdian Lake, a large shallow wetland in North China, using a time series of open water surface area data derived from 36 Landsat 8 multispectral images from 2013–2019 and in situ measured water level data. Water area classification was implemented using the Google Earth Engine. Six commonly used indexes for extracting water surface data from satellite images were compared and the best performing index was selected for the water classification. A composite hydrological connectivity index computed from open water area data derived from Landsat 8 images was developed based on several landscape pattern indices and applied to Baiyangdian Lake. The results show that, reflectance in the near-infrared band is the most accurate index for water classification with >98% overall accuracy because of its sensitivity to different land cover types. The slopes of the best-fit linear relationships between the computed hydrological connectivity and observed water level show high variability between years. In most years, hydrological connectivity generally increases when water levels increase, with an average R² of 0.88. The spatial distribution of emergent plants also varies year to year owing to interannual variations of the climate and hydrological regime. This presents a possible explanation for the variations in the annual relationship between hydrological connectivity and water level. For a given water level, the hydrological connectivity is generally higher in spring than summer and autumn. This can be explained by the fact that the drag force exerted by emergent plants, which reduces water flow, is smaller than that for summer and autumn owing to seasonal variations in the phenological characteristics of emergent plants. Our study reveals that both interannual and seasonal variations in the hydrological connectivity of Baiyangdian Lake are related to the growth of emergent plants, which occupy a large portion of the lake area. Proper vegetation management may therefore improve hydrological connectivity in this wetland.

1. Introduction

Wetlands are critical and vulnerable ecosystems that provide essential ecological services and account for 47% of the global ecosystem [1]. The service values of wetlands are strongly influenced by hydrological processes such as interception, infiltration, evapotranspiration, and surface water flow movement [2]. The interaction between a floodplain wetland and its adjacent river owing to flood pulses is the major control over floodplain biota [3]. Within a wetland, water-mediated connections contribute to the exchange of material, energy, and information with their surrounding environment, which plays an important role in maintaining the stability and functions of the ecosystem [4]. Hydrologic connectivity is defined as the water-mediated transport of matter, energy, and organisms within or between elements of the hydrologic cycle [5]. Hydrological connectivity has been identified as an important factor that influences water quality and provides suitable habitats for aquatic organisms in wetlands [6,7,8,9]. Climate change and human activities pose serious threats to the hydrological connectivity of wetlands, which can degrade wetland functions and lead to a series of environment problems [10,11]. An evaluation of hydrological connectivity changes is therefore important for understanding the driving mechanism of wetland ecosystem changes and formulating effective strategies to improve the health of wetland ecosystems [5]. However, there remains a lack of consensus regarding how to measure hydrological connectivity [12].

Tracking water surface dynamics is essential for understanding hydrological connectivity variations in wetlands. Satellite remote sensing technology has shown great potential in recent decades for tracing changes in terrestrial surface water from space [13,14,15]. Satellite images of optical bands have been widely used to trace changes in land surface features with the advantages of easy accessibility and diversity in terms of spatial and temporal resolutions [16,17,18]. They have been applied to evaluate the hydrologic connectivity of floodplains [19], freshwater deltas [20], and oxbow lakes [21]. Classification algorithms have been developed to discriminate between waterbody and non-waterbody pixels [22,23,24]. However, surface water areas beneath cloud and vegetation cover cannot be detected from optical images alone. Synthetic aperture radar imagery has been useful for detecting water surfaces in wetlands because of its strong penetration capability [25,26], whereas the image accuracy can be influenced by uneven topography and surface roughness. Light Detection and Ranging (LiDAR) technology also has potential to detect the extent of flood inundation under vegetation canopy [27,28,29]; however, high costs limit its routine use in studies that cover large spatial scales [30,31]. To more effectively monitor water surface dynamics in wetlands, the integration of optical images and ancillary data, such as terrain information from LiDAR or in situ water level data, has been shown to be effective [29,32,33]. Water level is also an important indicator to quantify hydrological connectivity [34,35].

Many shallow wetland lakes in arid or semiarid areas face the problem of reduced water resource storage owing to climate change and/or reservoir construction in upstream rivers, which can seriously damage wetland ecosystems and degrade their service functions [36,37]. Environmental water allocation to these lakes is an important approach to alleviate this damage and restore damaged wetland ecosystems. However, such water allocation changes the natural water level pulse, which is a primary factor that influences connectivity processes [38]. Exploring the dynamics of water surface area and hydrological connectivity on the landscape scale is important to better understand the influence of environmental water allocation on wetland ecosystems. Baiyangdian Lake (BYDL) is the largest natural freshwater wetland in the North China Plain, and is a typical wetland-type lake that highly relies on such water allocation. Previous studies have explored surface water changes and monitored landscape patterns in BYDL. Zhuang et al. [36] detected water body changes between 1974 and 2007 and analyzed the impacts of human activities in BYDL. Li et al. [37] used Landsat images from four periods to monitor BYDL landscape patterns. Liu et al. [39] simulated the water area under different water levels using the MIKE21 model and calculated landscape connectivity to determine the minimum ecological lake level in BYDL. Wang et al. [40] extracted BYDL water surface area data and analyzed water area changes in spring from 1984 to 2018. However, there remains a lack of knowledge regarding intra-annual or seasonal variations in surface water area as well as hydrological connectivity, which are vital for understanding the responses of wetland ecosystems to environmental water allocation processes that strongly influence water level dynamics. Numerous studies have shown intense interactions between hydrological connectivity and riparian vegetation [41,42]. In BYDL, previous studies have indicated that the spatial extent of open water area, which is the major water flow medium, is strongly affected by the spatial distribution of closed riparian canopies of emergent plants [43,44,45]. Dense roots in the water lead to low hydrological connectivity under these emergent plants (mainly reeds), whose growth is driven by their seasonal phenological characteristics. Detailed study of these seasonal variations is therefore required to better understand the influence of riparian vegetation on hydrological connectivity.

The objective of this study is to obtain knowledge regarding the interannual and seasonal variations of hydrological connectivity in BYDL and their driving mechanisms. Landsat 8 satellite images from 2013–2019 were used to detect the open water surface area dynamics of BYDL using the Google Earth Engine (GEE). Commonly used index methods for detecting water surface area from optical images were first applied and their accuracies were evaluated. The open water surface areas of BYDL were then extracted from all of the Landsat images using the most accurate method. A comprehensive index method based on several indices to describe landscape patterns was established to quantify hydrological connectivity in the wetland from water surface area data detected from the Landsat 8 images. The changes of water area and hydrological connectivity of the wetland in different seasons were analyzed, as well as their relationship with the in situ observed water level data. Previous studies in BYDL did not discuss the seasonal variation of hydrological connectivity in detail. Our study examines the interannual and seasonal variations of hydrological connectivity while considering the simultaneous influence of water level and vegetation growth. The results provide valuable insight for decision-making in environmental water allocation projects and ecosystem restoration for typical shallow lakes and other lakes in semiarid and arid regions worldwide.

2. Materials and Methods

2.1. Study Area and Data

BYDL is the largest shallow wetland in North China (Figure 1, 115°45′–116°05′E, 38°44′–38°59′N) with an area of approximately 320 km² and average water depth of 2.84 m. BYDL plays an important role in maintaining biological diversity and regulating climate in the Beijing-Tianjin-Hebei region. There are approximately 143 lakes of different sizes in the wetland and a few rural villages with more than 100,000 residents, which aggravates the plight of environmental water shortage and wetland fragmentation. The inflow of BYDL has gradually decreased owing to upstream reservoir operations and higher human water use, including agricultural irrigation and industrial and household water consumption. A series of environmental water allocation actions have been adopted in recent years to improve the health condition of the wetland ecosystem [46]. The wetland has also drawn increasing attention with the establishment of the Xiong’an New Area, which receives some noncapital functions from China’s capital, Beijing.

The region has a semihumid continental monsoon climate with a mean annual precipitation and evaporation of 563.9 and 1369 mm, respectively. More than half of the water area grows a variety of seasonal aquatic vegetation including emergent plants such as Phragmites australis, which cover the water area beneath [47]. In this study, the water body area observed by satellite images refers to the uncovered water surface and is hereafter called the open water area.

Landsat 8 multispectra images (spatial resolution: 30 m) from 2013 to 2019 were used to map the water surface. We queried and analyzed the surface reflectance products of Landsat 8 on the GEE, which has been atmospherically corrected. GEE is a cloud platform with a multi-petabyte catalog of satellite imagery and planetary-scale analysis capabilities, and has been widely used in time series to monitor and analyze terrestrial water [14,48,49]. As BYDL freezes in the winter season, our analysis in this study focuses on the three other seasons, i.e., spring, summer, and autumn. Using the GEE platform, 36 cloud-free Landsat 8 scenes of BYDL for the three seasons were selected and processed to extract the open water area. The daily surface water level data of the Shifangyuan Hydrological Station (Figure 1) operated by the Ministry of Water Resources was used in this study.

2.2. Methodology

2.2.1. Surface Open Water Mapping

Low reflectance in the near-infrared (NIR) band is a distinct characteristic of water surfaces compared with other land cover types. The detection of surface water from remote sensing data is most effective over NIR and visible wavelengths [50]. The use of water indexes computed from the reflectance of two or more spectral bands is an easy and effective method to identify water area from satellite images [16]. Five commonly used indexes were used and compared in this study. Water has higher green band reflectance within the visible band and lower NIR reflectance. The normalized difference of green and NIR reflectance (NDWI) has been established to enhance the presence of water [22]. The short-wave infrared band (SWIR) is less sensitive to sediments within the water than NIR, and the modified normalized difference water index MNDWI has been created for water detection [51]. Feyisa et al. [52] established a multiple-band index (AWEI_nsh) to accurately distinguish water from a shadowed surface with a stable threshold. Fisher et al. [53] created the WI2015 index that combines reflectance data from five bands, whose empirical parameters are determined by statistical analysis from more than 400,000 training pixels of different types of water and nonwater surfaces from Landsat 8 images.

We used the reflectance of six bands of Landsat 8 images (

Table 1. Bands of Landsat 8 image used in this study.

Band	Description	Wavelength (μm)
$B_2$	Blue Band	0.45–0.51
$B_3$	Green Band	0.53–0.59
$B_4$	Red Band	0.64–0.67
$B_5$	Near-infrared Band	0.85–0.88
$B_6$	Shortwave infrared Band 1	1.57–1.65
$B_7$	Shortwave infrared Band 2	2.11–2.29

) to compute the five indexes:

$$NDWI=\frac{B_3-B_5}{B_3+B_5}$$

(1)

$$MNDWI=\frac{B_3-B_6}{B_3+B_6}$$

(2)

$$AWE{I}_{nsh}=4\times \left(B_3-B_6\right)-\left(0.25\times B_5+2.75\times B_7\right)$$

(3)

$$AWE{I}_{sh}=B_2+2.5\times B_3-1.5\times \left(B_5+B_6\right)-0.25\times B_7$$

(4)

$$WI2015=1.7204+171\times B_3+3\times B_4-70\times B_5-45\times B_6-71\times B_7$$

(5)

where B₂, B₃, B₄, B₅, B₆, and B₇ are the surface reflectance of the blue band (wavelength: 0.45–0.51 μm), green band (wavelength: 0.53–0.59 μm), red band (wavelength: 0.64–0.67 μm), NIR band (wavelength: 0.85–0.88 μm), and shortwave infrared bands 1 (wavelength: 1.57–1.65 μm) and 2 (wavelength: 2.11–2.29 μm) of the Landsat 8 images, respectively.

The threshold of the above-mentioned water index methods must be determined to divide all of the pixels into water and nonwater areas. Two steps were carried out to determine the most feasible threshold. In the first step, box plots were drawn of the computed indexes for each land cover type. If the index ranges for open water and other landcover types did not overlap or only slightly overlap, water bodies were considered to be well distinguished. The land use of BYDL was divided into five types: open water; emergent plants; urban; paddy field; dryland. We chose 18 images of the different seasons from 2013 to 2019 as sampling images. For each land use type in each image, 70–100 sampling pixels were selected based on visual interpretation and index values were computed for each pixel. A similar box plot was also drawn for the grayscale images of the NIR band reflectance, because of its high sensitivity for water area detection [16]. A total of six box plots were drawn from the sampling images, five for the composite water indexes and one for the NIR band reflectance, which can also be considered as a water index. An assessment of the box plots allowed the index that most effectively discriminated water and nonwater areas to be determined based on the above-mentioned principle and that index was used for further study.

In the second step, the water area classification threshold was determined for each image considering the temporal variation of open water reflectance. If the histogram of the selected index for all pixels in an image showed a distinct bimodal distribution, the threshold was determined with the 2-mode method [54], which uses the index value with the lowest frequency between the two histogram peaks as the threshold. The Otsu method [55] was used in the other cases, in which we selected the index value that maximizes the variance between water pixels and nonwater pixels as the threshold. These methods were applied to all of the images to detect the water surface area.

To evaluate the classification accuracy, we chose seven images in three seasons (spring, summer and autumn) with different water levels, excluding the sampling images. Two hundred water and nonwater verification pixels from each of the seven images were identified and combined with high-resolution Google Earth images through visual interpretation. These pixels were nearly evenly distributed in the study area to maximize the degree of objectivity of the evaluation. The results of the classification using six indexes with a segmentation threshold derived from the box plot were compared. The 2-mode method and Otsu method were then applied to the seven images for further accuracy evaluation. The confusion metrics [56] and overall accuracy (OA) were calculated to assess the classification accuracy. The OA measures the proportion of the number of correctly classified pixels out to the total number of sampling pixels.

2.2.2. Assessment of Hydrological Connectivity in BYDL

The hydrological connectivity considered in this study is defined as the mobility of water flow and its carried materials in the wetland. Several indices describing landscape patterns were adopted to develop a composite index for quantifying the overall hydrological connectivity of BYDL based on water patches detected from the Landsat 8 images. The index was built from four aspects: shape of the water patches; distance between water patches; spatial distribution pattern of the water patches (aggregation); water patch fragmentation. A total of seven indexes were used to describe the four aspects of hydrological connectivity, as shown in

Table 2. Evaluation system for the wetland hydrological connectivity of Baiyangdian Lake (BYDL).

Aspect	Index	Interpretations
Shape	Related circumscribing circle (C1)	Movement efficiency within a water patch.
Distance	Euclidean nearest neighbor distance (C2)	Mean distance between water patches.
	Probability of connectivity (C3)	Possibility of connection between patches.
Aggregation	Shannon’s evenness index (C4)	Evenness of water patch distribution in the landscape.
	Aggregation index (C5)	Like adjacency among water patches
Fragmentation	Splitting index (C6)	Fragmentation degree of water patches
	Average area (C7)	Average water patch area.

, were computed as follows.

The related circumscribing circle index (C1) [57] was used to describe the closeness of the water patch shape to a circle. Patch shape influences the material transfer efficiency [58]. The shape of the circle processes the highest movement efficiency within a patch [59]. The index measures how close the water surface is to a circle and is calculated by the following equation:

$$\mathrm{C}1=\frac{1}{\mathrm{n}}{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}\left(1-\frac{{}_{{{\displaystyle \mathrm{a}}}_{\mathrm{i}}}}{{\mathrm{a}}_{\mathrm{i}}^{\mathrm{s}}}\right)}$$

(6)

where ${{\displaystyle \mathrm{a}}}_{\mathrm{i}}$ is the area of water patch i, ${\mathrm{a}}_{\mathrm{i}}^{\mathrm{s}}$ is the area of the smallest circumscribing circle around patch i, and n is the number of water patches. C1 is not influenced by patch size and equals 0 for circular patches and approaches 1 for linear patches.

The Euclidean nearest-neighbor distance index (C2) is defined using simple Euclidean geometry as the shortest distance between one water patch and its nearest neighboring water patch [60]. Lower values demonstrate stronger connectivity conditions. C2 can be calculated by the following equation:

$$\mathrm{C}2=\frac{1}{\mathrm{n}}{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{h}}_{\mathrm{i}}}$$

(7)

where ${\mathrm{h}}_{\mathrm{i}}$ is the distance from water patch i to its nearest neighboring water patch. C2 is larger than 0 and the minimum value is constrained by the cell size.

The probability of connectivity index (C3) [61] quantifies the possibility of connectivity between patches, which decreases with distance. C3 can be calculated by the following equation:

$$\mathrm{C}3=\frac{{\displaystyle {\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\displaystyle {\sum }_{\mathrm{j}=1}^{\mathrm{n}}{\mathrm{p}}_{\mathrm{i}\mathrm{j}}^{*}\cdot {\mathrm{a}}_{\mathrm{i}}\cdot {\mathrm{a}}_{\mathrm{j}}}}}{{\mathrm{A}}_{\mathrm{L}}^2}$$

(8)

$${\mathrm{p}}_{\mathrm{i}\mathrm{j}}^{}={\mathrm{e}}^{-\mathrm{k}\cdot {\mathrm{d}}^{\mathrm{i}\mathrm{j}}}$$

(9)

where ${\mathrm{p}}_{\mathrm{i},\mathrm{j}}^{*}$ is the maximum product probability of all possible paths between patch i and j, ${\mathrm{a}}_{\mathrm{i}}$ and ${\mathrm{a}}_{\mathrm{j}}$ are the area of water patch i and j, ${\mathrm{A}}_{\mathrm{L}}^{}$ is the total landscape area (both water and nonwater areas); ${\mathrm{d}}^{\mathrm{i}\mathrm{j}}$ is the edge-to-edge interpatch distance, and $\mathrm{k}$ is a constant based on the user-specified threshold distance. The distance threshold is set to 2000 m according to previous studies of BYDL [39]. C3’s value ranges from 0 to 1 and increases with increasing connectivity probability.

Shannon’s evenness index (C4) [62] describes the evenness of the water patch distribution in the landscape. An even distribution of water patches promotes the homogenization of physicochemical characteristics and increases the connectivity conditions in a wetland [63]. C4 is calculated as:

$$\mathrm{C}4=-\frac{{\displaystyle {\sum }_{\mathrm{i}=1}^2{\mathrm{P}}_{\mathrm{i}}}\cdot \mathrm{l}\mathrm{n}{\mathrm{P}}_{\mathrm{i}}}{\mathrm{l}\mathrm{n}2}$$

(10)

where ${\mathrm{P}}_{\mathrm{i}}$ is the proportion of the area of the ith patch class to the total lake area. BYDL is considered to be composed of water patches and nonwater patches when computing the index. C4 equals 1 when the distribution of area among different patch type is perfectly even and decreases to 0 as the distribution of different patch types becomes increasingly uneven.

The aggregation index (C5) [64] quantifies the aggregation level using the like adjacency among patches of the same type. Energy, nutrients and material exchanges are influenced by landscape aggregation [64,65]. C5 can be calculated by the following equation:

$$\mathrm{C}5=\frac{{\mathrm{g}}_{\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}}}{\mathrm{m}\mathrm{a}\mathrm{x}\to {\mathrm{g}}_{\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}}}$$

(11)

where ${\mathrm{g}}_{\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}}$ is the number of like adjacencies (joins) between water patch pixels and $\max \to {\mathrm{g}}_{\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}}$ is the maximum number of like adjacencies (joins) between water patch pixels. C5 is larger than 1, increases with increasing landscape aggregation, and equals 100 when the landscape aggregates into one patch.

The splitting index (C6) [66] measures the fragmentation degree of the same type of patch. It is sensitive to structural differences of landscape patterns and can be interpreted as the number of patches with a constant patch size by dividing the landscape into C6 number of patches. C6 can be calculated by the following equation:

$$\mathrm{C}6=\frac{{\mathrm{A}}_{\mathrm{L}}{}^2}{{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{a}}_{\mathrm{i}}^2}}$$

(12)

where ${\mathrm{a}}_{\mathrm{i}}^{}$ is the area of patch j and ${\mathrm{A}}_{\mathrm{L}}$ is the total landscape area. C6 equals 1 when the landscape contains one patch and increases with decreasing water patch area and increasing subdivision into smaller patches.

The average area (C7) is a fundamental measure of subdivision. It reflects the average conditions of patches and is sensitive to small patches. C7 is calculated by the following equation:

$$\mathrm{C}7=\frac{{\mathrm{A}}_{\mathrm{w}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{r}}}{\mathrm{n}}$$

(13)

where ${\mathrm{A}}_{\mathrm{water}}$ is the total area of the water patches and $\mathrm{n}$ is the total number of water patches.

The seven indexes were computed from the water area derived from Landsat 8 images. C1, C2, and C4–C7 were computed using FRAGSTATS [60] and C3 was calculated using Conefor software [67]. The Analytic Hierarchy Process (AHP) was applied to determine the weights of each indicator and compute a composite index that condenses all of the seven indicators of hydrological connectivity. We conducted pairwise comparisons of the importance of hydrological connectivity representation among the seven indexes to form a reciprocal decision matrix, which was solved to acquire the weight of each index. C3, C4, C5, and C7 show a positive relationship with hydrological connectivity, whereas C1, C2, and C6 show a negative relationship with hydrological connectivity. The values of all indexes were normalized to 0–1 by the following equation.

For indexes with positive relationship:

$$Z_i=\frac{X_i-X_{\min }}{X_{\max }-X_{\min }}$$

(14)

For indexes with negative relationship:

$$Z_i=\frac{X_{\max }-X_i}{X_{\max }-X_{\min }}$$

(15)

where Z_i is the normalized value of ith index and X_i, X_max, and X_min are the unnormalized, maximum, and minimum values of ith index, respectively.

The comprehensive index, called BYDLCO, was then computed for each image as follows:

$$\begin{array}{l}BYDL\mathrm{C}O=0.078\times \mathrm{C}1+0.13\times \mathrm{C}2+0.39\times \mathrm{C}3+0.1005\times \mathrm{C}4+0.1005\times \mathrm{C}5\\ +0.1005\times \mathrm{C}6+0.1005\times \mathrm{C}7\end{array}$$

(16)

3. Results

3.1. Comparison of Open Water Surface Extraction Methods

Box plots of six indices for different landcover types based on the chosen sample points are shown in Figure 2. All of the indices can distinguish open water surfaces from dryland, emergent plants, and urban. However, the ranges of water and paddy field overlap to some extent, except for the NIR band reflectance index. From the boxplots, the minimum values of NDWI, MNDWI, AWEI_nsh, AWEI_sh, and WI2015 and the maximum NIR band reflectance value for the sampled water area were applied as the thresholds to distinguish open water surface pixels from nonwater pixels to the seven images representing different water levels and seasons.

The OA values of the six indices for the seven images and their averages and standard deviations are listed in

Table 3. Overall accuracies of different open water surface classification methods.

Date	Water Level (m)	Water Indexes
		NDWI	MNDWI	AWEIsh	AWEInsh	WI2015	NIR	NIR	NIR
		Threshold for Open Water Surface
		>−0.049	>0.251	>0.005	>0.024	>3.219	<0.007	Based on Otsu Method	Based on 2-Mode Method
22 August 2015	7.68	0.9650	0.9775	0.9475	0.9775	0.9475	0.9925	0.9950	0.9875
28 September 2017	7.90	0.9775	0.9825	0.9650	0.9825	0.9600	0.9875	0.9900	0.9900
30 June 2019	8.21	0.9750	0.9700	0.9650	0.9700	0.9625	0.9925	0.9875	0.9925
4 March 2017	8.47	0.9725	0.9675	0.9700	0.9675	0.9700	0.9750	0.9750	0.9750
28 November 2016	8.50	0.9850	0.9725	0.9825	0.9775	0.9850	0.9825	0.9800	0.9900
3 October 2013	8.82	0.9825	0.9825	0.9850	0.9850	0.9850	0.9350	0.9775	0.9850
3 December 2018	8.69	0.9875	0.6000	0.8775	0.6075	0.8325	0.9950	0.9925	0.9925
Average		0.9779	0.9218	0.9561	0.9239	0.9489	0.9800	0.9854	0.9875
Standard deviation		0.0078	0.1420	0.0368	0.1397	0.0531	0.0210	0.0078	0.0061

. An OA value of 1 means that all of the sampled pixels of open water area were correctly detected. A low OA standard deviation among the different satellite images is also expected, which indicates the stability of the classification method under different hydrological conditions. The classification method that uses the NIR band reflectance demonstrates the highest average accuracy and lowest accuracy variation among the seven images, followed by NDWI, which uses the differences in green band and NIR band reflectance. The 2-mode and Otsu methods were then used for classification based on the NIR band reflectance. Compared with using a fixed threshold obtained from the box plots, the two methods perform better by adjusting the threshold for each image. We therefore used the NIR band reflectance to classify the water area in all of the satellite images. We adopted the 2-mode method when the histogram showed notable bimodal features, otherwise we used the threshold obtained from the Otsu method.

3.2. Variations in Open Water Surface with Water Level and Season

Figure 3a shows the correlation of the open water surface extracted from all of the Landsat 8 images from 2013–2019 and corresponding in situ observed water levels. A scattered and weakly positive relationship is observed. When the data are grouped by year, the positive correlation becomes considerably stronger in most years (Figure 3b–h). The R² is higher than 0.6 for all years except 2018, for which the number of satellite observations is very limited due to cloud cover. The slope of the best-fit linear relation varies from year to year. These features indicate large interannual variations in the relationship between open water surface area and water level.

The relationship between open water area and water level is also assessed in different seasons. Figure 4 shows that the positive correlations for spring, summer, and autumn are all higher than the correlation for all of the data (Figure 3a). This indicates that the pattern and magnitude of open water area variation with water level differ between spring, summer, and autumn. In all three seasons, the open water area generally increases with increasing water level, but the increase is not linear for all seasons and a certain degree of fluctuation is observed. The water area is highest when the water level reaches the upper limit of its range, whereas peaks of open water areas are also found when the water level is between 8.0 and 8.2 m. This implies, that the topographic features are unique in this elevation range, compared with other ranges, and further examination is required. For a given water level in most of the elevation ranges, the open water area is considerably higher in spring than in summer and autumn, the latter of which show a similar relationship between open water area and water level. This indicates that a certain type of landcover occupies a portion of the open water area in the spring season.

We examined the spatial distributions of the water area detected from satellite images in different seasons with the same water levels to constrain the seasonal variations in this relationship. Figure 5 depicts the open water areas at water levels of 7.95 and 8.50 m in spring, summer, and autumn. The water area in spring is substantially higher than that in summer and autumn, especially in the northern part of BYDL. For the water level, the spatial distribution of the open water area is similar for both summer and autumn. This can be explained by the fact that during spring, reeds are in the early growth stage and the biomass is low. After maturing in the summer, the area of reed land reaches a maximum and does not change with time for the remainder of the growth season. When the water level increases, the open water area increases significantly in summer and autumn. However, this is not the case in spring, when the water area only increases slightly when water level rises from 7.99 to 8.47 m. One possible reason is that the latter image was taken in early March when some part of the lake was still frozen. Lake ice has high reflectance in the NIR band, which causes the current classification method to underestimate open water area.

3.3. Temporal Variation in BYDL Hydrological Connectivity

The relationship between the hydrological connectivity composite index, BYDLCO, and water level for the entire period and each of the individual years is shown in Figure 6. Similar to the open water surface area, the relationship for the entire period shows a high level dispersion, whereas strong positive correlations are detected for each individual year, except 2018. The slope of the best-fit line varies from 0.1633 in 2015 to 1.5001 in 2014, which indicates that the relationship between BYDLCO and water level exhibits high interannual variability.

The BYDLCO values were also grouped by season and their relationship with water level is shown in Figure 7. For a similar given water level, the hydrological connectivity is higher in spring than the other two seasons for most of the water level range, similar to that observed for the open water surface area. There is a resemblance of a relationship for summer and autumn in the range that satellite observations are available in both seasons, with that of autumn being slightly higher than that in summer. For each season, the hydrological connectivity generally increases with increasing water level. A certain level of fluctuation occurs, which can be attributed to the interannual variation of the relationship.

The seasonal relationships between water level and each of the seven indexes composing BYDLCO were also analyzed (Figure 8). The index values are generally highest in the spring for a given water level that is the same, followed by the autumn value, except for C1 (Figure 8a). In most cases, each index value increases with increasing water level with different degrees of fluctuations. Interestingly, the fluctuation is highest for most indexes when the water level is in the range of 8.0–8.2 m, followed by a rapid rise and decent with further increasing water level.

The interpretations of each index provide insights regarding the variations of BYDLCO. The shape factor of C1 (Figure 8a) indicates how easily matter is exchanged within an individual water patch. The relationship of C1 and water level shows the highest fluctuation in autumn. The indexes related to mean distance (C2) and the possibility of a direct path between water patches (C3) are shown in Figure 8b,c, respectively. For most water level ranges, the two spring index values are higher than those for summer and autumn. However, their variations patterns with changing water level differ in spring. The properties of aggregation are quantified by C4 and C5, which focus on the distribution evenness and adjacency of the water patches, respectively. The seasonal relationships between these two indexes and water level are quite similar (Figure 8d,e): both C6 and C7 quantify the fragmentation level of the water patches. C6 puts more weight on relatively small water patches. Their relationships with water level (Figure 8f,g) differ substantially, indicating that the water patch size is uneven and that the existence of small patches cannot be ignored. The C6 value is close to 1 over the entire water level range. The degree of fragmentation level of water patches is lower in spring and does not change significantly with increasing water level.

4. Discussion

4.1. Accuracy of Extracting Open Water Area in BYDL from Landsat 8 Images

The water indexes combine reflectance values from different wavelengths and are formulated to maximize the reflectance differences between water and nonwater areas for classification. Our study shows that the accuracy of the water area classification is higher when using the NIR band reflectance index, and the variance is lower than the other five indexes using information from multiple bands. The index that performs second best, NDWI, also intensively uses NIR band reflectance. This can be explained by the fact that the reflectance differences between land cover types considered in this study are highest in the NIR band (Figure 9). For medium spatial resolution like Landsat 8, Fisher et al. [53] found that no single index performed better in all of the evaluated water bodies and depended on the composition of the land and water types. Huang et al. [16] also indicated that differences in sediments and other optical active constituents within the water can influence the applicability of certain water indexes to certain water bodies. Zhang et al. [68] compared the accuracy of commonly used water indexes in an open water body at five sites in North America, Australia, North Africa, and East Asia with different climates and found that NDWI, AWEI_nsh, and AWEI_sh have limitations when extracting turbid water. From three sites in the Yangtze River Basin, Yellow River Basin, and Haihe River Basin in China, Jiang et al. [23] concluded that AWEI_nsh and MNDWI are more vulnerable to shadow pixels. Our results show that for Landsat 8 images, the reflectance of the NIR band is the most effective water index in BYDL, which has a moderate level of eutrophication [69,70]. Han and Cui [71] found that water transparency has high spatial heterogeneity in BYDL that must be considered to better classify open water areas in future studies. Indeed, the reflectance characteristics of inland surface water bodies are complex and highly influenced by local conditions. A full evaluation of the commonly used water index methods, such as that performed here, is therefore necessary to determine the best index for a particular water body.

The water surface area and corresponding estimated hydrological connectivity in the water level range of 8.0–8.2 m show significantly high fluctuations for spring, summer, and autumn in all of the years, which implies high water surface classification errors over this range. Figure 10 shows the relationship between lake level and total lake water surface area in 10-cm increments (including open water area and surface of water under emergent plants) derived from a field survey [72]. The water area variation is approximately 10 km² when water level is lower than 8.0 m. The total water area increases significantly in the water level range of 8.0–8.3 m, which indicates that the lake topography changes remarkably and that the Landsat 8 images and water index methods used here do not trace the changes in open water areas with sufficient accuracy over this water level range. Many small dikes and sand bars exist in BYDL [73] due to aquaculture, rice farming, and the deposition of materials particles via inflow from upstream rivers. When the water level rises from 8.0 to 8.3 m, the above-mentioned landforms change from being exposed above the water surface to being submerged by the lake water. This explains the sharp increase of open water area over this range. A similar phenomenon [74] is found in some segments of the lower Mekong River, where sandbars exist in the river channel. The increment of water area variation decreases with increasing water level. The number of pixels that mix water and land is higher for the 8.0–8.3 m water level range in the Landsat 8 images, compared with other water level ranges, which increases the error in the water area classification. This error can be reduced by using pixel unmixing schemes (e.g., [75,76]) or by classifying water area based on satellite images with spatial resolutions higher than that in Landsat 8 images. High-resolution DEM data derived from a LiDAR system can also be introduced to quantify elevation variations under complex topographic conditions in water bodies to better the understand hydrological connectivity (e.g., Hudson et al. [77]).

4.2. Interannual and Seasonal Variations in Hydrological Connectivity of BYDL

Although the number of satellite observations is relatively low, our results clearly indicate that hydrological connectivity increases with increasing lake water level. This is intuitive because isolated water patches will connect with each other when the water level increases. The interannual and seasonal variations in the relationship between water level and hydrological connectivity are worthy of further attention.

The differences in extraction results may be affected from water surface classification errors, but more importantly, they are caused by differences in the spatial coverage of emergent plants, which occupy large areas of water surface in BYDL. The water area under the canopy of such plants cannot be detected from optical images, such as Landsat 8 images used in this study. The dense roots of such plants can obstruct water flow by increasing the local hydraulic roughness [78] and lowering the hydrological connectivity of the entire BDYL. This is the most important assumption made to assess hydrological connectivity using open water surface area data. Strong positive relationships between hydrological connectivity and water level were found in most years, but the slope of the best-fit linear equations vary between years, thus indicating that the spatial distribution and total area of emergent plants undergoes large interannual variations, which are likely influenced by the climate and lake hydrological regime [31,79].

For the same water level, hydrological connectivity is mostly higher in spring than summer and autumn, with minor differences between the latter two seasons. This can be attributed to the phenological characteristics of emergent plants. In spring, emergent plants are in the germination and leaf development stage. Their stems are thin and short compared with those in summer and autumn, and their canopy cannot be detected from satellite images. The drag force exerted by emergent plants, which reduces water flow [80,81], is therefore lower in spring than summer and autumn. Correspondingly, their influence on water flow is not apparent. After maturation in the summer season, their canopy and stems in the water do not significantly change before harvesting in early winter [82] and the obstruction to water flow is the strongest. This is consistent with the results of Liu et al. [73] who mapped hydrological connectivity using an interferometric synthetic aperture radar-based methodology and reported that 40% of permanent and 33% of conditional barriers to hydrological connectivity are related to reed land. Our study also reveals that the seasonal variation of hydrological connectivity may be connected with emergent plant phenological characteristics. The influence of vegetation growth in temperate wetlands where emergent plants occupy large portions of area must be examined to better understand the driving mechanisms of hydrological connectivity variations. Proper vegetation management may improve hydrological connectivity in such wetlands.

5. Conclusions

Hydrological connectivity is of great ecological importance to wetland systems. Although numerous studies have addressed hydrological connectivity in wetlands, the interannual and seasonal variations of hydrological connectivity and their driving mechanism remain poorly understood. This study established a composite index of hydrological connectivity based on the extent of open water derived from Landsat 8 images in Baiyangdian Lake and investigated its temporal dynamic changes with changing water levels. Reflectance in NIR band is the most accurate and robust index for detecting water compared with other indexes, including NDWI, MNDWI, AWEI_nsh, AWEI_sh, and WI2015. The relationship between extracted water area and water level demonstrates significant interannual variability, which may result from interannual growth differences of emergent plants. Hydrological connectivity generally increases with increasing water level. For a given water level, hydrological connectivity is mostly higher in spring than in summer and autumn largely because of seasonal variations in the phenological characteristics of emergent plants. These findings provide insight regarding the impact of water level regimes and phenological characteristics of emergent plants on the variation of hydrological connectivity in BYDL and guidance for decision-making in the management of this important wetland in North China.