Total Maximum Allocated Load of Chemical Oxygen Demand Near Qinhuangdao in Bohai Sea: Model and Field Observations

: Total maximum allocated load (TMAL) is the maximum sum total of all the pollutant loading a water body can carry without surpassing the water quality criterion, which is dependent on hydrodynamics and water quality conditions. A coupled hydrodynamic and water quality model combined with ﬁeld observation was used to study pollutant transport and TMAL for water environment management in Qinhuangdao (QHD) sea in the Bohai Sea in northeastern China for the ﬁrst time. Temporal and spatial variations of the chemical oxygen demand (COD) concentration were investigated based on MIKE suite (Danish Hydraulic Institute, Hørsholm, Denmark). A systematic optimization approach of adjusting the upstream pollutant emission load was used to calculate TMAL derived from the predicted COD concentration. The pollutant emission load, TMAL, and pollutant reduction of Luanhe River were the largest due to the massive runo ﬀ , which was identiﬁed as the most inﬂuential driving factor for water environmental capacity and total carrying capacity of COD. The correlation analysis and Spearman coe ﬃ cient indicate strong links between TMAL and forcing factors such as runo ﬀ , kinetic energy, and pollutant emission load. A comparison of total carrying capacity in 2011 and 2013 conﬁrms that the upstream pollutant control scheme is an e ﬀ ective strategy to improve water quality along the river and coast. Although, the present model results suggest that a monitoring system could provide more e ﬃ cient total capacity control. The outcome of this study establishes the theoretical foundation for coastal water environment management strategy in this region and worldwide.


Introduction
With growing industrialization and urbanization, effluent discharge and human activities impose increasing pressure on coastal ecosystems and natural resources, especially in densely populated coastal zones [1]. Coastal water is contaminated by discharge of land-based pollutant and marine aquaculture, which causes increasing water pollution for residents as well as risk of eutrophication. Red tide and harmful algal bloom due to eutrophication are common phenomena in coastal regions in linked with hydrodynamic variables and correlations between TMAL, runoff, kinetic energy, and emission load were obtained. The comparison of water environmental capacity and total capacity control of pollutants over the past years indicates that the environment management regulations and strategies efficiently improved QHD coastal water quality. The objective of the present study is to better understand the driving mechanisms of water environmental capacity in the coastal area of QHD. The outcome of this study should provide guidelines for robust coastal water environment management strategies in other coastal regions around the globe.

Field Site and Model Setup
QHD is located in the northeastern Hebei Province in China and adjacent to the Bohai Sea (39 • 24 ~40 • 37 N, 118 • 33 ~119 • 51 E). A large amount of recreational bathing beaches are distributed along the 163 km coastline of QHD, as shown in Figure 1. QHD is within the middle latitude zone with a humid and semi-humid continental monsoon climate of a warm temperate zone. QHD has an annual sunshine duration of 2700~2850 h, temperature of 8.8~11.3 • C, precipitation of 650~750 mm, and average evaporation of 1468.7 mm. Since the average temperature of QHD in winter is −3.7 • C, the rivers in QHD are frozen for about 108 days from 21 November to 31 March [35].
Water 2020, 12, x FOR PEER REVIEW 3 of 17 capacity in the coastal area of QHD. The outcome of this study should provide guidelines for robust coastal water environment management strategies in other coastal regions around the globe.

Field Site and Model Setup
QHD is located in the northeastern Hebei Province in China and adjacent to the Bohai Sea (39°24′ 40°37′ N, 118°33′~119°51′ E). A large amount of recreational bathing beaches are distributed along the 163 km coastline of QHD, as shown in Figure 1. QHD is within the middle latitude zone with a humid and semi-humid continental monsoon climate of a warm temperate zone. QHD has an annual sunshine duration of 2700~2850 h, temperature of 8.8~11.3 °C, precipitation of 650~750 mm, and average evaporation of 1468.7 mm. Since the average temperature of QHD in winter is −3.7 °C, the rivers in QHD are frozen for about 108 days from 21 November to 31 March [35]. QHD has a drainage basin over 30 km 2 and an annual average runoff volume of 12.6 × 10 8 m 3 . The interannual variation of runoff ranges from 2.5 × 10 8 to 41.8 × 10 8 m 3 . River flow in QHD area comes from a northwest to southeast direction before discharging into the Bohai Sea. Qilihai Lagoon collects water from four rivers and the Tanghe estuary from two rivers, which are approximated as one river. Thus there are 11 seagoing rivers including Luanhe River (LH), Qilihai Lagoon (QLH), Dapuhe River (DPH), Dongshahe River (DSH), Renzaohe River (RZH), Yanghe River (YH), Daihe River (DH), Xinhe River (XH), Tanghe River (TH), Xinkaihe River (XKH), and Shihe River (SH) from south to north, as numbered and indicated in Figure 1b. The rainfall mainly occurs during summer and autumn, and rarely during spring and winter, which causes significant seasonal variation of runoff. The annual mean runoff in 2013 ranged from 0.25 m 3 /s to 39.74 m 3 /s.
According to the marine function zoning regulation in Hebei Province (2011-2020), QHD is within the scenic tourism and recreational entertainment zone and boasts one of the best bathing beaches in China. In addition to tourism, mariculture is also important to the local economy. It follows that QHD has to meet the requirements of the criterion of second category water quality [36]. The COD has to be less than or equal to 3.0 mg/L based on the criterion of second category water quality according to the "National Marine Water Quality Criteria" (GB3097-1997) [37]. QHD has a drainage basin over 30 km 2 and an annual average runoff volume of 12.6 × 10 8 m 3 . The interannual variation of runoff ranges from 2.5 × 10 8 to 41.8 × 10 8 m 3 . River flow in QHD area comes from a northwest to southeast direction before discharging into the Bohai Sea. Qilihai Lagoon collects water from four rivers and the Tanghe estuary from two rivers, which are approximated as one river. Thus there are 11 seagoing rivers including Luanhe River (LH), Qilihai Lagoon (QLH), Dapuhe River (DPH), Dongshahe River (DSH), Renzaohe River (RZH), Yanghe River (YH), Daihe River (DH), Xinhe River (XH), Tanghe River (TH), Xinkaihe River (XKH), and Shihe River (SH) from south to north, as numbered and indicated in Figure 1b. The rainfall mainly occurs during summer and autumn, and rarely during spring and winter, which causes significant seasonal variation of runoff. The annual mean runoff in 2013 ranged from 0.25 m 3 /s to 39.74 m 3 /s.

Hydrodynamic and Water Quality Model
According to the marine function zoning regulation in Hebei Province (2011-2020), QHD is within the scenic tourism and recreational entertainment zone and boasts one of the best bathing beaches in China. In addition to tourism, mariculture is also important to the local economy. It follows that QHD has to meet the requirements of the criterion of second category water quality [36]. The COD has to be less than or equal to 3.0 mg/L based on the criterion of second category water quality according to the "National Marine Water Quality Criteria" (GB3097-1997) [37].

Hydrodynamic and Water Quality Model
TMAL was calculated based on hydrodynamic and water quality models. Field observations tended to be point measurements, whereas numerical simulation captured both spatial and temporal variations of hydrodynamics and contamination. MIKE 21 FM (flow model) is a two-dimensional modeling system in MIKE suite (Release 2014), which was applied in this study to simulate the hydrodynamic and pollutant transport. In MIKE 21 FM coupled with the hydrodynamic model, the water quality model was driven by tide and runoff. Using Boussinesq approximation, hydrostatic pressure, and shallow water conditions, the Reynolds averaged Navier-Stokes equations were reduced to generalized shallow water equations using the controlling volume method. The governing equations of the numerical model are demonstrated in the Appendix A [38].
The unstructured nested grid from the Bohai Sea to QHD was constructed based on World Geodetic System 1984 ( Figure 1b). The larger mesh of the Bohai Sea took the line connecting the tide stations of Dalian and Yantai as the offshore boundary, whose tidal level data were from the National Marine Data and Information Service of China [39]. Then the larger mesh domain provided a tidal level boundary for smaller mesh domain of the QHD sea. The mesh of QHD has 9442 nodes and 17,973 elements with WE span of 110 km and NS span of 120 km. The mesh size is 7 km in the open sea and 15 m near the coastline. The higher resolution near the coastline is necessary to resolve the complex coastal bathymetry and shoreline.
The boundary conditions at 11 river mouths were set by monthly average runoff and COD discharge in 2013. The flood season is from June to September and the largest runoff for each river occurs in July or August. December to March are the dry months with less runoff. LH is the largest river in QHD, with runoff much larger than other rivers. The mean discharge was over 1.72 × 10 9 t from July to September with a peak value of 4.28 × 10 9 t in August. Bottom stress was determined by a quadratic friction law with a bottom roughness from 0.01 to 0.015 s/m 1/3 through water depth and sediment diameter.
For the water quality model, non-conservative matter such as COD was not only influenced by flow convection and diffusion, but also by degradation processes such as microbial assimilation and metabolism. The initial COD concentration was set as the historical average value of 1.3 mg/L, and the decay rate for degradation was 0.03/d through model calibration [11]. The horizontal dispersion was determined by the eddy viscosity formulation dependent on the Prandtl number which was set as 1 by calibration in this study [40].

Field Observations
Predicted and observed water level and current were compared with each other to assess the performance of the hydrodynamic model. The hydrodynamic measurement work in May of 2013 was for other coastal projects by the Marine Environmental Monitoring Center of Hebei Province Qinhuangdao, China. However, measurements from July to August when the largest runoff occurs would be better for calibrating the model. Observations at 10 current stations (V1-V10) and 1 water level station (WL), as indicated in Figure 2a, provided measurements from 0:00 on 11 May to 24:00 on 12 May and from 0:00 on 16 May to 24:00 on 17 May 2013. The water level was measured at the tide station near XKH estuary based on the limnimeter (DCX-22, Keller, Zug, Switzerland). For current observation, 10 stations with an acoustic Doppler current profiler (ADCP) were set in QHD coastal area, and vertical mean current speed and direction were adopted. The observed COD concentrations given by the sampling laboratory report of 28 observation stations (C1-C28), as shown in Figure 2b, were used to validate the water quality model. The water samples were collected on 11 August 2013 and brought back to the laboratory for determination of the permanganate index (COD Mn ). The runoff was measured once a month in every seagoing river by a flowmeter in 2013; in the same manner, the water was sampled for determination of the dichromate index (COD Cr ) which was converted into COD Mn as the river boundary. In this study, COD Mn is treated as the major water quality research object and is uniformly called COD concentration in later sections.

Model Validations
A favorable linear correlation between measured and model predicted COD concentrations confirmed the good agreement between the model and data. R 2 was used to quantify the model-data agreement in the range of 0 to 1 (i.e., if R 2 was close to 1, the correlation was higher). Figure 3 indicates good agreements between the measured and computed water level, current speed, and current direction. The R 2 values were 0.80 and 0.87 for the computed and observed water level, 0.80 and 0.77 for current speed, and 0.94 and 0.92 for current direction on 11-12 May and 16-17 May, respectively. Except for the current speed on 16-17 May, all R 2 values were larger than 0.80, which indicates good performance of the hydrodynamic model, even though the hydrodynamic conditions in the QHD sea were complicated and variable due to the presence of an amphidromic point in this region. As shown in Figure 4, the computed and measured COD concentrations compare well. The differences between measured COD concentration once a day and daily average computed COD concentration were acceptable. The water quality model performance of COD concentration was evaluated by the root mean square error (RMSE) [41]. The RMSE of COD concentration was 0.17 mg/L, which indicated that an accurate COD concentration could be obtained from the validated TMAL model.

Model Validations
A favorable linear correlation between measured and model predicted COD concentrations confirmed the good agreement between the model and data. R 2 was used to quantify the model-data agreement in the range of 0 to 1 (i.e., if R 2 was close to 1, the correlation was higher). Figure 3 indicates good agreements between the measured and computed water level, current speed, and current direction. The R 2 values were 0.80 and 0.87 for the computed and observed water level, 0.80 and 0.77 for current speed, and 0.94 and 0.92 for current direction on 11-12 May and 16-17 May, respectively. Except for the current speed on 16-17 May, all R 2 values were larger than 0.80, which indicates good performance of the hydrodynamic model, even though the hydrodynamic conditions in the QHD sea were complicated and variable due to the presence of an amphidromic point in this region. As shown in Figure 4, the computed and measured COD concentrations compare well. The differences between measured COD concentration once a day and daily average computed COD concentration were acceptable. The water quality model performance of COD concentration was evaluated by the root mean square error (RMSE) [41]. The RMSE of COD concentration was 0.17 mg/L, which indicated that an accurate COD concentration could be obtained from the validated TMAL model. Water 2020, 12, x FOR PEER REVIEW 6 of 17

Model Results of Total Maximum Allocated Load (TMAL)
TMAL was determined by source of pollution and water quality distribution characteristics. Due to significant spatial variation of water quality in the QHD sea, it was necessary to execute the computation of the TMAL for each river. As shown in the computational procedure of the TMAL in Figure 5, COD was selected as the focus of this study, and it was necessary to impose an upper limit of COD concentration. The COD concentration control criterion of 2.0 mg/L was empirically determined by the water quality criteria. Based on the model simulation of COD concentration in the QHD sea, the area where the criterion exceeded in estuaries was identified. If no area exceeded the criterion, the TMAL was treated as the good practice pollutant emission load. Otherwise computation of TMAL needed to be determined through an optimization process by cutting down the pollutant emission load. The COD discharge at each river mouth boundary was determined through trial and error (i.e., with a reduction of −20%, −40%, −60%, −80%, and −100% for all rivers successively). If the reduction proportion could not satisfy the water quality criterion, a stronger cut of COD discharge would be adopted until that led to zero areas exceeding the criterion. Then an adjustment process of COD emission load for a single river could improve the accuracy of TMAL. The pollutant emission

Model Results of Total Maximum Allocated Load (TMAL)
TMAL was determined by source of pollution and water quality distribution characteristics. Due to significant spatial variation of water quality in the QHD sea, it was necessary to execute the computation of the TMAL for each river. As shown in the computational procedure of the TMAL in Figure 5, COD was selected as the focus of this study, and it was necessary to impose an upper limit of COD concentration. The COD concentration control criterion of 2.0 mg/L was empirically determined by the water quality criteria. Based on the model simulation of COD concentration in the QHD sea, the area where the criterion exceeded in estuaries was identified. If no area exceeded the criterion, the TMAL was treated as the good practice pollutant emission load. Otherwise computation of TMAL needed to be determined through an optimization process by cutting down the pollutant emission load. The COD discharge at each river mouth boundary was determined through trial and error (i.e., with a reduction of −20%, −40%, −60%, −80%, and −100% for all rivers successively). If the reduction proportion could not satisfy the water quality criterion, a stronger cut of COD discharge would be adopted until that led to zero areas exceeding the criterion. Then an adjustment process of COD emission load for a single river could improve the accuracy of TMAL. The pollutant emission

Model Results of Total Maximum Allocated Load (TMAL)
TMAL was determined by source of pollution and water quality distribution characteristics. Due to significant spatial variation of water quality in the QHD sea, it was necessary to execute the computation of the TMAL for each river. As shown in the computational procedure of the TMAL in Figure 5, COD was selected as the focus of this study, and it was necessary to impose an upper limit of COD concentration. The COD concentration control criterion of 2.0 mg/L was empirically determined by the water quality criteria. Based on the model simulation of COD concentration in the QHD sea, the area where the criterion exceeded in estuaries was identified. If no area exceeded the criterion, the TMAL was treated as the good practice pollutant emission load. Otherwise computation of TMAL needed to be determined through an optimization process by cutting down the pollutant emission load. The COD discharge at each river mouth boundary was determined through trial and error (i.e., with a reduction of −20%, −40%, −60%, −80%, and −100% for all rivers successively). If the reduction proportion could not satisfy the water quality criterion, a stronger cut of COD discharge would be adopted until that led to zero areas exceeding the criterion. Then an adjustment process of COD emission load for a single river could improve the accuracy of TMAL. The pollutant emission loads needed to be further decreased for rivers with relatively large pollutant discharge, meanwhile, the loads could be increased for less polluted rivers by ±5%, ±10%, ±15%, and ±20%. When the COD concentration started to exceed the control criterion, TMAL was the pollutant emission load after adjustments using the optimization scheme.
Water 2020, 12, x FOR PEER REVIEW 7 of 17 loads needed to be further decreased for rivers with relatively large pollutant discharge, meanwhile, the loads could be increased for less polluted rivers by ±5%, ±10%, ±15%, and ±20%. When the COD concentration started to exceed the control criterion, TMAL was the pollutant emission load after adjustments using the optimization scheme.

COD Transport in Qinhuangdao Sea
A coupled hydrodynamic and water quality numerical model was validated against field observations. The model was used to simulate hydrodynamics and COD concentrations of the QHD sea in 2013, which supplemented data for large-scale distribution of COD concentration. TMAL was calculated based on the COD concentration distribution which was dominated by tidal current forcing. The variations of water level and current speed near YH estuary (39°46′ N, 119°28′ E) are shown in Figure 6a. The regular semi-diurnal tide shows two peaks and two minima during a daily tidal cycle; however, the hydrodynamics were complicated at the QHD sea, due to approximation to the amphidromic point and the fact that the reversal of the tidal current lagged behind the reversal of the water level [42]. The four typical moments of tidal current were described as maximum flood, flood slack, maximum ebb, and ebb slack, and the flow fields of typical moments are shown in Figure  6b-e. The current was reciprocating flow parallel with coastline and the current direction was westsouthwest as flood tide and east-northeast as ebb tide. The current speed magnitude was relatively smaller with inconspicuous differences between flood tide and ebb tide. The average current speed was 0.13~0.3 m/s during flood tide and 0.14~0.32 m/s during ebb tide. At the moment of maximum flood and maximum ebb, the current speeds of the coastal area were smaller than those of open sea, and the current speeds near the LH estuary were relatively larger with 0.4~0.6 m/s. At the moment of flood slack, the current speed near YH and DH estuary was close to 0 m/s, but the current near LH estuary was still flood tide, and the current in the north of the QHD sea was ebb tide. At the moment of ebb slack, the current in the north of QHD sea was flood tide and near LH estuary was ebb tide. It meant that turns of tidal current were not at the same time in spatial distribution of QHD sea; the turns of the tidal current in the north area were in advance of the south area.

COD Transport in Qinhuangdao Sea
A coupled hydrodynamic and water quality numerical model was validated against field observations. The model was used to simulate hydrodynamics and COD concentrations of the QHD sea in 2013, which supplemented data for large-scale distribution of COD concentration. TMAL was calculated based on the COD concentration distribution which was dominated by tidal current forcing. The variations of water level and current speed near YH estuary (39 • 46 N, 119 • 28 E) are shown in Figure 6a. The regular semi-diurnal tide shows two peaks and two minima during a daily tidal cycle; however, the hydrodynamics were complicated at the QHD sea, due to approximation to the amphidromic point and the fact that the reversal of the tidal current lagged behind the reversal of the water level [42]. The four typical moments of tidal current were described as maximum flood, flood slack, maximum ebb, and ebb slack, and the flow fields of typical moments are shown in Figure 6b-e. The current was reciprocating flow parallel with coastline and the current direction was west-southwest as flood tide and east-northeast as ebb tide. The current speed magnitude was relatively smaller with inconspicuous differences between flood tide and ebb tide. The average current speed was 0.13~0.3 m/s during flood tide and 0.14~0.32 m/s during ebb tide. At the moment of maximum flood and maximum ebb, the current speeds of the coastal area were smaller than those of open sea, and the current speeds near the LH estuary were relatively larger with 0.4~0.6 m/s. At the moment of flood slack, the current speed near YH and DH estuary was close to 0 m/s, but the current near LH estuary was still flood tide, and the current in the north of the QHD sea was ebb tide. At the moment of ebb slack, the current in the north of QHD sea was flood tide and near LH estuary was ebb tide. It meant that turns of tidal current were not at the same time in spatial distribution of QHD sea; the turns of the tidal current in the north area were in advance of the south area. In the water quality model, the COD variation in QHD was combined with monthly hydrodynamic conditions and COD discharge at river boundaries. A favorable water exchange capacity helps to alleviate water pollutants and strong tidal stream energy indicates rapid water exchange [43]; therefore, the kinetic energy was derived from the current speed to analyze its influence on water quality distribution.
The spatial distributions of annual average COD concentration and annual average kinetic energy for 2013 are shown in Figure 7a,b. The COD concentration in the nearshore was higher than that in the offshore and it was higher in the southwest than northeast of QHD coastal water. It is evident from Figure 7a that isoline of 1.8 mg/L mainly appeared at the coast next to LH and from DPH to DH, while in other areas, the COD concentration was less than 1.8 mg/L. Serious pollution concentrated along the coast between LH and DH. The annual average COD concentrations for the 11 estuaries are shown in Figure 7c, and the COD concentrations of LH, YH, and DH estuaries were In the water quality model, the COD variation in QHD was combined with monthly hydrodynamic conditions and COD discharge at river boundaries. A favorable water exchange capacity helps to alleviate water pollutants and strong tidal stream energy indicates rapid water exchange [43]; therefore, the kinetic energy was derived from the current speed to analyze its influence on water quality distribution.
The spatial distributions of annual average COD concentration and annual average kinetic energy for 2013 are shown in Figure 7a,b. The COD concentration in the nearshore was higher than that in the offshore and it was higher in the southwest than northeast of QHD coastal water. It is evident from Figure 7a that isoline of 1.8 mg/L mainly appeared at the coast next to LH and from DPH to DH, while in other areas, the COD concentration was less than 1.8 mg/L. Serious pollution concentrated along the coast between LH and DH. The annual average COD concentrations for the 11 estuaries are shown in Figure 7c, and the COD concentrations of LH, YH, and DH estuaries were larger than 3.0 mg/L and exceeded the criterion of second category water quality. The COD concentration of LH estuary with about 6.4 mg/L in 2013 was the highest among the seagoing rivers in QHD. The COD concentrations at YH and DH were 5.7 mg/L and 3.2 mg/L, respectively. The water quality at DPH was close to the criterion of the second category with COD concentration of 2.8 mg/L. While the COD concentrations of RZH and DSH were lower than 2 mg/L, due to the shorter distance with YH, DH, and DPH, the coastal environment from DPH to DH was severely polluted. The COD concentration of XKH estuary was the lowest at 1.4 mg/L which showed fewer pollutant emissions, and the concentrations of TH, XH, and QLH were lower than 2 mg/L. COD in QLH coastal water was larger than 1.6 mg/L, which was mainly polluted by adjacent estuaries. From the perspective of mean kinetic energy distribution (Figure 7b), the kinetic energies in the southwest and northeast of QHD sea and open sea were relatively larger than those in the nearshore from LH to SH. Except for LH, QLH, and SH, the kinetic energies in adjacent water areas were lower than 0.01 J/kg, which meant a weaker coastal water exchange capacity in the QHD coastal area. Although SH had a larger COD concentration of 2.8 mg/L in estuary, a stronger kinetic energy with 0.025 J/kg led to cleaner water in the adjacent sea.

Water Environmental Capacity and Total Capacity Control of COD
QHD sea is a major tourist attraction for sea baths, therefore, it needs to meet the criterion of second category water quality, so that the COD concentration should be less than or equal to 3.0 mg/L. During red tide outbreaks in Xiamen, COD concentration may exceed 2.0 mg/L [44]. In previous studies of environmental capacity for Luoyuan Bay and Leqing Bay, the COD concentration of 2.0 mg/L was used as the environmental control criterion [45,46]. In QHD sea, before the red tide event in 2011, the COD concentrations of most field stations were lower than 2.0 mg/L, however, more than half of the field observed that COD concentrations became larger than 2.0 mg/L during red tide outbreaks. It follows that COD concentration of 2.0 mg/L was adopted as the indicating index for red tide outbreak [47]. Therefore, it was required to maintain COD concentration below 2.0 mg/L in QHD The three estuaries of LH, YH, and DH had COD concentrations over the second category water quality criterion and were chosen as the case study estuaries to investigate the monthly variation of COD concentrations, as shown in Figure 7d. Winter consists of December, January, and February in QHD, and as the temperature is lower than ice point, the rivers are frozen so that there is rarely any runoff discharge into the sea and pollutant is not liable to be carried by runoff. The COD concentrations of three estuaries in December were larger than those in January and February because rivers were just beginning to freeze in December. There were two peaks in the monthly variation for LH, YH, and DH. The first peak appeared in March and the second peak appeared from July to September. As spring was coming and temperature became warmer in March, accumulation of pollutant from melting ice discharged into the sea and COD concentration increased in each estuary. Although the COD concentration dropped down in April and May because the accumulative effect of pollutant was taped down during March. The precipitation during summer and autumn in QHD was larger than other seasons according to historical meteorological and hydrological data, so that there were larger amounts of water and agricultural pollution discharged during these seasons. Summer was from June to August and autumn was from September to November. The largest COD concentration occurred in July with a magnitude of 16.4 mg/L in LH, in September with a magnitude of 19.6 mg/L in YH, and in August with a magnitude of 5.5 mg/L in DH. Though the COD discharge of LH at 3.28 kg/s in July was larger than that of YH at 2.93 kg/s in September, the maximum COD concentration among the estuaries occurred in YH estuary in September, which was mainly due to the larger kinetic energy of LH estuary than that of YH estuary, as shown in Figure 7b. Therefore, the COD concentration in an estuary was not only controlled by emission load, but also by kinetic energy.

Water Environmental Capacity and Total Capacity Control of COD
QHD sea is a major tourist attraction for sea baths, therefore, it needs to meet the criterion of second category water quality, so that the COD concentration should be less than or equal to 3.0 mg/L. During red tide outbreaks in Xiamen, COD concentration may exceed 2.0 mg/L [44]. In previous studies of environmental capacity for Luoyuan Bay and Leqing Bay, the COD concentration of 2.0 mg/L was used as the environmental control criterion [45,46]. In QHD sea, before the red tide event in 2011, the COD concentrations of most field stations were lower than 2.0 mg/L, however, more than half of the field observed that COD concentrations became larger than 2.0 mg/L during red tide outbreaks. It follows that COD concentration of 2.0 mg/L was adopted as the indicating index for red tide outbreak [47]. Therefore, it was required to maintain COD concentration below 2.0 mg/L in QHD for sustainable socio-economy development and the natural environment, which satisfied the criterion of second category water quality at the same time. In case COD concentrations exceeded 2.0 mg/L, it was necessary to reduce the concentration by further decreasing the pollutant discharge at unqualified rivers while maintaining the emission at the same level for qualified rivers.
Through simulation of COD concentrations in estuaries and optimization processes, TMAL was calculated for all estuaries every month in 2013. The pollutant reduction was estimated as emission load minus the TMAL that corresponded to the desired total capacity control. The year-round emission load, TMAL, pollutant reduction, and runoff of each seagoing river in 2013 are demonstrated in Figure 8a. Annual mean runoff of LH was the largest at 39.74 m 3 /s; that of other rivers was lower than 5 m 3 /s, except for XKH and YH which had runoff of 9.53 m 3 /s and 5.27 m 3 /s respectively. It was shown that LH had the largest emission load of 37,817.4 t/a, TMAL of 22,221.3 t/a, and pollutant reduction of 17,993.3 t/a due to the largest runoff. These differences of emission load and TMAL followed the variation of runoff except for XKH which upstream had excellent water quality and little pollutant discharge. Five rivers, SH, XKH, XH, DSH, and QLH, did not need contaminant discharge reduction. Ranking from largest to smallest pollutant reduction was LH, YH, RZH, DPH, TH, and DH, and the effluent discharge of these rivers was necessary to be controlled.
Hydrological characteristics in QHD featured significant seasonal variations due to ephemeral rivers. TMAL, emission load, and pollutant reduction of QHD in 2013 were mainly controlled by runoff with strong seasonal variation. Since temperature of QHD in winter was below the frozen point, rivers were almost all covered by ice with little runoff and effluent discharge, so that January, February, and March were combined together as the average value to be considered. As shown in Figure 8b, variation of seasonal runoff characteristics revealed that small amounts of river discharge into the sea were observed for December to May. In summer and autumn, the runoff became larger. In August, the emission load became the largest with 1.9 × 10 4 t in 2013 with the abundant runoff of 279 m 3 /s, and the emission load as well as runoff in July occupied second place. In the other months, the emission loads were lower than 1.0 × 10 4 t, especially in January to May, which were lower than 1.0 × 10 3 t. It was remarkable that variation of emission load, TMAL, and runoff had the same tendency. In other words, runoff was the main dominant control factor for pollutant emission and TMAL. The total capacity control, during the rainy season from June to September had pollutant reduction larger than 200 t, and these months became the key months to regulate pollutant emissions.
rivers. TMAL, emission load, and pollutant reduction of QHD in 2013 were mainly controlled by runoff with strong seasonal variation. Since temperature of QHD in winter was below the frozen point, rivers were almost all covered by ice with little runoff and effluent discharge, so that January, February, and March were combined together as the average value to be considered. As shown in Figure 8b, variation of seasonal runoff characteristics revealed that small amounts of river discharge into the sea were observed for December to May. In summer and autumn, the runoff became larger. In August, the emission load became the largest with 1.9 × 10 4 t in 2013 with the abundant runoff of 279 m 3 /s, and the emission load as well as runoff in July occupied second place. In the other months, the emission loads were lower than 1.0 × 10 4 t, especially in January to May, which were lower than 1.0 × 10 3 t. It was remarkable that variation of emission load, TMAL, and runoff had the same tendency. In other words, runoff was the main dominant control factor for pollutant emission and TMAL. The total capacity control, during the rainy season from June to September had pollutant reduction larger than 200 t, and these months became the key months to regulate pollutant emissions.

Correlation of TMAL with Forcing Factors
TMAL was calculated by COD concentrations near the estuaries of QHD sea which were controlled by runoff, kinetic energy, and COD discharge. To conduct correlation analysis of various influence factors (runoff, kinetic energy, and emission load) with TMAL, variables including influence factors and TMAL were collected from 11 seagoing rivers every month in 2013. The linear

Correlation of TMAL with Forcing Factors
TMAL was calculated by COD concentrations near the estuaries of QHD sea which were controlled by runoff, kinetic energy, and COD discharge. To conduct correlation analysis of various influence factors (runoff, kinetic energy, and emission load) with TMAL, variables including influence factors and TMAL were collected from 11 seagoing rivers every month in 2013. The linear relationships were derived as shown in Figure 9a-c, which suggests that TMAL had approximately linear relationship with these key factors. Emission load featured a particularly high linear correlation with the R 2 of 0.92. The R 2 values of the correlation of TMAL with runoff and kinetic energy were 0.84 and 0.72, respectively.
Spearman and Pearson coefficients are frequently used to evaluate correlation, whereas the stabilization of Pearson coefficient is influenced by extremum [48]. From the box figure shown in Figure 9d, there were extremums in each factor, therefore the Spearman coefficient was used to evaluate the correlation instead of the Pearson coefficient and the Spearman values are shown in Table 1. The positive value showed the extent of positive correlation. The Spearman value between kinetic energy and runoff was the lowest at 0.84, due to that the current speed in QHD sea was mainly dominated by tide, terrain, and coastline instead of runoff. Similarly, the kinetic energy was mainly controlled by current speed rather than emission load, so that the Spearman value of kinetic energy with emission load was relatively small at 0.85. The Spearman value of 0.88 between kinetic energy and TMAL showed kinetic energy had a significant influence on TMAL. It was evident that runoff had a close correlation with emission load and TMAL with Spearman values at 0.96 and 0.92, respectively. Runoff contributed to pollutant discharge for emission load, meanwhile it could also accelerate pollutant dispersion and increase TMAL. Runoff became the key link between emission load and TMAL with Spearman value of 0.95. and 0.72, respectively.
Spearman and Pearson coefficients are frequently used to evaluate correlation, whereas the stabilization of Pearson coefficient is influenced by extremum [48]. From the box figure shown in Figure 9d, there were extremums in each factor, therefore the Spearman coefficient was used to evaluate the correlation instead of the Pearson coefficient and the Spearman values are shown in Table 1. The positive value showed the extent of positive correlation. The Spearman value between kinetic energy and runoff was the lowest at 0.84, due to that the current speed in QHD sea was mainly dominated by tide, terrain, and coastline instead of runoff. Similarly, the kinetic energy was mainly controlled by current speed rather than emission load, so that the Spearman value of kinetic energy with emission load was relatively small at 0.85. The Spearman value of 0.88 between kinetic energy and TMAL showed kinetic energy had a significant influence on TMAL. It was evident that runoff had a close correlation with emission load and TMAL with Spearman values at 0.96 and 0.92, respectively. Runoff contributed to pollutant discharge for emission load, meanwhile it could also accelerate pollutant dispersion and increase TMAL. Runoff became the key link between emission load and TMAL with Spearman value of 0.95.

Environmental Capacity Management
Comprehensive coastal environment and pollution treatment programs of seagoing rivers have been carried out by the QHD municipal government since 2012 [49]. Since March 2012, about 254 companies with illegal pollution discharge and 25 outdated production lines have been shut down. Among them, 167 companies were closed permanently, and 87 companies suspended productions [50]. To assess the effectiveness of management projects on seagoing rives, the emission load, TMAL, and pollutant reduction of COD in 2011 were calculated and compared with those in 2013, as shown in Figure 10. The emission loads of 11 rivers in 2013 were lower than those in 2011, which indicated that the pollutant discharge was mitigated to a certain degree. TMALs in 2013 differed from those in 2011 due to different hydrodynamic conditions in these years; however, TMAL of XKH in 2013 differed from that in 2011 due to pollutant discharge variation. The TMAL of LH in 2013 was larger than that in 2011, mainly due to a larger annual mean runoff of LH at 45.66 m 3 /s in 2013 than that at 39.74 m 3 /s in 2011. SH and DH no longer required pollutant reduction, and recommended pollutant reductions of DPH and TH decreased from 2011 to 2013. The total emission load dropped down to 64,272 t in 2013 which was about 83% of that in 2011. Overall, the water environment improved significantly through upstream pollutant control. In recent years, the upstream control program such as the River Governor System, which allocates the supervision and management of river water quality to a particular person in QHD, has improved water qualities in QHD sea and adjacent bathing beaches to satisfy the second category water quality criterion [51]. Reduction of pollutants discharged into the sea from rivers is the key for coastal pollution control and marine disaster prevention. that the pollutant discharge was mitigated to a certain degree. TMALs in 2013 differed from those in 2011 due to different hydrodynamic conditions in these years; however, TMAL of XKH in 2013 differed from that in 2011 due to pollutant discharge variation. The TMAL of LH in 2013 was larger than that in 2011, mainly due to a larger annual mean runoff of LH at 45.66 m 3 /s in 2013 than that at 39.74 m 3 /s in 2011. SH and DH no longer required pollutant reduction, and recommended pollutant reductions of DPH and TH decreased from 2011 to 2013. The total emission load dropped down to 64,272 t in 2013 which was about 83% of that in 2011. Overall, the water environment improved significantly through upstream pollutant control. In recent years, the upstream control program such as the River Governor System, which allocates the supervision and management of river water quality to a particular person in QHD, has improved water qualities in QHD sea and adjacent bathing beaches to satisfy the second category water quality criterion [51]. Reduction of pollutants discharged into the sea from rivers is the key for coastal pollution control and marine disaster prevention.
The environmental capacity and total capacity control of COD was investigated based on hydrologic and water quality field observational data collected once a month. During the model simulation setup, however, the data measured once a month were regarded as the monthly average data. Therefore, construction of monitoring stations was recommended to continuously monitor the discharge and water quality of seagoing rivers, which was adopted by the local government as a pilot project. This monitoring system could measure the water quality in the estuaries of SH, XKK, TH, DH, YH, RZH, DPH, and XKH, and provided sufficient and timely data for future total pollutant capacity control. The environmental capacity and total capacity control of COD was investigated based on hydrologic and water quality field observational data collected once a month. During the model simulation setup, however, the data measured once a month were regarded as the monthly average data. Therefore, construction of monitoring stations was recommended to continuously monitor the discharge and water quality of seagoing rivers, which was adopted by the local government as a pilot project. This monitoring system could measure the water quality in the estuaries of SH, XKK, TH, DH, YH, RZH, DPH, and XKH, and provided sufficient and timely data for future total pollutant capacity control.

Conclusions
Premium water quality at Qinhuangdao (QHD) sea in the Bohai Sea is necessary for recreational sea bathing activities, therefore, pollutant management projects have been implemented to improve environmental capacity in this region. In this study, a coupled hydrodynamic and water quality model was constructed to assess water environmental capacity in QHD sea, after model validation with field measurements of tide, current, and water quality. Spatial and monthly temporal variations of COD concentration from 11 rivers to QHD sea were investigated through model simulation for 2013 using MIKE. Total maximum allocated load (TMAL) was derived from the predicted COD concentration and optimized by adjusting the pollutant emission load.
It was found that runoff was the most influential driving factor for pollutant emission load, water environmental capacity, and total pollution capacity control. The seasonal variation of COD indicated that water quality was particularly problematic during the rainy season of summer and autumn when the runoff was at the peak, but less so during winter when rivers were frozen, and less water and pollutant were discharged into the sea. Spearman correlation coefficient results indicated a linear relationship between TMAL and forcing factors such as runoff, kinetic energy, and emission load. Comparisons of TMAL, emission load, and pollutant reduction at QHD Sea in 2011 and 2013 indicated that the upstream pollutant control management scheme such as River Governor System improved coastal water quality significantly. Similar practice may be applied to other coastal areas. Meanwhile, based on previous investigations of environmental capacity, a continuous monitoring system was adopted, and a pilot project is in operation to collect a large data set to achieve more accurate total capacity control.
The transport of pollutant variables that are present in the water column is simulated by advection-dispersion process based on water quality and the hydrodynamic model. The transport equation for water quality model is obtained as follows: where C is the depth average concentration; F C is the horizontal diffusion term; k p is the linear decay rate; C s is the concentration of source discharge.