Multidecadal changes in ocean transparency: Decrease in a coastal upwelling region and increase offshore

Detection of the effects of climate change on ocean ecosystems is often limited by the short duration of available time series. Here, we use ocean transparency measurements (the Secchi disk depth, ZSD) in the California Current Ecosystem since 1949 and combine them with satellite estimates. Historic in situ measurements of ZSD were irregular in space and time and are difficult to interpret in time series due to biases introduced by changing locations and timing. We normalize historic ZSD measurements with satellite‐derived mean climatology and create a merged in situ—satellite time series of ZSD for the last ~ 73 yr. Although interannual variability in ZSD is dominated by El Niño Southern Oscillation‐related variability (~ 50% of the total variance in many areas), a secular trend of decreasing transparency that is correlated with increasing productivity is detected 0–300 km from the coast in an area affected by coastal upwelling north of 27°N. In contrast, increasing transparency (correlated with decreasing productivity) is detected offshore (> 1000 km from the coast). In addition to those general trends, transparency is also increasing in coastal area off Baja California south of 27°N.

Oceanic phytoplankton play an important role in the global carbon budget and any changes in their productivity and composition are therefore of special interest for global climate change and biogeochemical cycling Falkowski et al. 2008). A claim by Boyce et al. (2010Boyce et al. ( , 2012Boyce et al. ( , 2014) that global phytoplankton biomass has drastically declined over the last 100 yr has been challenged (McQuatters-Gollop et al. 2011 and others) and is affected by uncertainties in the available data. Climate-driven trends in global phytoplankton are typically evaluated using satellite data but these time series are relatively short and open to different interpretation (Behrenfeld et al. 2006Martinez et al. 2009). Longer time series can be estimated using geochemical methods (Osman et al. 2019) or numerical models (Dutkiewicz et al. 2019;Cael et al. 2021).
Water transparency is an integrative variable that includes the effects of absorption by a wide variety of dissolved substances as well the effects of absorption and scattering by phytoplankton and non-algal particles (Preisendorfer 1986). Water transparency affects processes such as primary production, solar bleaching of organic material, and phytoplankton physiology (Kirk 1994), as well as the encounter rates of zooplankton with fish (Aksnes and Giske 1993; Aksnes and Utne 1997) and zooplankton diel vertical migration behavior (Ohman and Romagnan 2016). Also, the other way round, processes like primary productivity and bleaching affect transparency. Water transparency can be measured with optical instruments and estimated with satellite sensors but a simple method called the Secchi disk depth (Z SD , m) using the depth of disappearance of a white disk has been used for over a century (Tyler 1968;Preisendorfer 1986;Aksnes and Ohman 2009;Boyce et al. 2010;Pitarch 2020;Pitarch et al. 2021). Although the Secchi disk depth is a rather simple measurement, its advantage is the large number of measurements extending back many decades. Aksnes and Ohman (2009) documented a decrease in Z SD in the nearshore region of the southern California Current Ecosystem (CCE) over the period . However, the variability was too large to detect any Z SD trends offshore.
In oceanic waters, water transparency is closely correlated with phytoplankton concentration (typically measured as chlorophyll a [Chl a]; Falkowski and Wilson 1992) but is also sporadically reduced by river outflow plumes, other inputs of dissolved organic matter, and by resuspension of inorganic particles near the coast. Events such as phytoplankton blooms can markedly reduce Z SD . All these events typically last from days to weeks. In contrast, on longer time scales, interannual variability due to perturbations such as the El Niño Southern Oscillation (ENSO) and anomalous heatwaves (e.g., "the Blob," Bond et al. 2015) have strong effects on the whole marine ecosystem (Kahru and Mitchell 2000;Kahru et al. 2012Kahru et al. , 2018Jacox et al. 2018;Lilly et al. 2019) and can make it challenging to assess whether long-term secular changes have occurred in an ocean region. In the Baltic Sea, Dupont and Aksnes (2013) fitted an empirical relationship between Z SD and the distance to coast and bottom depth to remove the mean effect of those two variables. However, such effects are not uniform in space and time. Kahru et al. (2022) removed both the mean spatial and seasonal variability in Z SD using satellite-derived Z SD climatology to further investigate centennial trends in the Baltic.
Here, we adopt the same approach as Kahru et al. (2022) and normalize historic shipborne Z SD measurements by subtracting the mean 10-d seasonal Z SD value of the nearest 4.5-km satellite pixel from each Z SD measurement. We then compare in situ and satellite time series and merge them, in order to reconstruct a seven-decade record of water transparency in the CCE. The goal of this paper is to evaluate interannual variability and test for possible secular trends in water transparency in the CCE, using it as a proxy for primary productivity.

Satellite data
A number of models exist for calculating Z SD from optical data (Preisendorfer 1986), but most are empirical in nature (Prasad et al. 1998) and less reliable in variable conditions. The algorithm of Lee et al. (2015), on the other hand, is based on new theoretical interpretation of sighting a Secchi disk in water, which adapts to the change in the spectral composition of subsurface light. The algorithm and the model have been validated in different locations with different sensors (Lee et al. 2015(Lee et al. , 2018Liu et al. 2019;Kahru et al. 2022;Brewin et al. 2023).
We calculated net primary production (NPP, mg C m À2 d À1 ) using the well-known Vertically Generalized Productivity Model (VGPM) algorithm (Behrenfeld and Falkowski 1997) adapted to California Cooperative Oceanic Fisheries Investigations (CalCOFI) data (Kahru et al. 2009(Kahru et al. , 2020. As inputs, we used daily OC-CCI Chl a (mg m À3 ), photosynthetically active radiation (PAR, merged from multiple sensors processed by NASA OBPG) and optimally interpolated daily sea-surface temperature (Reynolds et al. 2007). NPP was also estimated using the more recent CAFE model (Silsbe et al. 2016). Monthly CAFE model NPP products based on SeaWiFS and MODISA data were downloaded from http://sites.science.oregonstate.edu/ocean. productivity/index.php and merged.
Monthly anomalies of Z SD and Kd 490 were calculated for each pixel by removing the climatological monthly mean pixel value. Monthly Chl a and NPP anomalies were calculated by using the ratio of monthly mean value to the mean climatological value of the respective month. The ratio anomaly was then expressed as percentage anomaly with 100 * (Anomaly -1).

In situ data
Aksnes and Ohman (2009) assembled in situ Z SD data in the California Current region for 1949-2007. Most of these data were from the CalCOFI cruises conducted approximately quarterly (Ohman and Venrick 2003). We added additional measurements from the Southern California Bight Study (SCBS, 1975-1987, Eppley et al. 1985 and from recent CalCOFI cruises (https:// calcofi.org/) until cruise 2021/11. Early Z SD measurements (until 1969) were made with different methods (including with artificial lighting at night and a disk with black and white quadrants by day) and these were transformed and standardized by Aksnes and Ohman (2009). The total number of in situ Z SD measurements assembled was 6557 for the time period 1949-2021 (Table 1). The locations of those in situ measurements changed over time ( Fig. 1), with more extensive coverage in earlier decades. During the period of satellite data availability (beginning September 1998), the areal coverage of in situ data was more limited ( Fig. 2A). We refer to that area as the main study area (MSA).
In situ Z SD measurements were converted to Z SD anomalies by subtracting the satellite-derived mean Z SD value Table 1. Number of in situ Secchi depth measurements per decade (see Fig. 1 for a map).

Normalization of Z SD data
Satellite measurements over 25 yr  were used to create Z SD and Kd 490 mean climatologies (i.e., mean 1950-1951 1949 1969 1970-1979 1980-1989 1990-1999 2000 seasonal cycles) with 10-d temporal and 4.5 km spatial resolution. Examples of Z SD and Kd 490 climatologies averaged over coastal, transitional, and offshore domains for the area of Fig. 2A are shown in the Supporting Information (Fig. S1). Daily datasets were composited over 10-d periods by averaging valid pixel values (i.e., the composited value is the mean of 1-10 valid values). During cloudy periods this procedure still resulted in some missing pixel values. Those missing pixels in 10-d composites were filled by linear interpolation between corresponding valid values of previous and following 10-d composites.
Corresponding 10-d periods over all available years were averaged to produce mean maps of satellite-derived Z SD and Kd 490 for each of the 10-d periods over a year (the last period consisted of either 5 or 6 d). These 37 datasets therefore consist of the mean pixel values for year days 1-10, 11-19, 20-29, and so on. The mean annual cycle for each pixel is made up of these 37 values. Anomalies of in situ Z SD were created by subtracting the respective climatology value of the nearest pixel in space and the 10-d period corresponding to the date of the measurement. This normalization of the numerous but spatially and temporally irregular historic samples made measurements in different areas (with different means) more comparable. The same procedure was applied to satellite data by subtracting the climatology value from the 10-d Z SD and Kd 490 composites. Z SD anomalies were averaged over monthly periods to produce monthly anomalies that were subsequently averaged annually to produce annual Z SD anomalies for both in situ and satellite values. Annual anomalies of each pixel were averaged over different spatial domains to produce time series for those domains.
As indices of ENSO, we used the Multivariate ENSO index (MEI.v2, Wolter and Timlin 2011; https://psl.noaa.gov/enso/ mei/) and the San Diego detrended Sea Level Anomaly (SDSLA; Lilly and Ohman 2018). The SDSLA timeseries consists of monthly average sea-level anomalies, from which the seasonal cycle and long-term trend were removed, and, while correlated with equatorial indices of ENSO, is a more appropriate index of ENSO expression in this mid-latitude region (Lilly and Ohman 2018).
Potential trends in time series and their statistical significance were evaluated using the nonparametric Sen slope estimator (Sen 1968;Gilbert 1987) with the nonparametric Mann-Kendall test to evaluate the statistical significance of the trend (Salmi et al. 2002).

Validation of satellite estimates
The Lee et al. (2015) Z SD model has been applied in a number of regions (Shang et al. 2016;Liu et al. 2020a,b;Kahru et al. 2022;Brewin et al. 2023). We validated the satellitederived Z SD datasets at daily and 4.5-km resolutions. We found a total of 443 same-day matchups with in situ measurements and used the average pixel value in a 3 Â 3 pixel window centered at the in situ measurement. Match-ups with high pixelto-pixel variability ([Max À Min]/Min > 0.5) and with less than five valid pixels out of nine were excluded. The relationship between in situ and satellite Z SD was evaluated using both the ordinary least squares linear regression and a Type II reduced major axis linear regression (York et al. 2004). We detected no significant bias between satellite and in situ data until Z SD of about 30 m. At in situ Z SD > 30 m, satellite-derived values underestimated the in situ values. To match satellite data product with in situ measurements for those deep Z SD waters, we applied an empirical conversion by fitting a power function (Y = 1.59 Â X 0.82 , where Y is satellite-derived Z SD and X is in situ Z SD ) and then inverting it. We applied the conversion only for Z SD > 15 m. The modified algorithm performed well (r 2 = 0.77, Fig. 2B) with a mean absolute percent difference $ 15% and median absolute percent difference $ 11%. Some of the scatter around the curve is expected to be due to the difference in scale between an in situ point measurement and the mean of a 3 Â 3 pixel area. The discrepancy at the clearest waters may be caused by problems with satellite estimates of R rs as the Lee et al. (2015) algorithm has been shown to work well also in clear waters when using in situ R rs data (Brewin et al. 2023). In conclusion, we assume that our converted satellite estimates of Z SD are reasonably accurate and representative of in situ Secchi depth measurements.

Time series of Z SD anomaly in the MSA
The MSA (Figs. 2A) was spatially well covered by both in situ and satellite measurements but the temporal coverage of the in situ data was restricted to approximately four cruises per year. The less frequent temporal coverage and the spatially less dense sampling probably explain why the monthly means of in situ Z SD were more variable than the satellite-derived monthly means (Fig. 3A vs. Fig. 3B). However, the annual means of the in situ and satellite time series (Fig. 3C) were in good agreement (r 2 = 0.78, p < 0.001). The annual in situ Z SD anomalies were slightly more variable than the satellitederived annual averages. Interestingly, the biggest discrepancy between the satellite and in situ values was during the 2014-2015 warm anomaly when the in situ Z SD values were higher than the satellite-derived values. Satellite-derived Kd 490 (Fig. 3D) behaved like an inverse of the satellite-derived Z SD with high correlation (r 2 = 0.80, p < 0.001). The minor differences may be due to the fact that Z SD is based on the most transparent band that is not necessarily the 490-nm band. The most transparent spectral band in the CCE is typically the 443-nm band for offshore waters that switches to longer wavelengths in coastal waters. Annual anomalies in Z SD are negatively correlated with anomalies of Chl a (Fig. 3E, r 2 = 0.98, p < 0.001) and NPP estimated with VGPM-CAL (r 2 = 0.96, p < 0.001) and CAFE (r 2 = 0.75, p < 0.001). It is notable that the coefficient of determination (r 2 ) is highest with Chl a, followed by VGPM-CAL, Kd 490 , and CAFE.
As satellite-derived Z SD estimates have better temporal and spatial coverage, we merge in situ and satellite data by using in situ values until 1997 and satellite-derived values from 1998 onward. The combined Z SD anomaly time series averaged over the entire MSA for the 1949-2021 time period (Fig. 4) had a small linear decreasing trend of À0.03 m yr À1 that would result in decreasing Z SD by 2.2 m over the 1949-2021 time period, but the trend was not statistically significant (p = 0.13) due to the large monthly and interannual variability. During the period of satellite data (1997-2021; Fig. 3B) the linear decreasing trend of monthly satellite-derived Z SD anomaly in MSA was significant at p < 0.05 but not at p < 0.01. Monthly Z SD means are correlated with ENSO indices such as MEI and SDSLA (Fig. 5). However, the strength and even the sign of the correlation is dependent on the location (Fig. 5B).
Time series of Z SD anomaly in 1 Â 1 spatial bins As temporal trends in Z SD can be either positive or negative dependent on location, the trend estimates averaged over large areas, such as MSA, may be misleading. Although satellite time series can be evaluated at high spatial resolution, in situ data are spatially sparse. We therefore used 1 Â 1 spatial bins to evaluate trends in Z SD time series. Due to the spatially and temporally irregular sampling pattern, particularly in the early decades, the start year and the total number of years with Z SD data were variable (Fig. 6). In most bins the measurements started around year 1950 (Fig. 6A) but were followed by a long gap (1952)(1953)(1954)(1955)(1956)(1957)(1958)(1959)(1960)(1961)(1962)(1963)(1964)(1965)(1966)(1967)(1968); Table 1; Fig. 4). We considered bins with at least 26 yr of Z SD data (i.e., 25 yr of satellite data and at least an additional 1 yr of in situ data). Most bins had at least 30 yr of data and the bins off southern California had 40-53 yr of data (Fig. 6b).
As Z SD anomalies are correlated with the ENSO cycle in this region (Fig. 5), the time series can be biased by sampling during a particular phase of the ENSO cycle. This was particularly relevant during the early period when the number of sampled years was small and sampling was not annual. We therefore removed the average effect of ENSO in each 1 Â 1 spatial bin by subtracting the expected Z SD value estimated from the ENSO index. We constructed linear regression models between an ENSO index (both MEI and SDSLA) and monthly Z SD for each 1 Â 1 bin (like Fig. 5A) using satellite data of 1997-2021. We then removed the mean ENSO effect by subtracting the predicted monthly Z SD value estimated from the monthly SDSLA index. The corrected monthly Z SD values were then averaged annually.
After applying the nonparametric Sen slope estimator and masking all 1 Â 1 bins where the trend was not significant at 95% ( p < 0.05), there were 40 bins with significant decreases and 10 bins with significant increases in Z SD (Fig. 7A). The ENSO correction removed some of the scatter and improved the detection of change. The bins with decreasing Z SD were primarily in coastal upwelling-affected areas along the coast of California and the bins with increasing Z SD were in a limited region off southern Baja California. In the majority of bins, the trend estimated as the Sen slope was not significant. The Sen slope is estimated as the median of all possible slope estimates. Therefore, for bins with 25 yr of satellite data and only a few years of in situ data the Sen slope can be significant only if the trend is also present during the satellite period. In offshore areas with only a few years of presatellite measurements, the Sen slope would not detect a significant change even if there was a big jump between the early period  and the satellite-period. To evaluate the change from the presatellite period 1949-1997 to the present 1998-2021 period we averaged the annual Z SD anomalies during those two periods and calculated the change in the mean Z SD (Fig. 7B).
Although the dominance of decreasing Z SD in coastal areas is similar to that detected by the Sen slope, its area is much wider. In contrast to the changes detected by the Sen slope, a number of offshore bins showed increasing Z SD . That was not evident in the Sen slope pattern due to the relatively small number of annual values during the presatellite period. This shows that oceanic offshore waters have probably experienced a shift to increased Z SD from $ 1950 to present but the change is statistically difficult to prove due to the availability of only a few annual values from the presatellite period. The trend of increasing transparency in offshore waters becomes apparent when the change in Z SD from the presatellite period to the present (1998-2021) is pooled for three regions based on distance from the coast (Fig. 8). The 25-75% percentile range of Z SD changes shows increasing Z SD offshore, decreasing Z SD in coastal (0-300 km) waters and intermediate trends in transition zone (300-1000 km).

Discussion
Water transparency is an important habitat characteristic for all aquatic organisms that are affected by light (Kirk 1994;Dupont and Aksnes 2013). In addition to in situ measurements, water transparency can be estimated with satellite ocean color data (Lee et al. 2015) with the advantage of consistent large-scale coverage at frequent intervals. High sampling frequency is essential to discover trends in areas of high variability. For example, around the Channel Islands in the Fig. 6. Combined in situ (1949-1997) and satellite (1997-2021 time series of the annual Z SD anomalies in 1 Â 1 bins: (A) year of the first in situ measurement after 1900; (B) total number of years with available data. Fig. 7. Combined in situ (1949-1997) and satellite (1997-2021 time series of the annual Z SD anomalies in 1 Â 1 bins. (A) Sen slope of the change (m yr À1 ) that is significant at p < 0.05. Slopes that are not significantly different from 0 are masked as white. (B) Change in the mean Z SD anomaly (m) from 1949-1997 to 1998-2021. Southern California Bight where transparency is highly variable and affected by both plumes of sediment-rich turbid waters and blooms of phytoplankton (Otero and Siegel 2004), trends were difficult to detect in spite of over 50 yr of measurements.
In oceanic waters transparency is highly correlated with phytoplankton biomass (Chl a) and primary productivity (Fig. 3E). As transparency measurements are robust and much easier to make than in situ measurements of primary productivity, transparency can be used as a proxy to estimate trends in primary productivity. The high correlation among Chl a, Z SD , NPP, and other satellite-derived variables is artificially increased when using satellite-derived estimates. All these products are derived from remote sensing reflectances that are inter-correlated and have few degrees of freedom (Cael et al. 2023). With these caveats in mind, it is not surprising that Z SD may be better correlated with a model of primary production than different models of primary production using different input data with each other.
The trend of decreasing water transparency in the California Current coastal waters was first shown by Aksnes and Ohman (2009) who estimated a reduction in the Secchi depth by 0.06-0.13 m yr À1 in the nearshore region of the southern California Current over the period 1969-2007. Our merged time series confirm the general decrease in transparency in the upwelling-affected coastal areas but also show a spatially complex pattern with areas of increasing transparency and areas of decreasing transparency. Decreasing transparency (increasing productivity) is evident in most upwelling-affected areas along the coast of California excluding a limited area off Southern Baja California. Offshore oceanic waters, on the contrary, have increased transparency (decreased productivity) when comparing the pre-satellite ($ 1950-1997) period with the modern (1998-2021) period. Due to the small number of measurements in the presatellite era, the statistical significance of that change can be problematic to prove.
Interannual variability associated with ENSO has a strong effect on water transparency and primary productivity and explains $ 50% of the total Z SD variability in many locations. Although the long-term trends in water transparency and productivity are relatively small compared to their interannual variability, they are likely to have significant long-term effects on the ecosystem in the California Current area. This underscores the importance of sustaining long records in order for underlying trends to emerge from background variability (Hameau et al. 2019).
Our results support the conclusions drawn from multiple previous studies that the effect of global climate change and warming of the upper ocean is leading to decrease of productivity in ocean gyres (Polovina et al. 2008;Kahru et al. 2009Kahru et al. , 2012Boyce et al. 2010Boyce et al. , 2014Behrenfeld et al. 2016;Osman et al. 2019;Meng et al. 2021;Leonelli et al. 2022) but increase in productivity in coastal zones (Kahru et al. 2009(Kahru et al. , 2012.

Conclusions
Time series of historic shipboard measurements of water transparency combined with modern satellite data show a long-term trend of decreasing transparency (increasing productivity) in the upwelling ecosystem of the California Current during 1949-2021 but an increasing transparency (decreasing productivity) in offshore waters. Although the secular changes are relatively small compared to the ENSO-type interannual variability, their long-term effects are likely to be significant on the functioning of the ecosystem.