3.1 Spatial changes of YRD wetland during 1986–2022
A total of 37 spatial distribution maps depicting the typical salt marsh vegetation in the YRD wetland for the period 1986–2022 were generated. The overall classification accuracy exceeded 96%, with a minimum Kappa coefficient value of 0.94, indicating that our classification results are reliable and highly accurate (Fig. 2c). Over the course of these 37 years, we present eight maps that are representative (Fig. 3). It is evident that in the past four decades, the YRD wetland has undergone substantial changes in terms of morphology, riverine pathways, and typical salt marsh vegetation.
In the early stages of the formation of the YRD, 1986, the wetland exhibited a relatively regular and open elliptical shape, with the typical salt marsh vegetation, PA, primarily distributed along the inland tidal flats on both sides of the Yellow River. Subsequently, as river erosion intensified, by 1990, the shape of the wetland became increasingly pointed, resembling a bird's beak.
Later, in 1996, a watershed year for the YRD wetland, the long-established southeastward outlet channel was abandoned, giving way to a new northeastward channel. Initially, the new channel had a smaller mouth, but after nearly 10 years of development, by 2006, the northeastward channel gradually expanded, and formed one of the two major branches of the YRD wetland, alongside the southeastward channel, making it a remarkable phenomenon among estuarine wetlands worldwide. Furthermore, in 2006, extensive seasonal water ponds started to emerge in the southern part of the YRD wetland.
In 2009, the fragmented patches of SA began to be observed, particularly close to the coastal tidal flat. Additionally, human activities such as aquaculture ponds, cultivated land, and salt fields gradually became evident during the same year. By 2012, the northeastern branch of the new flow path had further extended towards the ocean. And near the coastal areas, apparent patches of SA were observed, indicating its successful colonization. Simultaneously, larger seasonal water ponds appeared in the northern part of the YRD wetland.
Over the subsequent decade, the SA flourished and expanded its dominion, progressively extending its territory from the coastal tidal flats towards the inland areas. Meanwhile, the SA, itself has become another dominant species, soaring to prominence in this ecosystem. Furthermore, in 2022, it started to challenge the established dominant species in the region, the PA, which traditionally thrived along the Yellow River, that gradually encroached upon its stronghold. Moreover, the morphological transformation of the YRD wetland during this period witnessed the continuous strengthening of the northeastern flow path, while the southeastern flow path gradually declined.
3.2 Patch area changes of YRD wetland during 1986–2022
Through visualizing the temporal changes of patches area in the YRD wetland over the past 37 years, we identified three stages of landscape pattern evolution: rapid expansion stage, gradual decline stage, and bioinvasion stage (Fig. 4a). From 1986 to 1997 was the first stage of wetland expansion, during the 12 years, both wetland area and tidal flat area showed a rapid increase trend (R2 = 0.91 for wetland area and tidal flat area), with a 70.45% increase in the wetland area, reaching the peak of the entire time series in 1997 at 387.62 km2, which laid the foundation for the wetland pattern. While during this period, the dominant salt marsh vegetation, PA, was reduced by 25%. The second stage is the wetland decline stage, from 1997 to 2009, the wetland and tidal flat areas showed a slow decline trend (R2 = 0.33, Pearson's r = -0.57 for wetland; R2 = 0.58, Pearson's r = 0.76 for tidal flat area), with the area of wetland reduced by 21.33%, PA reduced by 15.96%. The third stage is the bioinvasion stage, SA demonstrated an unstoppable trend of rapid increase (R2 = 0.93), with an area increase of 68 times relative to 2009, expanding at an average rate of 3.44 km2 per year, approaching the area of PA by 2021 and becoming another dominant species in the YRD wetland. In the same term, the PA area increased by 17.62%, however, extensive tidal flat was being encroached upon by SA, resulting in its area being reduced by 23.02%.
Making a general observation of the transition matrix of patches area across the three stages (Fig. 4b), it can draw four findings: (1) the YRD wetland had undergone three different strengths of patch changes; (2) in the first stage, although there were significant variations in patch flows, the types of changes were relatively homogeneous; (3) in the second stage, there was a noticeable increase in human activities, leading to an increase in the transfer flows and the density of the flow network; (4) in the third stage, although the single flow of patch transfers decreased, the transfer network became more complex.
In the first stage, the area of tidal flats doubled in size, with 161.77 km2 contributed by wetland and 51.57 km2 contributed by PA, indicating a certain degree of degradation in the PA ecosystem. In the second stage, an area of 35.19 km2 of tidal flats transformed into wetlands, signifying increased seawater intrusion. Additionally, a small portion of the tidal flat was converted into a pond, salt pan, and cultivated land, indicating the manifestation of human activities. In the third stage, an area of 24.79 km2 of tidal flat and 19.27 km2 of water body were encroached upon by SA, while 38.38 km2 of tidal flat was occupied by water body. Furthermore, an area of 19.55 km2 of the water body transformed into a tidal flat, accounting for the tidal flat had experienced both erosion and deposition counterbalance in that time.
Compared to wetland patch area, the landscape indices can comprehensively reflect information about landscape composition and spatial configuration, representing a highly condensed characterization of landscape pattern features. The temporal variations of landscape indices from 1986 to 2022 displayed a dual-stage or tri-stage pattern, with notable differences observed in the timing and trends of various indices (Fig. 4c). LSI, SPLIT, LJI, and PD, these four landscape indices, which represent landscape shape complexity, fragmentation, interspersion and juxtaposition, and patch density, respectively, exhibiting two distinct stages with the year 2001 serving as the turning point. From 1986 to 2001, these indices displayed a decreasing trend (0.33 ≤ R² ≤ 0.57). However, after 2001, they exhibited a significant upward trend (0.37 ≤ R² ≤ 0.92). It is indicated that before 2001, the wetland tended to exhibit regularity and uniformity, while afterwards, it became increasingly fragmented with more complex shapes. AI reflects the connectivity between each patch type, in contrast, which exhibits a trend of initial increase followed by a decrease, with the year 2001 as a turning point.
LPI, representing the proportion of the largest patch in the wetland, showed a decreasing trend in 1997 (R² = 0.78) and remained relatively high thereafter. In the YRD wetland, the largest patch was the water body, which was replaced by the tidal flat during the wetland expansion stage. Subsequently, as the wetland gradually receded, the water body area increased, resulting in a consistent trend between LPI and the water body. SHEI exhibited a significant increasing trend before 1988 (R² = 0.79). A higher value of the index indicates a closer proportion of different patch types in the landscape and a higher level of evenness. This suggests that the evenness of the landscape continuously improved until 1998, after which it remained at a relatively low level. PAFRAC, reflecting the complexity of patch boundaries, gradually decreased before 1995 but experienced a sudden increase after 2006, maintaining a high level thereafter. Unlike AI, CONTAG reflects the overall connectivity of the landscape. Before 1998, the wetland exhibited poor connectivity with a declining trend (R² = 0.67). Subsequently, connectivity remained relatively high but also showed a decreasing trend.
3.3 Hydrological forces of wetland patch area
Stepwise linear regression results (Table 2) indicate that the dependent variables of wetland, tidal flat, and PA passed the model diagnostics (Pi<0.001, 1.35 < DW < 1.59, at the significance level of 0.05, n = 36 and k = 3; VIF3 − 4 < 10). SA is likely limited by sample size, with high VIF and overfitting indicated by adjusted R2, concluding that the model fails. Based on the model, the area of wetland, tidal flat, and PA in the YRD wetland are significantly influenced by SE + and RU+, explaining 61.5%, 75.7% and 63.8% of their variations respectively, while RU and SE are excluded from the model. Specifically, the wetland area and tidal flat area are positively correlated with the SE+ (β31 = 2.61, β32 = 4.62), while the RU + shows a weak negative effect (β41 = -0.023, β42 = -0.063). However, for PA, as a hydrophilic vegetation, the situation is the opposite as the RU + has a positive effect on its area (β43 = 0.008).
Table 2
Stepwise multiple linear regression results of YRD wetland patch area with hydrological variables. RU, SE+, RU + respectively represent runoff, cumulative runoff, and cumulative sediment transport. * means P ≤ 0.05, ** means p ≤ 0.01, and *** means P ≤ 0.001 (the same as below).
Patch area (Xi, i=1,2,3,4)
|
Adjusted Ri2
|
DWi
|
Fi
|
Constant
|
RU
|
SE+
|
RU+
|
\({{\beta }}_{0i}\)
|
\({{\beta }}_{1i}\)
|
VIF1
|
\({{\beta }}_{3i}\)
|
VIF3
|
\({{\beta }}_{4i}\)
|
VIF4
|
Wetland
|
0.615
|
1.463
|
28.944
|
223.984***
|
|
|
2.611***
|
7.367
|
-0.023***
|
7.367
|
Tidal flat
|
0.757
|
1.516
|
55.488
|
141.277***
|
|
|
4.620***
|
7.367
|
-0.063***
|
7.367
|
PA
|
0.638
|
1.374
|
31.788
|
87.259***
|
|
|
-0.868***
|
7.367
|
0.008***
|
7.367
|
SA
|
0.989
|
2.262
|
282.100
|
235.371**
|
− .025*
|
3.898
|
-6.738***
|
141.503
|
0.066***
|
119.275
|
3.4 Hydrological forces of landscape pattern
The correlation between RU + and SE + with the landscape indices demonstrated a higher overall significance compared to RU and SE (Fig. 5a1-5a2). Specifically, significant simple linear correlations (P < 0.001) were found between RU and AI, IJI, LSI, and PD, with Pearson correlation coefficients of -0.58, 0.61, 0.57, and 0.65, respectively (Fig. 5a1). Furthermore, RU + demonstrated significant (P < 0.001) linear correlations with AI, SPLIT, IJI, LSI, and PD, with Pearson correlation coefficients of -0.61, -0.53, 0.68, 0.63, and 0.68, respectively. Additionally, significant monotonic correlations (P < 0.001) were observed between RU + and SHEI, CONTAG, with Spearman correlation coefficients of -0.57 and 0.55, respectively (Fig. 5a2).
In contrast, the correlation between SE and the landscape pattern of YRD wetland was overall not substantial. Among the landscape indices, SHEI and CONTAG exhibited highly significant (P < 0.001) correlations with respective Pearson correlation coefficients of 0.47 and − 0.50. SE + displayed significant (P < 0.001) linear correlations with SPLIT, CONTAG, and LPI, with Pearson correlation coefficients of 0.58, 0.53, and − 0.57, respectively (Fig. 5a1). Furthermore, SE + exhibited highly significant (P < 0.001) monotonic correlations with AI, SHEI, IJI, LSI, and PD, as indicated by respective Spearman correlation coefficients of -0.58, -0.57, 0.56, 0.60, and 0.62 (Fig. 5a2).
The indices that satisfied Condition 1 included PD, LSI, CONTAG, IJI, SPLIT, and AI. LPI and SHEI satisfied Condition 2, while the index that satisfied Condition 3 was PAFRAC.
The patch density (PD) in the YRD wetland increases with an increase in RU and RU+ (Fig. 5b1). PD is more sensitive to RU, but the impact of RU + on its variation is slightly greater (R2 = 0.46). The two fitted curves of the landscape pattern index (LPI) are shown (Fig. 5b2). As RU + and SE + increase, the area of WA representing the largest proportion of landscape patches decreases, indicating an increase in the total wetland area. In other words, RU + and SE + promote the development of the YRD wetland, with both variables explaining 46% of the variation in LPI. The shape complexity, represented by the landscape shape index (LSI), significantly increases with the influence of RU and RU+ (Fig. 5b3). Similar to PD, LSI is more sensitive to RU, but the regression equation fitted with RU + performs better (R2 = 0.39).
Among all the landscape indices, PAFRAC is the only index that is not sensitive to any of the hydrological features. This is shown in the scatter plot (Fig. 5b4), where most of the points are concentrated in the upper left quadrant, indicating a positive correlation with SE and RU, but without statistical significance. CONTAG, which reflects landscape connectivity, decreases with an increase in SE (R2 = 0.30) (Fig. 5b5). However, it increases with an increase in SE+ (R2 = 0.28), and it is more sensitive to SE. Similar to PD and LSI, RU and RU + have a significant promoting effect on IJI (Fig. 5b6), with explanatory values of 37% and 47%, respectively. This indicates that increased runoff intensifies the fragmentation of wetland patches.
SPLIT is significantly influenced by both RU + and SE+ (Fig. 5b7), and it increases with an increase in both variables. SHEI is a special landscape index that reflects landscape uniformity. The results indicate that when RU + at the Lijin station is within 2040.63×108 m3 and SE + is within 53.94×108 t, SHEI increases significantly with RU + and SE+ (Fig. 5b8), with explanatory values of 82% and 79%, respectively. However, beyond these values, the scatter plot becomes disordered and remains at a lower level. In contrast to landscape fragmentation, AI, which reflects landscape spatial aggregation, decreases with an increase in RU (R2 = 0.34) and RU+ (R2 = 0.37) (Fig. 5b9).