Building for Nature: Preserving Threatened Bird Habitat in Port Design

: The fast economic development of the People’s Republic of China has created an increasing demand for usable land, resulting in large-scale land reclamations along the coastal zone. One of these regions is Tongzhou Bay (Jiangsu coast), a region characterized by large intertidal mudﬂats and deep tidal channels with potential for the development of agri-aquaculture and the construction of a deep-sea port. However, these intertidal mudﬂats also provide vital ecosystem services and support many wildlife species, including several endangered migratory shorebirds within the East Asian–Australasian Flyway. With increasing realization of the importance of maintaining such ecological values, a more integrated coastal development strategy is needed. This study aims to develop a sustainable integrated design for the Tongzhou Bay port, following a “Building with Nature” approach. We use a morphodynamic model to compute habitat suitability for two shorebird species (Great Knot Calidris tenuirostris and Bar-tailed Godwit Limosa lapponica ). Several port conﬁgurations were developed on the basis of three design criteria: (1) create area for future port development, whilst (2) preserving existing high-value ecotopes for shorebirds and (3) enhance the natural accretion rate of such ecotopes. Simulation results showed a clear di ﬀ erence in siltation patterns, preservation and enhancement of preferred ecotopes. This work therefore demonstrates the potential and importance of morphological and habitat suitability modelling when designing large-scale reclamations and port constructions, especially in dynamic areas such as Tongzhou Bay.


Introduction
Coastal habitats around the world have great ecological importance as well as economic potential. As a consequence of the increased human pressure, especially over the last few centuries, there have been huge losses to coastal wetlands [1]. China, with a net growth of 6.7% in GDP in 2016, is the most rapidly developing country in the world. Here, 43.3% of the total population lives in the coastal developed and simulated for the effect on morphology and habitat development. Finally, conclusions are drawn from these alternative designs for the Tongzhou Bay testcase, demonstrating the effectiveness of this systematic approach.
Methodologies to quantify such alternative designs are still largely under development. We expand on an approach by Muller et al. [28], who quantified the impact of designs on wetland habitat by constructing an ecotope map. In this paper, we develop an interactive ecotope map in which the morphodynamic response of the environment to an intervention (siltation, erosion) is also accounted for.

Methods for Sustainable Coastal Development
Water infrastructure works, such as port development, river deepening and land reclamations, are essential for economic development. However, such interventions also impact ecology and geomorphology (e.g., van Maren et al. [38]; Wang et al. [39]; Winterwerp et al. [40]). The traditional method to deal with such adverse impacts was to evaluate positive and negative effects on long and short term and intensify dredging legislation. Traditional port design primarily focuses on primary objectives, e.g., the accessibility for vessels and connection with the hinterland [41]. Impact on nature and ecosystems is often limited to an inventory of the risks and designating compensation. Moreover, this often conflicts with environmental legislation and results in delays [42]. This fairly passive approach toward infrastructural development is changing, however. In the last decade, significant steps have been taken in more integrated dredging and reclamation activities (e.g., port development, with contemporary cases in the Netherlands [43][44][45]) and elsewhere [46]. Here, economic development is combined with not only compensating, but also creating larger opportunities for wetland ecosystems. It is increasingly realised that water infrastructure development can also positively influence the environment (through, e.g., Building with Nature solutions; [19]).
Traditionally, designs were evaluated for their performance, focusing on operational objectives in terms of hydrology. However, a lack of knowledge regrading sediment transport and the relation of these processes to local and regional geomorphology resulted in negative effects on engineering performance (e.g., higher siltation rates in channels and harbour basins) and ecosystems (e.g., loss of habitat). Water infrastructure projects are part of a system (e.g., ecosystem) and the projects will both affect and be affected by the processes operating within that system. Integrating the concept of sustainability into our infrastructure projects will help us identify opportunities to cooperate and collaborate with natural processes, rather than seek to control them. The BwN program [19,47] implements this approach by allowing natural processes to contribute to objectives in the design process. This requires profound knowledge of the system and the parameters and key indicators on the state of a system that we want to change, i.e., a natural system, coral, mudflats or socio-economic state.
The original design for Tongzhou Bay port was the total reclamation of intertidal flat, resulting in a dramatic loss of high-value ecotopes and inevitable loss of important habitat. Moreover, the current wetland ecosystem is already under pressure due to declining wetland habitat from past reclamations activities [48]. In order to come up with more integrated alternative designs, we define the following criteria: (1) the design should provide a well exploitable area for port development (e.g., close to natural channels and over shallow easy-reclaimable shoals), and (2) simultaneously maintain crucial wetland habitat, providing vital ecosystem services and (3) enhance natural accumulation of these habitats to fully support the wetland ecosystem. For the port design, we were guided by various designs that were made (taking into account potential port layouts, economic viability assessments, port accessibility such as approach channels, turning basins and connectivity with the hinterland [26,49]).

The Jiangsu Coast and Tongzhou Bay
The Jiangsu coast is influenced by two distinct large-scale tidal circulations. These are the anti-clockwise rotary tidal wave in the northern part of the Southern Yellow Sea (SYS) and the progressive tidal wave in the SYS coming from the Pacific, converging at Jiangang [50][51][52][53] creating a "rope-skipping pattern" [54,55]. These patterns are mostly governed by the shape and dimensions of the Bohai, Yellow and East China Seas [56] and are negligibly influenced by the nearby Yangtze river [57]. Near Tongzhou Bay, the semi-diurnal tide has a tidal range of 3.9 m (measured at Lüsi station; [54,56]).
Confined by parallel intertidal ridges, tidal currents in the channels are nearly bilinear with average flow velocities of 1.0~1.2 m/s [55]. The mean velocities during flood and ebb current do not differ significantly [58]. Due to the complex topography between ridges and channels, waves have a limited fetch and only limitedly propagate into the tidal channels, and local wave heights are therefore small [59][60][61]. The wind and wave field varies seasonally, with larger wave heights during winter than in summer [62]. Although waves are important for re-entrainment and higher suspended sediment concentrations (SSC) in winter than in summer [55,63], they are too small to substantially influence long-term geomorphodynamics (as reflected in the tide-dominated large-scale morphology).
The Jiangsu coast can be divided into two distinct geomorphological zones: the northern Abandoned Yellow River Delta (AYD) and the fan-shaped Radial Sand Ridges (RSR) in the south ( Figure 1b). The AYD is a typical river-wave dominated delta that was formed by the large fine-grained sediment load carried by the Yellow river, until the river shifted its course towards the Bohai sea in 1855 AD [64]. Upon abandonment and continuing until the present, the former delta started to erode. The RSR consist of a large fan-shaped radiate system of mudflats and channels that vary from 2 m to 20 m deep [51]. Sediments were deposited here by the Paleo-Yangtze River Estuary [65][66][67]. As the Yangtze shifted south, the former river courses gradually became tidal channels and tidal currents reworked sandy deposits into long extensive shoals. The combination of abundant sediment supply, large-scale tidal currents create and maintain the typical sand ridge field at the Jiangsu coast [50,53]. The RSR system is in an overall trend of accretion, partly resulting from sediment supply from the AYD [68,69], maintaining a typical convex-up cross-shore profile [57,70,71].
Tongzhou Bay (Figure 1) is located in the southern part of the RSR. It is characterized by a deep tidal Dawanghong (DWH) and Wangcanghong (WCH) channel, which continue into the Xiaomiaohong (XMH) and Sanshahong (SSH) channel, respectively. The channels are intersected by narrow tidal ridges, extending into larger shoals such as the Yaosha shoal and Lengjiasha shoal further landward.
The grain size of sediments along the Jiangsu coast increases in the offshore direction with clayey silt or clay predominantly deposited near the Mean High Water (MHW) line, (sandy) silt along the coastline (3.9~62.5 µm) and fine sand more offshore (125~250 µm) [72][73][74][75]. The tidal flats in the Southern RSR and Tongzhou Bay are wider with finer sediments than in the Northern RSR, as most of the finer material originates for the Yangtze sub-aqueous delta instead of the AYD [74,76]. The southern deflection of the present-day Yangtze River plume results in a limited influence of fresh river-borne sediments at the Southern RSR [17,77]. . Study area with: (a) Regional setting of shoals and channels at Tongzhou Bay situated at the southern Jiangsu coast. Indicated are past reclamation activities, including the outline of the proposed Tongzhou Bay port (red dashed line) [49]. Measurement locations at Tongzhou Bay, where M-stations (blue) are waterlevel gauges and N-stations (red) measured current velocity and direction and sediment concentration among others; (b) Computational domains of the Jiangsu Regional Model (red) and Tongzhou Bay Model (blue); (c) A close-up of the grid domain shows the coupling and different resolutions between the two "online" coupled grids of the Jiangsu coast (JRM) and Tongzhou Bay (TBM).

Ecotope Classification
Since ecotopes describe the relation between the abiotic system and biotic elements (in this case: migratory birds), they are often used as tools in determining maintenance and habitat development policies [78]. Ecotopes are defined by basic hydrodynamic properties (flooding frequency and salinity) and sediment characteristics (mud content). These abiotic parameters vary in time and space, Figure 1. Study area with: (a) Regional setting of shoals and channels at Tongzhou Bay situated at the southern Jiangsu coast. Indicated are past reclamation activities, including the outline of the proposed Tongzhou Bay port (red dashed line) [49]. Measurement locations at Tongzhou Bay, where M-stations (blue) are waterlevel gauges and N-stations (red) measured current velocity and direction and sediment concentration among others; (b) Computational domains of the Jiangsu Regional Model (red) and Tongzhou Bay Model (blue); (c) A close-up of the grid domain shows the coupling and different resolutions between the two "online" coupled grids of the Jiangsu coast (JRM) and Tongzhou Bay (TBM).

Ecotope Classification
Since ecotopes describe the relation between the abiotic system and biotic elements (in this case: migratory birds), they are often used as tools in determining maintenance and habitat development policies [78]. Ecotopes are defined by basic hydrodynamic properties (flooding frequency and salinity) and sediment characteristics (mud content). These abiotic parameters vary in time and space, especially in response to an intervention. However, such modifications in abiotic parameters can be predicted with state-of-the-art geomorphological models. In order to predict the impact of interventions on ecotopes, an existing and well-validated model is used and extended with a morphological module (Section 4.1). We will first describe the ecotopes in more detail hereafter.
Ecotopes in the Tongzhou Bay system are quantified using the ZES.1 method [78], developed for saline open waters ecosystems in the Netherlands, but also applied for the complete Wadden Sea covering the coastline of Germany and Denmark as well [79]. The Bouma et al. [78] classification scheme is based on the flooding frequency, substrate, salinity, and flow velocity. However, the analysis of available field data [49,57] revealed only small variations in salinity as well as a homogenous soft substrate over the Tongzhou Bay domain, and therefore we simplify their classification scheme by omitting subdivisions based on salinity and substrate. Depth classification is based on tidal data (I), such as Mean High Water Spring (MHWS) and Mean Low Water Spring (MLWS). An additional depth sub classification (II) is based on a certain depth below MLWS, consistent to Bouma et al. [78]-see Figure 2. The flow velocity distinguishes low-energy from high-energy waters with a threshold value of 0.8 m/s. This value approximates conditions for the initiation of bed forms, signaling a change in ecotope distribution in the Wadden Sea [78,79]. The classification based on the dryfall period is largely consistent with ZES.1 method, although a lower limit of the high-range littoral of 60% was chosen to better match bird habitat field observations [32,80]. By assigning a color scheme to each ecotope ( Figure 2), an ecotope distribution map can be drawn.
Water 2020, 12, x FOR PEER REVIEW 6 of 22 especially in response to an intervention. However, such modifications in abiotic parameters can be predicted with state-of-the-art geomorphological models. In order to predict the impact of interventions on ecotopes, an existing and well-validated model is used and extended with a morphological module (Section 4.1). We will first describe the ecotopes in more detail hereafter. Ecotopes in the Tongzhou Bay system are quantified using the ZES.1 method [78], developed for saline open waters ecosystems in the Netherlands, but also applied for the complete Wadden Sea covering the coastline of Germany and Denmark as well [79]. The Bouma et al. [78] classification scheme is based on the flooding frequency, substrate, salinity, and flow velocity. However, the analysis of available field data [49,57] revealed only small variations in salinity as well as a homogenous soft substrate over the Tongzhou Bay domain, and therefore we simplify their classification scheme by omitting subdivisions based on salinity and substrate. Depth classification is based on tidal data (I), such as Mean High Water Spring (MHWS) and Mean Low Water Spring (MLWS). An additional depth sub classification (II) is based on a certain depth below MLWS, consistent to Bouma et. al. [78]-see Figure 2. The flow velocity distinguishes low-energy from highenergy waters with a threshold value of 0.8 m/s. This value approximates conditions for the initiation of bed forms, signaling a change in ecotope distribution in the Wadden Sea [78,79]. The classification based on the dryfall period is largely consistent with ZES.1 method, although a lower limit of the high-range littoral of 60% was chosen to better match bird habitat field observations [32,80]. By assigning a color scheme to each ecotope ( Figure 2), an ecotope distribution map can be drawn.

Shorebird Distribution
An ecotope only describes a potential niche for the possible occurrence of a certain habitat, and therefore the ecotope map needs to validated with field observations of certain species (Section 5.1). For this purpose, we use the distributional data of the two selected shorebird species, the Bar-tailed Godwits and Great Knots. These data were acquired through tracking of individual birds by satellitetags (5 g and 9.5 g solar Platform Terminal Transmitters, Microwave Telemetry, Maryland, USA) from their northwest Australian "non-breeding" areas, Roebuck Bay (17.98° S, 122.31° E and Eighty Mile Beach (19.40° S, 121.27° E) (for details on tag deployment, see [32]). We analyze 317 location fixes obtained between April 2015 and May 2018 within the study area, collected from the 13 tracked Bartailed Godwits and four tracked Great Knots which stopped there during their migration. We first retain the locations in the class 3, 2, 1 and A and remove locations of class B as those have an error radius that is too large for the size of our study area [81]. We then select locations collected during low and mid-tide, i.e., when water level is lower than the 70% quantile of a sample of predicted water level (every 10 min for a month) of the study area, assuming that this is the period when the birds are feeding during the tidal-cycle. From these locations we select those that are at least 15 min apart from one another. If there are more than one location within 15 min, we select the point with highest accuracy, or the earliest point in the case of ties. Next, we estimate Kernel Density home ranges from these locations for each species, using the "Heatmap" plugin in QGIS 2.18.11 [82]. The radius of each point is two times the published 68% percentile error radius [81] and is weighted by the inverse of

Shorebird Distribution
An ecotope only describes a potential niche for the possible occurrence of a certain habitat, and therefore the ecotope map needs to validated with field observations of certain species (Section 5.1). For this purpose, we use the distributional data of the two selected shorebird species, the Bar-tailed Godwits and Great Knots. These data were acquired through tracking of individual birds by satellite-tags (5 g and 9.5 g solar Platform Terminal Transmitters, Microwave Telemetry, MD, USA) from their northwest Australian "non-breeding" areas, Roebuck Bay (17.98 • S, 122.31 • E and Eighty Mile Beach (19.40 • S, 121.27 • E) (for details on tag deployment, see [32]). We analyze 317 location fixes obtained between April 2015 and May 2018 within the study area, collected from the 13 tracked Bar-tailed Godwits and four tracked Great Knots which stopped there during their migration. We first retain the locations in the class 3, 2, 1 and A and remove locations of class B as those have an error radius that is too large for the size of our study area [81]. We then select locations collected during low and mid-tide, i.e., when water level is lower than the 70% quantile of a sample of predicted water level (every 10 min for a month) of the study area, assuming that this is the period when the birds are feeding during the tidal-cycle. From these locations we select those that are at least 15 min apart from one another. If there are more than one location within 15 min, we select the point with highest accuracy, or the earliest point in the case of ties. Next, we estimate Kernel Density home ranges from these locations for each species, using the "Heatmap" plugin in QGIS 2.18.11 [82]. The radius of each point is two times the published 68% percentile error radius [81] and is weighted by the inverse of this radius. Thus, each point has the same "heat", but is more concentrated (for class 3 locations) or spread out (for the less precise class 2, 1 and A locations).

Setup
The applied hydrodynamic model is presented by Muller et al. [28], and for details on grid design and bed levels ( Figure 1) we refer thereto. We do briefly reiterate the hydrodynamic model results and subsequently extend the model with a morphodynamic module. The hydroand morphodynamics at Tongzhou Bay are quantified using the Delft3D numerical model [83]. Delft3D solves the Reynolds-averaged Navier-Stokes equations, under the shallow water and the Boussinesq assumption on a staggered grid by a finite difference scheme. The long-term morphologic development of the Jiangsu coast and Tongzhou bay is mainly driven by tidal flows [53,55]. Waves have only a minor effect, which is confirmed by sensitivity simulations with a SWAN-based wave model, but not reported here.
The whole model application includes the Jiangsu coast (JRM domain) to resolve large-scale tidal dynamics, and is refined near Tongzhou Bay (TBM domain) through domain decomposition ( Figure 1b).The JRM is driven by a series of astronomical tides at its two open boundaries. These were derived from a set of 13 tidal harmonic components from a large-scale tidal wave model of the Bohai, Yellow and East China Sea for a full morphological year time-series [56].
The JRM was set-up with four sediment fractions (16,45,90 and 180 µm, resembling fine silt, silt, very fine sand and fine sand respectively), following [84]. An extra cohesive fraction was added to account for finer suspended material prevalent in the shallow coastal waters of Jiangsu coast. Sediment transport is calculated using an adaptation of the transport formulas of van Rijn [85,86] by Yao et al. [84]. Cohesive sediment transport was modelled using the Partheniades [87] erosion equation and a permanent deposition flux.
At the two open boundaries of the JRM, a constant sediment concentration is prescribed: a constant sediment concentration of 0.5 mg/L and 20 mg/L at the eastern boundary and a linear decreasing concentration of 3 mg/L and 200 mg/L (from 0 > z b > −30 m, where z b is the bedlevel) to 0.5 mg/L and 20 mg/L (offshore) at the southern boundary for the 16 µm and cohesive fraction respectively. This was carried out to account for the higher concentrations of fine sediments near the coast that are normally caused by larger residual nearshore currents. Equilibrium sediment concentrations were prescribed for the coarser sand fractions. A 1.5 m thick sediment layer was prescribed at the bed: a thin layer (0.02) m of cohesive sediment and 25% of the remaining 1.48 m for each of the other four fractions. All fractions are vertically mixed in the bed.
The model was calibrated for the period 3 February to 3 March 2012. An initial sediment concentration was prescribed that resembled conditions similar to 2 months of spin-up: 200 mg/L for the cohesive fraction (and zero for sand and silt). During simulation of the different port configurations, the model is set to run for 20 morphological years using a MorFac acceleration technique [88].
The performance of the model was validated by comparing the water level, tidal current velocity, direction and the sediment concentration calculated and measured during spring and neap tidal conditions from 22 to 23 February 2012 (spring tide) and 29 February to 1 March 2012 (neap tide) at 14 mooring stations and six tidal gauges at Tongzhou Bay (Figure 1a).
The model performance was quantified by computing Model Efficiency (ME) and bias (PB). The Nash-Sutcliffe Model Efficiency coefficient [89] is defined as: Water 2020, 12, 2134 where m is the measured value; p is the predicted or calculated value by the model; m is the mean measured value and ∆m is the error in the measurements, as introduced by van Rijn et al. [90]. The model performance can be classified as#excellent (ME > 0.65), very good (0.5 > ME > 0.65), good (0.2 > ME > 0.5) or poor (ME < 0.2) [91,92]. The bias (expressed in percentages) is given by: where the numerator can either be ( The absolute value of the percentage bias is classified in four classes: excellent (|PB| < 10%), very good (10% <|PB| < 20%), good (20%< |PB| < 40%) and poor (|PB| > 40%) [91,92].

Validation
A large number of model runs were executed as part of a model calibration procedure not reported here. Instead, we only present the agreement between data and calibrated model run, focussing on (1) waterlevels, (2) flow velocity and SSC, and (3) ecotope distribution.

Validation
A large number of model runs were executed as part of a model calibration procedure not reported here. Instead, we only present the agreement between data and calibrated model run, focussing on (1) waterlevels, (2) flow velocity and SSC, and (3) ecotope distribution.
Waterlevels were collected at six tidal gauges. The model agreement with these data is excellent based on Model Efficiency scores (ME > 0.65) and good (Xiaomiaohong and Datang station) to very good (other stations) based on the model bias-see Figure 3. Current velocity and SSC were measured at 14 moorings for two consecutive tidal cycles during both spring and during neap. The agreement of the model with the current velocity, magnitude, and SSC is synthesized in Figure 4. The magnitude of the current velocity is slightly better resolved than its direction based on the Model Efficiency (nine stations excellent for magnitude, four for direction). The model bias, however, indicates the direction is better resolved (most stations classified as "very good") than the magnitude (mostly "good"). Stations N10, N12, and N13 provide a cluster of stations performing more poorly than the other stations (albeit still "good" according to both verification metrics). This could be related to errors in the bathymetric data (for which computed velocities are very sensitive) or complex local hydrodynamics. The flow velocity and direction is most poor reproduced at station N6 (scoring "poor"): here, significantly lower flow velocities were measured during ebb (peak flow of < 0.3 m/s) than computed. Station N6 is located in a relatively narrow channel, and therefore the flow velocity varies spatially. When the flow varies greatly on the scale of grid cell size, the model cannot accurately reproduce the observed flow velocity.
The accuracy of the computed SSC varies more among the various station than the hydrodynamics based for the Model Efficiency. At the upper reaches or shallower sections of the channels (i.e., N1, N2, N6 and N8) the model performance is "poor" (ME = −0.2-−0.5). In this case, the model overestimates the concentrations compared to the measured values. The model bias, however, suggests that the computed sediment concentration is at least as accurate as the computed hydrodynamics. Current velocity and SSC were measured at 14 moorings for two consecutive tidal cycles during both spring and during neap. The agreement of the model with the current velocity, magnitude, and SSC is synthesized in Figure 4. The magnitude of the current velocity is slightly better resolved than its direction based on the Model Efficiency (nine stations excellent for magnitude, four for direction). The model bias, however, indicates the direction is better resolved (most stations classified as "very good") than the magnitude (mostly "good"). Stations N10, N12, and N13 provide a cluster of stations performing more poorly than the other stations (albeit still "good" according to both verification metrics). This could be related to errors in the bathymetric data (for which computed velocities are very sensitive) or complex local hydrodynamics. The flow velocity and direction is most poor reproduced at station N6 (scoring "poor"): here, significantly lower flow velocities were measured during ebb (peak flow of < 0.3 m/s) than computed. Station N6 is located in a relatively narrow channel, and therefore the flow velocity varies spatially. When the flow varies greatly on the scale of grid cell size, the model cannot accurately reproduce the observed flow velocity.  The data-model agreement of a relatively poorly performing station (N1, XMH) and a wellperforming station (N7, SSH) is explored in more detail in the time domain ( Figure 5) for spring-neap and intra-tidal variability. The data show that tidal currents at both stations are (1) fairly symmetric (the ebb current velocity is comparable to the flood current velocity) and (2) the spring tidal current velocity is typically twice as high during spring tide than during neap tide. The modelled current velocity reveals the same patterns, which is important for modelling sediments. Although some deviations exist between the data and the model, there is no systematic ebb-flood asymmetry (as in the data).
The measured suspended sediment concentrations vary typically between 150 and 300 mg/L, which are typical values for the Jiangsu coast and Tongzhou Bay [55,59]. SSC varies fairly little over the tidal cycle-especially during neap tides. This suggests the presence of a slow-settling sediment fraction of which the sediment concentration is fairly independent of the hydrodynamics, providing a background sediment concentration. During spring tide, a tidally varying sediment concentration is superimposed on this background concentration, as more sediments are resuspended from the bed. This background concentration was the motivation for adding a cohesive sediment fraction to the model of Yao et al. [84]. The concentration modelled agrees better with the data during spring than during neap tide. At station N1, the measured concentration during neap tide is higher than during spring tide. The reasons for the higher neap tide concentrations are probably related to higher wind velocities during neap tide conditions (up to+m/s) than during spring tide conditions. The measured SSC at station N7 is higher during spring than during neap, and also better reproduced by the model (except for a slight underestimation of neap tide SSC). The model also suggests that the bulk of the suspended sediment is the fine cohesive sediment fraction and the fine silt fraction (spring tide only). The data-model agreement of a relatively poorly performing station (N1, XMH) and a well-performing station (N7, SSH) is explored in more detail in the time domain ( Figure 5) for spring-neap and intra-tidal variability. The data show that tidal currents at both stations are (1) fairly symmetric (the ebb current velocity is comparable to the flood current velocity) and (2) the spring tidal current velocity is typically twice as high during spring tide than during neap tide. The modelled current velocity reveals the same patterns, which is important for modelling sediments. Although some deviations exist between the data and the model, there is no systematic ebb-flood asymmetry (as in the data).
The measured suspended sediment concentrations vary typically between 150 and 300 mg/L, which are typical values for the Jiangsu coast and Tongzhou Bay [55,59]. SSC varies fairly little over the tidal cycle-especially during neap tides. This suggests the presence of a slow-settling sediment fraction of which the sediment concentration is fairly independent of the hydrodynamics, providing a background sediment concentration. During spring tide, a tidally varying sediment concentration is superimposed on this background concentration, as more sediments are resuspended from the bed. This background concentration was the motivation for adding a cohesive sediment fraction to the model of Yao et al. [84]. The concentration modelled agrees better with the data during spring than during neap tide. At station N1, the measured concentration during neap tide is higher than during spring tide. The reasons for the higher neap tide concentrations are probably related to higher wind velocities during neap tide conditions (up to+m/s) than during spring tide conditions. The measured SSC at station N7 is higher during spring than during neap, and also better reproduced by the model (except for a slight underestimation of neap tide SSC). The model also suggests that the bulk of the suspended sediment is the fine cohesive sediment fraction and the fine silt fraction (spring tide only). In absence of morphological data for calibration, we analyze the resulting sediment transport patterns phenomenologically. The total depth-averaged suspended sediment transport was retrieved for four consecutive spring-neap cycles with a 1-h interval ( Figure 6). The net bedload transport was several orders smaller than the net suspended transport and therefore not included in the analysis. Mean suspended transport is directed through XMH and SSH, showing a net import of sediment through the tidal channels. Inside SSH, the suspended transport becomes segregated. A large part of the suspended material crosses over the lower parts of the Lengjiasha shoal and a smaller part is transported up the Sanshahong channel. Although, smaller than the transport rates in the tidal In absence of morphological data for calibration, we analyze the resulting sediment transport patterns phenomenologically. The total depth-averaged suspended sediment transport was retrieved for four consecutive spring-neap cycles with a 1-h interval ( Figure 6). The net bedload transport was several orders smaller than the net suspended transport and therefore not included in the analysis. Mean suspended transport is directed through XMH and SSH, showing a net import of sediment through the tidal channels. Inside SSH, the suspended transport becomes segregated. A large part of the suspended material crosses over the lower parts of the Lengjiasha shoal and a smaller part is transported up the Sanshahong channel. Although, smaller than the transport rates in the tidal channel, a northward flux over Yaosha suggests most of the sediment supply comes from the Xiaomiaohong channel.
Water 2020, 12, x FOR PEER REVIEW  11 of 22 channel, a northward flux over Yaosha suggests most of the sediment supply comes from the Xiaomiaohong channel. Figure 6. Calculated mean total sediment transport at Yaosha for four consecutive spring-neap cycles.
To illustrate the significant mean transport over Yaosha, a lower transports rate limit is included of 7.0·10 −6 m 3 /m/s.

Initial Conditions
The model output is converted into an ecotope map (Figure 7) using the definitions summarized in Figure 2. The main color scale indicates the littoral zones, whereas the low and high hydrodynamic conditions are represented by lower or darker shading. The stronger hydrodynamic conditions roughly coincide with the overall bathymetry of the channels. The Yaosha shoal ranges from lowrange littoral at the outer edges up to high-range and supralittoral close towards the current seawall.  To illustrate the significant mean transport over Yaosha, a lower transports rate limit is included of 7.0 × 10 −6 m 3 /m/s.

Initial Conditions
The model output is converted into an ecotope map (Figure 7) using the definitions summarized in Figure 2. The main color scale indicates the littoral zones, whereas the low and high hydrodynamic conditions are represented by lower or darker shading. The stronger hydrodynamic conditions roughly coincide with the overall bathymetry of the channels. The Yaosha shoal ranges from low-range littoral at the outer edges up to high-range and supralittoral close towards the current seawall.
Water 2020, 12, x FOR PEER REVIEW 11 of 22 channel, a northward flux over Yaosha suggests most of the sediment supply comes from the Xiaomiaohong channel. Figure 6. Calculated mean total sediment transport at Yaosha for four consecutive spring-neap cycles.
To illustrate the significant mean transport over Yaosha, a lower transports rate limit is included of 7.0·10 −6 m 3 /m/s.

Initial Conditions
The model output is converted into an ecotope map (Figure 7) using the definitions summarized in Figure 2. The main color scale indicates the littoral zones, whereas the low and high hydrodynamic conditions are represented by lower or darker shading. The stronger hydrodynamic conditions roughly coincide with the overall bathymetry of the channels. The Yaosha shoal ranges from lowrange littoral at the outer edges up to high-range and supralittoral close towards the current seawall.  To determine ecotopes important for the two shorebird species (Bar-tailed Godwit and Great Knot), we analyse the satellite-tracked distributional data in two spatial scales. First, to show the areas where these species generally occurred within the entire study area, the 90% Kernel density contours, representing the home range of the species, are calculated and overlaid onto the ecotope classifications ( Figure 8). Bar-tailed Godwit and Great Knot occur mostly in the mid-range low dynamic littoral zone, with 30% and 35% of their 90% home range areas overlapping with this ecotope type (Figure 9).
Second, to show ecotopes that are more intensely-used within the 90% home range, we calculate, of the total "heat" (i.e., the sum of Kernel density in each output raster cell) within the 90% home range, the percentage that falls in each ecotope. The "heat" thus reflects the density of bird occurrence at that raster cell. For a specific ecotope, if the percentage of "heat" that falls within it equals the percentage of area, it means that the birds use the ecotope in proportion to its area. If the percentage of "heat" is higher than the percentage of area, it reflects that birds prefer to utilize those ecotopes, and vice versa. Put together, Bar-tailed Godwit occurs mostly in the mid-range low dynamic littoral zone and select for it (38% heat vs. 30% area). While 22% of their 90% home range overlap with the reclaimed zone, it was selected against with only 19% of the heat falls within that ecotope. The low-range low-dynamic littoral zone was used in proportion to its area (14%), and Godwits also show a slight preference for the high-range low-dynamic ecotope (7% heat vs. 6% area). For Great Knots, besides the mid-range ecotope, they show a preference for low-range (17% heat vs. 15% area) and shallow-subtidal range ecotopes (26% heat vs. 20% area).
This species difference resembles visual observations in areas along the Jiangsu coasts [5,25]. While the mid-range littoral zone is important for these species, the other zones occurring in its proximity are also important. The difference in intensity of use of ecotopes between the two species reflect that the higher (high and mid-range) littoral zones and the lower (low-range and shallow subtidal) littoral zones possibly supporting different sets of species that differ in prey preferences and foraging tactics. To determine ecotopes important for the two shorebird species (Bar-tailed Godwit and Great Knot), we analyse the satellite-tracked distributional data in two spatial scales. First, to show the areas where these species generally occurred within the entire study area, the 90% Kernel density contours, representing the home range of the species, are calculated and overlaid onto the ecotope classifications ( Figure 8). Bar-tailed Godwit and Great Knot occur mostly in the mid-range low dynamic littoral zone, with 30% and 35% of their 90% home range areas overlapping with this ecotope type ( Figure 9).
Second, to show ecotopes that are more intensely-used within the 90% home range, we calculate, of the total "heat" (i.e., the sum of Kernel density in each output raster cell) within the 90% home range, the percentage that falls in each ecotope. The "heat" thus reflects the density of bird occurrence at that raster cell. For a specific ecotope, if the percentage of "heat" that falls within it equals the percentage of area, it means that the birds use the ecotope in proportion to its area. If the percentage of "heat" is higher than the percentage of area, it reflects that birds prefer to utilize those ecotopes, and vice versa. Put together, Bar-tailed Godwit occurs mostly in the mid-range low dynamic littoral zone and select for it (38% heat vs. 30% area). While 22% of their 90% home range overlap with the reclaimed zone, it was selected against with only 19% of the heat falls within that ecotope. The lowrange low-dynamic littoral zone was used in proportion to its area (14%), and Godwits also show a slight preference for the high-range low-dynamic ecotope (7% heat vs. 6% area). For Great Knots, besides the mid-range ecotope, they show a preference for low-range (17% heat vs. 15% area) and shallow-subtidal range ecotopes (26% heat vs. 20% area).
This species difference resembles visual observations in areas along the Jiangsu coasts [5,25]. While the mid-range littoral zone is important for these species, the other zones occurring in its proximity are also important. The difference in intensity of use of ecotopes between the two species reflect that the higher (high and mid-range) littoral zones and the lower (low-range and shallow subtidal) littoral zones possibly supporting different sets of species that differ in prey preferences and foraging tactics.

Alternative Port Configurations
The originally proposed Tongzhou Bay port ( Figure 1) will result in the total reclamation of the earlier established high-value ecotopes and an inevitable loss of important habitat. To come to more sustainable alternatives, the BwN recipe is followed by applying the previously gained insight from initial ecotope distribution and driving natural mechanisms, e.g., model results on transport and siltation patterns, which yield the desired transition in ecotope distribution over time.
Three alternative designs for the Tongzhou Bay port are examined, these are alterations in overall shape, extent and use of offshore terminals based on the original port design. All configurations consist of the partial reclamation of the Yaosha shoal, where configuration v1 and v2 both reclaim the northside of Yaosha. The concept of the design is to capture the landward sediment transport though XMH over the southside of Yaosha and increase local siltation at the shoal. Configuration v2 is disconnected from the shoreline, creating a bypass for flow and transports along the reclamation. The third configuration v3, consists of the reclamation of the southern Yaosha shoal. Situated adjacent to the XMH channel, it provides good marine and hinterland accessibility for port development, as well as preserving the initial high-value ecotope at Yaosha. Conversely, this configuration would block the transport from XMH and yield less sedimentation at the shoal.
The impact of these designs on siltation patterns and ecotope distribution was simulated using the Tongzhou Bay model over a period of 20 years (Figure 10). A reference case was defined for the validated 2012 situation, in which no reclamation activity takes place (other than the pre-2012 reclamations). Each variant is assessed on their performance in enhancing the natural growth of the earlier established high-value ecotopes relative to the reference case. Although we do show results of all computed ecotope developments, we only discuss the littoral classification (i.e., supra-, high/mid/low-range littoral and shallow/sub littoral) as these are relevant for shorebirds. The initial and final area of each ecotope of each configuration and reference case are retrieved (Table 1), including the net gains or losses at the end of the simulation, due to direct loss to infrastructure or indirect loss due to temporal turnover of ecotopes. These net gains and losses of ecotope are further examined by gathering the net cumulative growth over time for each configuration ( Figure 10).
For the reference case, strong sedimentation takes place at the higher elevated intertidal areas and the southern edges of Yaosha and Lengjiasha and erosion in the adjacent tidal channels, making tidal flat edges progressively steeper (Figure 10a). Unfortunately, no historic data are available to setup nor validate a morphological hindcast. However, this steepening of the mudflat is also

Alternative Port Configurations
The originally proposed Tongzhou Bay port ( Figure 1) will result in the total reclamation of the earlier established high-value ecotopes and an inevitable loss of important habitat. To come to more sustainable alternatives, the BwN recipe is followed by applying the previously gained insight from initial ecotope distribution and driving natural mechanisms, e.g., model results on transport and siltation patterns, which yield the desired transition in ecotope distribution over time.
Three alternative designs for the Tongzhou Bay port are examined, these are alterations in overall shape, extent and use of offshore terminals based on the original port design. All configurations consist of the partial reclamation of the Yaosha shoal, where configuration v1 and v2 both reclaim the northside of Yaosha. The concept of the design is to capture the landward sediment transport though XMH over the southside of Yaosha and increase local siltation at the shoal. Configuration v2 is disconnected from the shoreline, creating a bypass for flow and transports along the reclamation. The third configuration v3, consists of the reclamation of the southern Yaosha shoal. Situated adjacent to the XMH channel, it provides good marine and hinterland accessibility for port development, as well as preserving the initial high-value ecotope at Yaosha. Conversely, this configuration would block the transport from XMH and yield less sedimentation at the shoal.
The impact of these designs on siltation patterns and ecotope distribution was simulated using the Tongzhou Bay model over a period of 20 years (Figure 10). A reference case was defined for the validated 2012 situation, in which no reclamation activity takes place (other than the pre-2012 reclamations). Each variant is assessed on their performance in enhancing the natural growth of the earlier established high-value ecotopes relative to the reference case. Although we do show results of all computed ecotope developments, we only discuss the littoral classification (i.e., supra-, high/mid/low-range littoral and shallow/sub littoral) as these are relevant for shorebirds. The initial and final area of each ecotope of each configuration and reference case are retrieved (Table 1), including the net gains or losses at the end of the simulation, due to direct loss to infrastructure or indirect loss due to temporal turnover of ecotopes. These net gains and losses of ecotope are further examined by gathering the net cumulative growth over time for each configuration ( Figure 10).
For the reference case, strong sedimentation takes place at the higher elevated intertidal areas and the southern edges of Yaosha and Lengjiasha and erosion in the adjacent tidal channels, making tidal flat edges progressively steeper (Figure 10a). Unfortunately, no historic data are available to setup nor validate a morphological hindcast. However, this steepening of the mudflat is also observed at similar mudflats along the Jiangsu coast, which is driven by (longshore) tidal currents [59,93]. Additionally, the Yaosha shoal-channel system indeed migrates southward [94] (similar to other channel-shoal systems along the Jiangsu coast [95,96]). The siltation at the center of Yaosha can also be seen in the increase of supralittoral at the center of Yaosha (Figure 10b). Comparing the growth rate, the increase of supralittoral (+2282 ha) and mid-range littoral (+991 ha) is much stronger than high-range littoral (+126 ha) (Table 1, Figure 11). As can be observed, the high-range littoral converts faster into supralittoral than it gains from mid-range littoral.
For configuration v1, the connected reclamation leads to higher siltation rates compared to the reference case, with most of the siltation shifting to the backside of Yoasha along the reclamation (Figure 10c). As a consequence, the mid-range littoral develops into high-range littoral and supralittoral and yields a larger growth of supralittoral (+2537 ha), relative to the reference case. The other surrounding intertidal areas do not show significant ecotope transition. However, strong erosion patterns can also be observed at the shoal itself (Figure 10d). This scouring of the shoal and along the edges of the reclamation is driven by a concentrated return flow during ebb. This can also be recognized in the development of high-dynamic sublittoral and an eventual smaller growth of mid-range littoral (+715 ha), compared to the reference case (Table 1, Figure 11).
Configuration v2 shows similar results compared to v1, including higher siltation and development of high-range and supralittoral at the center of Yaosha, as well as the development of ebb-tidal channels ( Figure 10e). However, due to the disconnected reclamation from the shoreline, the tide generates higher flow velocities over the shoal during high tide and local siltation rates are therefore lower. This effect is also notable by comparing the transition of mid-range littoral into high-range and supralittoral for v1 and v2. Due to the disconnection with the mainland, ecotope development tends to stabilize at high-range littoral ( Figure 10f). As such, v2 yields less growth in supralittoral (+2067 ha) and a larger area of high-range littoral (+222 ha) is gained (Table 1, Figure 11).
Although configuration v3 is similar to v2, the siltation rates of Yoasha shoal are lower and less of the mid-range littoral grows into high-range littoral (+55 ha) and supralittoral (+2012 ha), compared with the reference case and other configurations (Table 1, Figure 11). However, higher siltation is predicted at the southern banks of XMH and erosion rates within the channel are lower (Figure 10g). This can be readily explained with the reduction in the tidal prism in the mudflat-channel system and subsequent reductions in tidal flow velocities. The higher siltation rates at the southside are further recognizable by low-range littoral transitioned into mid-range littoral (Figure 10h), yielding a higher growth (+1110 ha) relative to the reference case (Table 1, Figure 11).

Discussion
In this study, we characterized the preferred ecotopes of the two shorebird species, the Great Knot and Bar-tailed Godwit, and we demonstrated the effect of different port designs on promoting the preservation and enhancement of their preferred ecotopes. However, we should be aware that ecotopes (as defined in this study by hydrodynamic conditions) could only characterize a main abiotic determinant of the habitat suitability for these species. Intertidal habitats serve as foraging sites for migratory shorebirds, and ultimately migratory shorebirds need a high density of good quality prey in these sites for fueling up for their migration [34,97]. Although habitats with the right abiotic features are crucial to support a high density of prey organisms, processes that determine prey densities, such as the settlement of the larvae of the prey species, could still have a large impact on prey densities. The relationship between hydrodynamic conditions and sediment properties and prey densities is complex and largely unknown. One example is in Yalu Jiang Estuary, Liaoning, China, where densities of the main prey of the Great Knot and the Bar-tailed Godwit, Potamocorbula laevis, have collapsed and one suspected cause is the changes in hydrodynamic conditions and sediment characteristics caused by a port development adjacent to the site [98]. Moreover, there are other factors determining habitat suitability, such as danger of predation, competition with other species, the level of disturbances in roosts and distances between suitable roosts and foraging sites [99]. In particular, the current preferred ecotopes (mid-tidal zone) happens to be close to the current seawall. Since shorebirds in the study area used aquacultural ponds with shallow water close to the seawall as roosting habitats during high tide, the current ecotope preferences might also be partly driven by their proximity to suitable roosting habitats, because shorebirds tend to forage at habitats close to their roosts, so they can conserve energy in daily commuting [99]. Therefore, our recommended port design would be only the first step towards the conservation of these species. Once the intertidal mudflats with the right abiotic conditions are created, settlement of larvae of the main prey of shorebirds, i.e., certain species of benthic macroinvertebrates (worms and shellfish) can be monitored. And, if necessary, artificial seeding and the translocation of shellfish stocks originating from other

Discussion
In this study, we characterized the preferred ecotopes of the two shorebird species, the Great Knot and Bar-tailed Godwit, and we demonstrated the effect of different port designs on promoting the preservation and enhancement of their preferred ecotopes. However, we should be aware that ecotopes (as defined in this study by hydrodynamic conditions) could only characterize a main abiotic determinant of the habitat suitability for these species. Intertidal habitats serve as foraging sites for migratory shorebirds, and ultimately migratory shorebirds need a high density of good quality prey in these sites for fueling up for their migration [34,97]. Although habitats with the right abiotic features are crucial to support a high density of prey organisms, processes that determine prey densities, such as the settlement of the larvae of the prey species, could still have a large impact on prey densities. The relationship between hydrodynamic conditions and sediment properties and prey densities is complex and largely unknown. One example is in Yalu Jiang Estuary, Liaoning, China, where densities of the main prey of the Great Knot and the Bar-tailed Godwit, Potamocorbula laevis, have collapsed and one suspected cause is the changes in hydrodynamic conditions and sediment characteristics caused by a port development adjacent to the site [98]. Moreover, there are other factors determining habitat suitability, such as danger of predation, competition with other species, the level of disturbances in roosts and distances between suitable roosts and foraging sites [99]. In particular, the current preferred ecotopes (mid-tidal zone) happens to be close to the current seawall. Since shorebirds in the study area used aquacultural ponds with shallow water close to the seawall as roosting habitats during high tide, the current ecotope preferences might also be partly driven by their proximity to suitable roosting habitats, because shorebirds tend to forage at habitats close to their roosts, so they can conserve energy in daily commuting [99]. Therefore, our recommended port design would be only the first step towards the conservation of these species. Once the intertidal mudflats with the right abiotic conditions are created, settlement of larvae of the main prey of shorebirds, i.e., certain species of benthic macroinvertebrates (worms and shellfish) can be monitored. And, if necessary, artificial seeding and the translocation of shellfish stocks originating from other mudflats, which are very common aquaculture practices in China [100], could be carried out. Moreover, reducing human disturbances and creating suitable roosting sites close to the foraging habitats could be the next steps to improve habitat quality for migratory shorebirds in the study area.
This contribution presents a new approach for the coastal management of the Jiangsu coast. For this large stretch of coast that has high ecological value, the current method would allow for a more sustainable port design with better insight into the combination of design requirements. Moreover, we see many opportunities for further improvement. The ecotope-species comparison could be expanded with more species of birds of variable ecological requirements, and even other taxa on lower trophic levels. A further improvement could be to use the bird distributional data to determine the boundary values of the parameters (hydrodynamics, depth) of the ecotopes, to create a species-specific "ecotope"/"species-tope". One way to achieve this is to integrate tools in the ecological literature, such as species distribution models [101], into our current approach. This allows for the direct calculation of how much preferable habitat of a species is lost/gained for each port design. Additionally, this contribution did not use all requirements for future potential port development reclamation. Including a broader set of requirements will lead to more robust solutions, and other options such as recreational needs can also be incorporated.

Conclusions
In this study, we propose and follow a novel approach to engineer new sustainable solutions for coastal infrastructure development and applied it to the Tongzhou Bay port case. System knowledge was acquired in the form of hydro-morphodynamical modelling, ecotope distribution and unique bird spatial distributional data. Gaining a better understanding of the natural dynamics and ecosystem services led to a potential solution that would not only allow for further port development, but also enhance the natural growth of valuable habitat. The following conclusions are made:

1.
Full reclamation of the main shoal (as originally planned) will result in the loss of high-value ecotopes and the populations of birds utilizing them. Partial reclamation of the shoal will reduce the loss of valuable ecotopes, especially when accounting for morphological feedbacks. Partially reclaiming the south side of the shoal, for instance, leads to a reduction in currents and transport fluxes whereas partial reclamation of the north side strengthens siltation rates (which is positive for habitat suitability).

2.
In order to increase the accuracy of the created ecotope map, it is recommended to increase the set of abiotic features and data sources for the ecotope classification and subsequent validation. More input of ecological and societal needs will increase the potential of the alternative design. Relevant and consistent measuring campaigns are important to achieve a better insight.

3.
The BwN approach allows for multi-value design, where economic-social and natural needs can be included at an early stage. Multidisciplinary cooperation between different specialties will lead to a more successful end solution in terms of sustainability and societal support.