Spatial Variation in Aragonite Saturation State and the Inﬂuencing Factors in Jiaozhou Bay, China

: Both natural processes and human activities a ﬀ ect seawater calcium carbonate saturation state ( Ω arag ), while the mechanisms are still far from being clearly understood. This study analysed the seawater surface Ω arag during summer and winter in Jiaozhou Bay (JZB), China, based on two cruises observations performed in January and June 2017. The ranges of Ω arag values were 1.55~2.92 in summer and 1.62~2.15 in winter. Regression analyses were conducted to identify the drivers of the change of Ω arag distribution, and then the relative contributions of temperature, mixing processes and biological processes to the spatial di ﬀ erences in Ω arag were evaluated by introducing the di ﬀ erence between total alkalinity (TA) and dissolved inorganic carbon (DIC) as a proxy for Ω arag . The results showed that biological processes were the main factor a ﬀ ecting the spatial di ﬀ erences in Ω arag , with relative contributions of 70% in summer and 50% in winter. The contributions of temperature (25% in summer and 20% in winter) and the mixing processes (5% in summer and 30% in winter) were lower. The increasing urbanization in o ﬀ shore areas can further worsen acidiﬁcation, therefore environmental protection in both o ﬀ shore and onshore is needed.


Introduction
Ocean acidification (OA) not only causes physical and chemical changes in the ocean but is also of widespread concern due to its serious impact on marine ecosystems and calcifying organisms [1][2][3][4]. Coastal waters are the most active areas for fishery production and shellfish farming, and studies related to OA mechanisms are a top priority in the region [5][6][7]. Unlike the open ocean, nearshore acidification is not only affected by global OA caused by increasing atmospheric CO 2 levels but also by the coupled effects of various processes, such as upwelling, organic matter degradation, primary production of phytoplankton and river input, due to the strong interaction between land and ocean, making the acidification mechanism in coastal waters more complicated [8][9][10][11]. In bays, which are highly affected by human disturbance, the mutual enhancement, interference and superimposition of the effects of the above processes are especially prominent. Therefore, clarifying the OA mechanism in this region will be helpful for understanding the OA process in coastal waters under the coupled influence of natural processes and human activities, and will be of great significance for predicting the developing trend and impact of OA and formulating appropriate mitigation measures in coastal waters. dominated by the East Asian Monsoon, with northerly wind in winter and southerly wind in summer. The tidal current in JZB is a regular semi-diurnal tide that induces strong vertical mixing, resulting in nearly homogeneous vertical profiles in term of seawater temperature and salinity [34]. Most of the rivers that enter JZB have no natural runoff under most conditions except for significant heavy rainfall. In the eastern area of the bay adjacent to the downtown of Qingdao, a wastewater treatment plant (WWTP) is located beside each of the major estuaries (Licun River, Haibo River and Loushan River). Therefore, the estuaries have become conduits for wastewater [27]. In the western area of the bay, which has an ecological wetland with an area over 50 km 2 , the seawater is relatively cleaner [35].
part of the bay. JZB covers a water area of 302.9 km 2 , and the average water depth is 7 m. The climate is dominated by the East Asian Monsoon, with northerly wind in winter and southerly wind in summer. The tidal current in JZB is a regular semi-diurnal tide that induces strong vertical mixing, resulting in nearly homogeneous vertical profiles in term of seawater temperature and salinity [34]. Most of the rivers that enter JZB have no natural runoff under most conditions except for significant heavy rainfall. In the eastern area of the bay adjacent to the downtown of Qingdao, a wastewater treatment plant (WWTP) is located beside each of the major estuaries (Licun River, Haibo River and Loushan River). Therefore, the estuaries have become conduits for wastewater [27]. In the western area of the bay, which has an ecological wetland with an area over 50 km 2 , the seawater is relatively cleaner [35].

Survey Stations and Sample Processing
The basic data included those collected at 28 stations on 28 June 2017 during summer [35] and on 17 January 2017 during winter ( Figure 1). All seawater samples were collected from the seawater surface at a depth of approximately 1.5 m because of the homogeneous vertical profiles in the JZB water column [34]. In addition to basic hydrological parameters, such as seawater temperature and salinity, dissolved oxygen saturation (DO%), dissolved inorganic carbon (DIC) and pH were also measured. Among the measured parameters, seawater temperature, salinity and DO% were automatically recorded each second and averaged over 1 min, while DIC was measured within one week after the seawater sample was transferred to the laboratory.
Surface temperature and salinity were measured using an SBE 45 Micro TSG (Sea-Bird, Inc., Bellevue, WA, USA), with a nominal precision of 0.002 °C for temperature and 0.005 for salinity. The DO% was measured with a YSI-5000 oxygen analyser (YSI Corporation, Yellow Spring, OH, USA), which was calibrated using the Winkler titration method (nominal precision: 0.1%). The DIC samples were directly collected from the Niskin bottles using a syringe and filtered through a 0.45 μm disposable syringe filter to avoid exchange with the air. The DIC samples were all poisoned with saturated mercury chloride (final concentration: c. 0.02% by volume) and measured by acid extraction using a DIC analyser (AS-C2, Apollo SciTech, Newark, DE, USA). In the process of DIC measurements, Certified Reference Materials (CRMs) from Scripps Institution of Oceanography were used for quality control, and the precision and accuracy level were 0.1%. pH was determined in situ using an Orion Star-A211 Plus pH Benchtop Meter with a Ross pH electrode (Thermo Fisher Scientific Inc., Beverly, MA, USA), calibrated according to standards from the National Bureau of Standards (NBS Standard) with values of 4.00, 6.86 and 9.18 at 25 °C. The precision of pH measurement was 0.01. Additionally, pH samples were also transferred to the laboratory and measured at a constant temperature of 25 °C for calibration.

Survey Stations and Sample Processing
The basic data included those collected at 28 stations on 28 June 2017 during summer [35] and on 17 January 2017 during winter ( Figure 1). All seawater samples were collected from the seawater surface at a depth of approximately 1.5 m because of the homogeneous vertical profiles in the JZB water column [34]. In addition to basic hydrological parameters, such as seawater temperature and salinity, dissolved oxygen saturation (DO%), dissolved inorganic carbon (DIC) and pH were also measured. Among the measured parameters, seawater temperature, salinity and DO% were automatically recorded each second and averaged over 1 min, while DIC was measured within one week after the seawater sample was transferred to the laboratory.
Surface temperature and salinity were measured using an SBE 45 Micro TSG (Sea-Bird, Inc., Bellevue, WA, USA), with a nominal precision of 0.002 • C for temperature and 0.005 for salinity. The DO% was measured with a YSI-5000 oxygen analyser (YSI Corporation, Yellow Spring, OH, USA), which was calibrated using the Winkler titration method (nominal precision: 0.1%). The DIC samples were directly collected from the Niskin bottles using a syringe and filtered through a 0.45 µm disposable syringe filter to avoid exchange with the air. The DIC samples were all poisoned with saturated mercury chloride (final concentration: c. 0.02% by volume) and measured by acid extraction using a DIC analyser (AS-C2, Apollo SciTech, Newark, DE, USA). In the process of DIC measurements, Certified Reference Materials (CRMs) from Scripps Institution of Oceanography were used for quality control, and the precision and accuracy level were 0.1%. pH was determined in situ using an Orion Star-A211 Plus pH Benchtop Meter with a Ross pH electrode (Thermo Fisher Scientific Inc., Beverly, MA, USA), calibrated according to standards from the National Bureau of Standards (NBS Standard) with values of 4.00, 6.86 and 9.18 at 25 • C. The precision of pH measurement was 0.01. Additionally, pH samples were also transferred to the laboratory and measured at a constant temperature of 25 • C for calibration.

Calculation of Ω arag and TA
Considering the possible effects of organic acids on TA, the calculation of Ω arag avoided using TA.
The Ω arag values at in situ temperature were calculated from the DIC, pH, in situ temperature and in situ salinity values using the CO 2 program [36] and the CO 2 system coefficients of Mehrbach et al. [37] as refitted by Dickson and Millero [38]. The K sp * values for aragonite were taken from Mucci [17], and the Ca 2+ concentrations were assumed to be proportional to salinity, as presented in Millero [39]. Additionally, TA was also obtained in the above calculation process.

Statistical Analyses
The correlations between Ω arag and environmental variables were analysed using Excel (Office 2016, Microsoft, Redmond, WA, USA). All statistical analyses were at significance level of 0.05.

Variations in Seawater Surface Temperature, Salinity and DO%
The seawater surface temperature was between 21.8 • C and 26.2 • C, with an average of 24.0 • C, in JZB during the summer, showing a decreasing trend from the upper end of the bay to the mouth of the bay. The highest temperature value appeared in the western area of the bay with the shallowest water depth. The temperature variation range in winter was 3.8~7.8 • C, and the average value was 5.8 • C. Due to the effects of land cooling, the temperature increased from the upper end of the bay to the mouth area, and the lowest value appeared in the western area. The ranges of seawater surface salinity were 31.0~32.1 in summer and 27.8~31.6 in winter, with average values of 31.8 and 30.9, respectively. The salinity variation gradient in winter was great and may have been related to the strong exchange of seawater between the Yellow Sea and the bay. In terms of distribution, the salinity showed an increasing trend from the northeastern area to the mouth area in the two periods. The salinity was the lowest and the gradient changed drastically in the northeastern area. Most of the rivers that enter JZB have no natural runoff at most of the time of a year and the estuaries have become conduits for wastewater discharge [27]. The total amount of wastewater treated by the three WWTPs located near the northeastern area of JZB is up to~510,000 tons per day ( Figure 1). Therefore, the direct input of treated wastewater was the main reason for the low salinity in this region [11,33]. Additionally, a high-salinity area with a value higher than 30.9 was found in the nearshore area of the western area in summer. Considering the shallow water depth (< 2 m) and the long water residence time (~100 d) in this region, it may be caused by strong evaporation. The range of seawater surface DO% was 81.1%~108.6% in summer, with an average of 98.7%. The DO values were unsaturated in the northeastern area and the area near the estuary of the Dagu River, while those in other areas were supersaturated. The lowest DO% value appeared in the northeastern area of the bay, and the highest value was observed in the southwestern area. The DO% values in winter were higher than those in summer, ranging from 96.5% to 113.8%, with an average of 103.9%. Except for the area near the estuary of the Haibo River, the DO in the entire bay was supersaturated, with the highest value (>110%) appearing in the northeastern area. In addition, the DO% distribution roughly showed a reverse pattern of the summer, and the DO% values decreased from the northeastern area to the mouth area ( Figure 2).

Variations in Seawater Surface Ωarag, DIC and pH
The seawater surface Ωarag was supersaturated in JZB during summer and winter. The Ωarag values in summer ranged from 1.55 to 2.92, with an average of 2.47, and the gradient variation was great. The lowest Ωarag value (< 2.10) was in the northeastern area, while values in most other areas exceeded 2.50, and the highest value (> 2.70) appeared in the nearshore area of the southwestern area. The Ωarag values in winter were between 1.62 and 2.15, with an average of 1.93. The Ωarag gradient variation in this period was obviously more gradual than that in summer, and the Ωarag value in most areas of JZB was approximately 1.90. The surface DIC ranged from 1979 μmol/kg to 2135 μmol/kg, with an average of 2057 μmol/kg, and the value showed a trend of increasing from the western area to the mouth area in summer. In winter, the DIC was in the range of 2027~2249 μmol/kg, with an average of 2122 μmol/kg, and generally increased from the northern area to the mouth area. The highest DIC value appeared in the area near the estuary of the Haibo River. The surface pH value ranged from 7.81 to 8.16, with an average of 8.06, and showed an increasing trend from the northeastern area to the mouth area in summer with high water temperature. In the winter with low water temperatures, the pH level was obviously higher, ranging from 8.12 to 8.37, with an average of

Variations in Seawater Surface Ω arag , DIC and pH
The seawater surface Ω arag was supersaturated in JZB during summer and winter. The Ω arag values in summer ranged from 1.55 to 2.92, with an average of 2.47, and the gradient variation was great. The lowest Ω arag value (< 2.10) was in the northeastern area, while values in most other areas exceeded 2.50, and the highest value (> 2.70) appeared in the nearshore area of the southwestern area.
The Ω arag values in winter were between 1.62 and 2.15, with an average of 1.93. The Ω arag gradient variation in this period was obviously more gradual than that in summer, and the Ω arag value in most areas of JZB was approximately 1.90. The surface DIC ranged from 1979 µmol/kg to 2135 µmol/kg, with an average of 2057 µmol/kg, and the value showed a trend of increasing from the western area to the mouth area in summer. In winter, the DIC was in the range of 2027~2249 µmol/kg, with an average of 2122 µmol/kg, and generally increased from the northern area to the mouth area. The highest DIC value appeared in the area near the estuary of the Haibo River. The surface pH value ranged from 7.81 to 8.16, with an average of 8.06, and showed an increasing trend from the northeastern area to the mouth area in summer with high water temperature. In the winter with low water temperatures, the pH level was obviously higher, ranging from 8.12 to 8.37, with an average of 8.24. Moreover, the pH distribution trend was opposite that in summer and showed a trend of decreasing from the northeastern area to the mouth area ( Figure 3).
Water 2020, 12, x FOR PEER REVIEW 6 of 17 8.24. Moreover, the pH distribution trend was opposite that in summer and showed a trend of decreasing from the northeastern area to the mouth area ( Figure 3).

Qualitative Identification of the Main Factors Controlling the Spatial Distribution of Surface Ωarag
In coastal waters, seawater surface Ω is mainly controlled by temperature, mixing processes (mixing between terrestrial input and seawater) and biological processes [8,20,40,41]. The influence of temperature on Ω is mainly manifested in two aspects. On the one hand, without changing TA and DIC, an increase in seawater temperature can increase the first and second carbonic acid dissociation constants (k1 and k2) of carbonic acid and decrease the Ksp* of calcium carbonate [17,37], thereby promoting an increase in Ω. The difference between highest and lowest seawater temperatures for each investigated cruise were small: 4.4 °C in summer and 4.0 °C in winter. By normalizing the Ωarag of each station at in situ seawater temperatures to the average seawater temperature (24.01 °C for the summer and 5.76 °C for the winter) of the corresponding cruises, it is found that the change of Ωarag before and after the normalization was very small: ≤ 0.04 in summer and ≤ 0.02 in winter. The variation in difference between highest and lowest Ωarag values for each investigated cruise did not exceed 4% after normalization. This indicated that this process had little effect on the spatial

Qualitative Identification of the Main Factors Controlling the Spatial Distribution of Surface Ω arag
In coastal waters, seawater surface Ω is mainly controlled by temperature, mixing processes (mixing between terrestrial input and seawater) and biological processes [8,20,40,41]. The influence of temperature on Ω is mainly manifested in two aspects. On the one hand, without changing TA and DIC, an increase in seawater temperature can increase the first and second carbonic acid dissociation constants (k1 and k2) of carbonic acid and decrease the K sp * of calcium carbonate [17,37], thereby promoting an increase in Ω. The difference between highest and lowest seawater temperatures for each investigated cruise were small: 4.4 • C in summer and 4.0 • C in winter. By normalizing the Ω arag of each station at in situ seawater temperatures to the average seawater temperature (24.01 • C for the summer and 5.76 • C for the winter) of the corresponding cruises, it is found that the change of Ω arag before and after the normalization was very small: ≤ 0.04 in summer and ≤ 0.02 in winter. The variation in difference between highest and lowest Ω arag values for each investigated cruise did not exceed 4% after normalization. This indicated that this process had little effect on the spatial distribution of Ω arag difference. On the other hand, an increase in seawater temperature can reduce the DIC level by reducing the solubility of CO 2 in seawater without changing TA [42], thereby increasing Ω. Since the change in CO 2 solubility in seawater often directly affects the CO 2 exchange process at the sea-air interface, this effect is usually strong [20,41]. However, the Ω arag had no obvious correlation with seawater temperature in JZB during winter and a somewhat negative correlation in summer ( Figure 4a). Obviously, temperature was not the main reason for the spatial distribution differences of Ω arag in these two cruises.
Terrestrial runoff usually has low pH and [CO 3 2− ] values, which directly lead to a decrease in Ω arag in the receiving water body [21,43]. The rivers entering JZB have no natural runoff, but the daily discharges of the Licun River WWTP, the Lushan River WWTP and the Haibo River WWTP, which are located at the three estuaries in the eastern area of the bay, are up to 250,000, 100,000 and 160,000 tons, respectively [11,44]. The direct input of treated wastewater has become the main form of terrestrial input, and unlike natural rivers affected by seasonal changes, treated wastewater input is a persistent process. The pH values of the treated wastewater from the three WWTPs were 7.16~7.33 in January and 7.44~7.50 in June [44]. The wastewater input with low pH values would directly reduce the Ω arag in the northeastern area of JZB. In summer, the Ω arag values in the northeastern area were just significantly lower than those in the central and mouth areas with high salinity values ( Figure 3c). However, the Ω arag in the northeastern area was relatively high in winter, indicating that the effect of treated wastewater input was masked by other factors (Figure 3b).
Water 2020, 12, x FOR PEER REVIEW 7 of 17 distribution of Ωarag difference. On the other hand, an increase in seawater temperature can reduce the DIC level by reducing the solubility of CO2 in seawater without changing TA [42], thereby increasing Ω. Since the change in CO2 solubility in seawater often directly affects the CO2 exchange process at the sea-air interface, this effect is usually strong [20,41]. However, the Ωarag had no obvious correlation with seawater temperature in JZB during winter and a somewhat negative correlation in summer ( Figure 4a). Obviously, temperature was not the main reason for the spatial distribution differences of Ωarag in these two cruises.
Terrestrial runoff usually has low pH and [CO3 2− ] values, which directly lead to a decrease in Ωarag in the receiving water body [21,43]. The rivers entering JZB have no natural runoff, but the daily discharges of the Licun River WWTP, the Lushan River WWTP and the Haibo River WWTP, which are located at the three estuaries in the eastern area of the bay, are up to 250,000, 100,000 and 160,000 tons, respectively [11,44]. The direct input of treated wastewater has become the main form of terrestrial input, and unlike natural rivers affected by seasonal changes, treated wastewater input is a persistent process. The pH values of the treated wastewater from the three WWTPs were 7.16~7.33 in January and 7.44~7.50 in June [44]. The wastewater input with low pH values would directly reduce the Ωarag in the northeastern area of JZB. In summer, the Ωarag values in the northeastern area were just significantly lower than those in the central and mouth areas with high salinity values ( Figure 3c). However, the Ωarag in the northeastern area was relatively high in winter, indicating that the effect of treated wastewater input was masked by other factors (Figure 3b). Biological processes mainly include primary production and aerobic respiration. Primary production produces oxygen and consumes CO2, thus increasing Ωarag [45]. Conversely, aerobic respiration consumes oxygen and reduces Ωarag. During the summer, Ωarag had an excellent positive correlation with DO% ( Figure 4c) (r = 0.94, p < 0.0001), suggesting that biological processes were the Biological processes mainly include primary production and aerobic respiration. Primary production produces oxygen and consumes CO 2 , thus increasing Ω arag [45]. Conversely, aerobic respiration consumes oxygen and reduces Ω arag . During the summer, Ω arag had an excellent positive correlation with DO% ( Figure 4c) (r = 0.94, p < 0.0001), suggesting that biological processes were the main factor affecting the distribution of Ω arag in JZB and might explain the main differences in the spatial distribution of Ω arag . According to the unsaturated DO levels, the dominance of aerobic respiration resulted in the lowest Ω arag in the northeastern area of JZB. However, in most other areas of JZB, the DO was supersaturated, and primary production was dominant, which led to higher Ω arag levels. During the winter, the oversaturated DO values in most areas of the entire bay meant that primary production was dominant, but the relationship between DO% and Ω arag only showed a relatively weaker negative correlation (Figure 4c) (r = 0.53, p < 0.01), indicating that the control of biological processes on the spatial distribution of Ω arag in winter was weaker than that in summer. Despite this, considering that Ω arag had no significant correlations with seawater temperature and salinity (Figure 4a,b), primary production may still be the main factor affecting the Ω arag distribution in winter, and the weak negative correlation between DO% and Ω arag may be the result of the effects of the mixing processes and regional temperature variation. For example, treated wastewater input with low pH values could decrease the Ω arag levels in the northeastern area where the DO% (> 107%) values were highest, the low seawater temperature could cause a decrease in Ω arag in the western area where the DO% values were approximately 105%, and the strong external seawater exchange could cause an increase in Ω arag in the mouth area where the DO% values were near 100% (Figure 2b,f and Figure 3b).
The regression analysis indicates that biological processes may be the main factor affecting the spatial difference of Ω arag in JZB during summer and winter, while the influence of temperature and mixing processes is relatively smaller. However, the specific contributions of the above three processes could not be calculated if only using the regression analysis, and the correlation between parameters may even obscure the important role of certain processes in some cases. For example, direct input of the treated wastewater had the most significant impact in the northeastern area with the lowest DO% values and the strongest degradation of organic matter during summer, and the excellent positive correlation between Ω arag and DO% (Figure 4c) partly included the contribution of the treated wastewater input. Therefore, a further quantification analysis involving the contributions of the above three processes to the spatial distribution of Ω arag is still necessary.

Processes Controlling the Spatial Distribution of Surface Ω arag Revealed by Quantitative Evaluation
In the carbonate system, the change in Ω arag with salinity is nonlinear and non-conservative, so Ω arag is not convenient for quantitative calculations. The Ω arag in seawater is essentially determined by the [CO 3 2− ] value, and the difference between TA and DIC ([TA-DIC]) can reflect the levels of [CO 3 2− ] well, thereby acting as a good proxy for Ω arag [20,46]. As shown in Figure 5a For the contribution of temperature variation, only the effect of temperature on the solubility of CO2 in seawater was evaluated in these two cruises, due to the negligible Ωarag changes caused by the variation of k1 and k2 of carbonic acid and Ksp* of calcium carbonate from temperature (Section 3.3). In the calculation process, taking the average TA and salinity of the cruise as the staring values (the average TA and salinity were 2260.1 μmol/kg and 31.76 in summer and 2291.4 μmol/kg and 30.93 in winter, respectively) and assuming seawater CO2 reached equilibrium with the atmosphere at in situ temperature of each station, this research calculated the corresponding DIC (DICeq) of each station at different temperatures and the DIC (DICavg) at the average temperature of the corresponding cruise at this equilibrium state using the CO2SYS program [36]. Then, the difference between DICavg and DICeq was the contribution from the temperature of the corresponding station to [TA-DIC]A ([TA-DIC]A-tem) (Equation (1)).  For the contribution of temperature variation, only the effect of temperature on the solubility of CO 2 in seawater was evaluated in these two cruises, due to the negligible Ω arag changes caused by the variation of k1 and k2 of carbonic acid and K sp * of calcium carbonate from temperature (Section 3.3). In the calculation process, taking the average TA and salinity of the cruise as the staring values (the average TA and salinity were 2260.1 µmol/kg and 31.76 in summer and 2291.4 µmol/kg and 30.93 in winter, respectively) and assuming seawater CO 2 reached equilibrium with the atmosphere at in situ temperature of each station, this research calculated the corresponding DIC (DIC eq ) of each station at different temperatures and the DIC (DIC avg ) at the average temperature of the corresponding cruise at this equilibrium state using the CO2SYS program [36]. Then, the difference between DIC avg and DIC eq was the contribution from the temperature of the corresponding station to [TA-DIC] A ([TA-DIC] A-tem ) (Equation (1)).
[TA-DIC] A-tem = DIC avg − DIC eq (1) where the DIC eq and DIC avg are the DIC value at in situ temperature of the station and at the average temperature of the corresponding cruise, respectively, when seawater CO 2 reached equilibrium with the atmosphere. Atmospheric CO 2 concentration data from flask measurements in the Tae  Meanwhile, compared with seawater, the treated wastewater had different temperatures, indicating that mixing term may include some bias from temperature changes. Temperature mixes linearly, so this component can be estimated if the proportion of treated wastewater in seawater can be obtained. According to mass balance, the sample salinity of a station can be calculated as follows [11]: where y is the proportion of wastewater in seawater and S i is the salinity of the station i. In January and June, the lowest salinity value in the eastern area of the bay was 27.8 and 31.0, respectively. According to the Equation (3), the y values were calculated as 0.13 and 0.03, respectively. The temperature difference between treated wastewater and seawater in the eastern area did not exceed 7.6 • C in January and 4.5 • C in June, indicating that the temperature changes caused by this mixing process were ≤ 1.0 • C in January and ≤ 0.14 • C in June. Therefore, the effect of this process on Ω arag could be ignored in the current study.
The [TA-DIC] A value after the deduction of the effects from temperature and mixing processes was the contribution of other processes ([TA-DIC] A-res ) (Equation (4)). The relative contributions of each process to [TA-DIC] A in summer are shown in Figure 6. The relative contribution of temperature to [TA-DIC] A was the largest in the northwestern area of the bay with the highest seawater temperature, exceeding 45%. However, in other areas, it was smaller, less than 30%. The mixing processes had the minimum relative contribution, and most of the regions featured values of less than 10%; however, the contribution reached 13% in the northeastern area, which had the strongest impact of treated wastewater input. The biological processes had the maximum relative contribution to [TA-DIC] A , which was more than 50% in most regions. In the northeastern area with the lowest DO% (<90%) and the southwestern area with the highest DO% (>105%), the relative contribution of biological processes exceeded 85%. Overall, the average contributions of temperature, mixing processes and biological processes to [TA-DIC] A were 25%, 5% and 70%, respectively. The biological processes were the main factor controlling the spatial distribution of Ω arag in JZB during summer. It is worth noting that the effects of temperature and mixing processes on [TA-DIC] A were mutually reinforced with biological processes in most areas in this period. In the northeastern area with the minimum value of Ω arag , the main effect was that treated wastewater input and aerobic respiration together caused the reduction of [TA-DIC] A, and the contributions of the two were up to 90%. In the southwestern area with the maximum value of Ω arag , primary production and the higher seawater temperature lead to an increase in [TA-DIC] A , and the combined contributions of the two reached approximately 90%. As such, there were marked difference in Ω arag among the various regions of JZB in summer.
Water 2020, 12, x FOR PEER REVIEW 11 of 17 contributions of temperature, mixing processes and biological processes to [TA-DIC]A were 25%, 5% and 70%, respectively. The biological processes were the main factor controlling the spatial distribution of Ωarag in JZB during summer. It is worth noting that the effects of temperature and mixing processes on [TA-DIC]A were mutually reinforced with biological processes in most areas in this period. In the northeastern area with the minimum value of Ωarag, the main effect was that treated wastewater input and aerobic respiration together caused the reduction of [TA-DIC]A, and the contributions of the two were up to 90%. In the southwestern area with the maximum value of Ωarag, primary production and the higher seawater temperature lead to an increase in [TA-DIC]A, and the combined contributions of the two reached approximately 90%. As such, there were marked difference in Ωarag among the various regions of JZB in summer. In winter (Figure 7), the relative contribution of temperature to [TA-DIC]A was largest and exceeded 45% in the western area where the seawater temperature was lowest; in the other regions, the contribution was less, below 25%. The relative contribution of the mixing processes to [TA-DIC]A was much higher in winter than in summer and was more than 30% in the northeastern area with the lowest salinity (< 30.0) and the mouth area with the highest salinity (> 31.5). This may be related to enhanced exchange of seawater between the bay and the Yellow Sea caused by the low seawater temperatures and high wind speeds in winter. The biological processes contributed the most to [TA-DIC]A, and the relative contribution was more than 60% in the northern area and the area near the estuary of the Haibo River. During this period, DO values in most areas of the entire bay were supersaturated, and primary production was dominant, which corresponds to the growth peak of phytoplankton in JZB in winter [25][26][27]. Overall, the average relative contributions of temperature, mixing processes and biological processes to [TA-DIC]A were 20%, 30% and 50%, respectively. The biological processes were also the main processes contributing to spatial differences in Ωarag in winter in JZB, but the total contribution of temperature and mixing processes was similar to that of the biological processes in this period. Furthermore, different from the situation in summer, the effects of the above three processes on [TA-DIC]A showed mutual cancellation in most areas. In the In winter (Figure 7), the relative contribution of temperature to [TA-DIC] A was largest and exceeded 45% in the western area where the seawater temperature was lowest; in the other regions, the contribution was less, below 25%. The relative contribution of the mixing processes to [TA-DIC] A was much higher in winter than in summer and was more than 30% in the northeastern area with the lowest salinity (<30.0) and the mouth area with the highest salinity (>31.5). This may be related to enhanced exchange of seawater between the bay and the Yellow Sea caused by the low seawater temperatures and high wind speeds in winter. The biological processes contributed the most to [TA-DIC] A , and the relative contribution was more than 60% in the northern area and the area near the estuary of the Haibo River. During this period, DO values in most areas of the entire bay were supersaturated, and primary production was dominant, which corresponds to the growth peak of phytoplankton in JZB in winter [25][26][27]. Overall, the average relative contributions of temperature, mixing processes and biological processes to [TA-DIC] A were 20%, 30% and 50%, respectively. The biological processes were also the main processes contributing to spatial differences in Ω arag in winter in JZB, but the total contribution of temperature and mixing processes was similar to that of the biological processes in this period. Furthermore, different from the situation in summer, the effects of the above three processes on [TA-DIC] A showed mutual cancellation in most areas. In the northeastern area where DO% values exceeded 109%, the low seawater temperature and treated wastewater input had a negative effect on [TA-DIC] A , while the strong primary production had a positive effect on [TA-DIC] A , resulting in a net increase in [TA-DIC] A of only approximately 10%. In the western area, the relative contribution of the mixing processes and primary production, corresponding to an increase in [TA-DIC] A , was approximately 60%, whereas the remaining 40% contribution corresponded to a decrease in [TA-DIC] A caused by low seawater temperatures. Therefore, the differences in Ω arag among the various regions in JZB were markedly reduced in winter.  In the open ocean, physical transport in the horizontal and vertical directions is a key factor causing differences in the distribution of Ω [47][48][49], while the increasing human disturbance has affected the biological processes more prominently in coastal waters. In ocean shelf regions, the effect of biological processes on Ω mainly appears in the summer. Due to the stratification of seawater, the surface water body often features strong primary production, which increases Ω under the appropriate illumination and nutrient supply conditions, while organic matter deposited from the upper water body often leads to strong aerobic respiration in the subsurface water body, resulting in low Ω and pH [8,50,51]. However, in nearshore bays where the impact of human activities is more significant, the situation is slightly different. The urbanization process affects biological processes in certain areas. For example, in the JZB, the eastern area of the bay, which is adjacent to downtown Qingdao, receives large amounts of terrestrial pollutants brought by wastewater all the year round and is dominated by aerobic respiration. This situation can weaken the living environment of calcified organisms in the region and does not change even after heavy rain [33,35]. In contrast, in the western area, which is far from the urban area and adjacent to the wetland, primary production produces an increase in Ω, and this region is an important area for shellfish farming in JZB [52]. A similar phenomenon was also observed in Fujian, Southeast China [53,54] and Guanabara Bay in Brazil [25]. Moreover, perennial eutrophication often causes algal blooms and produces strong primary production even under the low temperatures in winter in some waterbodies [55,56], offsetting the decrease in Ω caused by the low temperatures. In addition, bays often have a shallow In the open ocean, physical transport in the horizontal and vertical directions is a key factor causing differences in the distribution of Ω [47][48][49], while the increasing human disturbance has affected the biological processes more prominently in coastal waters. In ocean shelf regions, the effect of biological processes on Ω mainly appears in the summer. Due to the stratification of seawater, the surface water body often features strong primary production, which increases Ω under the appropriate illumination and nutrient supply conditions, while organic matter deposited from the upper water body often leads to strong aerobic respiration in the subsurface water body, resulting in low Ω and pH [8,50,51]. However, in nearshore bays where the impact of human activities is more significant, the situation is slightly different. The urbanization process affects biological processes in certain areas. For example, in the JZB, the eastern area of the bay, which is adjacent to downtown Qingdao, receives large amounts of terrestrial pollutants brought by wastewater all the year round and is dominated by aerobic respiration. This situation can weaken the living environment of calcified organisms in the region and does not change even after heavy rain [33,35]. In contrast, in the western area, which is far from the urban area and adjacent to the wetland, primary production produces an increase in Ω, and this region is an important area for shellfish farming in JZB [52]. A similar phenomenon was also observed in Fujian, Southeast China [53,54] and Guanabara Bay in Brazil [25]. Moreover, perennial eutrophication often causes algal blooms and produces strong primary production even under the low temperatures in winter in some waterbodies [55,56], offsetting the decrease in Ω caused by the low temperatures. In addition, bays often have a shallow water depth, and strong turbulence might cause full mixing of the seawater vertically, making the effect of biological activity on Ω extend through the entire water column [23]. With the rapid urbanization and worsening climate change in coastal areas, the effects of strong heterotrophic bacteria activity and vigorous plankton growth on Ω will become increasingly prominent. Thankfully, environmental protection has been received increasing attention in China. Central and local governments have accelerated efforts to tackle environmental pollution by raising environmental protection standards and revising laws [57,58]. However, strict enforcement of environmental protection laws is urgently needed to conserve the aquatic environment, particularly the effective treatment of wastewater which brought large amounts of organic matters and nutrients to coastal environment [59]. In the future, collaboration between government, industry and environmental non-governmental organizations (NGO) will be of great importance to mitigate ocean acidification and other problems in coastal areas.

Uncertainty Analysis of the Quantitative Estimate
Here the uncertainty associated with key terms in the quantitative evaluation method in Section 3.4 was analysed. First, the uncertainty of TA and Ω arag calculated from the pH and DIC were estimated to be~7 µmol/kg and~0.05, respectively, via the CO2SYS program (without considering the uncertainty induced by carbonate acid dissociation constants during calculations) [60], given the uncertainty of 2 µmol/kg for DIC and the uncertainty of 0.01 for pH measurement. The TA uncertainty will lead to an uncertainty of~7.3 µmol/kg, or (7 2 + 2 2 ) 0.5 , for [TA-DIC]. This uncertainty will propagate into [TA-DIC] A (~7.4 µmol/kg, or (7.3 2 + 7.3 2 /28) 0.5 ), Xi-tem and Xi-mix (< 5.5%) and Xi-res (< 8%). Second, [TA-DIC] has some limitations as a proxy for in situ Ω arag . In the current study, [TA-DIC] may not accurately reflect the influence of temperature. The effect of temperature only refers to its impact on thermodynamic constants or internal allocation of carbonate system species by temperature driven acid-base equilibrium under conditions of constant TA and DIC, which is called "internal temperature effect" [20]. For instance, for one unit increase in temperature, Ω arag will increase by~0.60% • C −1 based on our data. Fortunately, the thermal influence is relatively minor due to the small spatial variation in seawater temperature for each investigated cruise. Meanwhile, to calculate the effects of temperature on the solubility of CO 2 in seawater, it is assumed that seawater pCO 2 can reach an equilibrium state with the atmospheric CO 2 quickly. In fact, it is almost impossible to reach transient air-sea CO 2 equilibrium due to the slow air-sea CO 2 exchange [61], which will lead to overestimation or underestimation of DICeq. However, this influence on [TA-DIC] A-tem may be minor since this overestimated or underestimated part offset, to a large extent, each other after subtraction as shown from Equation (1). Third, to separate mixing effects, the slope values of [TA-DIC]-Salinity between the composite treated wastewater and seawater in January and June were estimated from end-members values. Treated wastewater end-member values were calculated by the [TA-DIC] values and discharge of three WWTPs in corresponding months which were reported by Liu et al. [44]. Though their studies showed that changes in DIC and TA were generally same in the two years, the uncertainties brought by the treated wastewater end-member in our cruises still exist in [TA-DIC] A-mix . Despite some uncertainties in our methods which are also widely applied in other studies, the current study provides some quantitative findings on the control mechanism of Ω arag in the bays under the influence of human activities.

Limitation and Future Research
Similar to many studies, some limitations exist in the current research. First, Ω arag and TA were calculated based on pH and DIC. The uncertainty of pH measurements made with electrodes calibrated on the NBS scale may cause a certain level of inaccuracy [62,63]. Thankfully, the spatiotemporal changes of pH in JZB were large enough, and the conclusions were valid based on these data. Second, this study focused the contributions of three factors (temperature, biological processes and mixing processes) to the distribution of Ω arag . However, the effects of other factors, for example, mineral precipitation, vertical mixing of water column, nitrification processes, may also exist, although their contributions may be minor in our research area. Additionally, the effects of urbanization on Ω arag include sewage input and other factors. Third, the seasonal changes of Ω arag is unavailable due to the data limitation. Obviously, more cruise investigations and related studies are still needed to improve our understanding of the complex ocean acidification under the impact of climate change and human activities.

Conclusions
Seawater surface aragonite saturation state (Ω arag ) during summer and winter in Jiaozhou Bay, China, was analysed in this study. The seawater surface Ω arag was supersaturated, and the ranges were 1.55~2.92 in summer and 1.62~2.15 in winter, respectively. In addition to the analysis of the effects of temperature, mixing processes and biological processes on the spatial difference of Ω arag based on regression analysis with surface temperature, salinity and DO%, the relative contributions of these three processes were quantified by introducing the difference between TA and DIC as a proxy for Ω arag . The results showed that the biological processes were the main processes affecting the spatial difference of Ω arag , and the relative contribution was 70% in summer and 50% in winter. The contributions of temperature (25% in summer and 20% in winter) and mixing processes (5% in summer and 30% in winter) were smaller. Moreover, the interaction among the three factors that influence Ω arag exhibited mutual promotion in summer and mutual cancellation in winter, causing the regional distribution of Ω arag to show marked difference in summer and a relatively uniform change in winter.

Conflicts of Interest:
The authors declare no conflict of interest.