Interannual Carbon and Nutrient Fluxes in Southeastern Taiwan Strait

: The Taiwan Strait (TS) is one of the main sources of phosphate that supports the large ﬁsh catches of the phosphate-limited East China Sea (ECS). The Penghu Channel is the deepest part of the TS, and most of the ﬂow of the TS towards the ECS is principally through this channel. Empirical equations that are based on measurements made during 19 cruises (2000–2011) were combined with water velocity, salinity, and temperature, which were modeled using HYCOM (the Hybrid Coordinate Ocean Model) to obtain the annual ﬂuxes for total alkalinity (TA), dissolved inorganic carbon (DIC), nitrate plus nitrite, phosphate, and silicate ﬂuxes. The TA and DIC are mainly transported in the top layer (0–55 m) because the current is much stronger there than in the bottom layer (55–125 m) whereas the TA and DIC concentrations in the top layer are only slightly smaller compared with the bottom layer. In contrast, the nitrate plus nitrite ﬂux is mainly transported in the bottom layer because the concentrations are much higher in the bottom layer. Generally, nutrient ﬂux increases with the Paciﬁc Decadal Oscillation (PDO) index, but TA and DIC ﬂuxes increase as the PDO index decreases.


Introduction
The Taiwan Strait (TS) directly connects the South China Sea (SCS) and the East China Sea (ECS), which are two bountiful marginal seas. The three major water masses in the TS are the Minzhe Coastal Water, the surface water of the SCS, and the surface water of the West Philippines Sea (WPS). Generally, the Minzhe Coastal Water has the highest nutrient concentrations among these three water masses, and the surface water of the WPS, which is brought into the TS through Kuroshio intruding the SCS, contains the lowest nutrient concentrations. The TS transports a mixture of SCS water and WPS water to the ECS mostly in the summer, while the Minzhe Coastal Water flows southward along the western side of the TS to the northern SCS in the winter. Through the exchange of water masses, the TS transports nutrients between these two seas [1,2].
The nutrient budgets in the ECS are influenced by fluvial transport, Kuroshio Intermediate Water, and TS seawater [3,4]. The nitrate and phosphate fluxes in the Yangtze River have increased dramatically with the population and consumption of fertilizer in the basin, but the silicon flux is reduced by sediment trapping in reservoirs [5]. The ecosystem in the Yangtze estuary has changed as the number of non-siliceous algae and the red tides of harmful algal blooms have increased because of the varied nutrient concentration ratios [6,7]. The biological uptake ratio between nitrogen and phosphate is around 16:1 in the ocean, and this is called the Redfile ratio.

Materials and Methods
To estimate continuous nutrient and inorganic carbon fluxes, the chemical concentrations were simulated based on the measured data. The daily salinity, temperature, and flow velocity in the PHC were obtained using the Hybrid Coordinate Ocean Model (HYCOM). The sampling sites are located between 22.5-23.5 • N and 119.5-120.1 • E (Figure 1), and most of them are close to 23 • N. A total of 555 bottle samples were collected during 19 cruises from 2000 to 2011 using a CTD/Rosette sampler, while salinity and temperature were recorded in situ. The nutrient concentrations were determined by published methods with flow injection analyzers. Nitrate plus nitrite (N) are the major species of dissolved inorganic nitrogen in the seawater, and the N concentration was obtained using the pink azo dye method [22], with a precision of approximately ±1% at 35 µmol kg −1 and ±3% at 1 µmol kg −1 . Phosphate (P) concentration was measured using the molybdenum blue method [23] with a precision of approximately ±0.5% at 2.5 µmol kg −1 and ±3% at 0.1 µmol kg −1 . Silicate (Si) concentration was measured using the silicon molybdenum blue method with a precision of around 0.6% at 150 µmol kg −1 and 2% at 5 µmol kg −1 . DIC concentration was determined using 10% phosphoric acid to acidify the samples, and then quantifying the produced CO 2 gas using an infrared gas analyzer (AS-C3 Apollo Scitech), a single operator multi-parameter metabolic analyzer (SOMMA), or a coulometric detector. TA was measured using the open-cell method of potentiometric titration at 25 • C with a PC-controlled titration cell [24,25]. The end points of titration were determined using the Gran Function with a precision of 0.1% [26]. To estimate the daily average chemical concentrations, simulation models were derived using in situ temperature, salinity, and measured chemical concentrations. A second-order polynomial regression model was selected using leave-one-out cross-validation because it minimized the mean square error of prediction [27]. (See Appendix A for more detail.) Table 1 presents the coefficients of regression fits, residual standard error, and adjusted coefficients of determination. Each adjusted determination coefficient between the seasonal chemical simulation result and the measured data exceeded 0.7, so the model captured more than 70% of the variability of the original data. The low probability values (p < 0.001) suggest that simulation models are compatible with the measured data (Table 1). To compare the measured and fitted results, the mentioned equations (Table 1) and the in situ temperature and salinity values that were collected with bottle samples were utilized to calculate the fitted chemical concentrations. Significant positive correlations existed between the approximately 500 measured and fitted values of TA, DIC, N, P, and Si concentrations, and the simulated performance values expressed more than 88% of measured data, according to the determination coefficients ( Figure 2).  To estimate the daily average chemical concentrations, simulation models were derived using in situ temperature, salinity, and measured chemical concentrations. A second-order polynomial regression model was selected using leave-one-out cross-validation because it minimized the mean square error of prediction [27]. (See Appendix A for more detail.) Table 1 presents the coefficients of regression fits, residual standard error, and adjusted coefficients of determination. Each adjusted determination coefficient between the seasonal chemical simulation result and the measured data exceeded 0.7, so the model captured more than 70% of the variability of the original data. The low probability values (p < 0.001) suggest that simulation models are compatible with the measured data (Table 1). To compare the measured and fitted results, the mentioned equations (Table 1) and the in situ temperature and salinity values that were collected with bottle samples were utilized to calculate the fitted chemical concentrations. Significant positive correlations existed between the approximately 500 measured and fitted values of TA, DIC, N, P, and Si concentrations, and the simulated performance values expressed more than 88% of measured data, according to the determination coefficients ( Figure 2).   Table 1 at the temperature and salinity of bottle samples.
To quantify the contribution from each predictor variable to simulation results, the R-squared values in different situations were calculated. The temperature is more important than the salinity as a predictor variable in the chemical concentration models except for TA fitted model during winter and spring (Table 2). Comparing in situ temperature and the daily temperature that was simulated using the HYCOM demonstrates that the simulation yielded underestimates in the surface layer (0-55 m) but overestimates in the bottom layer (55-125 m, Figure 3a). On the contrary, the HYCOM salinity and fitted N values were overestimated in the surface layer but were underestimated in the bottom (Figure 3b,c). Since temperature is the major controlling factor in the chemical model, the differences between the in situ temperature and the HYCOM model temperature caused simulated N concentrations to vary. There were larger variations in the bottom layers of temperature and N differences than in the surface layers, but the variation of salinity difference was lower in the bottom layer than in the surface layer ( Figure 3). The results also support the conclusion that the simulated N concentrations are less influenced by salinity than by temperature. The differences between simulated and measured data may have had the following few causes; (1) the variation between each single measurement and the daily average; (2) the assumption that biological production and consumption do not vary within a season; (3) the simplification of chemical concentrations using a fixed water mass mixing ratio between the SCS and WPS waters; and (4) the  Table 1 at the temperature and salinity of bottle samples.
To quantify the contribution from each predictor variable to simulation results, the R-squared values in different situations were calculated. The temperature is more important than the salinity as a predictor variable in the chemical concentration models except for TA fitted model during winter and spring (Table 2). Comparing in situ temperature and the daily temperature that was simulated using the HYCOM demonstrates that the simulation yielded underestimates in the surface layer (0-55 m) but overestimates in the bottom layer (55-125 m, Figure 3a). On the contrary, the HYCOM salinity and fitted N values were overestimated in the surface layer but were underestimated in the bottom (Figure 3b,c). Since temperature is the major controlling factor in the chemical model, the differences between the in situ temperature and the HYCOM model temperature caused simulated N concentrations to vary. There were larger variations in the bottom layers of temperature and N differences than in the surface layers, but the variation of salinity difference was lower in the bottom layer than in the surface layer ( Figure 3). The results also support the conclusion that the simulated N concentrations are less influenced by salinity than by temperature. The differences between simulated and measured data may have had the following few causes; (1) the variation between each single measurement and the daily average; (2) the assumption that biological production and consumption do not vary within a season; (3) the simplification of chemical concentrations using a fixed water mass mixing ratio between the SCS and WPS waters; and (4) the deviation of HYCOM results.  To estimate the daily average concentration of each chemical parameter, the double integral mean value theorem was used. By Chebyshev's inequality [28], the true values of daily salinity ( ) and temperature ( ) are 75% likely to be within two deviations of the daily average of salinity ( ̅ − 2 , ̅ + 2 ) and two deviations of the daily average temperature ( ̅ − 2 , ̅ + 2 ) if and are close to the true standard deviations. Let ̅ be the mean value of ( , ) over the above intervals under the model. According to the mean value theorem for a double integral [29], With properly defined and , ̅ can be regarded as an approximation to the true daily To estimate the daily average concentration of each chemical parameter, the double integral mean value theorem was used. By Chebyshev's inequality [28], the true values of daily salinity (x s ) and temperature (x t ) are 75% likely to be within two deviations of the daily average of salinity (x s − 2d s , x s + 2d s ) and two deviations of the daily average temperature (x t − 2d t , x t + 2d t ) if d s and d t are close to the true standard deviations. Let f be the mean value of f (x s , x t ) over the above intervals under the model. According to the mean value theorem for a double integral [29], With properly defined d s and d t , f can be regarded as an approximation to the true daily average of a chemical parameter in the model. Appendix A presents details of the model selection and the fitted daily average values of chemical parameters.
The daily salinity, temperature, and water velocity along 23.04 • N from 119.44 to 120.08 • E were obtained using HYCOM + NCODA (Navy Coupled Ocean Data Assimilation) Global 1/12 • Reanalysis from January 1993 to December 2012. Water flux was calculated using the daily model water velocity (v; m s −1 ), the depth layer (d i z; m), and the distance between two adjacent stations (dx; m): where i and n are the index of a station and the total number of stations, respectively; k is the layer number at the ith station; k i is the bottom layer at the ith station. The chemical fluxes were generated as: where C(x, z) was calculated from the aforementioned second-order polynomial regression models and the double integral mean value. The unit of the chemical flux estimated at the study section is mole s −1 . The positive value represents a northward transport chemical quantity through the profile, and negative value is a southward transport.

Seasonal Profiles
Since the PHC is directly and indirectly influenced by seasonal monsoons, the seawater mixing ratio and water transport velocity vary over time. The model salinity increased with depth from 33.9 to 34.8, but the temperature decreased with depth from 29 • C to 18 • C (Figure 4a-d). The salinity and temperature in winter and spring varied little between the surface and bottom layers, whereas those in summer and autumn varied more (Figure 4a,b). The differences in salinity and temperature between the surface and bottom layers were largest in the summer, and the highest and lowest salinities and temperatures were also found in this season (Figure 4c,d). The distribution of the P concentration followed the salinity distribution, and the P concentration is the highest in the bottom layer in summer. However, the distribution of the velocity differs from those of the aforementioned parameters. In autumn and winter, the velocity is highest in the middle layer of the deepest part of the PHC around 119.8 • E (Figure 4e,h), but in spring and summer, it is highest in the surface layer (Figure 4f,g).
The PHC is divided into top and bottom layers (red and blue lines in Figure 5
The PHC is divided into top and bottom layers (red and blue lines in Figure 5

Comparison between Measured and Simulated Data
To compare the data measured in situ in 19 cruises with the simulated data, the weighted-average gridding method in Ocean Data View [33] was used to interpolate and extrapolate data associated with a single profile along 23.04 • N from 119.4 to 120.08 • E. The average differences (simulated-measured in 19 cruises) in salinity (∆S) and temperature (∆T) are 0.063 ± 0.237 and 0.246 ± 1.086 • C, respectively. For TA, DIC, N, P, and Si fluxes, the average differences are 0.001 ± 0.017 × 10 6 mol C s −1 (∆TA), −0.010 ± 0.033 × 10 6 mol C s −1 (∆DIC), −409 ± 1064 mol N s −1 (∆N), −22 ± 68 mol P s −1 (∆P), and −588 ± 1664 mol Si s −1 (∆Si), respectively. Generally, the simulated data yield overestimates of salinity, temperature, and TA flux, but underestimates of DIC, N, P, and Si fluxes. The ∆TA and ∆DIC fluxes positively correlate with ∆S, but the differences ∆N, ∆P, and ∆Si fluxes are negatively correlated with ∆T ( Figure 6). The daily chemical concentration was calculated by averaging the maxima and minima of the estimated concentrations that were mentioned in Section 2. yield overestimates of salinity, temperature, and TA flux, but underestimates of DIC, N, P, and Si fluxes. The ΔTA and ΔDIC fluxes positively correlate with ΔS, but the differences ΔN, ΔP, and ΔSi fluxes are negatively correlated with ΔT ( Figure 6). The daily chemical concentration was calculated by averaging the maxima and minima of the estimated concentrations that were mentioned in Section 2.

Annual Variation
The two major factors that drive water transportation in the TS are wind stress and a northward pressure gradient force that arises from the fact that the surface of the south TS is higher than that of the north TS. Southwest monsoons accelerate northward flow during the summer, but northeast monsoons reduce it, even reversing it to southward flow during the period of strong northeasterly winds [10,15,34,35].
The monthly values of fluxes and parameters exhibit four patterns. First, the water, TA, and DIC fluxes (Figure 7a,d,e) are highest in July. The TA and DIC fluxes basically reflect the variation of water flux, since the seasonal percentage differences of TA and DIC concentrations are smaller than the variations of the nutrient concentrations ( Figure 7). Second, the salinity, TA, and DIC concentrations are lowest in September (Figure 7b,d,e). The third pattern contrasts with the second, as the temperature is highest in September (Figure 7c). During the summer and autumn, the main water mass in the PHC is the SCS water, which has a lower salinity than the WPS water at the same density in surface water. Since precipitation in the summer is heavy, the water in the top layer is warmer and fresher during the summer than the autumn. Yet, the strong summer water flow carries the colder and saltier deeper water into the PHC (Figure 4b,f). The northward current is weakened in autumn owing to the onset of the northeast monsoon, and reduces the coldness and saltiness of the seawater in the bottom. The area-weighted salinity is lowest in September, when the area-weighted temperature is highest. This phenomenon is consistent with the significant influence of the salty bottom water in the summer. Fourth, the N, P, and Si fluxes are highest in June (Figure 7f,g,h). The high water flux increases nutrient fluxes in the summer as the deeper water contains higher nutrient concentrations.

Annual Variation
The two major factors that drive water transportation in the TS are wind stress and a northward pressure gradient force that arises from the fact that the surface of the south TS is higher than that of the north TS. Southwest monsoons accelerate northward flow during the summer, but northeast monsoons reduce it, even reversing it to southward flow during the period of strong northeasterly winds [10,15,34,35].
The monthly values of fluxes and parameters exhibit four patterns. First, the water, TA, and DIC fluxes (Figure 7a,d,e) are highest in July. The TA and DIC fluxes basically reflect the variation of water flux, since the seasonal percentage differences of TA and DIC concentrations are smaller than the variations of the nutrient concentrations ( Figure 7). Second, the salinity, TA, and DIC concentrations are lowest in September (Figure 7b,d,e). The third pattern contrasts with the second, as the temperature is highest in September (Figure 7c). During the summer and autumn, the main water mass in the PHC is the SCS water, which has a lower salinity than the WPS water at the same density in surface water. Since precipitation in the summer is heavy, the water in the top layer is warmer and fresher during the summer than the autumn. Yet, the strong summer water flow carries the colder and saltier deeper water into the PHC (Figure 4b,f). The northward current is weakened in autumn owing to the onset of the northeast monsoon, and reduces the coldness and saltiness of the seawater in the bottom. The area-weighted salinity is lowest in September, when the area-weighted temperature is highest. This phenomenon is consistent with the significant influence of the salty bottom water in the summer. Fourth, the N, P, and Si fluxes are highest in June (Figure 7f,g,h). The high water flux increases nutrient fluxes in the summer as the deeper water contains higher nutrient concentrations.

Interannual Variations
To determine simplified patterns of the anomalous flux of TA, DIC, N, P, and Si, 30-month moving averages were used to eliminate any seasonal pattern (Figure 8

Interannual Variations
To determine simplified patterns of the anomalous flux of TA, DIC, N, P, and Si, 30-month moving averages were used to eliminate any seasonal pattern (Figure 8  Generally, the variation of the 30-month moving average PDO index tends to oppose that of water volume in the top layer, but the values in both layers increase slightly with time (Figure 8b). For the top layer, the increased PDO index is related to the reduced water flux and vice versa. The anomalous salinity patterns in both layers mirror the pattern of the PDO index, and decrease slightly over time (Figure 8b). The salinity abatement may be an indirect effect of the reduced Kuroshio Generally, the variation of the 30-month moving average PDO index tends to oppose that of water volume in the top layer, but the values in both layers increase slightly with time ( Figure 8b). For the top layer, the increased PDO index is related to the reduced water flux and vice versa. The anomalous salinity patterns in both layers mirror the pattern of the PDO index, and decrease slightly over time (Figure 8b). The salinity abatement may be an indirect effect of the reduced Kuroshio intrusion, and a simulated low salinity in 2001 has been reported [36]. If the decreased salinity is caused by the reduced Kuroshio intrusion, then the temperature will decrease with the low salinity event in 2001. However, the low salinity event in 2008 is associated with a high temperature signal, suggesting that Kuroshio intrusion may not be the only factor that controls water mass mixing in the PHC (Figure 8c,d). The anomalous TA and DIC patterns are similar to the anomalous water flux, and also increase with time (Figure 8e,f). Generally, the TA and DIC fluxes rise as the PDO index decreases. However, the nutrient fluxes vary with the PDO index, indicating that the nutrient concentrations are the main controlling factor. The anomalous nutrient fluxes also increase slightly with time (Figure 8g-i).

Conclusions
The simulated chemical fluxes of TA, DIC, N, P, and Si are calculated from estimated chemical concentrations using individual second-order polynomial regression equations and water flux that is obtained using the HYCOM. Equations for TA, DIC, N, P, and Si concentrations were derived using chemical concentrations that were measured from bottle samples, in situ salinity, and temperature. All the chemical models presented herein explained more than 70% of the variability of the original data. The individual daily chemical concentrations in various locations and at different water depths were calculated using the derived equations and the daily salinity and temperature values that were obtained using the HYCOM. The error ranges of the simulated chemical concentrations were estimated using the double integral mean value theorem. The estimated annual northward fluxes of TA, DIC, N, P, and Si are, 1.87 ± 0.79 × 10 6 mol C s −1 , 1.65 ± 0.71 × 10 6 mol C s −1 , 1203 ± 994 mol N s −1 , 132 ± 85 mol P s −1 , and 3062 ± 1828 mol Si s −1 , respectively. The TA, DIC, P, and Si are transported principally in the top layer, but the N flux is mainly transported in the bottom layer. The high water flux in the summer is mostly responsible for the highest chemical fluxes in the whole year. There are two kinds of chemical flux patterns. For the anomalous 30-month moving average TA and DIC, the patterns follow that of the water flux. On the contrary, the patterns of anomalous 30-month moving average nutrient fluxes are oppose to that of the water flux, but follow that of the PDO index. The largest increases of N, P, and Si fluxes are around 10 to 20% in the period of rising PDO index, especially during 2004-2007. set containing the remaining observations. By repeating the prediction of the validation set using the model fitted by the training set, MSE was obtained by computing the average of summation of squared residuals, which are the differences between observed and predicted values for each order polynomial function. The second-order polynomial regression model with interaction term was found with the smallest MSE of prediction on all chemical parameters and hence was chosen for fitting TA, DIC, N, P and Si, respectively.