Patterns in the temporal complexity of global chlorophyll concentration

Satellite ocean color products, particularly chlorophyll-a ([chl-a]), have been used for decades to estimate primary productivity, develop ecological indicators, monitor long-term trends, and evaluate Earth system models. Despite widespread use of [chl-a] time series, the role of their temporal complexity in explaining region-specific differences has been unclear. Phytoplankton dynamics can be chaotic and nonlinear, suggesting time series complexity could capture large-scale structuring of phytoplankton communities. Various complexity measures exist (fractal/Hausdorff dimension, permutation entropy, Lyapunov exponents), but global analyses for [chl-a] have been lacking. The authors assembled ~25 years (1998–2022) of daily, 25 km resolution [chl-a] from a merged GlobColour GSM product (SeaWiFS, MODIS, MERIS, VIIRS, OLCI-A/B) and also tested single-mission (MODIS) data. They introduce two complementary metrics—an elasticity-based measure of regularity in daily changes and the fractal dimension—to quantify and map global patterns in [chl-a] time series complexity and assess decadal changes and sensitivity to methodological choices.

The study builds on extensive literature using satellite ocean color to assess global trends in productivity, runoff impacts, sea ice changes, harmful algal blooms, and bloom phenology. Prior work has examined [chl-a] variability, bloom timing, and climate-driven components, and used diverse complexity metrics (e.g., fractal dimension, permutation entropy, Lyapunov exponents) in other contexts or at limited locations. Merged ocean color products (e.g., GSM) are widely used but can exhibit biases among sensors; prior studies highlight algorithmic and uncertainty considerations in ocean-color climate applications. The authors position temporal complexity metrics as complementary to existing approaches and potentially informative for redefining or refining marine ecological provinces.

Data: Daily [chl-a] at 25 km resolution, 1998–2022, from the GlobColour Hermes portal using the Garver-Siegel-Maritorena (GSM) merged product (SeaWiFS, MODIS-Aqua, MERIS, VIIRS, OLCI-A/B). For comparison/sensitivity, MODIS-only (2002–2022) analyses were also performed. Filtering and coverage: Time series with >20% of days sampled over 25 years were retained (~≥1800 days). This balances data availability and spatial coverage. Elasticity metric (regularity of daily change): - Construct a 1-day-lagged time series; compute daily change as original minus lagged [chl-a]. Missing samples propagate to change series. - Retain time series with ≥400 samples of daily change. - Create two uniform threshold sets T1 and T2 (1000 values each) spanning 1–50% of the median [chl-a] for that pixel. - For each threshold set, count the number of days where daily [chl-a] change exceeds the threshold to obtain X1 and X2. - Compute ε = X2/X1 − 1 and N = X2 − X1. - Fit linear model in log–log space: log N = −β1 log ε + β2. Define elasticity as |β1|. Larger elasticity implies stronger negative log–log relationship and greater sensitivity of exceedance counts to threshold (i.e., rougher, less regular day-to-day changes). No outlier removal or formal significance testing was applied for elasticity. Fractal dimension (FD; roughness/anomalous spikes): - Clean each [chl-a] time series by removing outliers >3 SD from the mean (typically ~1–2 days per year across locations, ~0.004% of all pixel time series flagged as outliers). - Fill gaps via bootstrap resampling: randomly sample (with replacement) from existing values within the same time series to create complete series per pixel. - Estimate fractal dimension using the box-counting method via the R package fractaldim. Annual metrics and change maps: - For each year (1998–2022), compute mean [chl-a] and FD if ≥100 sampled days; compute elasticity if ≥100 consecutive sampled days. Report global means with 95% CIs (1.96 × S.E.). - Compute percent change between 2003 and 2022 for each metric per pixel: % change = (ρ2022 − ρ2003)/ρ2003 × 100. Sensitivity analyses: - Compare merged vs MODIS-only results (2002–2022). - Assess impact of temporal aggregation (daily/weekly/monthly) on metrics. - Estimate uncertainty by recomputing metrics using random 4000-day windows (30 trials) to obtain standard errors. - Examine relationships between metrics and traditional variability measures (mean [chl-a], SD, relative SD, seasonal amplitude), and between elasticity and FD themselves. Software: Analyses in R using raster, maps, ggplot2, cowplot, fractaldim, MASS, viridis, ncdf4.

- Elasticity patterns: Oligotrophic gyres exhibited the highest elasticity (rougher, less regular daily changes). Regions with higher [chl-a] and pronounced seasonal cycles (e.g., Patagonian Shelf, Baltic Sea, west coast of Central America, southwest Pacific) showed lower elasticity, consistent with strong periodic forcing. The North Atlantic gyre had lower-than-average elasticity compared with Pacific gyres, likely due to stronger seasonal structure (e.g., spring bloom). If all daily changes were equal, elasticity would be 0. - Fractal dimension (FD) patterns: Global FD values were high and spatially consistent, typically >1.85, with a reported global mean of ~1.95. Lower-than-average FD occurred in the subtropical North and South Pacific, tropical North Atlantic, Amazon plume, and off eastern Madagascar. These spatial differences likely reflect physical processes (currents, eddies, freshwater discharge, tropical instability waves) and the episodic nature of elevated [chl-a], which can manifest as outliers. After removing outliers before FD estimation, remaining FD patterns indicate the likelihood of anomalous spikes (long tails) within time series; regions with fewer outliers show higher FD. - Sensitivity and robustness: MODIS-only analyses yielded broadly similar global patterns to the merged dataset, with minor magnitude differences in oligotrophic gyres. Temporal aggregation produced smoother series (lower FD); elasticity also changed with coarser resolution while retaining geographically similar patterns. From random 4000-day window resampling (30 trials), elasticity S.E. up to ±0.20 in oligotrophic regions; FD S.E. up to ±0.04. FD showed a nonlinear positive relationship with mean [chl-a], SD, relative SD, and seasonal amplitude; elasticity showed nonlinear negative relationships with these metrics. Elasticity and FD had no strong mutual correlation. - Interannual trends: From 1998–2022, both FD and elasticity exhibited a marked shift around 2003 followed by a gradual return toward 1998 levels. Mean [chl-a] generally decreased until 2019. A second transition around 2019 featured rising mean [chl-a], a sharp increase in FD, and a decrease in elasticity. Larger errors prior to 2003 reflect fewer observations (SeaWiFS-only). Post-2003 changes are likely influenced by expanded sensor coverage (merged record) and sensor transitions (SeaWiFS end in 2010, MERIS in 2012, VIIRS start in 2012). - Spatial changes 2003–2022: Changes in mean [chl-a] were regionally inconsistent (e.g., increases in the Black Sea, parts of the western equatorial Pacific, subtropical Indian Ocean), whereas changes in elasticity and FD were more spatially consistent. This indicates the complexity metrics are relatively insensitive to [chl-a] magnitude and can be independent of processes setting [chl-a] levels. Overall, an increase in FD with a decrease in elasticity implies fewer anomalous events and stronger, repeatable cycles in many regions, potentially reflecting both observing system changes and underlying dynamics.

The study addresses the gap in understanding the intrinsic temporal variability of satellite-derived [chl-a] by formalizing and mapping two complementary complexity metrics globally. Elasticity captures the regularity and sensitivity of day-to-day [chl-a] changes to thresholding, while the fractal dimension quantifies roughness and the propensity for anomalous spikes. The metrics reveal that regions differing in [chl-a] magnitude can exhibit similar temporal complexity, indicating that time series structure provides information distinct from concentration levels. Spatial patterns suggest climatic and seasonal forcing strength is reflected in elasticity, and physical processes and episodic events shape FD. Interannual changes in complexity (notably around 2003 and 2019) imply that observing system characteristics (sensor transitions, merged product coverage/quality) and potentially true changes in ocean dynamics both influence time series structure. The consistency of complexity changes across space, contrasted with heterogeneous [chl-a] magnitude trends, supports the notion that complexity is comparatively independent of concentration drivers. These insights have implications for model evaluation (matching climatic forcing strength to observed elasticity), sampling design (greater precision and sampling frequency in more dynamic regions), and for augmenting classifications of marine ecological provinces using temporal structure of [chl-a]. Future improvements in satellite capabilities and uncertainty characterization should enhance the interpretability and predictive utility of complexity metrics.

This work introduces and applies two global metrics of [chl-a] time series complexity—elasticity of daily changes and fractal dimension—using 25 years of daily, 25 km satellite data. The authors map robust spatial patterns, show that oligotrophic gyres have high elasticity while many seasonally forced regions have low elasticity, and identify regions with lower fractal dimensions associated with episodic dynamics and anomalous spikes. They document interannual shifts in complexity (notably around 2003 and 2019) and demonstrate that complexity changes are more spatially consistent than changes in [chl-a] magnitude, suggesting relative independence from concentration levels. These metrics provide new avenues for monitoring, model benchmarking, and potentially redefining ecological provinces based on temporal structure. Future research should: (i) incorporate explicit uncertainty from ocean color algorithms and sensor-specific errors; (ii) explore spectral drivers of temporal complexity; (iii) relate complexity to environmental and ecological processes (e.g., mixed layer dynamics, eddies, temperature, plankton growth/grazing); (iv) perform formal statistical assessments (e.g., false discovery rates, threshold selection sensitivity, linearity assumptions) to solidify inference; and (v) evaluate impacts of data resolution and coverage to optimize observing strategies.

- Data and coverage: Analyses required >20% temporal coverage over 25 years, excluding many pixels (notably high latitudes). Annual metrics required ≥100 sampled days (and ≥100 consecutive days for elasticity), potentially biasing spatial representativeness. - Sensor/algorithm dependencies: The merged GSM product can exhibit biases among sensors and over time. Sensor transitions (SeaWiFS/MERIS end; VIIRS start) and changes in spatial/temporal coverage likely influenced complexity metrics, especially around 2003 and 2010–2012. - Outlier handling and methods: FD via box-counting is sensitive to outliers; authors removed values >3 SD and filled gaps through bootstrap resampling, which may affect fine-scale properties. Elasticity did not remove outliers and assumes a linear relationship in log–log space without formal significance testing. - Resolution effects: Coarser temporal aggregation smooths time series (lower FD) and alters elasticity. Reported patterns may vary with chosen resolution. - Uncertainty not propagated: The study did not incorporate per-pixel [chl-a] retrieval uncertainties or algorithm-specific error models, limiting quantitative attribution of changes in complexity to environmental versus measurement effects.