Impact of Anthropogenic Activities and Sea Level Rise on a Lagoon System: Model and Field Observations

: The long-term geomorphological evolution of a coastal lagoon is driven by hydrodynamic forcing and is inﬂuenced by climate changes and human activities. In this study, a numerical model of the Qilihai lagoon (QL) system was established based on ﬁeld measurements, previous hydrology data and satellite remote sensing measurements, to simulate the geomorphological evolution of QL from 1900 to 2018. The inﬂuences of sea level rise, runoff and human activities on the evolution of geomorphology were investigated. The results of the model show that the construction projects including the tide gate, the bridge, reclamation and the straightening or widening of the tidal channel increased the net deposition within the QL system. Furthermore, the spatial distribution of tidal asymmetry during the natural time period was similar to that of the change in bed thickness. However, bed erosion or deposition was not only dependent on tidal asymmetry but it was also affected by the external sediment supply and the discharge of upstream rivers. Moreover, sea level rise had a signiﬁcant effect on the tidal asymmetry; therefore, it enhanced the accumulation of sediments in the QL system, while runoff had little effect on the tidal asymmetry or geomorphological changes in the system.


Introduction
Coastal lagoons are highly productive coastal systems in the transitional area between ocean and land that provide a range of ecosystem services and habitats for many species that society values [1][2][3][4]. The tidal inlet provides the major connection from the lagoon to the open sea, and the barrier means that the lagoon becomes a semi-enclosed water [5,6]. A stable hydrodynamic condition and a shallower bathymetry of the lagoon provide a biologically suitable natural environment [7][8][9]. Geomorphological evolution of lagoon is determined by terrestrial and marine sediment transport driven by hydrodynamic forces which can be wave-dominated and tide-dominated [10,11]. The spatial distribution of sediment is mainly dictated by the asymmetrical structure of flow induced by runoff, tide and wave. Because coastal areas are densely populated and vital ecosystems are at the same time important for biodiversity [12], extensive human interventions, including reclamation and tidal gate and bridge construction have occurred in the lagoon and other coastal zones. These human activities as well as climate changes including sea level rises may seriously threaten the physical and ecological environment of the coastal lagoon. Understanding and predicting geomorphological evolution is critical for assessing the influence of human activities and climatic variations on the ecological system of lagoons [13].
Geomorphological evolution mechanisms of many lagoons around the world have been studied, for example, the Venice Lagoon [14], the Ria de Aveiro Lagoon, Portugal [15] and the Laolonggou Lagoon, China [16]. Tidal inlets and barriers around lagoons are the major controlling factors for the hydrodynamics in the lagoon and have attracted lots of attention recently, especially if the tidal inlet acts as the key factor for the dynamic equilibrium in the lagoon system [17]. At the same time, the morphological development of tidal creeks and shoals inside the lagoon have mutual interactions with the hydrodynamics within the lagoon system [18,19]. Constructions of seawalls, reclamations and dams decrease the tidal prism, which is the primary driving force for the narrowing of the tidal inlet and swing [20][21][22]. The tidal amplitude and the water circulation inside the lagoon increase with the size of the tidal prism [12]. Tidal creeks and shoals make up tidal networks in lagoons habituated by variety of species, and a large number of tidal creeks develop simultaneously and interlace with each other from time to time [23]. The networks are shaped by tidal currents mainly but may also be affected by waves and river discharge [24]. Physical experiments and numerical simulations have been used to study the evolution of tidal networks over flat or sloping beds. Initial bathymetry, friction coefficients, bed slopes and tidal ranges are deemed to be determining factors for geomorphological changes and final bathymetry [25][26][27]. Stefanon et al. [28] conducted laboratory experiments and found that the synthetic tidal networks develop via headward growth driven by the exceedance of a critical bottom shear stress and reproduced the statistical network characteristics of geomorphic relevance.
It is well known that tidal asymmetry is the primary driving force for net sediment transport and geomorphological evolution in a tide-dominated system. River discharge, current speed and water level have a great influence on the propagation and deformation of the tide [29][30][31]. Flood and ebb dominances contribute to net sediment import and export, respectively. Geomorphological responses to human activities are more obvious in semi-enclosed waters than open seas [32,33]. The hydrodynamic variation from 1935 to 2008 in Jiaozhou Bay, China indicates that deposition appears with increased flood dominance caused by reclamation, and similar behavior was observed in Marennes-Oléron Bay, France, Murray Mouth, Australia and Xiangshan Bay, China [20,[34][35][36]. The change from a flood to an ebb dominance situation will probably reduce deposition in the Nador Lagoon (Morocco), namely the cost of dredging the tidal inlet [12]. The potential loss or accumulation of sediment in the lagoon should be monitored and modelled.
The Qilihai Lagoon (QL) is located in north of the Luanhe Estuary in the Bohai Sea, northeastern China. It is close to Qinhuangdao (QHD) city which has been voted as one of the top ten ecological cities in China and the most well recognized coastal resort. QL is of great ecological value for this region, and an investigation into the geomorphological evolution of QL is beneficial to ecological protection, planning and construction. The tidal channel and terrain of QL has formed naturally and has been transformed by human activities in the past decades. Our previous study evaluated the short-term water exchange capacity of QL influenced by runoff, water level and ice, using the water flushing and residence time [37]. The main objective of this study was to investigate the influence of sea level rise and anthropogenic activities on the geomorphological evolution of QL in a medium-term and long-term period from 1950 to 2018 using a coupled hydrodynamic and sediment transport model.

Study Area
As the conservation area of QHD city, Hebei province, China, QL is well known for its abundant natural resources, and is located in the north of the Luanhe Estuary in the Bohai Sea (39 • 35 N, 119 • 15 E) ( Figure 1). Frequent human activities such as aquaculture and reclamation have decreased the surface water area and aggravated water pollution. The lagoon basin covers an area of 16 km 2 whereas the water area is only 2.3 km 2 , and there are four ephemeral streams (Zhaojiaganggou, Nijinggou, Liutuogou and Daozigou) discharging into the lagoon [38]. An irregular semi-diurnal tide in QL and adjacent sea is actuated by M 2 tidal constituent, and the low average tidal range of 0.7 m is caused by the amphidromic point near QL [39]. A shallow lagoon with an average depth of 0.5 m makes the tidal flat a transfer station for migratory birds, and a great area of tidal flat emerges during the ebb tide. A long tidal channel of Xinkaikou is the unique link between the adjacent sea and the internal lagoon with a longitudinal length of 1800 m, a width of 80~300 m and an average depth of 3 m, which is used as a fishing port nowadays. The sheltered water is formed by a long tidal channel with a pair of breakwaters, so that the influence of waves on the hydrodynamics within the lagoon system is negligible. J. Mar. Sci. Eng. 2021, 9, x FOR PEER REVIEW 3 of 26 The lagoon basin covers an area of 16 km 2 whereas the water area is only 2.3 km 2 , and there are four ephemeral streams (Zhaojiaganggou, Nijinggou, Liutuogou and Daozigou) discharging into the lagoon [38]. An irregular semi-diurnal tide in QL and adjacent sea is actuated by M2 tidal constituent, and the low average tidal range of 0.7 m is caused by the amphidromic point near QL [39]. A shallow lagoon with an average depth of 0.5 m makes the tidal flat a transfer station for migratory birds, and a great area of tidal flat emerges during the ebb tide. A long tidal channel of Xinkaikou is the unique link between the adjacent sea and the internal lagoon with a longitudinal length of 1800 m, a width of 80~300 m and an average depth of 3 m, which is used as a fishing port nowadays. The sheltered water is formed by a long tidal channel with a pair of breakwaters, so that the influence of waves on the hydrodynamics within the lagoon system is negligible. Similar to the previous study of the Venice lagoon by Carniello et al., the present study is based on a simulation of the geomorphological evolution of the QL system in combination with field observations and satellite imagery [40,41]. The early maps recorded in ancient books are diagrammatic drawings without accurate shorelines. Since satellites were launched in 1950s, intuitive maps became available. We have collected satellite images since the image of QL was first recorded by the Keyhole satellite in 1970. The other remote-sensing images were recorded from Landsat, Airbus, NASA and Maxar Technologies and they were downloaded from Google Earth. The earliest Google Earth image is from 1984; the resolution of images was low before 2007 and became high from 2008 to 2018. The photos of QL are shown in Figure 2. Partial shoals in the south-west lagoon were embanked while the tidal inlet kept its nature without transformation in 1970 (Figure 2a), and the water area was 6.2 km 2 . On account of the missing images from 1971 to 1983, the development process could only be deduced from the literature. The rivers into lagoon were regulated for floodway and irrigation from 1970 to 1972. A project of land reclamation with dikes over 30km and a tide gate in tidal channel was executed from 1975 to 1979, which attempted to change the lagoon into a freshwater reservoir. However, the project was halted in 1979 due to ecological concerns and the retained tide gate has remained open since then [42]. The obvious differences between 1970 and 1984 ( Figure 2b) were concentrated in the tidal channel. The tide gate was built across the mouth adjacent to the lagoon and impeded the flow. The alongshore current and weakened lagoon discharge moved the mouth connection to the open sea and the main channel to the south slightly. Similar to the previous study of the Venice lagoon by Carniello et al., the present study is based on a simulation of the geomorphological evolution of the QL system in combination with field observations and satellite imagery [40,41]. The early maps recorded in ancient books are diagrammatic drawings without accurate shorelines. Since satellites were launched in 1950s, intuitive maps became available. We have collected satellite images since the image of QL was first recorded by the Keyhole satellite in 1970. The other remote-sensing images were recorded from Landsat, Airbus, NASA and Maxar Technologies and they were downloaded from Google Earth. The earliest Google Earth image is from 1984; the resolution of images was low before 2007 and became high from 2008 to 2018. The photos of QL are shown in Figure 2. Partial shoals in the south-west lagoon were embanked while the tidal inlet kept its nature without transformation in 1970 (Figure 2a), and the water area was 6.2 km 2 . On account of the missing images from 1971 to 1983, the development process could only be deduced from the literature. The rivers into lagoon were regulated for floodway and irrigation from 1970 to 1972. A project of land reclamation with dikes over 30km and a tide gate in tidal channel was executed from 1975 to 1979, which attempted to change the lagoon into a freshwater reservoir. However, the project was halted in 1979 due to ecological concerns and the retained tide gate has remained open since then [42]. The obvious differences between 1970 and 1984 ( Figure 2b) were concentrated in the tidal channel. The tide gate was built across the mouth adjacent to the lagoon and impeded the flow. The alongshore current and weakened lagoon discharge moved the mouth connection to the open sea and the main channel to the south slightly.
From 1984 to 1990, the reclamation was mainly used for aquaculture and the water area decreased from 6.2 km 2 to 5.2 km 2 . A channel regulation project in the tidal channel, as shown in the image from 1991 (Figure 2c), straightened the tidal channel by dredging and blocking off the original channel. Meanwhile, a pair of breakwaters were built outside the tidal inlet. The water area decreased to 2.6 km 2 from 1991 to 2000 by expanding aquaculture ponds (Figure 2d), but the water area reduction slowed down in 2001 to 2008. The upper half part of tide gate was dismantled in 2008 and the concrete pier under the water was retained with a top elevation of −1.5 m. At the same time, an overflow bridge was built near the tide gate (Figure 2e). After 2008, the tidal channel was widened and the average width increased from 100 m to 150 m, which enhanced the water exchange capacity. The protruding jetties were built to prevent the beach from eroding in 2013, and the latest remote-sensing image was recorded in 2018 ( Figure 2f). From 1984 to 1990, the reclamation was mainly used for aquaculture and the water area decreased from 6.2 km 2 to 5.2 km 2 . A channel regulation project in the tidal channel, as shown in the image from 1991 (Figure 2c), straightened the tidal channel by dredging and blocking off the original channel. Meanwhile, a pair of breakwaters were built outside the tidal inlet. The water area decreased to 2.6 km 2 from 1991 to 2000 by expanding aquaculture ponds (Figure 2d), but the water area reduction slowed down in 2001 to 2008. The upper half part of tide gate was dismantled in 2008 and the concrete pier under the water was retained with a top elevation of −1.5m. At the same time, an overflow bridge was built near the tide gate (Figure 2e). After 2008, the tidal channel was widened and the average width increased from 100m to 150m, which enhanced the water exchange capacity. The protruding jetties were built to prevent the beach from eroding in 2013, and the latest remote-sensing image was recorded in 2018 ( Figure 2f). The bathymetric mapping of QL was executed by the Tianjin North China Geological Exploration Bureau in 2016. The horizontal control survey is based on VRS (virtual reference station) technique [43], and digitized terrain measurement is through RTK (real time kinematic) positioning and SOKKIA (a brand name in Japan) total station. The projected coordinate system is Beijing 1954 (120 °E), and the horizontal measurement step is 5~15 m. The tidal channel features the deeper portion with average elevation of −3 m. The deepest area is near the bridge and tide gate with water depth around 5 m. The tidal creeks are significant inside the lagoon, whose widths become narrow and whose depths become shallower closer to river mouth. A major tidal creek has formed at the entrance to lagoon and divides into three bypass tidal creeks. The shallowest area is facing the tidal channel The bathymetric mapping of QL was executed by the Tianjin North China Geological Exploration Bureau in 2016. The horizontal control survey is based on VRS (virtual reference station) technique [43], and digitized terrain measurement is through RTK (real time kinematic) positioning and SOKKIA (a brand name in Japan) total station. The projected coordinate system is Beijing 1954 (120 • E), and the horizontal measurement step is 5~15 m. The tidal channel features the deeper portion with average elevation of −3 m. The deepest area is near the bridge and tide gate with water depth around 5 m. The tidal creeks are significant inside the lagoon, whose widths become narrow and whose depths become shallower closer to river mouth. A major tidal creek has formed at the entrance to lagoon and divides into three bypass tidal creeks. The shallowest area is facing the tidal channel with an average elevation of 0.5 m, which is submerged under water at high tide. Most areas of the lagoon have a shallow terrain except the tidal creeks, and the average water depth is 0.8 m.
The sediment characteristics play a vital role in bed roughness and evolution. The sediment in QL is unevenly distributed and the sediment diameter ranges from 0.001 mm to 0.314 mm. The area of silty clay with diameters below 0.063 mm occupies 62% of QL and is mainly located within lagoon and river. However, the sandy area occupies 38% of QL, with sand diameters greater than 0.2 mm mainly located within the tidal channel near the open sea, and fine sands with diameters of 0.063~0.2 mm are mainly located in the transition region between the tidal channel and lagoon.

Methodology
For shallow water with a water depth of~0.5 m at low tide, the shallow water equations and the associated depth-averaged velocities may capture the hydrodynamics in QL adequately. The MIKE (Michael Barry Abbott, the name of an important founder in water numerical modeling industry) mathematical model has been widely used in the simulation and study of hydrodynamic, sediment and pollutant transport, and geomorphological evolution at multiple time scales in coastal waters, estuarine deltas and lake areas [44][45][46][47]. The two-dimensional coupled hydrodynamic and sediment transport models in the MIKE suite (https://www.mikepoweredbydhi.com/, accessed on 3 December 2021) were used to examine the morphodynamics at QL in the past century.

Model Setup
Assuming Boussinesq hydrostatic pressure and shallow water, the incompressible Reynolds averaged Navier-Stokes equation was reduced to the shallow water equations. The hydrodynamic model was used to simulate ocean circulation driven by tide, wind and runoff. To achieve a stable model simulation and the Courant Friedrich Levy (CFL) limited at 0.8, the adaptive time step ranging from 0.001 to 60 s was used in the model runs. The wet-dry dynamic boundary processing technique was used to resolve the submerged surface under water. The bottom stress τ b = (τ x , τ y ) was calculated by the quadratic friction law: where ρ 0 is the reference density of water; c b is the drag coefficient; u s = (u, v) is the depth-average flow velocity; g is the gravitational acceleration; h is the water depth; n b is the Manning roughness coefficient. The n b of Bohai Sea and QHD Sea were set as 0.0135 according to a previous study [39,48], while the n b in QL was chosen as 0.0167, determined by the average sediment size at the site. A three-level nested unstructured finite element mesh was constructed to cover QL and Bohai Sea and QHD Sea as shown in Figure 1, which could resolve complicated terrain and irregular coastline accurately. The simulation was run on a high-performance computer with 24 cores. It took about 15 days to run each of the 5 cases in this study. Over the simulation period, varying meshes were adopted as the shoreline evolved, the coarse mesh (Level 1) covered the entire Bohai Sea to ensure a correct tidal flow field whose boundary was located through the Bohai Strait. The intermediate mesh (Level 2) covered the QHD Sea and tidal inlet of QL, whose boundary conditions were provided from the model output of the Level 1 domain. The fine mesh (Level 3) covered QL and its four adjacent rivers and had a mesh size of 3~20 m which was sufficient to obtain a high-resolution flow field and terrain variation.
The aquaculture ponds and reclamation areas were treated as impermeable zones through closed boundaries in the model system. The open boundaries at upstream rivers and tidal inlet are driven by runoff and tide, respectively. The sea level rise recorded by the tidal gauge at the nearby Tanggu station [49] and the China Sea Level Bulletin (Figure 3a) was adopted to examine the influence of sea level rise in the model simulation. Because the sea level rise at Dalian and Yantai station at the ends of Bohai Strait is the same as that at the Tanggu station [39], the water level at the offshore boundary of Bohai Sea was set with the tidal table from National Marine Data and Information Service of China in 2016 [50] plus the sea level rise at the Tanggu station. The hourly sea level variation was obtained by linear interpolation within a year, which reflected the process of the rise in sea level. To avoid the instabilities caused by using the water level boundary condition, at the exchange boundary between the computation domains of Level 2 and Level 1 or Level 3, we adopted the Flather condition [51] specified by the velocities in the xand y-direction and water level, which are extracted from Level 1 and Level 2, respectively. The interface boundaries between the lagoon and the four adjacent rivers are controlled by the monthly discharge, but the runoff discharge was measured only in 2017. The four rivers belong to the Luanhe River system, so that their runoff discharge could be estimated by the ratio of discharge among them and based on the total runoff discharge in Luanhe River recorded at Luanxian Station from 1950 to 2009 [52] by National Hydrological Statistics Annual Report ( Figure 3b). Sediment supply from rivers mainly comes from the watershed surface and riverbed erosion. The rivers entering the QL are tributaries formed by the ancient Luanhe River inlet channel, which are narrow and shallow; therefore, hardly any sediment in the upstream of the Luanhe River eventually enters the lagoon. Thus, the sediment supply from river tributaries to the lagoon was neglected at the river-lagoon boundary in the model set-up. Table 1 shows the ratio between the runoff discharge in upstream rivers and that in Luanhe River. The initial terrain, dynamic conditions (water level, current velocity) and sediment conditions (concentration, bed thickness) at each stage were based on the simulation results at the end of the previous stage, and other parameters were kept the same.
The aquaculture ponds and reclamation areas were treated as impermeable zones through closed boundaries in the model system. The open boundaries at upstream rivers and tidal inlet are driven by runoff and tide, respectively. The sea level rise recorded by the tidal gauge at the nearby Tanggu station [49] and the China Sea Level Bulletin ( Figure  3a) was adopted to examine the influence of sea level rise in the model simulation. Because the sea level rise at Dalian and Yantai station at the ends of Bohai Strait is the same as that at the Tanggu station [39], the water level at the offshore boundary of Bohai Sea was set with the tidal table from National Marine Data and Information Service of China in 2016 [50] plus the sea level rise at the Tanggu station. The hourly sea level variation was obtained by linear interpolation within a year, which reflected the process of the rise in sea level. To avoid the instabilities caused by using the water level boundary condition, at the exchange boundary between the computation domains of Level 2 and Level 1 or Level 3, we adopted the Flather condition [51] specified by the velocities in the x-and y-direction and water level, which are extracted from Level 1 and Level 2, respectively. The interface boundaries between the lagoon and the four adjacent rivers are controlled by the monthly discharge, but the runoff discharge was measured only in 2017. The four rivers belong to the Luanhe River system, so that their runoff discharge could be estimated by the ratio of discharge among them and based on the total runoff discharge in Luanhe River recorded at Luanxian Station from 1950 to 2009 [52] by National Hydrological Statistics Annual Report (Figure 3b). Sediment supply from rivers mainly comes from the watershed surface and riverbed erosion. The rivers entering the QL are tributaries formed by the ancient Luanhe River inlet channel, which are narrow and shallow; therefore, hardly any sediment in the upstream of the Luanhe River eventually enters the lagoon. Thus, the sediment supply from river tributaries to the lagoon was neglected at the river-lagoon boundary in the model set-up. Table 1 shows the ratio between the runoff discharge in upstream rivers and that in Luanhe River. The initial terrain, dynamic conditions (water level, current velocity) and sediment conditions (concentration, bed thickness) at each stage were based on the simulation results at the end of the previous stage, and other parameters were kept the same.   20 1950 1954 1958 1962 1966 1970 1974 1978 1982 1986 1990 1994 1998 1950 1954 1958 1962 1966 1970 1974 1978 1982 1986 1990 1994 1998  The sediment transport simulation was based on the advection-dispersion equation. The sediment settling velocity was determined by particle diameter and calculated by the Stokes law [53]. Different sediment characteristics were captured by variable settling velocity. Influenced by organic matter content and salinity, flocculation gathers multiple sediment particles into a group with size up to 10 times the particle size while the sediment is silty clay [54]. The average sediment size of silty clay is 0.003 mm so that the reference size with flocculation is 0.03 mm with the corresponding settling velocity of 0.001 m/s. The average sediment size of sand is 0.239 mm with the settling velocity of 0.051 m/s. The critical erosion and deposition shear stresses for sediment were the key controllers for the geomorphological simulation. The critical bed shear stress τ ce is related to the sediment incipient velocity u ce through [55]: The incipient velocity u ce is calculated by [56]: where ρ s is the sediment dry density; d is the sediment grain size; υ is the coefficient of viscosity. The erosion of different sediment fractions in the model are driven by the same erosion shear stress, so that the critical shear stresses of both erosion and deposition are determined by the average sediment diameter of 0.093 mm. The associated sediment incipient velocity is 0.097 m/s, then the critical erosion stress is 0.12 N/m 2 . The empirical critical deposition stress is assumed to be four-ninths of the critical erosion stress, so the critical deposition stress is 0.05 N/m 2 .

Model Verification
The current, water level and suspended sediment concentration were measured hourly to validate the model results at the field sites shown in Figure 4a. The water level was collected at the WL station in tidal channel from 23 September 2016 to 23 October 2016 by a limnimeter (DCX-22 series, Keller, Zurich, Switzerland), while the current and suspended sediment concentrations were measured at V1~V4 stations on 23 September 2016 by ADCP (Acoustic Doppler Current Profiler, Nortek, Vangkroken, Norway). V1 and V2 are located in the tidal channel but V3 and V4 are in the inner QL. The variation of the current magnitude was evident near and inside tidal channel and the average current speed inside the lagoon was lower than 0.2 m/s. The current and suspended sediment between V1 or V2 inside tidal channel and V3 or V4 near the mouth to lagoon differed from each other; therefore, they were better suited for model verifications. Figure 4b-e show great linear relationships between the measured and computed data of current speed and direction, water level and suspended sediment concentration. The percent bias (PBIAS) index values were 12.83% (current speed), 2.36% (current direction), 11.08% (water level) and 11.59% (suspended sediment concentration), respectively. These values were lower than 20% and satisfied the accuracy requirement of the simulation [57,58]. These values were lower than 20% and satisfied the accuracy requirement of the simulation [57,58].

Geomorphological Evolution Process
The historical data of sea level and runoff discharge have been collected since 1950. The geomorphological simulation with a speedup technique for QL during 1900~1950 was based on an initial plane bathymetry of −0.5 m, which was the average bed elevation of QL in 2016. The speedup technique can effectively reduce the simulation time by the speedup factor. If the factor is set to 100, the model result reflects the geomorphological evolution in 100 years by setting a period of 1 year. The first satellite image available in 1970 was used to obtain the coastline for simulation before 1970. It was found that the terrain of QL could have reached a relatively equilibrium state within about 50 years [59], which was used for the initial terrain of 1950 ( Figure 5). The predicted tidal networks were similar to those indicated in the satellite image in 1970 (Figure 2a), which shows that the model could predict the basic geomorphology characteristics.
The simulation contained process-based computations considering variations of sea level, runoff and human activities. The baseline map was used to determine the shoreline in the simulation and shoreline variation over a few years was generalized into several types of scenarios and the corresponding maps for different typical periods are illustrated in Table 2. The geomorphological simulation was executed for five phases, i.e., 1950~1969, 1970~1979, 1980~1990, 1991~2008 and 2009~2018. The water area and tidal channel of QL kept natural evolution conditions from 1950 to 1969. The lack of satellite images of QL from 1970 to 1979 meant that the simulation of geomorphological evolution for this time period was based on the satellite image of 1984. The tide gate at the entrance of the lagoon played a double role in both water storage and flood water discharging during flooding and an ebb tide. The tide gate acted as an overflow weir, whose crest level was 0.5 m above mean water level. The coastline from 1980 to 1989 based on the satellite image of 1990 was used for the model simulation of this time period. The reclamation and tide gate opening and close were also incorporated in the model. The tide gate opening was taken as an overflow weir with the crest level changed to −1.5 m. The coastline from 1990 to 2007 was

Geomorphological Evolution Process
The historical data of sea level and runoff discharge have been collected since 1950. The geomorphological simulation with a speedup technique for QL during 1900~1950 was based on an initial plane bathymetry of −0.5 m, which was the average bed elevation of QL in 2016. The speedup technique can effectively reduce the simulation time by the speedup factor. If the factor is set to 100, the model result reflects the geomorphological evolution in 100 years by setting a period of 1 year. The first satellite image available in 1970 was used to obtain the coastline for simulation before 1970. It was found that the terrain of QL could have reached a relatively equilibrium state within about 50 years [59], which was used for the initial terrain of 1950 ( Figure 5). The predicted tidal networks were similar to those indicated in the satellite image in 1970 (Figure 2a), which shows that the model could predict the basic geomorphology characteristics.
The simulation contained process-based computations considering variations of sea level, runoff and human activities. The baseline map was used to determine the shoreline in the simulation and shoreline variation over a few years was generalized into several types of scenarios and the corresponding maps for different typical periods are illustrated in Table 2. The geomorphological simulation was executed for five phases, i.e., 1950~1969, 1970~1979, 1980~1990, 1991~2008 and 2009~2018. The water area and tidal channel of QL kept natural evolution conditions from 1950 to 1969. The lack of satellite images of QL from 1970 to 1979 meant that the simulation of geomorphological evolution for this time period was based on the satellite image of 1984. The tide gate at the entrance of the lagoon played a double role in both water storage and flood water discharging during flooding and an ebb tide. The tide gate acted as an overflow weir, whose crest level was 0.5 m above mean water level. The coastline from 1980 to 1989 based on the satellite image of 1990 was used for the model simulation of this time period. The reclamation and tide gate opening and close were also incorporated in the model. The tide gate opening was taken as an overflow weir with the crest level changed to −1.

1950~1969
Based on the initial terrain which was relatively stable, the geomorphological evolution of QL from 1950 to 1969 was simulated considering the effects of sea level and runoff discharge. The terrain on 1 January 1970 and the change in bed thickness from 1950 to 1969 were obtained, as shown in Figure 6. Compared with the initial terrain of 1950, a narrower width of the tidal creek and a higher elevation of the tidal flat were observed in 1970. Tidal creek B extended deeper to Daozigou and tidal creek D took a more obvious shape. However, there was little change in the geomorphological characteristics of the whole tidal network. Erosion mainly occurred in the tidal channel and tidal creeks, while deposition appeared in the shallowest water area. The area of deposition was larger than that of erosion, and the main deposition area was concentrated in the southeastern tidal inlet and both sides of the tidal creek.

1950~1969
Based on the initial terrain which was relatively stable, the geomorphological evolution of QL from 1950 to 1969 was simulated considering the effects of sea level and runoff discharge. The terrain on 1 January 1970 and the change in bed thickness from 1950 to 1969 were obtained, as shown in Figure 6. Compared with the initial terrain of 1950, a narrower width of the tidal creek and a higher elevation of the tidal flat were observed in 1970. Tidal creek B extended deeper to Daozigou and tidal creek D took a more obvious shape. However, there was little change in the geomorphological characteristics of the whole tidal network. Erosion mainly occurred in the tidal channel and tidal creeks, while deposition appeared in the shallowest water area. The area of deposition was larger than that of erosion, and the main deposition area was concentrated in the southeastern tidal inlet and both sides of the tidal creek.
Six observation points in Figure 7 were used to quantitatively study the change process of bed thickness. It can be seen from T1 that the bed thickness decreased continuously, indicating that the tidal channel was in a state of continuous erosion. By 1970, the erosion thickness of T1 was 1.139 m, which was also the largest because the average velocity of T1 was 0.297 m/s, the largest among the six points. T2 was located near the entrance to the lagoon, whose erosion thickness increased with the decreasing rate from 0.039 m/y to 0.022 m/y and reached 0.608 m in 1970. T3 was located at the tidal flat on the northern tidal creek A, whose deposition thickness increased slowly with the reducing rate and the average deposition rate was 0.006 m/y. T4 was located in tidal creek D, and the bed thickness change fluctuated. The erosion thickness reached the maximum of 0.375 m in 1961. After 1968, the erosion thickness was basically maintained at 0.248 m. T5 and T6 represent the changes of tidal creek C and B, respectively, whose terrains were both under continuous erosion. After 1967, the change rate in the erosion thickness was less than 0.010 m/y at T5, and the erosion thickness was 0.652 m in 1970. The average annual erosion thickness at T6 in tidal creek B was 0.030 m/y, and the maximum change rate was 0.072 m/y in 1961. It can also be seen from T4 and T5 that the change rate in 1961 had a sudden increase with 0.171 m/y and 0.046 m/y, respectively, because the runoff discharge in 1960 was the maximum value compared with that in any other years. The runoff discharge of the Luanxian Station was the maximum of 300 m 3 /s in 1960, indicating that the increase in the runoff discharge intensified the erosion in the tidal creeks. Six observation points in Figure 7 were used to quantitatively study the change process of bed thickness. It can be seen from T1 that the bed thickness decreased continuously, indicating that the tidal channel was in a state of continuous erosion. By 1970, the erosion thickness of T1 was 1.139 m, which was also the largest because the average velocity of T1 was 0.297 m/s, the largest among the six points. T2 was located near the entrance to the lagoon, whose erosion thickness increased with the decreasing rate from 0.039 m/y to 0.022 m/y and reached 0.608 m in 1970. T3 was located at the tidal flat on the northern tidal creek A, whose deposition thickness increased slowly with the reducing rate and the average deposition rate was 0.006 m/y. T4 was located in tidal creek D, and the bed thickness change fluctuated. The erosion thickness reached the maximum of 0.375 m in 1961. After 1968, the erosion thickness was basically maintained at 0.248 m. T5 and T6 represent the changes of tidal creek C and B, respectively, whose terrains were both under continuous erosion. After 1967, the change rate in the erosion thickness was less than 0.010 m/y at T5, and the erosion thickness was 0.652 m in 1970. The average annual erosion thickness at T6 in tidal creek B was 0.030 m/y, and the maximum change rate was 0.072 m/y in 1961. It can also be seen from T4 and T5 that the change rate in 1961 had a sudden increase with 0.171 m/y and 0.046 m/y, respectively, because the runoff discharge in 1960 was the maximum value compared with that in any other years. The runoff discharge of the Luanxian Station was the maximum of 300 m 3 /s in 1960, indicating that the increase in the runoff discharge intensified the erosion in the tidal creeks.  Six observation points in Figure 7 were used to quantitatively study the change process of bed thickness. It can be seen from T1 that the bed thickness decreased continuously, indicating that the tidal channel was in a state of continuous erosion. By 1970, the erosion thickness of T1 was 1.139 m, which was also the largest because the average velocity of T1 was 0.297 m/s, the largest among the six points. T2 was located near the entrance to the lagoon, whose erosion thickness increased with the decreasing rate from 0.039 m/y to 0.022 m/y and reached 0.608 m in 1970. T3 was located at the tidal flat on the northern tidal creek A, whose deposition thickness increased slowly with the reducing rate and the average deposition rate was 0.006 m/y. T4 was located in tidal creek D, and the bed thickness change fluctuated. The erosion thickness reached the maximum of 0.375 m in 1961. After 1968, the erosion thickness was basically maintained at 0.248 m. T5 and T6 represent the changes of tidal creek C and B, respectively, whose terrains were both under continuous erosion. After 1967, the change rate in the erosion thickness was less than 0.010 m/y at T5, and the erosion thickness was 0.652 m in 1970. The average annual erosion thickness at T6 in tidal creek B was 0.030 m/y, and the maximum change rate was 0.072 m/y in 1961. It can also be seen from T4 and T5 that the change rate in 1961 had a sudden increase with 0.171 m/y and 0.046 m/y, respectively, because the runoff discharge in 1960 was the maximum value compared with that in any other years. The runoff discharge of the Luanxian Station was the maximum of 300 m 3 /s in 1960, indicating that the increase in the runoff discharge intensified the erosion in the tidal creeks. Due to changes in sea level and runoff, the lagoon could not reach a complete equilibrium state and the average bed thickness kept rising. As can be seen from Figure 8, both Due to changes in sea level and runoff, the lagoon could not reach a complete equilibrium state and the average bed thickness kept rising. As can be seen from Figure 8, both erosion and deposition increased with time continuously with a similar magnitude. The average rate of change of erosion was 2.07 kg/m 2 /y, and that of deposition was 2.29 kg/m 2 /y, so the net deposition remained positive, which indicates that the lagoon maintained at a deposition state likely due to the ample sediment supply from the open sea. The change rate of net deposition fluctuated around an average increasing rate of 0.22 kg/m 2 /y. Due to changes in sea level and runoff, the lagoon could not reach a complete equilibrium state and the average bed thickness kept rising. As can be seen from Figure 8, both erosion and deposition increased with time continuously with a similar magnitude. The average rate of change of erosion was 2.07 kg/m 2 /y, and that of deposition was 2.29 kg/m 2 /y, so the net deposition remained positive, which indicates that the lagoon maintained at a deposition state likely due to the ample sediment supply from the open sea. The change rate of net deposition fluctuated around an average increasing rate of 0.22 kg/m 2 /y. In natural conditions, with changes in the influences of sea level rise and runoff discharge, the lagoon stayed in a slight deposition state, while the tidal creeks were eroded slowly and the increasing river discharge intensified the erosion in the tidal creeks.

1970~1979
During 1970~1979, a tide gate was built at the entrance to the lagoon. The terrain of the QL in 1980 (1 January 1980) and changes in the bed thickness from 1970 to 1979 are shown in Figure 9. Overall, there was little change in the geomorphological characteristics of the tidal network from 1970 to 1979. The bed thickness changed slightly with an average In natural conditions, with changes in the influences of sea level rise and runoff discharge, the lagoon stayed in a slight deposition state, while the tidal creeks were eroded slowly and the increasing river discharge intensified the erosion in the tidal creeks.

1970~1979
During 1970~1979, a tide gate was built at the entrance to the lagoon. The terrain of the QL in 1980 (1 January 1980) and changes in the bed thickness from 1970 to 1979 are shown in Figure 9. Overall, there was little change in the geomorphological characteristics of the tidal network from 1970 to 1979. The bed thickness changed slightly with an average bed thickness variation of 0.026 m. Deposition occurred in 74% of the water area with an average deposition thickness of 0.031 m, while erosion occurred in 7% of the water area with an average erosion thickness of 0.020 m. The average current speed decreased to 0.01 m/s from 0.06 m/s of 1969, so that the bed shear stress was weakened significantly, which led to a litter change in the geomorphology. The most obvious change in the bed thickness occurred in the tidal inlet, where obvious deposition occurred with the maximum bed change thickness of 0.892 m. Due to the existence of the tide gate, the dynamic condition was weakened during ebb tide, causing the offshore sediment to deposit in the tidal inlet. There was a hotspot of erosion at the seaward side of tide gate, which may have been the result of a vortex induced by the large hydraulic head drop when the water level at tidal inlet was lower than 0.5 m. The maximum erosion thickness reached 1.080 m at that location.
There was little change in the bed elevation in the lagoon, so that only erosion and deposition are considered in this section. As shown in Figure 10, the erosion and deposition were significantly less than those in 1950~1969. The erosion, deposition and the net deposition all grew with time and the erosion was much less than the deposition so that the deposition was approximately same as the net deposition, as shown in Figure 10. Due to the construction of the tide gate, the lagoon dynamic was weakened considerably, along with the bottom shear stress and erosion ability. The net deposition amount increased continuously, indicating that the lagoon was in a continuous deposition state due to the sustained supply of offshore sediment. Generally speaking, the construction of the tide gate helped the lagoon to reach an equilibrium state, while greatly influencing the erosion and deposition downstream of the tide gate.
with an average erosion thickness of 0.020 m. The average current speed decreased to 0.01 m/s from 0.06m/s of 1969, so that the bed shear stress was weakened significantly, which led to a litter change in the geomorphology. The most obvious change in the bed thickness occurred in the tidal inlet, where obvious deposition occurred with the maximum bed change thickness of 0.892 m. Due to the existence of the tide gate, the dynamic condition was weakened during ebb tide, causing the offshore sediment to deposit in the tidal inlet. There was a hotspot of erosion at the seaward side of tide gate, which may have been the result of a vortex induced by the large hydraulic head drop when the water level at tidal inlet was lower than 0.5 m. The maximum erosion thickness reached 1.080 m at that location. There was little change in the bed elevation in the lagoon, so that only erosion and deposition are considered in this section. As shown in Figure 10, the erosion and deposition were significantly less than those in 1950~1969. The erosion, deposition and the net deposition all grew with time and the erosion was much less than the deposition so that the deposition was approximately same as the net deposition, as shown in Figure 10. Due to the construction of the tide gate, the lagoon dynamic was weakened considerably, along with the bottom shear stress and erosion ability. The net deposition amount increased continuously, indicating that the lagoon was in a continuous deposition state due to the sustained supply of offshore sediment. Generally speaking, the construction of the tide gate helped the lagoon to reach an equilibrium state, while greatly influencing the erosion and deposition downstream of the tide gate.

1980~1989
During 1980~1989, the 1.5 km 2 tidal flat in the northwest of QL was reclaimed, and the tide gate became open with its crest elevation set at −1.5 m. The opening tide gate connected the lagoon with offshore waters, which enhanced the dynamic condition significantly. The average velocity of QL increased from 0.01 m/s in 1970-1979 to 0.07 m/s. The width of the entrance to the lagoon decreased from 300 m in 1970 to 85 m in 1980, resulting in a significant increase in the velocity and stronger erosion at the tide gate. As can be seen from Figure 11a, the bed elevation at tide gate reached −5.908 m. Compared with the terrain on 1 January 1980, the change in creek A~C was not obvious, and creek D migrated slightly to the west. From Figure 11b, the tidal channel and tidal creeks were mainly eroded with the maximum erosion thickness of 3.582 m. The main deposition areas were at the tidal inlet offshore and on both sides of the tidal creek and the maximum deposition thickness was 2.099 m.

1980~1989
During 1980~1989, the 1.5 km 2 tidal flat in the northwest of QL was reclaimed, and the tide gate became open with its crest elevation set at −1.5 m. The opening tide gate connected the lagoon with offshore waters, which enhanced the dynamic condition significantly. The average velocity of QL increased from 0.01 m/s in 1970-1979 to 0.07 m/s. The width of the entrance to the lagoon decreased from 300 m in 1970 to 85 m in 1980, resulting in a significant increase in the velocity and stronger erosion at the tide gate. As can be seen from Figure 11a, the bed elevation at tide gate reached −5.908 m. Compared with the terrain on 1 January 1980, the change in creek A~C was not obvious, and creek D migrated slightly to the west. From Figure 11b, the tidal channel and tidal creeks were mainly eroded with the maximum erosion thickness of 3.582 m. The main deposition areas were at the tidal inlet offshore and on both sides of the tidal creek and the maximum deposition thickness was 2.099 m.
Due to the tide gate construction, the observation station T2 in Figure 7 was moved to be at the entrance of the lagoon, and the observation stations T1~T6 as shown in Figure 12 were selected to study the bed thickness change over time. The erosion thickness was 1.090 m after 10 years at T1 with a decreasing erosion rate from 0.241 m/y in the first year to below 0.050 m/y after 1986. T2 was located at the entrance to lagoon and had an erosion thickness of 3.788 m after 10 years. Nonetheless, in 1987, the deposition became dominant with the rate of 0.031 m/y in the following years. As indicated by the observation at station T3, the geomorphic evolution process of the tidal flat was in a continuous deposition state and the change rate of the thickness decreased to less than 0.010 m/y after 1988. T4, T5 and T6 represent the geomorphic changes of tidal creek D, C and B, respectively. They were all in a continuous erosion state with a decreasing rate of change with time, which was less than 0.010 m/y after 1987. In summary, the construction project resulted in considerable geomorphological evolution during the initial period and the change rate decreased gradually with time. It had the most significant effect on the tidal channel, especially that at the entrance to lagoon but had a relatively weak effect inside the lagoon.

1980~1989
During 1980~1989, the 1.5 km 2 tidal flat in the northwest of QL was reclaimed, and the tide gate became open with its crest elevation set at −1.5 m. The opening tide gate connected the lagoon with offshore waters, which enhanced the dynamic condition significantly. The average velocity of QL increased from 0.01 m/s in 1970-1979 to 0.07 m/s. The width of the entrance to the lagoon decreased from 300 m in 1970 to 85 m in 1980, resulting in a significant increase in the velocity and stronger erosion at the tide gate. As can be seen from Figure 11a, the bed elevation at tide gate reached −5.908 m. Compared with the terrain on 1 January 1980, the change in creek A~C was not obvious, and creek D migrated slightly to the west. From Figure 11b, the tidal channel and tidal creeks were mainly eroded with the maximum erosion thickness of 3.582 m. The main deposition areas were at the tidal inlet offshore and on both sides of the tidal creek and the maximum deposition thickness was 2.099 m. Due to the tide gate construction, the observation station T2 in Figure 7 was moved to be at the entrance of the lagoon, and the observation stations T1~T6 as shown in Figure  12 were selected to study the bed thickness change over time. The erosion thickness was 1.090 m after 10 years at T1 with a decreasing erosion rate from 0.241 m/y in the first year to below 0.050 m/y after 1986. T2 was located at the entrance to lagoon and had an erosion thickness of 3.788 m after 10 years. Nonetheless, in 1987, the deposition became dominant with the rate of 0.031 m/y in the following years. As indicated by the observation at station T3, the geomorphic evolution process of the tidal flat was in a continuous deposition state and the change rate of the thickness decreased to less than 0.010 m/y after 1988. T4, T5 and T6 represent the geomorphic changes of tidal creek D, C and B, respectively. They were all in a continuous erosion state with a decreasing rate of change with time, which was less than 0.010 m/y after 1987. In summary, the construction project resulted in considerable geomorphological evolution during the initial period and the change rate decreased gradually with time. It had the most significant effect on the tidal channel, especially that at the entrance to lagoon but had a relatively weak effect inside the lagoon. The changes in the amounts of erosion, deposition and net deposition from 1980 to 1989 are shown in Figure 13. The deposition amount increased with fluctuation, and the average annual deposition amount was 6.74 kg/m 2 /y. The amount of erosion fluctuated  1980 1981 1982 1983 1984 1985 1986 1987 1988 1989 1990 Bed thickness change (m) Time (year)   T1  T2  T3  T4  T5  The changes in the amounts of erosion, deposition and net deposition from 1980 to 1989 are shown in Figure 13. The deposition amount increased with fluctuation, and the average annual deposition amount was 6.74 kg/m 2 /y. The amount of erosion fluctuated and increased before 1986, reaching the maximum value of 26.92 kg/m 2 in 1986, and then fluctuated in the range of 25.00~27.00 kg/m 2 . The dynamic condition induced by the opening tide gate intensified the erosion in the lagoon, which gradually decreased over time. The net deposition increased almost linearly with time, and the average annual deposition amount reached 4.43 kg/m 2 /y during 1980~1989. In summary, the shrinkage of the lagoon area caused by reclamation and the narrowed width of the entrance to lagoon resulted in an increased deposition of the lagoon system, and the rate was higher than that in the previous years, as shown in Figure 10.

1990~2007
During 1990~2007, QL was further reclaimed, the water area shrunk from 5.2 km 2 to 2.3 km 2 , and the tidal channel shape changed from curved to straight with an average width of 110 m. The average bed elevation rose by 0.240 m by 1 January 2008. Creek B and C inside the lagoon remained similar with their initial state in 1990, while creek D became inconspicuous after significant deposition accumulation (Figure 14a). According to the bed thickness change (Figure 14b), the major erosion area was near the tidal inlet where the maximum erosion thickness was 2.103 m. This was because of a small cross section width near the tidal inlet, and the minimum width was only 60 m, which led to a fast flow velocity of 0.22 m/s. The most deposition was concentrated at the mouth of the lagoon where the maximum deposition thickness reached 3.878 m. The average bed thickness change in the lagoon was 0.240 m, indicating that the sediment supplement from open sea was more favorable to the deposition of QL after the reclamation and the tidal channel straightening.

1990~2007
During 1990~2007, QL was further reclaimed, the water area shrunk from 5.2 km 2 to 2.3 km 2 , and the tidal channel shape changed from curved to straight with an average width of 110 m. The average bed elevation rose by 0.240 m by 1 January 2008. Creek B and C inside the lagoon remained similar with their initial state in 1990, while creek D became inconspicuous after significant deposition accumulation (Figure 14a). According to the bed thickness change (Figure 14b), the major erosion area was near the tidal inlet where the maximum erosion thickness was 2.103 m. This was because of a small cross section width near the tidal inlet, and the minimum width was only 60 m, which led to a fast flow velocity of 0.22 m/s. The most deposition was concentrated at the mouth of the lagoon where the maximum deposition thickness reached 3.878 m. The average bed thickness change in the lagoon was 0.240 m, indicating that the sediment supplement from open sea was more favorable to the deposition of QL after the reclamation and the tidal channel straightening.
The observation stations in Figure 12a were selected to examine the evolution of bed thickness from 1990 to 2008, as shown in Figure 15. The tidal flat (T3) remained in a continual deposition state with a change rate less than 0.010 m/y after 2004. The tidal creek (T4, T5, T6) and tidal channel (T1) were firstly in a deposition state and in a slight erosion state later on. In 2001, the deposition thickness at the entrance of the lagoon (T2), reached the maximum of 2.040 m, and then decreased continuously. These results suggest that the reclamation and strengthened tidal channel gave rise to extensive deposition in the lagoon system.
The amounts of erosion, deposition and net deposition from 1990 to 2007 are shown in Figure 16. The deposition was much greater than the erosion. The amounts of deposition, erosion and net deposition were123.80 kg/m 2 , 30.34 kg/m 2 , and 93.46 kg/m 2 , respectively in 2007, indicating that QL experienced serious deposition during this period, and the aver-age annual net deposition was 5.19 kg/m 2 /y. The change rates tended to rise linearly with time after 1995 with constant rates of 0.88 kg/m 2 /y (erosion), 4.59 kg/m 2 /y (deposition) and 3.71 kg/m 2 /y (net deposition), respectively. C inside the lagoon remained similar with their initial state in 1990, while creek D became inconspicuous after significant deposition accumulation (Figure 14a). According to the bed thickness change (Figure 14b), the major erosion area was near the tidal inlet where the maximum erosion thickness was 2.103 m. This was because of a small cross section width near the tidal inlet, and the minimum width was only 60 m, which led to a fast flow velocity of 0.22 m/s. The most deposition was concentrated at the mouth of the lagoon where the maximum deposition thickness reached 3.878 m. The average bed thickness change in the lagoon was 0.240 m, indicating that the sediment supplement from open sea was more favorable to the deposition of QL after the reclamation and the tidal channel straightening. The observation stations in Figure 12a were selected to examine the evolution of bed thickness from 1990 to 2008, as shown in Figure 15. The tidal flat (T3) remained in a continual deposition state with a change rate less than 0.010 m/y after 2004. The tidal creek (T4, T5, T6) and tidal channel (T1) were firstly in a deposition state and in a slight erosion state later on. In 2001, the deposition thickness at the entrance of the lagoon (T2), reached the maximum of 2.040 m, and then decreased continuously. These results suggest that the reclamation and strengthened tidal channel gave rise to extensive deposition in the lagoon system. The amounts of erosion, deposition and net deposition from 1990 to 2007 are shown in Figure 16. The deposition was much greater than the erosion. The amounts of deposition, erosion and net deposition were123.80 kg/m 2 , 30.34 kg/m 2 , and 93.46 kg/m 2 , respectively in 2007, indicating that QL experienced serious deposition during this period, and the average annual net deposition was 5.19 kg/m 2 /y. The change rates tended to rise linearly with time after 1995 with constant rates of 0.88 kg/m 2 /y (erosion), 4.59 kg/m 2 /y (deposition) and 3.71 kg/m 2 /y (net deposition), respectively. The amounts of erosion, deposition and net deposition from 1990 to 2007 are sh in Figure 16. The deposition was much greater than the erosion. The amounts of dep tion, erosion and net deposition were123.80 kg/m 2 , 30.34 kg/m 2 , and 93.46 kg/m 2 , res tively in 2007, indicating that QL experienced serious deposition during this period, the average annual net deposition was 5.19 kg/m 2 /y. The change rates tended to rise arly with time after 1995 with constant rates of 0.88 kg/m 2 /y (erosion), 4.59 kg/m 2 /y ( osition) and 3.71 kg/m 2 /y (net deposition), respectively.

2008~2018
During 2008~2018, the reclamation of QL attained a saturation situation, and the terfront of the lagoon changed little compared to 1990~2007. A bridge was built at

2008~2018
During 2008~2018, the reclamation of QL attained a saturation situation, and the waterfront of the lagoon changed little compared to 1990~2007. A bridge was built at the entrance to lagoon to connect the two sides of the tidal channel. The top part of the tide gate was removed and the tidal channel was widened near the tidal inlet. The average width of the tidal channel increased from 110 m to 150 m. As is shown in Figure 17, the bathymetry of QL did not change significantly, but the overall bed elevation rose by 0.106 m. The erosion occurred mainly in the tidal creeks and the northern side of tidal channel with the maximum erosion thickness of 0.692 m, and various degrees of deposition occurred in other areas. The deposition along the southern side of tidal channel was larger with the maximum thickness of 3.691 m. Due to the presence of the bridge piers, the cross section was narrowed, resulting in erosion of the shoreline at the northeast side of the bridge piers with the maximum erosion thickness of 3.513 m. The eroded sediment was deposited at the southwest side of bridge piers and the deposition at the lagoon-side of the bridge piers was larger.  The field stations T1~T6 shown in Figure 12a were selected to analyze the terrain changes ( Figure 18). T1 was in a deposition state with a maximum thickness of 0.014 m in 2010 and then remained in a continuous erosion state. T2 was located near the piers and was in a state of sediment accumulation with a decreasing rate. After 2012, the change tended to be linear with an average rate of 0.023 m/y. T3 was located at tidal flat which stayed at a deposition state before 2015, then the deposition thickness remained at a constant of 0.086 m. T4 was located in tidal creek D, which is close to T3, so the change trend was similar to T3, and the deposition thickness remained at a constant of 0.036 m. T5 (tidal creek C) and T6 (tidal creek B) were in an erosion state with the average annual erosion rates of 0.055 m/y and 0.019 m/y, respectively. The highest erosion rates of T5 and T6 were 0.348 m/y and 0.122 m/y in the first year, respectively. Slight deposition occurred at T6 from 2013 to 2016, and the average deposition rate was 0.004 m/y. In general, the construction projects of the bridge construction and tidal channel widening made the tidal creeks enter an equilibrium state rapidly and had a little effect on them, while the tidal flat had a linear deposition tendency. The field stations T1~T6 shown in Figure 12a were selected to analyze the terrain changes ( Figure 18). T1 was in a deposition state with a maximum thickness of 0.014 m in 2010 and then remained in a continuous erosion state. T2 was located near the piers and was in a state of sediment accumulation with a decreasing rate. After 2012, the change tended to be linear with an average rate of 0.023 m/y. T3 was located at tidal flat which stayed at a deposition state before 2015, then the deposition thickness remained at a constant of 0.086 m. T4 was located in tidal creek D, which is close to T3, so the change trend was similar to T3, and the deposition thickness remained at a constant of 0.036 m. T5 (tidal creek C) and T6 (tidal creek B) were in an erosion state with the average annual erosion rates of 0.055 m/y and 0.019 m/y, respectively. The highest erosion rates of T5 and T6 were 0.348 m/y and 0.122 m/y in the first year, respectively. Slight deposition occurred at T6 from 2013 to 2016, and the average deposition rate was 0.004 m/y. In general, the construction projects of the bridge construction and tidal channel widening made the tidal creeks enter an equilibrium state rapidly and had a little effect on them, while the tidal flat had a linear deposition tendency.
As is shown in Figure 19, the erosion, deposition and the net deposition all grew with time and the erosion was much less than the deposition so that the deposition was almost same as the net deposition., and the average annual erosion and deposition rates were 2.57 kg/m 2 /y and 9.30 kg/m 2 /y, respectively, which were higher than the rates from 1990 to 2007 (1.68 kg/m 2 /y and 6.87 kg/m 2 /y). This indicates that the widening of the tidal channel increased the sediment supplement from the open sea. The annual average amount of erosion rate decreased to less than 1.00 kg/m 2 /y after 2013, and that of deposition reached the maximum of 15.78 kg/m 2 /y in 2008, then fluctuated in the range of 6.50~10.50 kg/m 2 /y. From 2008 to 2018, the annual average net deposition rate was 6.73 kg/m 2 /y. The change in the net deposition amount tended to be linear after 2011. From 2011 to 2018, the annual net deposition rate was 7.42 kg/m 2 /y, indicating that an accumulation of sediments stayed within the lagoon system. As is shown in Figure 19, the erosion, deposition and the net deposition all grew with time and the erosion was much less than the deposition so that the deposition was almost same as the net deposition., and the average annual erosion and deposition rates were 2.57 kg/m 2 /y and 9.30 kg/m 2 /y, respectively, which were higher than the rates from 1990 to 2007 (1.68 kg/m 2 /y and 6.87 kg/m 2 /y). This indicates that the widening of the tidal channel increased the sediment supplement from the open sea. The annual average amount of erosion rate decreased to less than 1.00 kg/m 2 /y after 2013, and that of deposition reached the maximum of 15.78 kg/m 2 /y in 2008, then fluctuated in the range of 6.50~10.50 kg/m 2 /y. From 2008 to 2018, the annual average net deposition rate was 6.73 kg/m 2 /y. The change in the net deposition amount tended to be linear after 2011. From 2011 to 2018, the annual net deposition rate was 7.42 kg/m 2 /y, indicating that an accumulation of sediments stayed within the lagoon system.   As is shown in Figure 19, the erosion, deposition and the net deposition all grew with time and the erosion was much less than the deposition so that the deposition was almost same as the net deposition., and the average annual erosion and deposition rates were 2.57 kg/m 2 /y and 9.30 kg/m 2 /y, respectively, which were higher than the rates from 1990 to 2007 (1.68 kg/m 2 /y and 6.87 kg/m 2 /y). This indicates that the widening of the tidal channel increased the sediment supplement from the open sea. The annual average amount of erosion rate decreased to less than 1.00 kg/m 2 /y after 2013, and that of deposition reached the maximum of 15.78 kg/m 2 /y in 2008, then fluctuated in the range of 6.50~10.50 kg/m 2 /y. From 2008 to 2018, the annual average net deposition rate was 6.73 kg/m 2 /y. The change in the net deposition amount tended to be linear after 2011. From 2011 to 2018, the annual net deposition rate was 7.42 kg/m 2 /y, indicating that an accumulation of sediments stayed within the lagoon system.

Discussions
In tide-dominated coastal systems, tidal asymmetry is one of the primary controlling factors of sediment transport and geomorphological evolution. Tidal asymmetry controls net sediment transport, which fundamentally determines the long-term evolution of geomorphological systems [60]. Tidal asymmetry is affected by human activities together with hydrodynamic conditions, and it was found that human activities have more noteworthy influences on tidal asymmetry in a semi-closed sea than that in an open sea [32,33]. Therefore, tidal asymmetry was critical for our study of geomorphological evolution. It is caused by the tidal deformation due to the nonlinear interaction of tidal components, so the analysis of the interaction through the phase difference and amplitude ratio among two or more tidal components can be used to quantify the tidal asymmetry [61]. M 2 is the main tidal component in the QHD Sea; therefore, the tidal asymmetry is determined by M 2 and the shallow water tidal harmonics M 4 . In this section, the following formula is used to calculate tidal asymmetry γ [62].
where the amplitudes of M 2 and M 4 are a M2 and a M4 , and the phases are φ M2 and φ M4 , respectively. The variation range of γ is from −1 to 1. Positive γ indicates the flood dominance, namely fast flood tide and slow ebb tide, while negative γ indicates the ebb dominance, namely fast ebb tide and slow flood tide.

Spatial Distribution of Tide Asymmetry
The tidal asymmetry γ at each grid mesh was calculated, and the distribution of tidal asymmetry of QL during the five phase time periods was obtained by Kriging interpolation [63], as shown in Figure 20.
It can be seen that the tidal asymmetry distribution during 1950~1969 followed the bathymetry characteristics (Figure 20a). The negative γ in the tidal creeks within the range of −0.35~−0.25, indicates the ebb dominance and erosion. The positive γ at the shallow area on both sides of tidal creeks in the range of 0.25~0.45, indicates flood dominance and deposition. In the tidal channel, there is an area of flood dominance with γ greater than 0.45 in the southwest near the tidal inlet, and obvious deposition occurs in this area. Most parts of the tidal channel are dominated by an ebb tide with a negative γ within the range from −0.25 to −0.20. The tidal asymmetry distribution during this period is similar to the distribution of bed thickness change (Figure 6b).
In 1970~1979, a tide gate was built at the entrance to the lagoon. The width of the entrance decreased from 300 m to 85 m, which significantly altered the flow condition in the lagoon. It can be seen that γ fluctuated between 0.70 and 0.80 indicating a flood dominance inside the lagoon (Figure 20b). Ebb dominance occurred in the tidal channel with γ ranging from −0.30 to −0.20. However, the ebb dominance was weaker near the tidal inlet than the side close to the tide gate. Despite the flood dominance and sediment accumulation inside the lagoon, the deposition was not obvious. This may have been because the sediment input was not considered at the boundary of upstream rivers in the simulation. An ebb dominance occurred in the tidal channel, but it can be seen from Figure 9b that the tidal channel was in a state of deposition, and the deposition was larger near the tidal inlet with the maximum deposition thickness of 0.892 m. Therefore, bed erosion or deposition is not only dependent on tidal asymmetry, but also on the external sediment supply.
The reclamation area of the lagoon increased by 1.5 km 2 and the tide gate stayed open during 1980~1989. The distribution of tidal asymmetry (Figure 20c) was similar to that during 1950~1969, and ebb dominance presented in the tidal creeks with γ of −0.35~−0.20. However, inside the lagoon, the closer to the tide gate, the stronger the ebb dominance was, which is opposite with the period of 1950~1969. This was mainly because the average discharge into lagoon decreased from 6.60 m 3 /s during 1950~1969, to 3.12 m 3 /s during 1980~1989. The decreasing discharge reduced the runoff influence on ebb dominance. At the same time, the velocity near the tide gate increased significantly, which resulted in a stronger erosion effect. Therefore, the position with a greater ebb dominance during 1980~1989 was located near the tide gate. The tidal flat on both sides of the tidal creek experienced flood dominance with γ ranging from 0.35 to 0.45, which is greater than that during 1950~1969, but the area of flood dominance decreased. The tidal asymmetry γ of tidal channel ranged from −0.20 to −0.10 and the closer the channel was to the open sea, the weaker the ebb dominance is, but it was weaker than that during 1950~1969. The reclamation area of the lagoon increased by 1.5 km 2 and the tide gate stayed open during 1980~1989. The distribution of tidal asymmetry (Figure 20c) was similar to that during 1950~1969, and ebb dominance presented in the tidal creeks with γ of −0.35~−0.20. However, inside the lagoon, the closer to the tide gate, the stronger the ebb dominance was, which is opposite with the period of 1950~1969. This was mainly because the average discharge into lagoon decreased from 6.60 m 3 /s during 1950~1969, to 3.12 m 3 /s during 1980~1989. The decreasing discharge reduced the runoff influence on ebb dominance. At the same time, the velocity near the tide gate increased significantly, which resulted in a stronger erosion effect. Therefore, the position with a greater ebb dominance during 1980~1989 was located near the tide gate. The tidal flat on both sides of the tidal creek During 1990~2007, the reclamation area in the lagoon increased by 4.4 km 2 , and the tidal channel was straightened. Compared with that during 1980~1989, the flood dominance of the tidal flat increased and the tidal asymmetry γ ranged from 0.40 to 0.60; the ebb dominance of the tidal creek increased and the tidal asymmetry γ ranged from −0.50 to −0.45 (Figure 20d). In the tidal channel, the tidal asymmetry γ ranged from −0.45 to 0.18, and ebb dominance near the entrance to lagoon was the largest, which decreased gradually and shifted to flood dominance near the tidal inlet. Both flood and ebb dominance were strong in QL, but the area of flood dominance was larger, resulting in a significant increase in the net deposition (Figure 14b). During 2008~2018, the tidal asymmetry γ of the tidal flat ranged from 0.06 to 0.63, and the area with γ greater than 0.50 (Figure 20e) reduced significantly compared with that during 1990~2007, indicating that the flood dominance of the tidal flat was weakened. Ebb dominance of the tidal creek, with γ ranging from −0.32 to −0.22, weakened compared with that during 1990~2007. The variation of γ near tide gate was discontinuous, and the ebb dominance at lagoon side (average value γ = −0.28) was stronger than that at channel side (average value γ = −0.15). In the main axis of tidal channel, the ebb tide was dominant, and the tidal asymmetry γ was from −0.15 to −0.09. The ebb dominance became weaker towards the open sea. In the southwest bank of tidal channel, the flood tide was dominant, and the tidal asymmetry γ ranged from 0.05 to 0.60.

Effects of Sea Level and Runoff
The variation of γ in QL from 1950 to 2018 is shown in Figure 21. Under the natural evolution process, the tidal asymmetry γ ranged from −0.37 to −0.02 (ebb dominance) during 1950~1969, indicating that the lagoon was prone to erosion. Due to the ample supply of offshore sediment, the net deposition increased continuously (Figure 8), but the average annual net deposition was small at only 0.22 kg/m 2 /y. During 1970~1979, the tidal asymmetry γ increased to the range of 0.50~0.78, that is to say, ebb dominance changed to flood dominance, which gave arise to sediment accumulation in the lagoon, but the average annual net deposition increased to only 0.56kg/m 2 /y due to an insufficient sediment supply. During 1980~1989, the average annual tide asymmetry γ changed from 0.61 (flood dominance) in 1970 to −0.11 (ebb dominance), the average annual net deposition reached 4.43 kg/m 2 /y ( Figure 13). Although the flood dominance during 1970~1979 was obvious, the average annual net deposition during 1980~1989 was larger than that during 1970~1979 because the capacity of carrying sediment into the lagoon was weakened by the operation of the tide gate during 1970~1979. During 1990~2007, the reclamation area increased by 2.90 km 2 and the tidal channel was straightened. In this period, the tidal asymmetry γ in the QL was positive. The ebb dominance changed to flood dominance, the tidal asymmetry γ increased to 0.01, and the average annual net deposition rate increased to 5.19 kg/m 2 /y. During 2008~2018, the tidal inlet offshore was widened and a bridge was built near the entrance to lagoon. The flood dominance was further enhanced, the average annual tidal asymmetry γ increased to 0.17 and the annual net deposition rate increased to 6.73 kg/m 2 /y. It can be seen from Figure 21, under the same construction project condition, the variation of γ was similar to that of sea level, but different from that of runoff. The pearson correlation coefficient [64] was used to evaluate the influence of sea level and runoff on tidal asymmetry, as shown in Table 3. The correlation between sea level and tidal asymmetry γ was strong, with correlation coefficients of 0.716~0.916, while the correlation between runoff and tidal asymmetry γ was weak, with correlation coefficients of 0.003~0.507, indicating that sea level has a more significant impact on tidal asymmetry than runoff. When sea level rises, the tidal asymmetry γ increases, which strengthens the flood dominance and weakens the ebb dominance. In other words, sea level rise intensifies sediment accumulation in the QL.
From 1950 to 1969, the hydrodynamics in QL were hardly affected by construction projects, so the tidal asymmetry was mainly affected by runoff and sea level. Moreover, the correlation between runoff and tidal asymmetry γ was weak, so the tidal asymmetry determined by sea level was selected as the representation without the influence of construction projects. The relationship between sea level and tidal asymmetric γ from 1950 to 1969 is shown in Figure 22. By the linear fitting of sea level (SL) and tidal asymmetry γ from 1950 to 1969, the following relationship was obtained: The determination coefficient R 2 of the linear relationship reaches 0.86 with the applicable fitting effect. Therefore, the equation can be used to estimate tidal asymmetry without construction projects according to sea level, then, the geomorphological evolution can be predicted preliminarily. It can be seen from Figure 21, under the same construction project condition, the variation of γ was similar to that of sea level, but different from that of runoff. The pearson correlation coefficient [64] was used to evaluate the influence of sea level and runoff on tidal asymmetry, as shown in Table 3. The correlation between sea level and tidal asymmetry γ was strong, with correlation coefficients of 0.716~0.916, while the correlation between runoff and tidal asymmetry γ was weak, with correlation coefficients of 0.003~0.507, indicating that sea level has a more significant impact on tidal asymmetry than runoff. When sea level rises, the tidal asymmetry γ increases, which strengthens the flood dominance and weakens the ebb dominance. In other words, sea level rise intensifies sediment accumulation in the QL. From 1950 to 1969, the hydrodynamics in QL were hardly affected by construction projects, so the tidal asymmetry was mainly affected by runoff and sea level. Moreover, the correlation between runoff and tidal asymmetry γ was weak, so the tidal asymmetry determined by sea level was selected as the representation without the influence of construction projects. The relationship between sea level and tidal asymmetric γ from 1950 to 1969 is shown in Figure 22. By the linear fitting of sea level (SL) and tidal asymmetry γ from 1950 to 1969, the following relationship was obtained: 1.49 0.06 SL γ = + (6) The determination coefficient R 2 of the linear relationship reaches 0.86 with the applicable fitting effect. Therefore, the equation can be used to estimate tidal asymmetry without construction projects according to sea level, then, the geomorphological evolution can be predicted preliminarily.

Effects of Construction Projects
In order to analyze the effect of construction projects on tidal asymmetry, the variations of tidal asymmetry with and without construction projects are shown in Figure 23. From 1950 to 1969, the tidal asymmetry with and without construction were almost identical.
During 1970~1979, the tidal asymmetry γ increased greatly by 0.77 to a positive value due to the tide gate, indicating that the construction made QL more susceptible to deposition, but then the tidal asymmetry decreased sharply by almost the same amount due to the opening of the tidal gate around 1979. This may have been caused by the dramatic decrease in the water exchange capacity inside lagoon with the tide gate construction and then the improved capacity with the opening of the gate. During 1980~1989, due to the reclamation and the opening of the tide gate, the average γ increased by 0.05 compared with the natural condition, which increased the flood dominance. During 1990~2007, the average value of γ increased by 0.08 due to the reclamation and tidal channel straightening, so the flood dominance increased further. This is consistent with the finding of Liang et al. [32,33] in the study of Laolonggou lagoon in Caofeidian that tidal flat reclamation increases the flood dominance. During 2008~2018, the average tidal asymmetry γ increased to 0.17 due to tidal channel widening and bridge construction, and it was 0.10 larger than that without construction project. In summary, the construction projects enhanced the flood dominance and therefore the deposition.

Effects of Construction Projects
In order to analyze the effect of construction projects on tidal asymmetry, the variations of tidal asymmetry with and without construction projects are shown in Figure 23. From 1950 to 1969, the tidal asymmetry with and without construction were almost identical. During 1970~1979, the tidal asymmetry γ increased greatly by 0.77 to a positive value due to the tide gate, indicating that the construction made QL more susceptible to deposition, but then the tidal asymmetry decreased sharply by almost the same amount due to the opening of the tidal gate around 1979. This may have been caused by the dramatic decrease in the water exchange capacity inside lagoon with the tide gate construction and then the improved capacity with the opening of the gate. During 1980~1989, due to the reclamation and the opening of the tide gate, the average γ increased by 0.05 compared with the natural condition, which increased the flood dominance. During 1990~2007, the average value of γ increased by 0.08 due to the reclamation and tidal channel straightening, so the flood dominance increased further. This is consistent with the finding of Liang et al. [32,33] in the study of Laolonggou lagoon in Caofeidian that tidal flat reclamation increases the flood dominance. During 2008~2018, the average tidal asymmetry γ increased to 0.17 due to tidal channel widening and bridge construction, and it was 0.10 larger than that without construction project. In summary, the construction projects enhanced the flood dominance and therefore the deposition.

Conclusions
The geomorphological evolution of QL from 1950 to 2018 was simulated to assess the impact of sea level rise, runoff and human activities and then the mechanism of erosion -0.5 0.0 0.5 1.0 1950 1955 1960 1965 1970 1975 1980 1985 1990 1995

Conclusions
The geomorphological evolution of QL from 1950 to 2018 was simulated to assess the impact of sea level rise, runoff and human activities and then the mechanism of erosion and deposition was examined in relation to tidal asymmetry, which provides a reference for the sediment transport mechanism and geomorphological evolution of other lagoons in the global region.
In natural conditions, under the influences of sea level and runoff discharge, the lagoon stayed in a slight deposition state, while the tidal creeks were eroded slowly and the increased river discharge could intensify the erosion in tidal creeks. The construction of tide gate drove the lagoon into an equilibrium state and had a great influence on the erosion and deposition downstream. Noteworthy geomorphological evolution occurred during the initial period after the tide gate construction and the rate of change decreased gradually with time. The shrinkage of the lagoon area by reclamation and the narrowed width of the entrance to lagoon increased the deposition in the lagoon system. Further reclamation and tidal channel strengthening resulted in an even more extensive deposition of the lagoon system. The bridge construction and tidal channel widening cause tidal creeks to attain an equilibrium state rapidly and had a little effect on tidal creeks afterwards, while the tidal flat has a deposition rate that increases with time almost linearly.
Spatial distribution of tidal asymmetry during natural period without human interventions was similar to that of the change in the bed thickness. The construction of the tide gate led to a sharp increase in the tidal asymmetry resulting from weak water exchange capacity inside lagoon, but it did not cause serious deposition due to little sediment supply. So, bed erosion or deposition is not only dependent on tidal asymmetry but is also related to the external sediment supply and the discharge of upstream rivers. Sea level rise can intensify sediments accumulation through tidal asymmetry, and the tidal asymmetry (γ) can be estimated from sea level (SL) before construction projects using the relationship of γ = 1.49SL + 0.06, then the geomorphological evolution can be estimated. There is little correlation between runoff and tidal asymmetry. Human activities such as construction projects enhance the flood dominance and thus intensify the deposition to certain degrees.