Impacts of Climate Change on a Coastal Wetland from Model Simulation Combining Satellite and Gauge Observations: A Case Study of Jiangsu, China

Abstract

Coastal wetlands are affected by both natural processes and human activities. In the present study, the impacts of natural processes on wetland area variations along the Jiangsu coast in the East China Sea were investigated using a long-term high-resolution numerical model (55 years from 1955 to 2010) validated by satellite and gauge observations. In our 55-year simulation, a 1.33 km² yr⁻¹ decreasing trend of wetland area as a result of global warming was identified. It was found that the wetland area varied depending on the following temporal scales: intra-seasonal, seasonal, interannual, and decadal. Tides and steric sea level changes are responsible for the intra-seasonal and seasonal wetland variations, respectively. The long-term variations are attributable to the El Niño-Southern Oscillation and Pacific Decadal Oscillation. These basin-scale phenomena cause changes in the local wind patterns along the Jiangsu coast and impact the wetland variation on long-term scales.

1. Introduction

Wetlands are one of the most productive environments in the world, and are of significance for biological diversity [1]. As a distinct ecosystem, the coastal wetland links the ocean and land and plays an important role in the local natural environment and social economy. It has been reported that up to one-third of the world’s population lives in the limited coastal areas (~4% of the total land area [2]) and more than 10% of the population lives in coastal areas with terrain elevations of less than 10 m [3,4]. Hence, as an important component of coastal areas, coastal wetlands have become a significant natural resource for humans. Coastal wetland not only provides habitats for different types of wildlife, but also affects sedimentation, and controls storms and floods [5,6,7,8,9]. However, a trend of reduced coastal wetland has been reported globally, decreasing at a rate of 0.949% yr⁻¹ during the 20th and early 21st century [10].

Different mechanisms, including both natural factors and human activities, regulate the variability of the coastal wetland area. In terms of natural factors, the impact of the rise in sea levels induced by global warming has drawn much attention [11,12,13]. A study by Nicholls and Cazenave [14] shows that the low-lying coastal areas of countries in Africa, and South, Southeast and East Asia are vulnerable to flooding as the sea level rises, implying the potential erosion of coastal wetland. In addition, dynamical oceanic and atmospheric processes including tides, typhoons and monsoons can induce land flooding or exposure, which can impact low-lying coastal wetlands as well [8,11]. Human exploration and dredging have eroded a number of river deltas and islands and adversely affected coastal wetlands [15,16,17,18].

Coastal wetlands accounts for about 3% of the total wetland in China according to Niu et al. [19]. The coastline of China is low-lying, and 1.3% of China’s land area (140 million population) has an elevation of below 10 m accounts for. The shoreline and wetlands are thus sensitive to sea level changes and are susceptible to ocean flooding, which requires study due to the dense population. The coastal wetland is demonstrated to have decreased by about 174 km² yr⁻¹ during 1978~2008 [19]. As one of the coastal provinces in China, Jiangsu Province has shallow coastal seas and a very low-lying land elevation (Figure 1), which is the province that has the third largest wetland area in China (accounting for 16% of the total [19]). Jiangsu Province is adjacent to the East China Sea (ECS), and the ECS is subject to complicated and dynamic processes which deeply modulate the coastal wetland variability. Global warming contributed to a rising-sea-level rate of 0.24 cm yr⁻¹ in the ECS during 1980~2002 [20]. Seas along the Jiangsu coast are dominated by semidiurnal tides with tidal amplitudes of up to 2 m [21]. In addition to tides, the sea level variability along the Jiangsu coast is also modulated by the steric effect (such as seasonal solar radiation and Pacific Decadal Oscillation (PDO)) and water mass exchanges with the Kuroshio current and the Taiwan Strait [20]. Meanwhile, the local wind is also greatly affected by the El Niño-Southern Oscillation (ENSO), inducing interannual changes of the sea level in the ECS [22].

The impacts of human activities on wetland area variability are highly complicated and unpredictable. In the present study, we focused on the impacts of natural factors on the wetland area variability in the Jiangsu coastal area in the ECS. The Jiangsu coastal area, which is northwest-oriented, is located north of the Yangzi River estuary. A numerical model of the Jiangsu coast was utilized to investigate the multi-scale temporal variations of the coastal wetland area and to identify the mechanisms regulating the wetland area variability. This work focused on the impacts of physical dynamical processes such as tides and winds, while other factors including biochemical processes, sedimentation, and anthropogenic factors are not considered. The rest of the paper is organized as follows: Section 2 introduces the model used and the validation of the model results based on observations. The numerical results are analyzed in Section 3 and Section 4 presents a summary of the study.

2. Data and Methods

2.1. Model Setup

The simulation was conducted using the Regional Ocean Modeling System (ROMS) [23]. Since this work focuses on coastal wetland area variability, we applied a two-dimensional barotropic model with a wet–dry module [24]. The model domain was defined over the coast of Jiangsu Province with a horizontal resolution of 1 km (Figure 1a). The topography of the model domain was derived using a combination of the following two datasets: the ocean depth was obtained from bathymetry charts, while the land elevation was downloaded from the Chinese Academy of Sciences data center (http://www.resdc.cn/; accessed on 6 September 2020). Both datasets have original resolutions of 1 km.

The model was driven by atmospheric and tidal forces. The atmospheric forces considered here were 6-hour winds from the European Centre for Medium-Range Weather Forecasts with a spatial resolution of 0.125°. The tidal force condition included eight tidal constituents (M₂, K₁, O₁, S₂, N₂, P₁, K₂, Q₁) from TPXO8 with a spatial resolution of 1/30° near the shore. The flux and sea level data along the lateral boundaries of the model domain were obtained from the Simple Ocean Data Assimilation (SODA v2.2.4) dataset, with a spatial resolution of 1/2° (https://iridl.ldeo.columbia.edu/SOURCES/.CARTON-GIESE/.SODA/.v2p2p4/?Set-Language=en; accessed on 29 September 2020). The climatological monthly mean mass transport of the river runoff of the Yangtze River was provided by the Global River Flow and Continental Discharge Dataset [25]. The model was run for 55 years, from 1955 to 2010.

2.2. Model Validation

Three tidal gauges for tidal observations exist along the coast of Jiangsu Province, namely, Binhai, Dafeng, and Lvsi. They are denoted by red dots in Figure 1b. To validate the model simulation, the tidal harmonic constants from the simulation were derived at the gauge locations and compared with those from the observations.

A comparison of the amplitudes and phase lags of the eight main tidal constituents at the three tidal gauges revealed that the model reproduced the tidal fluctuations well (Figure 2). The mean absolute errors (MAEs) for amplitudes were 7.04, 4.65, 4.32, 1.94, 2.14, 1.48, 1.56, and 0.58 cm for M₂, S₂, N₂, K₂, K₁, O₁, P₁, and Q₁, respectively, while the MAEs of the phase lags were 5.87°, 10.77°, 17.32°, 9.34°, 8.77°, 4.50°, 3.18°, and 16.81°, respectively. Due to the location of one of the nodal points of Q₁ near Lvsi Station, the phase lag error appeared to increase for Q₁ at Lvsi Station. Hence, the phase lag error of Q₁ at Lvsi Station was excluded from the MAE calculation. In addition to the evaluation of the tides, the simulated monthly averaged sea surface height was compared with observations at Lvsi Station as well (Figure 3). The simulated sea surface height demonstrated strong seasonal and obvious long-term variations, which were, in general, in accordance with the observations. The MAEs of the tidal harmonic constants and the comparison of monthly averaged sea surface heights indicated the ability of the model to simulate the sea level variation.

Moreover, a qualitative evaluation of the model’s simulation ability was conducted based on satellite images. As observed in Figure 1b, there are shoals along the coast, especially between Dafeng and Lvsi Station. These shoals are very sensitive to sea level changes due to the shallow water depth of the surrounding waters (Figure 4). The changes of the shoals can be directly seen in the satellite images. According to the satellite image captured at 2:15 UTC on 13 January 2009, the shoals were clearly defined, especially between Dafeng and Lvsi Station (Figure 4a). It can be inferred that they were taken at low tide, thus the shoals emerged because of the low sea level. In contrast, as can be seen in the image captured at 2:18 UTC on 6 June 2009, these same shoals were submerged and not identifiable at high tide. In addition, the coastline moves seaward at low tide and landward at high tide, indicating the sensitivity of the coastline to sea level changes.

A similar scenario occurred in the simulation (Figure 5). The upper subplots replicate the variation of land area with time on 13 January 2009. At 2:00 UTC (roughly at the same time as in Figure 4a), the land area had a relatively large value, implying that the images were taken around low tide with a generally lower sea level along the coast. The distribution of the sea surface height at low tide shows distinct shoals, exhibiting a similar pattern to the satellite image. In comparison, at high tide (6:00 UTC 13 January 2009), the sea surface height distribution shows that the shoals were submerged (Figure 5b). The similar pattern between the model results and satellite images demonstrates the capability of the model in simulating wetland area variability.

2.3. Definition of Wetland Variability

Variations of the shoals and the coastline are always accompanied by wetland area variability. From a physical perspective, the coastal wetland can be defined as the area where an exchange between land and ocean occurs. Hence, in this study, the coastal wetland area was calculated as the difference between the minimum land area and the land area at a certain time. The wetland area variability can thus be represented by the land area anomaly from the mean value. In the following sections, the land area anomaly (hereafter called wetland area) was used as an indicator of the wetland variability. To simplify the problem, it is noteworthy that this definition is solely from a physical perspective. Traditionally, a wetland is always tightly related to biochemical processes, which will not be considered in this work.

3. Results

The power density spectrum of the wetland area was calculated to identify the significant temporal scales of the wetland area variability (Figure 6). The results show several significant peaks for different time periods, ranging from hours to up to one year. These peaks suggest that the wetland area variability has a multi-scale temporal variation feature. In general, these peaks can be classified into several bands, namely, high-frequency (less than a half-day), semidiurnal, diurnal, semimonthly, monthly, semiannual, and annual. The interannual (3~10 years) and decadal (>10 years) variations are discussed as well, although no significant temporal scales are observable in these two bands from the power spectrum result.

3.1. Intra-Seasonal Variation

As the dominant temporal scale, the semidiurnal component contributes to more than 60% of the total wetland area variability. The time series in January 2009 shows that the variation amplitude of the semidiurnal component was able to exceed 4000 km², with obvious semimonthly envelopes (red line in Figure 7a). Compared with the superposition of tidal elevations due to the primary semidiurnal tidal constitutes (black line in Figure 7a), the wetland area variation is considered as its opposite, but with synchronous semimonthly envelope variations. As the semidiurnal tides arrive, the land becomes submerged and as the tides ebb, the land emerges. The wetland variation amplitude decreases when the sea surface height variation decreases at the neap tide, and increases when the sea surface height variation increases at the spring tide. Hence, the semidiurnal wetland area variability mainly results from the primary semidiurnal tides.

This is also applicable to diurnal tides, which act as a primary factor in diurnal wetland variation. The amplitude of the diurnal wetland area variability component is one order of magnitude smaller than the semidiurnal component (around 400 km²) (red line in Figure 7b). Similarly, the sum of the tidal elevations of the main primary diurnal constituents has an anti-phase variation with the wetland variation (black line in Figure 7b). The coastal wetland variation observed in the satellite images (Figure 3) is believed to be mainly induced by tides. The area of wetland uncovered and covered by tides is termed as the “intertidal zone”. Intertidal zones exist along coastal areas all over the world. They have been proven to be a significant part of wetlands and are thus a hot topic for oceanographers [26,27].

In shallow and coastal waters, bottom friction induces significant nonlinearity to tides. Higher and lower frequency harmonics are generated, according to trigonometric function formulas. These harmonics are revealed to be important constituents in coastal areas with shallow water [28,29,30]. Since the model domain has shallow water, it can be inferred that different harmonics exist in the model result due to the nonlinear interaction, although only eight primary tides were used as the forcing conditions in the simulation.

The potential periods of harmonics caused by nonlinear interactions between different tides are listed in

Table 1. Periods of harmonics caused by nonlinear interaction of primary tidal constituents.

Tides	M2 (12.4 h)	K1 (23.9 h)	O1 (25.8 h)	S2 (12.0 h)	N2 (12.7 h)	P1 (24.1 h)	K1 (12.0 h)	Q1 (26.9 h)
M2 (12.4 h)	6.2 h, NAN	\	\	\	\	\	\	\
K1 (23.9 h)	8.2 h, 1.1 d	12.0 h, NAN	\	\	\	\	\	\
O1 (25.8 h)	8.4 h, 1.0 d	12.4 h, 13.7 d	12.9 h, NAN	\	\	\	\	\
S2 (12.0 h)	6.1 h, 14.8 d	8.0 h, 1.0 d	8.2 h, 0.9 d	6.0 h, NAN	\	\	\	\
N2 (12.7 h)	6.3 h, 27.6 d	8.3 h, 1.1 d	8.5 h, 0.4 d	6.2 h, 9.6 d	6.3 h, NAN	\	\	\
P1 (24.1 h)	6.2 h, 1.1 d	12.0 h, 182.7 d	12.5 h, 14.8 d	8.0 h, 1.0 d	8.3 h, 1.1 d	12.0 h, NAN	\	\
K1 (12.0 h)	6.1 h, 13.7 d	8.0 h, 1.0 d	6.1 h, 0.5 d	6.0 h, 182.4 d	6.2 h, 9.1 h	8.0 h, 1.0 d	6.0 h, NAN	\
Q1 (26.9 h)	8.5 h, 1.0 d	12.7 h, 9.1 d	8.5 h, 1.2 d	8.3 h, 0.9 d	8.6 h, 1.0 h	12.7 h, 9.6 d	8.3 h, 0.9 d	13.4 h, NAN

. A number of harmonic frequencies are generated by nonlinear interactions between primary constituents, with periods ranging from hours to semiannual cycles. Taking the semidiurnal tides as an example, the multiplication of the two tides, M₂ and S₂, can be divided into the sum of two harmonics with periods of 6.1 h and 14.8 days. The nonlinear interaction between the constituents, K_1, and P₁, can even generate a harmonic period of around 180 days. As a result, not only are higher frequency harmonics with periods of a few hours produced during the interaction, but also lower frequency harmonics with periods including semimonthly, monthly and semiannual cycles. Comparing the harmonic periods in Table 1 and the power spectrum result in Figure 6, the temporal scales of the wetland variations from the power spectrum are in accordance with the harmonic periods from the table. Hence, for periods of the wetland variations including hours, semidiurnal, diurnal, semimonthly, monthly, and semiannual, nonlinear interactions between different primary constituents are suggested to contribute significantly.

It should be noted that the linear superposition of two tidal constituents can lead to an envelope variation with a lower cycle. For example, the linear superposition of M₂ and S₂ can be considered as a multiplication of two trigonometric functions with periods of 0.5 and 29.5 days (Figure 7). However, the envelopes resulting from linear superposition are not reflected in the power spectrum.

3.2. Seasonal Variation

According to the power spectrum in Figure 6, the seasonal frequency is the only remaining significant peak for wetland variations, in addition to short-term wetland variations. The amplitude of the seasonal wetland variation is about 400 km², which is comparable to the diurnal component (Figure 8b). The seasonal wetland area reaches its height in March and decreases to its minimum in September. Undoubtedly, this wetland variation is tightly related to the seasonal variation of the domain-averaged sea surface height (red line in Figure 8a). The sea surface height anomaly is lowest in March (around −0.1 m), with the maximum wetland area and height (around 0.1 m) in September resulting in a minimum wetland area.

The seasonal variations of sea surface height in the ECS are modulated by different factors, including the steric effect due to ocean density change, the inverse barometric effect due to atmospheric pressure change, and Ekman transport due to monsoons [20,31]. According to Figure 3, the model results reproduced the sea surface height variation well compared to the observations, especially for the seasonal variation, despite the fact that the SODA dataset does not consider the inverse barometric effect. Hence, the inverse barometric effect is minimal in the region along the Jiangsu coast. On the other hand, the steric effect due to density change can be decomposed into thermosteric and halosteric components, which can be calculated as,

$${\eta }_{steric}={\eta }_{thermo}+{\eta }_{halo}={{\displaystyle \int }}_{-H}^0(-\frac{1}{\rho }\frac{\partial \rho }{\partial T}\Delta T-\frac{1}{\rho }\frac{\partial \rho }{\partial S}\Delta S)dz$$

(1)

where $\rho$ is the water density, H is water depth, and $\Delta T$ and $\Delta S$ are the temperature and salinity differences relative to the temporally averaged values. The partial derivatives of the density due to temperature and salinity imply the occurrence of thermal expansion and salt compression. Using the SODA dataset, the steric sea surface height can be estimated. The local steric effect is usually very small since the study area has a very shallow water depth. Here, the steric height was estimated as the domain average over 122°E~126°E, 30°N~36°N, assuming that the steric effect in deep water can be transmitted to the coastal area [31,32].

A comparison between the seasonal components of the total sea surface height and steric height is provided (Figure 8a). It can be seen that the total sea surface height is highly consistent with the steric height, not only at the seasonal scale but also at long-term time scales. This consistency indicates that the seasonal variation of the sea surface height can be mostly explained by the steric effect, despite occasional slight biases. The sea surface height increases in summer due to water expansion, leading to a decrease in the wetland area, but decreases in winter due to water contraction causing an increase in the wetland area.

In addition to the steric effect, the role of the strong, along-coast monsoon has been emphasized due to Ekman transport [33]. The along-coast wind blows from northwest to southeast in winter and reverses in summer. As a result, the sea surface height increases due to shoreward Ekman transport in winter but decreases due to offshore Ekman transport. However, the northwesterly wind is usually strongest around January and the southeasterly wind is strongest around July, implying a lag in the wind-induced surface height variation compared with steric height. Given the good consistency between the total and steric heights, it can be inferred that the contribution of monsoons is limited in the study region. A quantitative analysis is provided in the last section.

3.3. Interannual and Decadal Variations

Long-term variations at interannual and decadal scales are important components of wetland variations since these variations can deeply influence the wetland ecosystem. Using the Butterworth bandpass filter, the interannual and decadal components of the wetland variation were extracted and analyzed. Meanwhile, considering the link between the local wind and wetland variability, the annual and interannual along-coast winds were also extracted using the same bandpass filtering method.

For the interannual variation of wetland, no regular or significant periods were observed from the time series of the power spectrum. Similar variational patterns were noted for the interannual wind stress where no significant regularity was observed. However, the correlation coefficient of 0.55, which was calculated between the along-coast wind stress and wetland area, suggests a tight relationship between the two variables, that is, the Ekman transport induced by the interannual along-coast wind was found to be an important factor in modulating the interannual wetland variation (Figure 9a).

As an interannual dynamical phenomenon, the ENSO is primarily defined as the interannual sea surface temperature anomaly in the tropical Pacific Ocean. However, the impact of the ENSO is global rather than local. To determine further mechanisms associated with interannual wetland area variability, a 30-year moving correlation analysis between the wetland area and ENSO index was conducted (Figure 9b). Regions with values exceeding the 95% significant level are represented by white lines. The moving correlation analysis displays a positive correlation coefficient of up to 0.39 with a lag of several months, implying a potential relationship between the wetland variation and the ENSO.

ENSO reportedly induces anticyclonic (cyclonic) wind anomalies during El Niño (La Niña) events over the Philippine Sea, although this is a phenomenon that was observed over the tropical Pacific Ocean [34,35]. Considering the relative location of the ECS, this anticyclonic (cyclonic) wind anomaly leads to a southerly (northerly) wind over the ECS during El Niño (La Niña). As a result, an along-coast wind component is induced, which finally leads to offshore (shoreward) Ekman transport and an increase (a decrease) in the wetland area during El Niño (La Niña). The specific process by which ENSO induces wetland variability via wind anomaly is reflected in the positive correlation coefficient calculated above.

For the decadal wetland variation, a more pronounced regularity was found. Generally, the wetland variation was positive before the 1970s and negative after the 1970s, indicating the wetland area was larger before the 1970s and subsequently shrunk (Figure 9c). Consistent with the wetland variation, the decadal wind anomaly blew southeasterly before the 1970s and northwesterly after that. The correlation coefficient calculated between wind and wetland area was 0.62, implying a significant role of the decadal along-coast wind in modulating decadal wetland variability via Ekman transport.

The PDO, which is an important climate variability mode in the Pacific Ocean, occurs over the mid-latitude Pacific Ocean, with a decadal temperature anomaly dipole over the west and east Pacific. The temperature anomaly dipole is accompanied by wind anomalies, which can influence the local wind spatial structure over the ECS [36]. According to the PDO index, one can see that the PDO index was negative before the 1970s and positive after the 1970s, indicating a PDO phase shift from a cool phase to a warm phase before and after the 1970s, respectively.

Compared to the PDO index in Figure 9c, the decadal wetland variation demonstrated an opposite trend to the PDO index with a correlation coefficient of −0.42. During the PDO cool phase, the west Pacific water became warmer and an anomalous east wind was generated. As a result, the along-coast wind anomaly component blew southeasterly and induced offshore transport near the Jiangsu coast. An increase in wetland area thus occurred. In contrast, during the warm phase of the PDO, the west wind strengthened and the along-coast wind anomaly component blew northwesterly, leading to a decrease in the wetland area.

3.4. Long-Term Variation Trend

In addition to the periodic variations, the wetland area shows a linearly decreasing trend according to the long-term (>1 year) time series (Figure 10b). The linearly fitted decreasing trend was approximately 1.33 km² yr⁻¹, indicating a decrease in wetland area over the five decades. The wetland area decreased by more than 70 km² from 1955 to 2010. This decreasing trend implies an increase in the sea surface height. The linearly fitted increasing trend of the domain-averaged sea surface height was 0.56 mm yr⁻¹ (red line in Figure 10a). For comparison, the linear trend of the steric height was 0.47 mm yr⁻¹ (green line in Figure 10a). Thus, the variation trend of the simulated sea surface height is mainly caused by the steric effect. The upward trend of the steric height is suggested to be a result of global warming.

However, an obvious bias can be found between the variation trends of the simulated and observed sea surface heights, despite the good consistency illustrated in Figure 3. The observed sea surface height had a sharper increasing trend, compared with the simulated one (Figure 10a). The calculated increasing trend of the observed sea surface height was 1.84 mm yr⁻¹, which is more than 3 times larger than the simulated result (black line in Figure 10a). The effects of global warming on rising sea levels occur via two aspects, including mass redistribution via modulating ocean currents and melting sea ice, and water expansion due to higher water temperatures. It has been demonstrated that the ocean mass redistribution is the main factor contributing to rising sea levels in the ECS [37]. However, the SODA dataset only reproduces rising sea levels caused by the steric effect, but greatly underestimates the part played by mass redistribution. So, it can be inferred that the decreasing trend of the wetland area was underestimated in the model result as well.

4. Discussion

4.1. Wind Contribution to Diurnal and Seasonal Wetland Area Variabilities

The contribution of wind is discussed in terms of interannual and decadal wetland variations. Two obvious energy peaks at diurnal and seasonal periods can be seen from the wind spectrum (Figure 11a). These two peaks contribute to the diurnal and seasonal wetland variations, respectively. The results provided above demonstrated the significant roles of the diurnal tides and the steric effect based on a correlation analysis, but the extent of wind contribution requires further clarification.

Under a constant wind drag at the initial condition, the ocean reaches a steady state due to the balance of the wind stress, Coriolis force, and pressure gradient force. Considering the shallow water of the study area, the balance can be described by the following equation,

$$f\overrightarrow{k}\times \overrightarrow{V}=gH \nabla \eta +\frac{\overrightarrow{\tau }}{\rho }$$

(2)

where f is the Coriolis parameter, g = 9.8 m s⁻² is the gravity acceleration, $\rho =$ 1024 kg m⁻³ is the seawater density, $\overrightarrow{V}$ is the depth-averaged horizontal velocity, H is the water depth, $\nabla \eta$ is the spatial gradient of the sea surface height and $\overrightarrow{\tau }$ is the wind stress. The equation is not solved here. The sea surface height due to the along-coast wind can be estimated using a scale analysis. It can be inferred that the mass transport perpendicular to the coast should be balanced by the transport parallel to the coast at a steady state. Hence, these terms should have the same magnitude, which can be expressed as,

$$\frac{\Delta \eta }{\Delta x}~\frac{\tau }{gH\rho }$$

(3)

Finally, the magnitude of the sea surface height induced by wind drag can be estimated by,

$$\frac{\Delta \eta }{\Delta x}~\frac{\tau }{gH\rho }$$

(4)

The diurnal wind component was extracted from the wind speed spectra using the Butterworth bandpass filter (frequency band: 0.5~1.5 days). Climatological averages of the zonal and meridional wind speeds during both daytime and nighttime were then calculated (blue vectors in Figure 11b,c). The wind distribution over the ECS during daytime and nighttime revealed an anticyclonic and cyclonic pattern, respectively. The corresponding Ekman transport was calculated based on the finite-depth Ekman transport equation (red vectors in Figure 11b,c) [38]. The anticyclonic wind during daytime is roughly along-coast with a southeast orientation. The induced offshore Ekman transport results in a sea level decrease and expansion of the wetland area. In contrast, a shoreward Ekman transport is induced by the northwest wind which leads to an increase in the sea level and wetland area shrinking. Based on Equation (4), the sea surface height induced by sea–land breeze can be estimated. The magnitude of the sea–land breeze stress was about 5 × 10⁻⁴ N m⁻². For the ECS shelf, the mean water depth was H~50 m and the zonal range was $\Delta x ~$500 km. As a result, the magnitude of the sea surface height caused by the sea–land breeze was approximately 5 × 10⁻⁴ m. This magnitude is three orders smaller than that from diurnal tides. So, compared to the diurnal tides, the contribution of the diurnal sea–land breeze to the wetland area variation is negligible.

Similarly, the sea surface height anomaly caused by seasonal monsoon can be estimated. A magnitude of 0.03 m was obtained, given that the magnitude of the stress was around 0.03 N m⁻². Compared with the total seasonal sea surface height, this magnitude accounts for about 30% of the variation. Hence, the seasonal monsoon plays a non-negligible contribution to the wetland seasonal variation. In winter, the East Asian winter monsoon (EAWM) blows from land to the ocean. As a result, a northwesterly (negative value) along-coast wind dominates, inducing shoreward Ekman transport (Figure 11d). The shoreward transport leads to a rise in the sea level near the shore, causing the wetland area to decrease. On the other hand, the East Asian summer monsoon (EASM) blows from the ocean to land in summer, turning the along-coast wind into southeasterly, thus inducing offshore transport (Figure 11e). As a result, the sea level decreases, and the wetland area increases.

4.2. The Contribution of Sea Level Rise

As has been mentioned, the model results underestimate the effect of sea level rise on wetland variation due to the bias of the SODA data. Based on the estimated sea level rise at the Lvsi tidal gauge, the contribution of sea level rise can be corrected. Taking the bathymetry of the model domain as the mean state during the simulation, the sea surface height increases by adding the sea level rise trend, and the land is submerged. The wetland area can be instantaneously calculated according to the submerged land area. The estimation indicates that the observed sea level rise of 1.8 mm yr⁻¹ can cause a decreasing trend of 6.68 km² yr⁻¹ in the wetland area, which is 5 times larger than that for the model simulation. The decreasing trend implies that the wetland along the Jiangsu coast lost an area of more than 335 km² over the last 50 years. The ecosystem should be greatly affected by this decreasing wetland area. In addition to the direct contribution of rising sea levels associated with global warming, an implicit effect from global warming should also be noted. As the sea level rises, the amplitudes of the multi-scale wetland variabilities are potentially modulated since the wetland area is deeply dependent on the topography. Flat-land topography behind the wetland indicates a larger wetland area when the sea level rises and vice versa. However, this implicit effect is nonlinear and difficult to quantify, and was therefore not analyzed in this study.

Meanwhile, the wetland area demonstrated a linearly decreasing trend of about 1.33 km² yr⁻¹. The analysis indicates that this decreasing trend of the simulated sea surface height was mainly caused by the steric effect, which underestimates the decreasing trend due to the sea level rise. Based on the estimated sea level rise trend at the Lvsi tidal gauge, the contribution of rising sea levels is estimated as 6.68 km² yr⁻¹ in the wetland area, which is 5 times larger than that from the model simulation.

4.3. Limitations of the Model

Despite the fact that the model reproduces the multi-scale variabilities of the wetland along the Jiangsu coast, several limitations of the model need to be noted. First, the model is a simple 2D model and the effect of vertical stratification was not included in the simulation. Even though the study region has relatively shallow water depths, the stratification in summer is strong [39], which potentially modulates the local circulation and sea level. As a result, the 2D assumption may introduce biases in the simulation of wetland variability, especially in summer. Second, the model uses invariant topography during the 55-year simulation period. However, it should be mentioned that the coastline and the coastal topography of Jiangsu experienced remarkable changes due to land reclamation and erosion/deposition. These changes potentially affect the wetland area and its variabilities. Lastly, the spatial resolution of the model is another aspect limiting the accuracy of the model. Both the ocean bathymetry and land elevation datasets have spatial resolutions of 1 km, and include publicly available data offering the highest spatial resolution. For the purpose of consistency, the model resolution was taken as 1 km. One can speculate that the higher the model resolution is, the more precise the simulated wetland will be.

5. Conclusions

Based on a 55-year simulation from a barotropic model validated by satellite and gauge observations, the multi-scale variability of the wetland area along the Jiangsu coast was investigated. According to the results of the analysis, the wetland along the Jiangsu coast shows multi-scale temporal variabilities which are modulated by several different dynamical processes. Although other factors such as biochemical processes, sedimentation, and anthropogenic factors also contribute to the wetland area variability, only the dynamic process effects are considered in this work. As a more predictable factor, studying the dynamic process effects on the wetland area variability is of crucial importance to improving the wetland area prediction.

For short-term variations (from hours to half a year), tides are the dominant contributors to wetland variability, which account for more than 65% of wetland area variability. The eight primary tidal constituents considered in the model mainly modulate the semi-diurnal and diurnal variations of the wetland area. Meanwhile, the harmonics generated from nonlinear interactions of primary constituents due to bottom friction give rise to short-term wetland variations.

For seasonal variation, the sea level change due to the steric effect is the main factor influencing the wetland area. The local steric effect leads to a sea level variation of about 0.2 m, which causes the maximum wetland area in March and the minimum area in September. For long-term variations, wind plays a significant role in modulating the wetland area variability on interannual and decadal scales. According to the analysis presented, the long-term wetland area variability is tightly regulated by the ENSO and the PDO. The correlation coefficients were between 0.39 and −0.42 for wetland area variability and ENSO, wetland area variability, and PDO, respectively. Given that the wetland area variability is closely related to the Ekman transport induced by the along-coast wind, these dynamical processes can influence the strength and direction of the along-coast wind which then modulate the long-term wetland area variability.

As the protection of the coastal environment and wetland has attracted increased attention globally, the impacts of anthropogenic factors on wetland variability are expected to decrease in the future. Hence, wetland variability will be more dependent on natural factors and become more predictable in the future. The results in this work demonstrate the feasibility of using regional ocean models to simulate wetland variations due to ocean dynamic processes. To establish a more accurate model for wetland prediction, the coupling of different models to include biochemical processes and sedimentation processes needs to be studied.