Introduction

Tidal flats, defined as unvegetated mud or sand flats (distinct from mangroves and salt marshes) with regular tidal inundation1, are one of the most extensive coastal ecosystems worldwide, providing essential ecosystem services such as storm surge protection2,3, carbon sequestration4,5,6, and serving as critical habitats for large populations of waterfowls5,7. With the recent advancements in remote sensing, global-scale tidal flat dynamics have been mapped in detail1,8. Strikingly, approximately 16% of the global tidal flat area was lost during 1984–2016. 66% of the area losses were attributed to indirect drivers, rather than direct human impact through conversion to other land uses8, although in some countries such conversion was conducted at large scales (e.g., 240–400 km2 per year during 1950–2010 in China9). Indirect drivers include both natural coastal processes such as wave erosion, and environmental changes influenced by human activities such as altered sediment availability. They often operate at large scales and originate far from sites with emerging dynamics.

Previous studies have demonstrated that tidal range, wave height, and sediment availability can all influence tidal flat morphology at a local scale10,11,12,13,14, but which factors would dominate over a large spatial scale and how they lead to temporal changes remain unclear. For instance, sediment availability is widely regarded as a key factor influencing tidal flat morphology7,15, and there have been dramatic alterations in global suspended sediment flux to the coast16,17. In the global hydrologic north (north of ~20 °N), a 49% reduction in river sediment flux has been observed due to damming, while in the global hydrologic south (south of ~20°S), river sediment flux has increased by approximately 41% since the 1980s, primarily due to mining in the upstream areas17. These global changes in sediment flux are expected to have a profound impact on large-scale tidal flat ecosystems, but their impact and the underlying mechanisms are yet to be revealed to provide guidance for management.

China possesses the 2nd largest tidal flat area in the world, with a variety of tidal ranges, wave heights, and sediment availability1,18,19. Additionally, China’s tidal flat ecosystem has experienced some of the most rapid loss globally (at 2.6% per year)1,7,8. These facts make China’s tidal flats an ideal model system to study the large-scale tidal flat morphology. In the current study, we track the spatio-temporal variation of tidal flat slope and width20, which are available from global datasets1,21. In total, we included 2538 tidal flat transects along China’s coastline for width analysis from the global tidal flat maps and 1620 transects for slope analysis from the Foreshore Assessment using Space Technology (FAST) database. These transects are selected at 1 km resolution, covering open-coast and estuaries, but those affected by direct human impact or where the width or slope information cannot be extracted were excluded from analysis (see “Methods”).

In this work, by combining tidal flat morphological data with environmental factors (i.e., hydrodynamic forcing and sediment availability), we identify the main drivers influencing spatio-temporal variation of tidal flat morphology. We further find different responses of sediment-rich and sediment-starved tidal flats to the same level of sediment availability change, driven by their tendency to restore morphological equilibrium. To assess the universality of these findings, we extended our analysis to three other countries with extensive tidal flat areas but significantly lower sediment availability compared to China, i.e., Indonesia, the USA, and Australia. Finally, the observed diverse morphological responses were explained by a model based on dynamic equilibrium theory10,11,20,22,23,24.

Results

Spatial variation and drivers in tidal flat morphology

Remote sensing observations show that tidal flat width ranges from 30 m to 2.6 × 104 m in 2016 (Fig. 1a). The median value of tidal flat width is 536 m. A large group of the tidal flats (34.3%) has a width of 30–300 m, whereas another large group of the tidal flats (32.4%) is wider than 1000 m (Supplementary Fig. 1). The widest tidal flats are found in Jiangsu (JS) province with a median width of 3997 m, while the smallest median width is observed in the northern Shandong Peninsula, being 65 m (Fig. 1a). To compare the tidal flats in estuaries and on open coasts, we selected 5 major estuaries in China, namely the Changjiang River Estuary, the Yellow River Estuary, the Pearl River Estuary, the Minjiang River Estuary, and the Liaohe River Estuary (Supplementary Fig. 2)25 for analysis. Notably, the median flat width in these estuaries (807 m, N = 275) is greater than those on open coasts (512 m, N = 2263), which is likely attributed to the higher median SSC for estuarine tidal flats (71.0 mg/L) than their open coast counterparts (44.0 mg/L).

Fig. 1: Distribution of tidal flats width and slope along the coast of China.
figure 1

a Tidal flat width distribution. b Tidal flat slope distribution. Coastal provinces of China are listed in the maps, namely Liaoning (LN), Hebei (HB), Shandong (SD), Jiangsu (JS), Zhejiang (ZJ), Fujian (FJ), Guangdong (GD), Guangxi (GX), Hainan (HN) and Taiwan (TW). Red boxes on the maps indicate regions with large widths or slopes, and blue boxes indicate regions with minimal values. The values in the figure represent the median values of data in the box. The variation of tidal flat width and slope along latitude is displayed on their right panel. The shaded areas represent the 25th and the 75th percentile values. c The relationship between tidal flat width and median slope in each group. The number N represents the number of data points. The upper whisker extends to the 75th percentile plus 1.5 times the interquartile range, and the lower whisker extends to the 25th percentile minus 1.5 times the interquartile range. Source data are provided as a Source Data file.

Full size image

The slope of tidal flats in China ranges from 3.6 × 10−5 to 1.2 × 10−1 (Fig. 1b), and the median slope is 7.2 × 10−3. The median slope of tidal flats in the five major estuaries (4.1 × 10−3, N = 153) is gentler than that on open coasts (7.5 × 10−3, N = 1467). The tidal flats in southern Fujian have the highest median slope in the country, i.e., 1.5 × 10−2 (Fig. 1b), which is influenced by the large tidal range (3.6 m) and a low SSC level (29 mg/L). In contrast, those in southern Jiangsu show the lowest slope i.e., 6.7 × 10−4, influenced by the exceptionally high SSC (110 mg/L) and a moderate tidal range (2.93 m). Furthermore, based on the national-scale tidal flat width and slope data, it is clear that the tidal flat slope becomes gentler as the tidal flat width increases (Fig. 1c).

To identify the main drivers shaping the tidal flat morphology, we analyze the distribution of width (Fig. 2a, c, e) and slope (Fig. 2b, d, f) in conjunction with local wave height, tidal range, and suspended sediment concentration (SSC) (details in Methods). The results show that tidal flat width exhibits a positive correlation with SSC (R = 0.71, P < 0.01) and tidal range (R = 0.44, P < 0.01) but a weaker negative correlation with wave height (R = −0.20, P < 0.05, Fig. 2e). A stepwise linear regression analysis indicates that the combination of SSC (80% relative importance) and tidal range (20% relative importance) best explains the distribution of tidal flat width (R2 = 0.61) (Fig. 2g, Eq. 3), in which the effect of spatial autocorrelation is excluded (See “Methods”). Furthermore, the width of tidal flats with flood dominance is markedly greater than that of those with ebb dominance (Supplementary Fig. 4), which is in line with previous studies20.

Fig. 2: Relations between the spatial distribution of tidal flat width and slope with environmental drivers.
figure 2

a, b Tidal flat width and slope versus SSC. c, d Tidal flat width and slope versus tidal range. e, f Tidal flat width and slope versus wave height. The number displayed at the base of each bar represents the number of data points of tidal flat width (average of 2014−2016) and slope (average of 1997-2017). Source data are provided as a Source Data file. g, h Linear regression models (solid lines) that best explain the distribution of width and slope with 95% confidence intervals (shaded areas), in which data are transformed to ensure normal distribution. The effect of spatial autocorrelation is excluded through the application of a spatial error model (details of data processing are included in Methods).

Full size image

The tidal flat slope exhibits a positive correlation with the tidal range (R = 0.34, P < 0.01), but a negative correlation with SSC (R = −0.26, P < 0.01). Waves are known to have an impact on tidal flat slope10,20, but it is not significant in the current dataset (R = 0.13, P > 0.05; Fig. 2f). The combination of tidal range (58% of relative importance) and SSC (42% of relative importance) best explains the distribution of tidal flat slope (R2 = 0.30) (Fig. 2h, Eq. 4).

The high SSC leads to a large sediment settling flux during each tide. Hence, there must also be more energetic hydrodynamic forcing to re-erode the additional sediment to maintain the morphological equilibrium at such high SSC. In general, wider tidal flats have higher bed shear stress from tidal current and waves (τ90 in Eq. (12)) compared to narrower flats10 (Supplementary Fig. 6). Thus, it explains the correlation between high SSC and larger tidal flat width in space10,20. Additionally, tidal range also plays a secondary role in tidal flat width, as naturally larger tidal ranges increase the vertical extent of intertidal area, which contribute to wider tidal flats20. Lastly, tidal flat width exhibits a weak negative correlation with wave height, which may be caused by wave-induced mobilization of sediments10. The relation between the slope of tidal flats and tidal range is related to the spatial distribution of bed shear stress11,20. On tidal flats, the representative bed shear stress from tidal currents and wind waves is generally higher at the seaward boundary and decreases towards the landward boundary, leading to an erosion trend at the lower tidal flats but a deposition trend at the higher tidal flats. Greater tidal ranges can increase such differences and lead to steeper tidal flat profiles (see also Supplementary Fig. 7).

Temporal dynamics in tidal flat width and the underlying elastic responses

Tidal flat width shows a clear trend of overall erosion (Fig. 3a, c). Alarmingly, the nationwide median tidal flat width was reduced by 10% (596 m to 536 m) from 2002 to 2016 (Fig. 3c). 30.3% of tidal flats experienced a width reduction rate greater than 1% per year. Such a trend is particularly prominent on Hainan Island and in Changjiang River Estuary, where the interannual change rate in tidal flat width was −1.89% and −1.90% per year, respectively (Fig. 3a). In contrast, there were sites (28.3%) that exhibited an expansion trend of > 1% per year. Tidal flats on the southern Shandong (SD) Peninsula and the coast of Fujian Province increased their width at rates of 0.76% and 0.56% per year, respectively (Fig. 3a).

Fig. 3: Temporal variation of tidal flat width and suspended sediment concentration along the coast of China.
figure 3

a Inter-annual trend of tidal flat width (2002-2016) and (b) suspended sediment concentration (2003-2011). The values near the box represent the median values of the area. 51.3% of the tidal flat transects show a decreasing trend in width, while 54.9% of the transects experience a decrease in SSC. The insets represent the probability distribution of the data points, revealing that a majority of the instances for tidal flat width and SSC are experiencing a decrease. c Temporal changes in the median value of tidal flat width (linear regression, P = 0.12) and (d) SSC along the coast of China (linear regression, P = 0.03). Source data are provided as a Source Data file.

Full size image

Concurrently, there was a nationwide declining trend in sediment availability at -0.4% per year as observed in interannual SSC from 2002 to 2012 (Fig. 3d). To avoid spatial autocorrelation, we conducted a K-means clustering analysis, a measure for partitioning data points into non-overlapping groups. It reveals a strong correlation between the interannual SSC reduction and tidal flat width decline (R = 0.64, P < 0.01) (Supplementary Fig. 8). Additionally, 56% of all observation transects exhibit concurrent declines or increases of tidal flat width with sediment availability changes (Fig. 3a, b), whereas the rest of transects do not share this pattern as 86% of them do not bear a clear SSC change trend. These statistics suggest a close linkage between altered sediment availability and changes in tidal flat morphology.

Beyond the overall connection between SSC reduction and tidal flat width loss, we further discover diverse responses of sediment-starved (SSC < 50 mg/L, average SSC of China’s tidal flats) and sediment-rich (SSC > 75 mg/L, average SSC in five major Chinese estuaries, Supplementary Fig. 2; the Changjiang River Estuary, the Yellow River Estuary, the Pearl River Estuary, the Minjiang River Estuary, and the Liaohe River Estuary) sites to altered sediment availability. Notably, sediment-rich sites undergo 2.1 times faster width reduction than sediment-starved sites (0.27% vs. 0.13% width loss per year) when SSC decreases between -1.2% to -4.8% per year (Fig. 4a). Conversely, sediment-starved systems exhibit 2.1 times more rapid width expansion (0.62% vs. 0.30% width gain per year) when SSC increases between 1.2% to 4.8% per year.

Fig. 4: Elastic response of tidal flats width changes to altered sediment availability.
figure 4

Flat width changes with sediment concentration change in sediment-rich and sediment-starved areas based on observation (a) and explanatory modeling (b). The numbers at the base of each bar represent the number of observations. Model results are based on an exemplified sediment-rich site and a sediment-starved site (see “Methods”). Source data are provided as a Source Data file. c Diagram showing a sediment-rich tidal flat is similar to a stretched spring that tends to return to its morphological equilibrium (dashed line) and is resistant to being stretched further. d Diagram showing a sediment-starved tidal flat is similar to a compressed spring that tends to restore its morphological equilibrium (dashed line) and resists being compressed further. The convex shape of sediment-rich flats and the concave shape of the sediment-starved flats correspond to the ideal profiles from the dynamic equilibrium theory14,20.

Full size image

To elucidate these observed divergent responses, we conducted explanatory simulations of a sediment-rich and a sediment-starved tidal flat using the DET-ESTMORF model (Fig. 4b and Supplementary Fig. 9). The modeling results show similar qualitative conditional responses to sediment availability changes, but have different response rates in width change compared to observations. This difference is likely caused by the fact that these simulations are based on two exemplified sites and do not account for all the observed variations across the country. Following a previous study26, we further assess the relaxation time needed for new stable profiles. When external sediment availability decreases by 2%, the sediment-rich tidal flat reaches a new stable profile more quickly (12 years) than the sediment-starved tidal flat (24 years). However, when sediment availability increases by 2%, the sediment-starved tidal flat reaches a stable state faster (19 years) than the sediment-rich site (26 years) (Supplementary Fig. 10). These modeling results are also clearly in line with the observed pattern of tidal flat width changes.

These non-linear dynamics of tidal flat morphology depending on the current sediment input level can also be explained by an analogy to the elastic response of springs to compressional and stretching force. Thus, tidal flats in sediment-rich environments can be regarded as stretched springs due to the above-average sediment inputs (Fig. 4c). These systems have a tendency to retract to their long-term morphological equilibriums. When facing sediment input reduction, they tend to shrink quickly to restore equilibriums, but when sediment input rises, they are resistant to being stretched further from their equilibriums. By contrast, tidal flats in sediment-starved environments can be regarded as compressed springs (Fig. 4d). They have the tendency to bounce back to their equilibriums: expanding their width fast when sediment availability increases, while reducing their width slowly when sediment input is further reduced.

Finally, we conducted a synthetic modeling experiment to explore the elastic response of tidal flat morphology to sediment availability changes on a national scale (Fig. 5a). The results delineate that there are two ‘hot’ zones and a ‘neutral’ zone in tidal flat width dynamics. The two ‘hot’ zones are located at both extremes, i.e., increased sediment inputs in sediment-starved systems and decreased sediment inputs in sediment-rich systems, where the tidal flats are expected to experience rapid widening or narrowing ( > 0.3% per year). The ‘neutral’ zone, situated between the two hot zones, covers the space of sediment reduction in sediment-starved systems, and also sediment increase in the sediment-rich systems, where the change in tidal flat width is expected to be slow. The mechanism of tidal flat elastic responses in the ‘hot’ and ‘neutral’ zones is included in Supplementary Text. This synthetic modeling aligns with the observed data points along China’s coastline. Sites with higher rates of widening or narrowing mostly fall in the ‘hot’ zones, while sites with lower rates generally fall in the ‘neutral’ zone. It is noteworthy that there are no data points with a rapid SSC increase ( > 1% per year) in sediment-rich systems (Fig. 5a), reflecting the fact that China’s tidal flat ecosystem is currently experiencing a nationwide reduction in sediment availability.

Fig. 5: Modeled and observed tidal flat width response to sediment availability changes.
figure 5

a The contour map represents the tidal flat width change corresponding to different variation rates in SSC from DET-ESTMORF modeling, whereas bubbles indicate the representative data points for tidal flat width changes and the number indicates the location along China’s coast (see “methods” for data selection). Red bubbles indicate tidal flats expansions, and green bubbles indicate tidal flats retreats. The size of the bubbles indicates the variation rate. The upper left and lower right corners of the diagram are two ‘hot’ zones, i.e., sediment-starved sites with an increase of SSC and sediment-rich sites with SSC drop. The area between the two solid black lines indicates a ‘neutral’ zone where the interannual variation of tidal flat widths is within -0.05% and 0.3%. There are no data points in the upper right corner, i.e., rapid SSC increase ( > 1% per year) in sediment-rich systems. b The rate of tidal flat width changes in relation to the initial SSC level when SSC increases by 3 ~ 9 %/year in all four countries (p < 0.01), indicating slower expansion of sediment-rich areas if SSC increases. c The rate of tidal flat width changes in relation to the initial SSC level when sediment decreases by 3 ~ 9 %/year in the four countries (p < 0.01), indicating faster retreat of sediment-rich areas if SSC decreases. The brown line is the linear regression line with the p-value. Data from these countries are randomly selected based on their initial SSC level and SSC change rate (see “Methods”). ‘N’ means the number of data points selected from each country.

Full size image

In our analysis, we focus on understanding how tidal flats respond to external disturbances, since losses due to direct human interventions (such as conversion to other land use) are generally less predictable. We identified 9 transects that are converted at the beginning of our study period (e.g., 2002, see “Method”), and they are observed to generally accumulate at fast rates regardless of the SSC level (Supplementary Fig. 11). Due to the different response pattern, we excluded transects affected by direct conversion, accounting for 17% of the total transects (Supplementary Fig. 12). In addition to direct conversion, we conducted analysis excluding transects in 9 estuaries or bays, as well as those adjacent to converted transects along the open coast, where the far-field effect of conversion could be important (see “methods”). We show that the tidal flat width response to SSC change still hold when the data of these areas are excluded (Supplementary Fig. 13).

The observed elastic response of tidal flat morphology in China also has global relevance, as similar dynamics are observed from three other countries with substantial tidal flat areas, i.e., Indonesia, Australia, and the United States. These three countries are ranked 1st, 3rd, and 4th place globally in tidal flat areas (with China being the 2nd place). Data from these four countries collectively show: 1) tidal flat width grows if SSC increases and shrinks if SSC declines; 2) tidal flats with higher SSC have lower expansion rates as SSC increases, but higher retreat rates if SSC reduces (Fig. 5b, c). Furthermore, based on the data from China’s coast, we found that in sediment-rich areas tidal flat width losses would increase rapidly once the SSC is reduced by 24% from 2002 (Supplementary Fig. 14). This finding can serve as an important rule of thumb for coastal management.

Discussion

Tidal flat systems are highly valued for their ecosystem services, but their morphological evolution on a grand scale remains unclear. In this study, we use China’s tidal flat ecosystem as an example and elucidate that tidal range predominantly influences the slope of tidal flats, compared to incident waves or SSC. In our system, a higher tidal range typically leads to a steeper slope. This is complementary to the observations from microtidal systems like the Venice Lagoon (Italy) or the Virginia Coast Reserve (USA), where the tidal flats are exceptionally level due to the low tidal ranges10,27. Different from the slope, the spatial distribution of tidal flat width is mainly controlled by SSC. The widespread reduction in nearshore SSC is likely the main cause of the national-scale tidal flat narrowing. These insights may provide global implications. For instance, the observed large intertidal area retreat in other countries in the temperate Northern Pacific region (e.g., the United States)28,29 is likely also caused by the reduction in sediment availability.

Interestingly, our observation, combined with synthetic modeling, unveils a non-linear response of sediment-rich and sediment-starved tidal flats to altered sediment availability across four major countries with extensive tidal flat areas and diverse sediment supply levels (i.e., China, Indonesia, Australia, and the United States). The conditional responses of tidal flat morphology can be explained by the balance between sediment availability and hydrodynamic forcing. Sediment-rich tidal flats typically have flatter and shallower profiles, which are subject to greater wave and current induced bed shear stresses10. On these flats, surplus sediment is less likely to settle, hindering tidal flat expansion, but a sediment availability deficit can easily lead to rapid erosion due to the primitive high level of bed shear stress. Conversely, sediment-starved flats have relatively lower bed shear stresses due to the steeper and deeper profiles10, making them able to accommodate surplus sediment and resistant to further erosion during sediment deficit. These conditional responses lead to the observed spring-like behaviors of tidal flat morphodynamics that tend to move towards the equilibrium state rather than being driven further away from it (Fig. 4). This mechanism can be well explained from a modeling perspective, provided in Supplementary Text.

Sediment-starved tidal flats (often found on the open coast) are in a compressed state compared to the equilibrium, which is more sensitive to sediment availability increase, but less sensitive to further reduction in sediment availability. Thus, these systems may be more resilient than previously thought, whereas management efforts that increase sediment availability, e.g., controlled river diversion30 and nourishment operations31, would effectively extend tidal flat width in these systems. In contrast, sediment-rich tidal flats are more susceptible to sediment availability declines. Therefore, SSC in these systems (typically in estuaries) should be cautiously monitored and maintained (if possible) to prevent unexpectedly fast intertidal area losses. Moreover, considering the widespread SSC reduction due to river damming, it is essential to preserve at least 76% of the long-term stable SSC to prevent severe tidal flat area loss (Supplementary Fig. 14). On the other hand, actions to increase sediment inputs are not that efficient in sediment-rich systems, as tidal flats in these systems resist being stretched further away from their equilibrium and would not accommodate more sediment. The redundant sediment is better allocated to sediment-starved systems to avoid diminished light penetration to adjacent coastal ecosystems, such as macroalgae32, seagrasses33, and coral reefs34,35, which is vital to their survival.

Due to the decadal time scale (2002–2016), sea-level rise does not have a clear impact on tidal flat morphology based on data from 14 Chinese long-term water level observation sites (p = 0.89, see “Methods”). However, over longer time scales (e.g., a few decades), its impact can be pronounced, especially with concurrent SSC reduction24,36,37,38. Sediment-starved tidal flats face greater threats under accelerated sea-level rise scenarios39, whereas sediment-rich tidal flats must also maintain their sediment availability to mitigate these risks36. The applied global data sources (e.g., FAST, ERA-5, tidal flat map) in the current studies may contain errors in specific sites. For example, the tidal flat map data may misidentify highly turbid waters1, which may lead to an overestimation of the tidal flat widths. Furthermore, the 30 m spatial resolution of the tidal flat map affects the accuracy of relatively narrow tidal flat widths, however, the overall impact is limited given the median flat width being 536 m in China. Despite these potential errors, our analysis aligns with dynamic equilibrium modeling and is consistent with numerous previous studies10,11,20.

In light of global sediment flux alterations16,17 and widespread tidal flat distribution1,8, our findings underline the imperative of developing sustainable sediment allocation strategies, in which sediment availability history and its changing trends should be considered synthetically. Thus, in the global hydrologic north, where sediment availability is decreasing, efforts should focus on preserving sediment levels in sediment-rich tidal flats to prevent rapid loss of intertidal zones. Conversely, in the hydrological south, where sediment availability is increasing, excess sediment should be redirected from sediment-rich tidal flats to sediment-starved ones, promoting tidal flat expansion in areas more responsive to increased sediment flux. These strategies would pave the way towards practical guidelines to fine-tune the sediment fluxes to our coasts, enabling effective tidal flat conservation and restoration.

Methods

Tidal flat width and slope data

The width of the tidal flat is defined as the horizontal distance from mean high water to low water (Supplementary Fig. 15). The data of tidal flat area along China’s coastline were obtained based on the global tidal flat maps1. These maps used over 7.0 × 105 satellite images to record tidal flat extent, and the data explicitly exclude vegetation-dominated intertidal ecosystems (such as mangroves and vegetated marshes). The resultant global map of tidal flats exhibits a horizontal resolution of 30 m, delineating the spatiotemporal variation of tidal flats on a global scale. To obtain the width of China’s tidal flat, we measured the distance of cross-shore transects along the shoreline in Open Street Map (OSM)40 every 1 km. The overall accuracy of the tidal flat area map is 82.2%41. Misclassification of tidal flats (pixels) can occur in areas of highly turbid water, polluted areas, aquaculture ponds, and areas subject to strong wave action. Despite the potential errors, to our knowledge, this map is the only freely accessible tidal flat map that covers large spatial (national) and temporal (decadal) scales. The accuracy of this map is sufficient ( > 80%)41 to reveal the overall dynamics of tidal flat areas and to support the main findings of the current study.

The slope of the tidal flats was ascertained from the Foreshore Assessment using Space Technology (FAST) database21, which provides a high-resolution intertidal bathymetry/elevation dataset. The bed level data has a 20 m horizontal resolution and typically a 30-50 cm vertical accuracy. To obtain the tidal flat slope, we extracted the bed level data of cross-shore transects along the shoreline. The transect endpoints are based on the range of the tidal flats. The slope of the tidal flat was then calculated by the cross-shore horizontal length L of the transect and vertical height difference ∆H between the endpoints (i.e., slope = ∆H/L, Supplementary Fig. 16). Transects with discontinuous elevation were excluded from analysis (917 of 2538 transects). There are mainly two reasons that cause the discontinuous elevation: 1) gaps within the continuous elevation data (51.4%); 2) none or too few elevation data in the FAST database (48.6%).

To identify the key factors that sculpt the spatial morphology of tidal flats, we compared tidal flat width data (2014–2016) and slope data (1997–2017) with SSC (2016), tidal range (2010), and wave height (1979–2021). To reveal the temporal changes in tidal flat width, we traced the tidal flat width change data from 2002–2016 and SSC change data from 2003-2011, constrained by data availability. Changes in tidal flat width were obtained from tidal flats maps every three years (2002–2004, 2005–2007, 2008–2010, 2011–2013, and 2014–2016)1. The data on the annual variation in tidal flat width shows no significant spatial autocorrelation (Moran’s I < 0.01, p = 0.44).

Hydrodynamic data

We utilized the TPXO9-atlas dataset42 to obtain the mean tidal range for each cross-shore tidal flat transect along China’s coast. The spatial resolution of the dataset is 1/30°. The accuracy of the TPXO9 on China’s coast can be shown by an average Root Sum Squares (RSS) value of 19.92 cm43. Employing the tide model driver (TMD) tool, we processed the TPXO dataset to calculate the water levels at each transect for the year 2010. This was done by incorporating the amplitudes and phases of eight major tidal constituents: K1, K2, M2, N2, O1, P1, Q1, and S2. The average tidal range data was then interpolated from the mean high water level and the mean low water level. Tidal asymmetry is calculated from the phase difference between the M2 and M4 tidal constituents \((2{\theta }_{M2}-{\theta }_{M4})\). A phase difference within the range of 0–180° results in a flood-dominant tide. Conversely, a phase difference between 180° and 360° indicates an ebb-dominant tide44,45. Under the same phase difference, the amplitude ratio of M4 and M2 tidal constituents is used to indicate the magnitude of the tidal asymmetry. The tidal flats along China’s coastline exhibit a tidal range that varies from a minimum of 0.59 m near the Yellow River mouth to a maximum of 7.33 m on the southern Shandong Peninsula, with an average of 2.78 m.

We collected SLR data from 14 coastal sites in China46 and calculated the correlation between the SLR data and changes in local tidal flat width. However, the results indicate no significant correlation (p = 0.89). The annual mean sea level data of the 14 coastal sites are collected from the Permanent Service for Mean Seal Level47,48 (see Supplementary Table 1). To test the possible effect tidal range has on tidal flat width, we compared the tidal range change from the 12 UHSLC tide gauge sites49 with the tidal width changes. The results indicated no significant correlation between the tidal range changes and the nearby tidal flat width changes (p = 0.32), which is likely due to the relatively small changes in tidal range (averaged 1.1 cm per year).

The offshore wave conditions were derived from the ERA5 reanalysis database50, which encompasses comprehensive wave height and period data from 2002–2016. The database has been proven suitable for China’s coastal areas51, with data bias estimated at less than 0.1 m in deep waters. These data were firstly reprojected to align with the offshore positions of the tidal flat transects. Subsequently, the offshore wave height values were transferred to the wave height situated 3 km off the coast, equating to 78% of the offshore wave height52. In the absence of detailed nearshore bathymetry data, we assumed linear profiles from 3 km off the coast to the mean low water level53, while the intertidal bathymetry (low water to high water level) data are available from FAST21. The offshore depth is obtained from the ERA5 database. This approach allowed for the estimation of wave heights in the immediate vicinity of the tidal flats at mean low water level, accounting for wave shallowing and bottom friction54:

$${D}_{w}=\frac{2}{3\pi }\rho {f}_{e}{{U}_{\delta }}^{3}$$
(1)

with

$${U}_{\delta }=\frac{\pi {H}_{s}}{{T} \, {sinh} \, \, {kh}}$$
(2)

where Dw is the wave height dissipation, ρ is the water density, fe is the wave friction factor, which is assigned a value of 0.05555, Uδ is the wave orbital velocity, Hs is the significant wave height, T is the wave period, k is the wave number, and h is the water depth. Wave refraction and diffraction processes are not accounted due to the lack of nearshore bathymetry data.

To assess if a change in waves might affect tidal flat width, we calculated the change rate of significant wave height at 105 sites along the coast of China (with a resolution of 0.5 degrees) from ERA5 data. It shows that 95% of tidal flat areas have experienced a minimal variation in significant wave height below 0.1 cm/year from 2002 to 2016. Thus, there is no significant correlation between the wave height change and tidal flat width dynamics (p = 0.42).

To validate the applicability of the ERA5 data, we collected wave data from tidal flats at various locations along China’s coast, including Chongming Island56,57 and Leizhou Bay58 for both wave height and period, Taiwan Longdong59 and the Yellow River mouth60 for wave height, and Jiangsu61 for wave period. The agreement between our estimated wave heights and periods with in-situ observations (Supplementary Tables 2 and 3) confirms good applicability of ERA5 wave data. In the current study, the impact of alongshore currents is not explicitly considered, since there is a lack of nationwide dataset on alongshore currents. Yet, as alongshore currents transport sediment over large scales62,63, their impact is implicitly reflected by the SSC data we analyzed.

Suspended sediment concentration (SSC) data

SSC is considered representative of the local sediment availability36,37,38,64, and a key external driver on tidal flat morphodynamics. The SSC data near the monitoring tidal flat transects were obtained from the GlobColour project (https://www.globcolour.info). The data product of total suspended matter concentration (TSM, mg/L) produced by the project is used to represent SSC (mg/L) because in turbid coastal waters most suspended matters are due to suspended sediments. It is estimated from satellite-derived particle backscattering coefficient in the blue wavelength by assuming a backscattering/scattering ratio of 0.015, and such a product has been previously used and found reliable for large scale analyzes64,65. The spatial patterns of such derived TSM are consistent with those from independent studies using other algorithms66. Although satellite-derived TSM (or SSC) represents the particle concentration in surface waters (down to an optical depth), for both well-mixed water column or highly stratified waters it can be used as a proxy for suspended particles in the water column66. This is because, for the former case it is obvious that TSM reflect the overall sediment concentration and, for the latter case, the signals from the bottom, particle-rich layer will only be attenuated slightly by the surface, particle-poor layer. The applicability of this approach has been validated in terms of both spatial distribution67 and temporal variability68, and the method has been widely applied in previous studies (e.g., sediment availability data in ref. 64 and turbidity for estuaries in ref. 69).

The SSC data product has a horizontal resolution of 1/24 degree (generally 2.2 km off the coast) and a temporal resolution of one month. To investigate the drivers behind the spatial distribution of tidal flat morphology, we analyzed the monthly average SSC data obtained through the inversion of OLCI satellite data for the year 2016, synchronized with the tidal flat width (2002 - 2016) and slope data (1997 - 2017). The highest annual average SSC in the Pearl River estuary is 100 mg/L, and within the Changjiang River Delta, particularly in the Hangzhou Bay area, the maximum SSC reaches up to 200 mg/L, which is in line with previous field observations70,71. To assign a local sediment concentration level to each tidal flat data point, the TSM data point closest to the target point was used.

We define sediment-starved systems in China as the areas with an annual mean SSC below the national mean SSC level (50 mg/L). This value is obtained by calculating the average SSC of pixels containing tidal flats (4 km resolution). Correspondingly, areas with an annual average SSC higher than 75 mg/L are defined as sediment-rich systems, which are the annual mean SSC in 5 major estuaries (the Changjiang River estuary, the Yellow River estuary, the Pearl River estuary, the Minjiang River estuary, and the Liaohe River estuary). These thresholds can be considerably higher than other countries17,69, due to the overall higher SSC on China’s coast. Thus, sediment-starved and sediment-rich systems should be defined differently in other regions of the world.

Human impact through conversion

Conversion of tidal flats to other land use directly and indirectly influence the evolution of tidal flats. To investigate the morphological evolution of large-scale tidal flat systems with environmental factors (i.e., hydrodynamic forcing and sediment availability), we excluded transects affected by direct conversion. First, we identified tidal flat transects with direct conversion by analyzing satellite imagery (see Supplementary Fig. 12). The tidal flats are mainly converted to agriculture or aquaculture land use, and also (about 1.2%) to artificial sandy beaches in areas like Xiamen city72. 750 transects are affected by direct conversion, which are excluded from our analysis.

To compare the effects of direct conversion and external drivers, we identified nine transects that were converted at the beginning of our study period and have SSC change data, so that we can track their development. The results show that eight out of nine transects exhibited a seaward progradation trend after the conversion (Supplementary Fig. 11). The rate of tidal flat width change exceeds the rate of tidal flat evolution driven by natural drivers, highlighting the important difference imposed by direct conversion and the necessity to exclude those transects in our study.

As conversion can also affect nearby hydrodynamic conditions, thereby influencing adjacent remaining tidal flat evolution by far-field effects. Therefore, we further excluded transects adjacent to those directly impacted by conversion. We identified nine estuaries or bays where the conversion level is high (over 40% of the transects were directly converted), including Qinjiang Estuary, Yingzai River Estuary, Zhenhai Bay, Luojiang Estuary, Wenzhou Bay, Wulei Island Bay, Liaohe Estuary, the southern coast of Hangzhou Bay, and Qixing Reef. The remaining transects were identified as being influenced by far-field effects of conversion. Additionally, those adjacent to converted transects along the open coast are also deemed influenced by far-field effects. After excluding all these 156 transects, our primary conclusions regarding spatial and temporal patterns remain unchanged (See Supplementary Fig. 13).

Statistical analysis

Multiple Linear Models (MLM) were used to identify the key drivers that shape the spatial distribution of tidal flat morphology. The statistical analysis were performed using the ‘R’ software. MLMs were created for tidal flat width and tidal flat slope, respectively. Given the low resolution of offshore wave height data (0.5°), we consolidated data points that corresponded to the same offshore wave height and computed the mean values for tidal flat width, significant wave height (Hs), suspended sediment concentration (SSC), and tidal range (TR) within each offshore wave height pixel.

First, we used stepwise linear regression to select the best combination of independent variables that explain the spatial distribution of the tidal flat morphology. This process was guided by the selection of the model yielding the lowest Akaike Information Criterion (AIC) scores, a widely accepted criterion for model selection73,74,75. For tidal flat width and slope, SSC and tidal range were included in the most streamlined final model. Subsequently, we established a multiple linear regression model for the selected variables (Eq. 3 and Eq. 4) and calculated their relative importance. The relative importance is calculated by assessing their contribution to the model’s determination coefficient (R2).

$${{width}}^{\frac{1}{3}} \sim {{SSC}}^{\frac{1}{3}}+{{TR}}^{\frac{1}{3}}$$
(3)
$$\log \left({slope}\right) \sim {{SSC}}^{\frac{1}{3}}+{{TR}}^{\frac{1}{3}}$$
(4)

To obtain a normal distribution of the data for proper statistical analysis, we applied transformations to the tidal flat width, tidal range, and SSC data by taking their cube roots, and log-transformed the tidal flat slope data prior to their inclusion in the statistical models. Considering the potential impact of spatial autocorrelation, we calculated Moran’s I for the residuals in our regression models for width (0.35, p < 0.01) and slope (0.11, p = 0.17). We then used a spatial error model to eliminate the effects of spatial autocorrelation for the width model, achieving a better fit than the multiple linear regression (R2 = 0.61). P-values and R-squared values were calculated using the ‘summary’ package in R software, and the relative importance of each parameter in the linear model was calculated using the ‘relaimpo’ package75.

DET-ESTMORF model

To explore the relationship between sediment availability variations and the corresponding temporal shifts in tidal flat width, we utilized a one-dimensional tidal flats morphological model, DET-ESTMORF (Fig. 4 and Fig. 5)11. This model is a hybrid approach that integrates morphological equilibrium20,22 with hydrodynamic and sediment transport processes, in contrast to fully process-based models12. This model has been validated previously11, and applied in intertidal biogeomorphology studies76.

Tidal flats in morphological equilibrium are characterized by a consistent distribution of bottom bed shear stress and SSC throughout their profile. When the actual local bed shear stress surpasses the equilibrium shear stress, there is a tendency of erosion, requiring a higher SSC to counteract the tendency and maintain equilibrium, i.e., a higher local equilibrium concentration than average. Hence, the magnitude of bottom bed shear stress determines the tendency of tidal flats’ geomorphic change, while the difference between actual SSC and local equilibrium SSC indicates the occurrence of geomorphic change, which is similar to the process-based models, e.g., Delft3D.

In DET-ESTMORF, we considered both the effects of wave and tidal shear stresses on the bottom bed. To obtain the tidal-induced bottom bed shear stress, we adapted the method based on mass conservation from a previous study55. Concerning the wave propagation modeling, we applied the wave energy flux conservation equation considering bottom friction54:

$$\frac{d}{{dx}}\left(E{c}_{g}\right)=\frac{d}{{dx}}\left(\frac{1}{8}\rho g{{H}_{s}}^{2}{c}_{g}\right)=-{D}_{w}$$
(5)

where E is the wave energy [kg m−2 s−2)], cg is the group velocity (s2 m−1), and g is the acceleration of gravity (m s−2). The incident Hs is obtained by a random draw from the Rayleigh distribution around a given mean value54.

The attenuation of wave energy by friction is given by Eq. 1 and Eq. 2. We also considered the effect of wave breaking in the wave propagation process, where the maximum wave height Hs is associated with the local water depth and the breaking height coefficient:

$${H}_{s}=\min \left({H}_{s},{\gamma }_{b}h\right)$$
(6)

where γb is a constant value (0.5). The bed stress induced by waves τw is quantified as77:

$${\tau }_{w}=\frac{1}{4}{\rho {fw}{U}_{\delta }}^{2}$$
(7)

where fw is a friction factor estimated as:

$${f}_{w}=1.39{\left(\frac{{A}_{\delta }}{{k}_{b}}\right)}^{-0.52}$$
(8)

Where the wave orbital semi-excursion at the bed \({A}_{\delta }={U}_{\delta }T\), bed roughness \({k}_{b}=2\pi {d}_{50}/12\), and d50 is the median grain size of sediment particles. The mean bed shear stress under combined waves and tidal currents during a wave cycle is calculated as78

$${\tau }_{m}={\tau }_{{cur}}\left[1+1.2{\left(\frac{{\tau }_{{wave}}}{{\tau }_{{cur}}+{\tau }_{{wave}}}\right)}^{3.2}\right]$$
(9)

Thus, the maximum bed shear stress during a wave cycle is calculated as

$${\tau }_{\max }={\left[{\left({\tau }_{m}+{\tau }_{{wave}}\left|\cos \theta \right|\right)}^{2}+{({\tau }_{{wave}}\left|\sin \theta \right|)}^{2}\right]}^{0.5}$$
(10)

where θ is the angle between current direction and direction of wave propagation. In our study, θ = 0 as both currents and wave propagation are in the cross-shore direction. τmax is quantified every 30 minutes during a tidal cycle. Subsequently, τ90 is derived as the 90th percentile of τmax over a tidal cycle, which is considered as the hydrodynamic forcing that drives morphological changes.

In morphological equilibrium, the balance between erosion and accretion in the vertical direction is achieved, which can be expressed as follows79,80:

$${m}_{e}\left(\frac{{\tau }_{E}}{{\tau }_{{cr}}}-1\right)={c}_{E}{w}_{s}$$
(11)

Here, me represents the erosion coefficient [kg m−2 s−1)], τcr denotes the critical shear stress (Pa), and ws represents the sedimentation rate (m/s), cE is the long-term average SSC value (mg/L) at the seaward boundary. Via this equation, τE is the uniform bed shear stress (Pa) that can be derived from the known variables. Subsequently, the local equilibrium concentration ce based on the ratio of τ90 to τE can be derived as:

$${c}_{e}={c}_{E}{\left(\frac{{\tau }_{90}}{{\tau }_{E}}\right)}^{n}$$
(12)

The local equilibrium concentration ce can be converted to volume concentration by dividing the sediment density (ρs = 2650 kg/m3)81. The actual spatial and temporal distributions of local SSC are obtained by solving the diffusion equation:

$$\frac{\partial ({hc})}{\partial t}={w}_{s}\left({c}_{e}-c\right)+\frac{\partial }{\partial x}\left({Dh}\frac{\partial c}{\partial x}\right)$$
(13)

where c is the volume sediment concentration (m3/m3), and D represents the diffusion coefficient (m2 s−1). The change in tidal flats morphology following each tidal cycle is calculated based on the difference between c and ce:

$$\frac{\partial z}{\partial t}=\frac{1}{1-p}{w}_{s}\left(c-{c}_{e}\right)$$
(14)

where z is the bed elevation (m), and p is the bed porosity. The key parameters are detailed in supplementary Table 4, and the values of these parameters are based on previous modeling work (Supplementary Table 5)10,12,82.

To reveal the nonlinear response of tidal flat width to SSC change, we selected two sites for explanatory modeling, namely the Western Scheldt site in the Netherlands (sediment-starved)76 and Jiangsu in China (sediment-rich)63. We did not select a Chinese site for the case study of sediment-starved area, as it was not possible to find detailed field data for model validation in China, while the Western Scheldt has a unique detailed dataset available. At the Western Scheldt site (SSC = 58 mg/L), the DET-ESTMORF model has been applied and validated in previous studies11,76. At the Jiangsu site, we obtained the equilibrium profile by performing a 100-year modeling from the initial linear profile based on the tide level, while wave height and offshore SSC data (134 mg/L) are based on the collected datasets in the current study. We validated the model’s accuracy by comparing the obtained equilibrium profile with the observed profile reported in a previous study (Supplementary Fig. 17)63. Building on this validation, we proceeded to calculate equilibrium profiles for both sites under varying interannual SSC change rates, including −1.2%, −0.8%, 0.8%, and 1.2%. This exercise was designed to ascertain the interannual rates of tidal flat width change over a 13-year period, in line with the observational timeframe spanning from 2003 to 2015.

Analysis of temporal changes in tidal flats width

To analyze the temporal changes in tidal flat width, we integrated measured data (the bubble chart of Fig. 5a) with experimental results (colormap of Fig. 5a). We extracted tidal width change rate of tidal flats with SSC above 50 mg/L and those below 50 mg/L across different SSC change rate groups (−3% to −1.5%/year, −1.5% to 0%/year, 0% to 1.5%/year, 1.5% to 3%/year). Within each group, we chose median tidal width change (mean of the 48%–52% percentiles) as the representative value (Fig. 5a). To minimize the influence of anthropogenic factors on the geomorphology of tidal flats, data points exhibiting interannual tidal width trends exceeding 20% per year were omitted.

We then employed the DET-ESTMORF model to simulate the variation of tidal flat width under different SSC (20 ~ 120 mg/L) and varying interannual rates of change of sediment concentration (−3% ~ 3%/year). For a consistent basis of comparison, we initialized the model with a wider profile for tidal flats with higher sediment availability, maintaining uniform parameters across all scenarios except for the initial SSC level and the SSC variation rates. For a given SSC level, we first obtained the equilibrium profile from the initial linear profile after 100 years of simulation, after which the change in the tidal flat morphology is negligible76. Based on the obtained equilibrium profile, we simulated the tidal flat morphological variations with changing SSC over 13 years (2003 to 2015, as a reflection of the observation), and finally compared the simulated flat width changes resulting from various SSC changing rates.

To test the observed elastic response of tidal flat morphology in China, we further analyzed the data from three other countries with substantial tidal flat areas, i.e., Indonesia, Australia, and the United States. There are a large number of observation transects in these countries and many of them are affected by direct human intervention. To select the representative transects, we categorized tidal flat transects from these three countries into four groups for each country, based on initial SSC levels (0-20 mg/L and > 20 mg/L) and SSC change rates (-9 to -3 %/year and 3 to 9 %/year). Within each group, we randomly selected 50 tidal flat transects, and subsequently excluded the transects with direct human impact based on the satellite base map from Google Earth Engine. This resulted in 34 transects in Indonesia, 49 in Australia, and 33 in the USA. For China, we grouped SSC into sediment-starved and sediment-rich areas as mentioned in the previous section.