Eutrophication trends in the coastal region of the Great Tokyo area based on long-term trends of Secchi depth

Database of the ocean-going monthly observations

Monthly ocean-going observations have been conducted since April 1963 in Tokyo and Sagami Bays (Fig. 1), using vessels belonging to Kanagawa Prefectural Fisheries Technology Center. There were total 29,460 stations with unique ID based on the observation year, month, and geographic position. The data observed during the 20^th century were primarily collected from field books or paper documents issued by Kanagawa Prefecture; resulting in some missing information such as equipment types. Observations of the target month were always conducted between the end of the preceding month and the beginning of the target month. For example, when the target month is February, the observation was always conducted between the 25th of January and the 10th of February.

Various items described in the documents included date; location; Forel-Ule scale water color; Secchi depth; and vertical profiles of temperature, salinity, nitrate, nitrite, ammonium, and phosphate concentrations. Only the year and target month were recorded as dates. Temperature and salinity were recently measured using a conductivity-temperature-salinity (CTD) sensor (SBE911plus, Sea-Bird Electronics, Inc., Bellevue, WA, USA); however, there was no information on when the CTD observations started and what equipment was used until then. Temperature and salinity were documented at the 14 depths (0, 10, 20, 30, 50, 75, 100, 150, 200, 250, 300, 400, 500, and 600 m). Samples for nutrient analysis were collected at 0, 10, 30, 50, 100, and 250 m. The sampling layer was unchanged when the bottom depth of the station was shallow. Water samples for nutrient analyses were kept frozen until on-land measurements during these 30 years, but details before then are not described. Nutrient concentrations were manually measured from April 1967 to June 2010, after April 2012, they were measured using an autoanalyzer (QuAAtro2-HR, BLTEC) based on standard methods. Nutrient data from July 2010 to March 2012 were unavailable.

Other database

The monthly variations of precipitation data at the Tokyo observatory (35°41.5′N, 139°45.0′E) were downloaded from Japan Meteorological Agency. The annual means of nitrate plus nitrite, total nitrogen (TN), and total phosphorus (TP) concentrations in the Sagami River (Mairi Bridge, at 35°19.9′N 139°22.0′E) which flows into Sagami Bay were observed from 1970, 1984 and 1980, respectively, and downloaded from the water information system, Ministry of Land, Infrastructure, Transport and Tourism (http://www1.river.go.jp/; accessed April 9^th 2023). Those monthly concentrations could not be obtained. The Kuroshio axis position was provided from Marine Information Research Center, Japan Hydrographic Association. The Kuroshio axis position is the fee-charging content, and thus we did not show it as figures and did not put it in the data file.

Data quality check and trimming

Because the methodology for historical ocean-going observations was unknown, data quality checks and validations were conducted on our database. First, we removed data whose temperature was >30 °C or whose salinity was >36. When the data exceeded these thresholds, the value was one to two orders of magnitude higher than the usual value. Data with Forel-Ule scale watercolor ≥22 and the Secchi depth ≥100 m was also removed. A Forel-Ule scale watercolor of ≥12 was only observed in 1973 and 1983 before 1993, while Secchi disk depth was sometimes recorded ≤2 m before 1993. This indicates that Ule-scale observations were not frequently conducted before 1993; thus, the Forel-Ule scale watercolor data were not used for further trend analysis. While Ule-scale observations were partly limited, the relationship between Secchi depth and Forel-Ule scale water color was significantly negative (Fig. S1), indicating that the Secchi depth is a good indicator of ocean surface color. Surface chlorophyll a concentration was not measured in the monitoring program by Kanagawa Prefecture, but Shibata & Aruga (1982) reported that the surface chlorophyll a concentration showed a strong agreement with Secchi depth even in Tokyo Bay. Thus, the Secchi depth could be treated as the index of phytoplankton biomass in our study.

Regarding nutrient concentrations, variations in nitrate and ammonium concentrations indicated possible artificial errors. At a depth of 100 m, where nutrient samples were only collected in Sagami Bay after 1993, the yearly median values of nitrate concentration were usually ~10 µM; however, those in 2006 and 2007 were ~2 µM, and those in 2012 and 2013 were >40 µM. Forty-micromolar nitrate at a depth of 100 m has never been observed downstream of the Kuroshio and Sagami Bays (Aoki et al., 2022; Hashihama et al., 2008; Kodama et al., 2014b). The reliable ammonium measurement was not enough, but rarely >0.5 µM in this area (Kodama et al., 2014a; Yasunaka et al., 2018). However, the yearly median ammonium concentration at 100 m depth was usually >0.5 µM before 2010. Thus, the nitrate and ammonium concentrations during all periods were not used in the downstream analysis.

Using the same criterion, phosphate and nitrite concentrations were considered reliable because manual standard methods were simpler than those for nitrate and ammonium. For the phosphate concentration, further validation was conducted using its relationship with temperature. A negative correlation between phosphate and temperature was expected based on Kodama et al. (2014b), using Eq. (1):

(1) $$\left[P{O}_4\right]\sim a\times T+b$$where, [PO₄] and T represent the phosphate concentration and temperature, respectively. We calculated the slope (a) and intercept (b) for each year using data at depths ≥20 m. The mean slope and intercept with standard deviation (SD) values were −0.0673 ± 0.185 µM/°C and 1.69 ± 3.36 µM, respectively. The outlier slope and intercept values, defined as those without mean ± 1 × SD, were observed in 1967–1969 and 2003, and thus the phosphate concentration in these years was not considered. No outliers were found in the nitrite data.

Observations were conducted at over 70 unique stations during this observation period, including stations observed for only 1 year (Fig. 1). The frequency of observations indicated that 35 stations have been repeatedly observed over the last 50 years (Fig. 1).

Division to subareas

At some of 35 stations showed similar changes of parameters, and the observational errors are not completely removed by the above quality check and trimming processes. Thus, we classified 35 stations into several subareas using Ward’s hierarchical clustering analysis with Euclidean distance of temperature and salinity at a depth of 10 m. Consequently, the 35 stations were classified in four clusters (Fig. 2), and the cluster distributions agreed with the geographic distributions (Sagami Bay, Sagami Nada, Uraga Channel, and Tokyo Bay). Cluster I, which contained 13 stations, was observed in Sagami Bay. Cluster II, containing 13 stations, was primarily observed in offshore area of Sagami Bay, named as “Sagami Nada.” Cluster III, which contained only three stations, was observed in offshore of Tokyo Bay, named as “the Uraga Channel.” Cluster IV, which contained six stations, was observed in the Inner Tokyo Bay. A downstream analysis was conducted based on the median for each subarea calculated based on routine 35 stations.

Statistical analysis for long-term trend

Using the median of each subarea and month, all statistical calculations were conducted using R (R Core Team, 2023). The significance of seasonal and long-term trends in Secchi depth, temperature, salinity, phosphate concentration, and nitrite concentration at 10 m depth were investigated using a linear model with a Prais-Winsten estimator. A Mann-Kendall analysis is usually used to detect a trend (Kubo et al., 2019), but this approach can only show significantly positive, significantly negative, or insignificant trends, and we could not determine the change rate. Therefore, we calculated the trend using a linear model with the Prais-Winsten estimator. These trend analyses were conducted for the medians of each month and subarea using Eqs. (2) and (3).

(2) $$Y_{i,m}\sim a\times year+b$$

(3) $$Y_i\sim a\times year+f\left(month\right)$$where, Y, i, m, a and b denote the response variables (Secchi depth, temperature, salinity, phosphate, nitrite, precipitation amount at Tokyo), subarea, target month, coefficient and intercept, respectively. In Eq. (3), the f value indicates that the month was transformed into categorical values to allow for nonlinear seasonal variations as intercepts. Thus, the interannual trend was identified regardless of the month.

To identify the factors controlling the variation in Secchi depth, a linear model was developed using temperature, salinity, nitrite, phosphate, precipitation at Tokyo, Kuroshio axis position, and year as response variables. Before setting the model, seasonality was removed from all the data as follows: (1) the monthly median (monthly climatological values) was calculated for each area and parameter and (2) the difference from the monthly climatological values was calculated. These differences were not seasonal (analysis of variance (ANOVA): df = 11, p > 0.05). Most of these differences were autocorrelated (Durbin-Watson test, threshold p = 0.05); autocorrelation was not identified for the Secchi depth in Tokyo Bay (n = 539, DW = 1.8652, p = 0.05355) or phosphate concentration in Sagami Nada (n = 527, DW = 1.9102, p = 0.1408). Thus, the Prais-Winsten estimator was used in the following linear model (Eq. (4)): (4) $$\begin{array}{rl}& \mathrm{\Delta }{Z}_{SD}\sim a_1\times \mathrm{\Delta }T+a_2\times \mathrm{\Delta }S+a_3\times \mathrm{\Delta }N{O}_2+a_4\times \mathrm{\Delta }P{O}_4+a_5\times \mathrm{\Delta }rain+a_6\times \mathrm{\Delta }K\\ & +a_7\times year+b\end{array}$$where, ∆Z_SD, ∆T, ∆S, ∆NO₂, ∆PO₄, ∆rain, ∆K are the differences from the monthly climatological values of Secchi depth, temperature, salinity, nitrite, phosphate of the ocean, precipitation at Tokyo, and mean Kuroshio position between 139°E and 140°E, respectively. The a₁ to a₇ were the coefficients and b was the intercept. The most suitable model was defined as the one with the lowest Akaike information criterion.

We used 10 m depth values in this study because the environments at 0 m depth may be impacted by short-term phenomena, such as rainfall and river plumes, and are inappropriate for long-term analysis. The comparison between the medians at 0 and 10 m depth was mostly similar, except for salinity (Fig. S2). The median salinity at a depth of 10 m was usually higher than that at 0 m for every subarea (Fig. S2C). In particular, the salinity at 0 m in Tokyo Bay was occasionally <25, indicating that the effects of freshwater must be considered in this subarea.

Seasonality of subareas

Monthly variations in temperature at a depth of 10 m showed similar variations among the subareas (Fig. 3A); the temperature was lowest in March (February in Tokyo Bay) and highest in September. The highest temperature was not largely different (~24 °C in Tokyo Bay and the Uraga Channel, and ~25 °C in Sagami Nada and Sagami Bay), but the lowest temperature was spatially different (medians: 10.0 °C in Tokyo Bay, 12.4 °C in the Uraga Channel, 14.6 °C in Sagami Bay, and 15.1 °C in Sagami Nada).

For salinity, the spatial variation was larger than the seasonal variation (Fig. 3B). In all areas, the lowest salinity was observed in September and the highest salinity was observed in March. The monthly median salinities in Tokyo Bay were 31.1–32.6, those in the Uraga Channel were 32.6–33.9, those in Sagami Bay were 33.6–34.6, and those in Sagami Nada were 33.9–34.6. In Tokyo Bay, interannual salinity variations were large in September and October (interquartile ranges: 2 and 1.6).

Variations in nitrite and phosphate levels also differed among the subareas (Figs. 3C and 3D). In Tokyo Bay and the Uraga Channel, nitrite and phosphate concentrations were not depleted ≤0.1 µM, and nitrite concentration was lowest in August, but phosphate concentration was lowest in May (Tokyo Bay) and June (the Uraga Channel). Between June and September, both nitrite and phosphate were depleted in both Sagami Nada and Sagami Bay. The phosphate concentration was highest in February in these two subareas, and the nitrite concentration was similar from December to April.

The monthly variation pattern of Secchi depth differed among the subareas (Fig. 3E). Secchi depth was deepest in January for all the subareas. In Tokyo Bay, the median Secchi depth was the shallowest, varying from 2 m (June–September) to 6 m (January). These varied from 3 m (July) to 12 m (January) in the Uraga Channel, 7 m (July) to 20 m (January) in Sagami Bay, and 12 m (May) to 23 m (January and February) in Sagami Nada.

Long-term variations

The monthly and annual rates of change in temperature at 10 m depth showed a warming trend when it was significant (Figs. 4A and 5): the threshold p-value of significance was 0.05. The annual trends based on Eq. (3) were significantly elevated in all subareas (Fig. 4A). The trends were 0.022 ± 0.0028 °C year⁻¹, 0.014 ± 0.0029 °C year⁻¹, 0.011 ± 0.0033 °C year⁻¹, and 0.019 ± 0.0030 °C year⁻¹ in Tokyo Bay, the Uraga Channel, Sagami Nada, and Sagami Bay, respectively. The detailed results of the statistical analyses (degrees of freedom, t-value, and p-value) are shown in Table S1. The monthly temperature was significantly elevated in Tokyo Bay, except during August and September (Figs. 4A and 5A). The significant monthly warming trend in Tokyo Bay with standard error was ≥0.010 ± 0.004 °C year⁻¹, most of which were ≥0.023 °C year⁻¹. In the Uraga Channel, Sagami Nada, and Sagami Bay, a warming trend in temperature was observed at 5, 3, and 7 months, respectively (Figs. 4A, 5B and 5C).

The annual salinity trend indicated that salinity at 10 m depth significantly (p < 0.05) decreased in Tokyo Bay (−0.0040 ± 0.0017 year⁻¹), the Uraga Channel (−0.0048 ± 0.0013 year⁻¹), and Sagami Nada (−0.0014 ± 0.0006 year⁻¹), and insignificantly in Sagami Bay (−0.00091 ± 0.00090 year⁻¹, Fig. 4B, Table S2). The significantly decreasing monthly salinity trends were temporally limited in every subarea (Figs. 4B and 6, Table S2): October in Tokyo Bay and the Uraga Channel, January and December in Sagami Nada, and November in Sagami Bay.

The annual trends of both nitrite and phosphate concentrations decreased significantly in all subareas (p < 0.001, Figs. 4C and 4D, Tables S3 and S4). The annual decreasing trends of these two nutrient concentrations were higher in Tokyo Bay (nitrite: −0.0077 ± 0.0024 µM year⁻¹, phosphate: −0.0084 ± 0.0018 µM year⁻¹) and the Uraga Channel (nitrite: −0.0081 ± 0.0020 µM year⁻¹, phosphate: −0.0055 ± 0.0010 µM year⁻¹), and lower in Sagami Nada (nitrite: −0.0029 ± 0.00088 µM year⁻¹, phosphate: −0.0022 ± 0.00059 µM year⁻¹), and Sagami Bay (nitrite: −0.0018 ± 0.00058 µM year⁻¹, phosphate: −0.0017 ± 0.00036 µM year⁻¹). When we see the details, the peak of nitrite concentration in Tokyo Bay and the Uraga Channel was observed around 1990, in particular, from February to September (Fig. 7). This peak was also observed from June to September at Sagami Nada and Sagami Bay (Fig. 7). However, such peaks were not observed for phosphate concentration (Fig. 8). The phosphate decreasing trends were significant from November to April in Tokyo Bay and Uraga Channel (Fig. 8). A significant decrease in phosphate was observed only in February in Sagami Nada and in February, April, August, October, and December in Sagami Bay (Fig. 8).

The annual Secchi depth trend showed significant shallowing, except in Tokyo Bay (p < 0.05, Fig. 4E, Table S5). The annual shallowing trend was similar in Uraga Chanel (−0.076 ± 0.012 m year⁻¹) and Sagami Nada (−0.075 ± 0.011 m year⁻¹), and was slower in Sagami Bay (−0.046 ± 0.011 m year⁻¹). At the monthly level, a shallowing trend was usually observed in winter (January–March, Figs. 4E and 9). In the Uraga Channel, shallow trends were observed from December to June, with the exception of May. Sagami Nada was observed from November to June, with the exception of January and March, whereas Sagami Bay was observed from January to March. In Tokyo Bay, a shallowing trend was only observed in April, while a deepening trend was observed in August with Secchi depth ≤3 m (Fig. 9A). No peak in Secchi depth was observed in the Uraga Channel, Sagami Nada, or Sagami Bay (Fig. 9).

The monthly variations were shown in Figs. 5–9, but the annual mean values calculated using the least squared mean technique were shown in the supplemental document (Fig. S3).

The annual precipitation amount at the Tokyo observatory was significantly increased (p = 0.0003, Fig. 10). Focusing on monthly values, the significant increasing trends were observed in August and September (p = 0.048 and 0.032, respectively, Fig. 10).

The nitrate plus nitrite concentration in Sagami River water was significantly increased from 1970 to 2020 (coefficient ± SE: 0.021755 ± 0.005785, p = 0.000425), while the peak was observed 1990–2000 and gradually decreased after 2000 (Fig. 11A). In fact, nitrate plus nitrite concentration in Sagami River water significantly decreased after 1995 (coefficient ± SE: −0.006199 ± 0.001120, p < 10⁻⁵). The TN and TP concentrations were significantly decreased (p < 0.005). As well as nitrate plus nitrite concentration, the TN and TP gradually decreased after 2000 (Figs. 11B and 11C).

Contributions of Secchi depth anomaly

Based on Eq. (4) and the AIC values, the contributions of temperature, salinity, nitrite, phosphate, and Kuroshio axis position to the Secchi depth anomaly variation were estimated throughout the year (Fig. 12, Table S6). The precipitation anomaly at Tokyo has never been selected in all the subareas. In Tokyo Bay, only salinity anomalies remained in the final (least-AIC) model to explain the Secchi depth anomaly, and the relationship between the Secchi depth and salinity anomalies was significantly positive (i.e., an increase in salinity deepened the Secchi depth). This significant positive relationship between Secchi depth and salinity anomalies was the same as that in the other subareas (Fig. 12B). In the Uraga Channel, temperature, salinity, and phosphate and Kuroshio axis anomalies remained in the final model, and the coefficient indicated that the Secchi depth increased with increases in the temperature and salinity anomalies (Fig. 12). The effects of phosphate and Kuroshio axis on the Uraga Channel remained in the least-AIC model, but insignificant (p > 0.07). In the Sagami Nada, the variables of the full model except precipitation and Kuroshio axis anomaly were selected. Elevations in temperature were related to the deepening of the Secchi, but the effects of phosphate were not significant. In Sagami Bay, this trend was not selected as the final model, and warming and phosphate decrease deepened the Secchi depth; the coefficient of the nitrite anomaly was found to be insignificant (p > 0.05).

Since the Secchi depth of Tokyo Bay and Uraga Channel was several meters, the 10 m depth was not in the euphotic layer, 1% of surface radiance, and the contributions of parameters would be different when we used the 0 m variables. However, we noted that when we used the 0 m variables instead of the 10 m depth variables in Tokyo Bay, the phosphate anomaly was selected in the final model as well as the salinity anomaly, but it was insignificant (p = 0.319). The coefficient of salinity anomaly was significantly positive (0.279 ± 0.117, p = 0.017), which was similar to that when using 10 m depth values. In Uraga Channel, the Kuroshio axis was not in the final model compared to that using 10 m depth values, and the other remaining explanatory variables were the same: the relationships with temperature (0.6585 ± 0.1373, p < 10⁻⁵) and salinity (2.3044 ± 0.2829, p < 10⁻²⁰) were significantly positive as well as that using 10 m depth values, and that of phosphate was insignificant (0.4816 ± 0.2781, p = 0.084).