This section introduces information on remote sensing images and historical charts of JZB. Coastlines were extracted from remote sensing images or historical charts while water depth was derived from historical charts.

Historical charts of JZB, including charts from 1863, 1935, 1966, 1986, 1992, 1996, and 2005 (Table 1), were collected and then registered, as were the Landsat-5 Thematic Mapper (TM) images with resolution of 30 m in 1988, 1900, 1997, 2000, 2003, and 2005 and SPOT-5 satellite images in 2008 and 2020 with resolution of 2.5 m.

The coastline of the charts was modified according to the remote sensing map of nearby years if there were any. The extracted coastlines were divided into artificial coastlines (such as port, salt or aquaculture pond, and other artificial coastlines) and natural ones (estuary, island, sand or gravel beach, muddy tidal flats, and rocky). The length of each coastline type from each remote sensing image or historical chart was characterized and measured under the same principle.

The original depth values were derived from each chart including scattered water depth points and depth contours. Due to the large range of time, the projections, depth units, and chart datum of JZB charts are not the same. All the Mercator projections in the six charts collected were converted into Gauss-Kruger.

For the same sea area, the depth datum of various charts published in different historical periods is not always the same, because the determination of the depth datum is related to the position of the tide station, the length of the tide data, and the calculation method. The time interval and length of the tide data used in the chart of JZB in different periods are different, as are the positions of the depth datum (Table 1), thus, it is necessary to unify the depth datum. Here all the depth datum were corrected respective to the datum of the chart from 2005, by referencing the relationship of mean sea level and tide tables listed in charts, since the mean sea level is relatively stable. For example, according to the tide information table of 1996 chart, mean sea level is 2.4 m higher than the datum (Lowest Normal Low Water), while tide information table of 2005s chart indicated that mean sea level is 2.3 m higher than the chart datum (Lowest Normal Low Water). The mean sea level of JZB stayed the same from 1996 to 2005, thus, the Lowest Normal Low Water in 1996 is 0.1 m lower than that in 2005. Thus, the water depth of 1996s chart is uniformly added to by 0.1 m to be consistent with datum of 2005 chart. The standard deviation of depth correction in this way is less than ± 7 cm, which is higher than the current accuracy of depth sounding, thus the accuracy of water depth is guaranteed by correcting in this way. The above unified water depth points are grid-interpolated with the same grid distance (100 × 100 m) for each chart, thus the area and volume below or between certain contour values can be calculated accordingly based on the mesh water depth for each period.

MIKE 21, a very popular numerical model for simulating hydrodynamics, especially in bays and coastal zones, was applied here to simulate the tidal dynamics and pollution transport for each historical scenario. MIKE 21 is based on the numerical solution of Navier-Stokes equations, using unstructured triangular mesh and dry-wet boundary (Danish Hydraulic Institute, 2012), which can adapt well to the complicated coastlines of JZB.

The continuity equation is as Equation 1,

where h is the total water depth, t is time, x and y are the Cartesian coordinates, u¯ and v¯ are depth average value of velocity components in x and y direction, and S is the magnitude of the discharge due to point sources.

The two horizontal momentum equations for the x and y components are given in Equations 2 and 3, respectively:

where f is the Coriolis parameter, g is gravitational acceleration, η is surface elevation, ρ₀ is reference density of water, P_a is atmospheric pressure, ρ is the density of water, τ_sx and τ_sy are the x and y components of surface wind and τ_bx and τ_by are the components of bottom stress, u_s, v_s is the velocity by which the water is discharged into the ambient water, S_xx, S_xy, S_yx and S_yy are components of the wave radiation stress tensor. T_xy, T_x, and T_yy includes viscous friction, turbulent friction, and differential advection estimated using eddy viscosity formulation.

A Eulerian tracer-tracking module (Equation 4) was coupled to the Mike21 hydrodynamic model to simulate tracer transport for each scenario.

where C¯ is the depth-averaged concentration of the tracer, D_h is the horizontal diffusion coefficient, k_p is the linear decay rate of the tracer, and Cs is the concentration of tracer at the source.

The bay entrance of JZB is very narrow and the wave in JZB is generally very weak. Thus, the wave was not taken into account in this paper. As there is a large area of tidal flats in JZB, the flow of which sometimes dries out and sometimes is flooded, the dry-wet boundary is used in the model. When the water depth is less than the wetting depth at 0.1 m, the problem is reformulated, and only if the water depth is less than the drying depth of 0.005 m is the element removed from the calculation. The flooding depth, set as 0.05 m, is used to determine when an element is flooded (i.e., reentered into the computations). The reformulation is made by setting the momentum fluxes to zero and only taking the mass fluxes into consideration. In this way, continuity is fully preserved during the flooding and drying process as the water depths at the elements which are dried out are saved and then reused when the elements become flooded again. The bottom friction is set as Manning coefficient of 58 m^1/3/s. Eight major tidal constituents (M2, S2, K1, O1, N2, K2, P1, and Q1) are used to generate the tidal water levels at the open boundary located in the Yellow Sea.

The model in August 2012 was used for verification. Its coastline was derived from satellite in 2012 while topography was extracted mainly from 2005 nautical chart with some corrections where the change of coastline partially affected the water depth, because the depth of JZB has not changed obviously since 2005 (Li et al., 2014). The mesh of 2012 has 36,094 triangular elements with mesh size 60–100 m in the JZB (Figure 2A). The effect of 864 bridge square piers, with a side length of 6 m, is modeled as sub-grid structures using a simple drag-law to capture the increasing resistance imposed by the piers as the flow speed increases (Figure 2B).

The measured tidal level from the long-term tide station DG and current from station C1 and station C2, as shown in Figure 2A, were used for model validation. Figure 3 shows that the modeled tide levels match the observation data very well, with an average absolute error of 0.02 m. It is shown from Figure 4 that the model results of tidal velocity and direction agree with the observation data fairly well. Independence tests of the mesh size were then carried out to exclude the possibility that the results were affected by the intensification of mesh.

The coastline change in JZB mainly happened in the tidal flats. Thus, change of tidal flats area was inevitable, as was the change of tidal prism. Tidal prism is an important index for the environmental assessment of a bay and it reflects the self-purification capacity of a bay and determines the exchange intensity between the bay and the open sea, which is of great significance to the environment, ecology, erosion, and deposition of the bay (Shi et al., 2011).

Tidal prism P_m for any tidal period is calculated as the product of the mean of the high- and low-water surface areas of the bay behind the bay entrance and the tidal range in each segment, shown as Equation 5:

where H₁ is low-water level, H₂ is high-water level, and S(l) refers to the water area corresponding to the water level H and water boundary l. According to the definition, tidal prism can be calculated as the volume between the spring high-water and the spring low-water based on the detailed bathymetric data in section “Coastline and Water Depth Extraction.” The volume between the spring high-level and the spring low-water are taken as the tidal prism in this paper. Similarly, the area between the spring high-water and the spring low-water are taken as the area of tidal flats.

Water-exchange of a semi-closed bay determines the transport and redistribution of pollutants, nutrients, and suspended sediments, which is essential for the ecosystem of the bay. Residence time is one of the most frequently used timescales to reflect the water-exchange capacity of a bay.

Eulerian definition of residence time is applied in this paper, which has been widely used in numerous shallow coastal areas such as the Archachon Bay (Plus et al., 2009), Bay of Biscay (Grifoll et al., 2013), Jervis Bay (Sun et al., 2018), and Jiaozhou Bay (Yuan et al., 2021). A Eulerian passive tracer with a concentration of 1 was released at each grid cell inside the JZB on 1st February (time t₀). No further tracer was added after t₀, and the volume of each grid cell V_i in the bay remained constant over time. The concentration of tracer at each grid cell C_i evolved with time (t) as the bay exchanged water with the open sea. The relative tracer mass (RTM) is given as Equation 6:

where the subscript i labels the grid points within the bay. Here the C_i(t₀) = 1 as the initial conditions. The numerator and denominator of Equation 6 are the total tracer mass inside the bay at time t and time t₀, respectively.

The basin residence time is taken as the time at which the RTM inside the bay first fell below e^–1 (37%). Similarly, the local residence time at a grid cell was obtained as the time required for the tracer concentration at that grid cell to fell below 37% of the initial concentration.

Three runs to investigate the tidal current pattern, residence time, and pollutant transport, respectively, were then carried out for each year in 1863, 1935, 1966, 1986, 2005, and 2020.

As shown in Figure 1, the coastline barely changed from 1863 to 1935, only the northwest coastline of JZB moved seaward slightly due to the construction of salt ponds. From 1935 to 1966, the coastline in the northwest and northeast of JZB changed significantly, especially the isolated Hong Island in the northeast of JZB which merged into mainland due to the great expansion of salt and aquaculture ponds. From 1966 to 1986, the salt and aquaculture ponds in the northwest and northeast of JZB were further expanded toward the sea. During this period, another isolated island, Huang Island, in the southwest of JZB was incorporated into the mainland. From 1986 to 2020, the coastline in the southwest and east of JZB changed greatly, and the JZB Bridge across the northern bay, with a length of 31.63 km, was built in 2011. During this period, the coastline generally moved toward the sea and tended to be straight, due to port construction and urban expansion. Overall, the area of JZB shrank significantly by 42% (567.95 km²) from 1863 to 2020.

According to the coastline attribute indicated by charts and images, the length of each coastline type in JZB was measured.

Figure 5 shows that the length of nature coastline barely changed from 1863 to 1935, then the natural coastline decreased sharply from 1935, especially in the 50 years from 1935 to 1986, during which the natural coastline had been reduced by 137.9 km (70%). In 1986, the tidal flat coastline in a large area had basically disappeared, most of which had been replaced by salt and aquaculture ponds and port, and the rocky coastline had also decreased by 83%. After 1986, the decrease rate of natural coastline became slower.

The significant coastline change corresponded with rapid economic development and the coastline type change was mainly affected by human intervention. JZB had experienced salt and aquaculture pond construction in 1950s, reclamation of tidal flats from the mid-1960s to the 1970s, and reclamation for agriculture, port development, highway construction, and port adjacent projects since the 1980s. With the intensification of human development activities, the coastline moved toward the sea gradually while the natural coastline was lessening, and was gradually replaced by artificial coastlines.

As presented in section “Tidal Prism Calculation,” the water volume and area below spring high-water and the spring low-water in JZB can be calculated and the resulting tidal prism and area of tidal flats are listed in Table 2. According to Table 2 and Figure 6, the area of tidal flats decreased only by 6% (15.67 km²) in 72 years from 1863 to 1935 due to little human influence. Thereafter, the area of tidal flats decreased sharply by 76% (182.22 km²) in 85 years from 1935 to 2020 due to intense human activities. Large areas of natural tidal flats were replaced by salt/aquaculture ponds and artificial reclamation. The tidal prism decreased by 0.002 km³ during the 72 years from 1863 to 1935, and decreased by 0.290 km³ during the next 85 years from 1935 to 2020.

Generally, JZB was predominated by natural evolution from 1863 to 1935, during which the tidal prism barely changed. Thereafter, intensified human intervention resulted in more natural tidal flats being occupied by salt ponds, aquaculture area, and reclamation very quickly. Consequently, tidal prism, area of tidal flats, and area of JZB decreased sharply by 0.290 km³, 182 km², and 223 km², respectively, from 1935 to 2020, corresponding to shrinking rates of 123, 10, and 12 times, respectively, as that of before 1935.

Tidal dynamics have a significant effect on the water exchange capacity and pollutant transport in a bay. In this section, the evolution of mean tidal velocity and tidal residual currents of JZB are discussed.

The mean tidal velocity during 1 month (August) is calculated for each year and is shown in Figure 7. Generally, the average velocity in the area surrounded by Huangdao island, Tuandao, and the bay entrance was the largest, up to 0.80 m/s, and gradually decreased northward from the bay entrance. Along with the land reclamation, the mean tidal velocity in JZB decreased and the tidal dynamic weakened. The area with mean velocity greater than 0.2 m/s shrank considerably, especially in the northern part with only sporadic distribution in the deep channel in 2020.

Tidal residual currents, though much smaller than tidal currents, are essential for long-term mass-transport and evaluating the water exchange capacity (Chen et al., 2012; Li et al., 2014). The most calculated residual current in studies on tidal dynamics is the Euler residual current due to its simplicity. Eulerian residual current u_Euler is defined as the averaged velocity at a fixed location, over one or many tidal cycles (Leonardi et al., 2015), calculated as Equation 7 (Rodríguez et al., 2017).

where u is the tidal current at time t, T is the period of a tidal cycle, and n is the number of tidal cycles. In this section, the Euler residual current of JZB in each year is calculated on the basis of 1 month’s model results.

u_Euler is decomposed into U_E and V_E in x and y direction,

where t₀ is the initial time for calculation and u¯⁢(x0,t) and v¯⁢(x0,t) are depth-average value of velocity components in x and y direction at a fixed location x₀.

Figure 8 shows that residual currents are mostly less than 0.10 m/s in most scenarios, and larger (0.3 m/s) near the bay entrance. There are about 7–8 residual current eddies in JZB with blue character denoting clockwise eddy and red representing counterclockwise one. The strongest eddies of residual currents a, b and c (one anticyclonic between two cyclonic) are near the bay entrance and Huang Island in all scenarios.

Circulation pattern of residual current in JZB also varied due to the change of coastline and bathymetry. Residual currents eddy a was basically stable over the 150 years, but became a bit stronger in 2020, which is mainly related to the smaller cross-sectional area near the bay entrance due to large land reclamation in the southwest of JZB during 2005–2020. The clockwise residual current eddy b had been basically stable, a small counterclockwise residual current eddy b’ was formed at the southwest of b since 2005 (Figures 8E,F), mainly due to great coastline changes in the southwest of JZB. Increased reclamation projects on the north side of Huang Island had gradually narrowed the scope of counterclockwise residual current eddy c. The clockwise residual current eddy d tended to move westward and its intensity became weaker over time, which was probably related to the position change of Cangkou Channel and reclamation seaward of east coastline. The counterclockwise residual current eddy e barely changed before 1986, but its scope decreased and intensity became weaker and weaker thereafter due to large land reclamation in the northeast of JZB.

The residual current was very weak in the northeast of JZB before 1935 because of large areas of tidal flats. Since 1966, a counterclockwise residual current eddy f formed in the northeast of JZB (Figure 8C). The counterclockwise residual current eddy g between Dagu estuary and Hong Island tends to move to the southwest along with land reclamation, while the counterclockwise residual current eddy h became weaker and weaker, and finally disappeared in 2005.

Overall, most residual current eddies became weaker and smaller with geometry change of JZB. Some local residual current eddies disappeared while new eddies appeared at the places with the greatest coastline change. Change of residual flow field will directly affect the water exchange and transport of pollutants in JZB.

The coastline changed a lot due to land reclamation over the last 150 years as discussed in section “Coastline and Tidal Prism Change of Jiaozhou Bay,” which leads to some differences in grid resolution, especially near the coastline change area. According to the sensitivity analysis of horizontal diffusion coefficient conducted by Yuan et al. (2021), the horizontal diffusion coefficient D_h in Equation 4 is set at a constant of 3 m²/s, in order to eliminate the effects of grid resolution on the residence time between different scenarios. Thus, the differences in the residence time were affected only by reclamation.

Spatial distribution of residence time in each scenario shown in Figure 9 indicated that the residence time increased northward from the bay entrance in general, with the smallest in the southern bay entrance (less than 5 days) and largest in the northeast (130–160 days).

Distribution of residence time in JZB is mainly controlled by three factors: the distance to the bay entrance, the pattern of the tidal current field, and size of the inter-tidal flat. The southern bay entrance is the only channel connecting JZB and the open water area of Yellow Sea. Therefore, the water near the bay entrance is preferentially diluted by the open sea water, which makes the water exchange time lower in the southern bay entrance and higher in the north bay head. In addition, the water near the bay entrance easily exchanged with the open sea water due to strong eddies of residual current of a and b as shown in Figure 8. The gradient of residence time nearby the north of Huang Island is relatively large, mainly due to the existence of a counterclockwise residual current eddy c in this area, making the tracers entering this eddy difficult to move outside the bay. Most of the tracers in the northwest of JZB are transported southeastward through the residual current eddy d, then southward along the east coast of JZB, and finally transported outside the bay by the two strong eddies of residual current a and b. The northeast of JZB is mainly controlled by the counter-clockwise eddy of residual current e, which makes it very difficult for the tracers to be transported to the open sea. Thus, the residence time in the northeast is the largest.

Most of the tidal flats, located at the northern bay head and the west of Huang Island, were far from the bay entrance, resulting in longer exchange time. In addition, the tidal current on the tidal flats is very weak, making it especially difficult for the tracer to exchange with the nearby waters. Thus, the residence time in tidal flats is relatively large and area of tidal flats also affects the basin residence time of JZB.

Basin residence time, corresponding to the TM of 37% calculated by Equation 6, was shown in Figure 10. Basin residence time slightly decreased during 1863–1935, because the area of tidal flats slightly decreased and the residual current stayed the same. Basin residence time then increased from 15.5 days (1935) to 17.6 days (2020), which indicated that the basin water-exchange capacity deteriorated significantly since 1935. The increased rate of basin residence time was smaller since 2005, because the land reclamation was on a smaller scale.

Along with the quick economic development and population growth of Qingdao, the wastewater discharge increased significantly from 147 million tons per year in 1987 to 633 million tons in 2019 (Qingdao Municipal Statistics Bureau, 2020). In addition, with large areas of land reclamation, the ecosystem was deteriorating. JZB is now a eutrophic ecosystem with high nitrogen and phosphorus concentrations, resulting in both the recession of marine resources and the environmental carrying capacity (Zhang et al., 2017). Currently, the main pollutants are dissolved inorganic nitrogen (DIN), phosphorus, and petroleum pollutants. The concentration of DIN has increased considerably and is regarded as the dominant pollutant. According to field measurements by Qingdao Environment Monitoring Centre, the average concentration of JZB is 286 μg/L in 2014, corresponding to second grade of the Chinese national seawater quality standard. The water quality in the northeast of JZB is even worse than the fourth grade. Therefore, the DIN was chosen as a pollutant indicator to investigate the effect of long-term land reclamation on the pollutant transport in this section. The DIN discharge measured in 2007 at nine sewage outlets and rivers in JZB was set in each year (Yuan et al., 2016). Positions of sources were adjusted according to the coastline evolution. All the scenarios have the same initial condition and boundary forcing.

The spatial distribution of spring-neap tidal averaged DIN in each year is presented in Figure 11. The area of seawater with concentrations of DIN larger than 5 μg/L, which can reflect the diffusion range of DIN from the discharge sources, decreased and the concentration gradient became larger from 1986 to 2020, which indicated worse water exchange after intense human intervention of JZB. The relatively large offshore residual current near the Dagu River estuary shown in Figures 8A–C favored the pollutant spreading to offshore area. However, this residual current became very weak since 1986 (Figure 8D), so the spread area narrowed and the contour lines of DIN concentration were becoming almost parallel to the coastlines. Similarly, situations happened in waters near Yang River estuary, the Lianwan sewage, and the eastern bay.

According to the Chinese national seawater quality standard, seawaters with concentrations of DIN less than 200, 300, 400, and 500 μg/L are classified as the first, second, third, and fourth grade respectively, where the fourth grade is regarded as poor quality. Thus, seawater with DIN concentrations larger than 500 μg/L represent serious deterioration of seawater quality and urgently needs to be treated and improved. The water area with DIN concentrations larger than 500 μg/L increased along with intense land reclamation, especially in the northeastern area of JZB. Decision-makers in the local government should pay more attention to improving the water quality of this area.