Statistical combination of different types of chlorofyll-a measurements in the Dutch North Sea

Both the EO and in-situ datasets were acquired through OSPAR, whereby the in-situ data can also downloaded through the ICES data portal. The chl-a indicator is used to address eutrophication, and thus only reflects concentrations in the upper 10 meters of the water column. While this depth limitation is inherent to the EO data, any in-situ measurements taken at greater depths are excluded from the analysis. Additionally, the indicator relies on growing season means, with the growing season, as defined by OSPAR, extending from March to September. The focus of this analysis is on these seven months, although algal blooms may occur outside this window, especially as climate change influences water temperatures in the North Sea^{1) 2)}. The analysis covers the period from 1998 to 2020, as both in-situ and EO data are available for these years.

2.1 In-situ data

The number of in-situ measurements and chl-a concentrations in the Southern North Sea (SNS) assessment area are summarized in Table 1. The sample size fluctuated over the study period, with a peak in 2019, featuring approximately 131 samples during the growing season. The spatial distribution of these measurements is shown in Figure 1, which reveals that sampling sites are unevenly distributed, with a higher concentration of samples taken near the British, Belgian, and Dutch coasts. Only a few sites on the Dutch Continental Shelf were sampled throughout all months of the growing season, whereas most sites were sampled only during a single month (although not all in the same month, see distribution over months sampled per year in Table 1).

2019 was the most sampled year, and the development of quantile chl-a concentrations during that year’s growing season is depicted in Figure 2. Although no data is available for March, the results show elevated chl-a concentrations during the spring bloom, particularly near the Dutch coast. However, the sparse spatial distribution of samples makes it difficult to discern broader patterns across the entire SNS area.

2.2 Earth Observation data

EO data is extremely well sampled throughout the SNS assessment area (Table 2), with millions of measurements that span the entire duration of the growing season. All months, and almost all days, are sampled during each year in contrast to the in-situ data (Table 1). For more background on the acquisition of EO data, see Van der Zande et al. (2019)³⁾ and Lavigne et al. (2021)⁴⁾.

Table 2 highlights the temporal coverage of the Earth Observation (EO) data, while Figure 3 presents the spatial distribution of this data for the year 2019. To appropriately visualize the data, chl-a concentrations were aggregated to a 5x5 km grid (the original resolution of the EO data being 1x1 km). These aggregated values were then categorized into quantiles. The results demonstrate that the EO data effectively captures both spatial and temporal patterns of chl-a concentrations throughout the growing season. Notably, higher concentrations are observed during early phytoplankton blooms from March to May, with levels being generally elevated near the coast.

In the later years of the monitoring period, the EO data provides excellent coverage of the SNS area. However, coverage was less comprehensive in the earlier years (see differences in n Table 2; and Figure 4), with some months showing gaps in 5x5 km grid data. Despite these early gaps, the EO data coverage of the SNS area remains nearly complete.

2.3 Conclusion data exploration

Large differences exist between the spatial and temporal distribution of the in-situ and EO datasets, with the two datasets not being directly comparable, as illustrated in Table 3. The number of in-situ samples is substantially smaller than that of the EO data, both in terms of temporal coverage and spatial distribution. In-situ data is not representative of the entire SNS assessment area, and in many years, there are gaps where not all months of the growing season were sampled.

Given these discrepancies, a strong recommendation is to treat in-situ samples as individual data points alongside the EO data. Aggregating both datasets before calculating growing season averages may provide a more accurate representation of the overall conditions. If and how the current weighting method may introduce biases is explored in the next chapter.

3.1 Confidence rating

In the current approach for combining in-situ and EO data to determine chl-a growing season means, the means of both datasets are first calculated separately. A weighted average is then computed based on the confidence rating of each dataset, following criteria that consider both spatial and temporal confidence aspects. These criteria are detailed in Annex 13 of the OSPAR eutrophication status assessment procedure . If the in-situ data has high confidence, a weighting of 50:50 (in-situ/EO) is applied; for moderate confidence, the ratio is 30:70 (in-situ/EO), and for low confidence, it shifts to 10:90 (in-situ/EO).

Even when the confidence of the in-situ data is low, assigning it a 10% weight remains significant, especially considering the large discrepancy in the representativeness of in-situ versus EO data across the growing season (see chapter 2), with the in-situ data representing at best 0,0039% of all combined measurements taken in 1999.

To demonstrate the impact of different weightings on growing season means, Figure 5 illustrates all weighting scenarios considered in the present study, including the individual means of in-situ and EO data, as well as the weighting used in OSPAR. With the exception of 2002, a 50:50 weighting was applied in all years. However, these weighting factors were determined somewhat subjectively. There is no statistical backing for boundaries at which the confidence classes and the corresponding weights are set⁵⁾, but as a result, high in-situ values can disproportionately affect the calculated growing season means.

Since in-situ data is collected sporadically compared to EO data, there is a possibility that extremely high chl-a concentrations may be captured due to the timing of in-situ sampling during periods of high phytoplankton density, such as algal blooms, and/or the higher variation in sample sites. For instance, in 1999, the area around the British shore was extensively sampled in contrast to other years (Figure 1), and samples were only taken during four months of the growing season (Table 1).

An argument often cited to support the use of in-situ data is the claimed higher precision of chl-a measurements compared to EO data. However, comparisons between in-situ and EO chl-a measurements in the Baltic Sea reveal that the uncertainties associated with both methods fall within the same range⁶⁾. Furthermore, these comparisons were made using chl-a products derived from MERIS and OLCI satellites. In contrast, the Sentinel satellites, which have been in use since 2016 and incorporate the JMP EUNOSAT quality control, offer even greater precision.

3.2 Conclusion current method

The current method of integrating in-situ and EO data is highly sensitive to temporal and spatial biases in the in-situ measurements. As a result, the elevated chl-a concentrations observed at the start of the time series are likely to be artifacts of the methodology. Investigating alternative approaches for data integration could provide more robust results and is recommended.

In this chapter, we present an alternative approach for determining chl-a growing season means and their associated confidence ratings, addressing the imbalances between in-situ and EO measurements. This method considers the variation in spatio-temporal distribution during the growing season and accounts for missing data.

Given the discrepancies between in-situ and EO measurements, an aggregation of these datasets is recommended before performing any further calculations, as the weighting method currently overemphasizes the influence of in-situ samples (chapter 3). Moreover, variations in the spatio-temporal distribution of both EO and in-situ measurements necessitate aggregation prior to calculating growing season means (Figure 6). Before combining in-situ and EO data into a single growing season mean, separate means are calculated. However, these means are derived by simply averaging the data over the entire growing season without accounting for temporal and spatial variations (Figure 6). In certain months and regions, more measurements may be available due to favorable weather conditions, leading to overrepresentation of specific areas or periods (for both EO and in-situ data). This approach can result in a skewed representation of chl-a concentrations, particularly as coastal regions and certain months tend to exhibit higher concentrations (Figure 3 and Figure 4).

Our approach will thus treat in-situ measurements as individual data points alongside EO data, aggregating both on a year-month-grid basis before proceeding with analyses. The focus will be on 1) determining the optimal grid size for aggregation, 2) imputing missing values, 3) calculating growing season means, and 4) introducing a new method for confidence ratings.

4.1 Optimal grid size

As stated above, a grid-based approach can effectively cope with temporal and spatial variation in chl-a measurements. The optimal grid size may be determined by identifying the sample size at which the variation in measurements becomes less pronounced. Figure 7 illustrates this relationship by plotting the standard deviation against sample size for four random 10x10 km grids across different months in 2008. While the sample size at which the standard deviation stabilizes varies between grids, a sample size of approximately 100 appears to provide a reliable estimate of chl-a. When dividing the SNS into 5x5 km grids, 12% of the grids have fewer than 100 samples, whereas only 6% of grids fall below this threshold when using a 10x10 km grid.

Grids are assigned to assessment areas based on their initial shapes, and may thus overlap between two neighboring areas. However, only the corresponding datapoints that lie within the initial shape of the assessment area are considered for further calculations. Dividing the study area into grids introduces edge effects, particularly where only part of the geographical area falls within the boundary of a grid cell (Figure 6). This results in smaller sample sizes at the edges. While in a large area like the SNS, these effects do not significantly impact the overall confidence rating, they become more pronounced in smaller coastal OSPAR assessment areas, where edge areas make up a relatively larger portion (Figure 8). Over the entire study period, using a 5x5 km grid results in 46% and 40% of grids containing fewer than 100 samples in RHPM and MPM, respectively. In contrast, using a larger 10x10 km grid reduces these percentages to 28% and 25%, thereby improving the confidence rating. Also in later years, with higher resolution EO data, edge effects are more prominent in smaller assessment areas.

Another important consideration when determining grid size is its impact on the mean growing season values. Utilizing a smaller grid size helps minimize the loss of resolution during the aggregation process. Table 4 illustrates the differences in growing season means for various grid sizes. Across the whole study period, the mean growing season difference between the 5x5 and 10x10 grids is 0.05 μg L^-1 chl-a, while the difference between the 5x5 and 25x25 grids is 0.2 μg L^-1 chl-a. Although increasing the sample size up to approximately 250 samples can improve the reliability of chl-a estimates, selecting a grid size that is too coarse beyond this threshold may lead to a significant loss of resolution. Do note that mean values in the coarser grid size are consistently higher than the smaller grid size, due to the relatively higher contribution of grids with higher mean chl-a values (high chl-a concentrations can sometimes be 100 times higher than background concentrations, which impacts mean concentrations).

A 10x10 km grid may represent the optimal choice, as it offers a sufficient sample size to provide reliable chl-a estimates on a year – month – grid basis while minimizing edge effects and resolution losses. Consequently, further analyses will be conducted using this grid size.

4.2 Missing values

Particularly in the beginning of the study period, certain year–month–grid combinations are missing (Table 5; SNS area). Since only a small portion of the data is missing, the remaining data can be used to impute these missing values. One approach is the use of a random forest algorithm, implemented through the R package missForest⁷⁾. This method preserves the relationships and distributions of the observed data, is resistant to overfitting, and is capable of handling high-dimensional data.

The accuracy of the imputation using missForest was evaluated using the Normalized Root Mean Square Error (NRMSE) metric. An exceptionally low NRMSE of 3.4 × 10^-7 was obtained, indicating that the imputed values closely align with the observed data and demonstrate strong model performance with minimal prediction error.

Given the small proportion of missing data, the impact of imputing these values on the overall mean per growing season is minimal (Table 6), with a maximum difference of only 0.02 μg L^-1 chl-a in 1998. For the SNS area, imputing missing values may therefore not be necessary, as the effect on the overall mean is negligible, and the process introduces an additional step in data handling. However, this may not be the case for other OSPAR assessment areas, where imputing missing values may be more relevant for a comprehensive analysis. For instance, the proportion of missing year–month–grid data in the MPM and RHPM regions combined varies from 11 to 17% between 1998 and 2001 (Table 5; MPM & RHPM), although missing values do not pose an issue in later years.

4.3 Trends

To illustrate the effect of the new method to calculate chl-a growing season means as compared to the current weighting method, we calculated the trend over the years using TrendSpotter⁸⁾. This approach detects and estimates non-linear trends in environmental time series, and calculates standard errors and concomitant confidence intervals based on the Kalman-filter framework. In addition, the difference between the modelled trend value (the red line in figure 10) in each year and the model value in the last year is calculated. This enables the assessment of multiplicative trends (the percentage yearly change) and confidence intervals⁹⁾, which are the base for a classification in increasing, stable, decreasing or uncertain trends (Figure 9).

By giving the in-situ data 50% weight in most years, the high in-situ chl-a concentrations, especially at the start of the monitoring period, influence the overall mean and this results in an overall moderate decrease in chl-a between 1998 and 2020, while the trend is stable in the new method (Figure 10). In the new method, in-situ and EO data were first aggregated on the 10x10 km scale, and missing values were imputed (Table 6).

4.4 Confidence rating – new method

In the current COMPEAT assessment the confidence rating per OSPAR assessment area for the chl-a measurements for both in-situ and EO data depends on the same criteria in Annex 13 of the OSPAR eutrophication status procedure6. These criteria, originally based on in-situ measurements, assign weights of 50%, 30%, or 10% to the in-situ data, though these weights lack a formal statistical foundation. Here, we propose a method that assigns a confidence rating to each year–month–grid combination based on sample size, which can then be aggregated to derive an overall confidence rating for the chl-a growing season mean each year.

To determine the confidence classes with sample size for each OSPAR assessment area, we will use the relative margin of error (MOE). The MOE quantifies the amount of random sampling error and represents the radius (= half) of the confidence interval. A 95% confidence interval reflects the 95% likelihood that the true value lies within that interval and is bounded by one MOE below and one MOE above the estimated mean. By using the relative MOE (the MOE divided by the mean), we adjust for variations in mean chl-a values across year–month–grid combinations. As shown in Figure 11, the relative MOE decreases as the sample size increases.

The sample size at which the relative MOE falls below a specified threshold can be used to define the boundaries for ‘low’, ‘moderate,’ or ‘high’ confidence ratings. The choice of relative MOE thresholds for confidence ratings depends on the acceptable error range. Given the typically large sample sizes per year–month–grid combination, we suggest to set the ‘moderate confidence’ boundary at a relative MOE of 0.1 and the ‘high confidence’ boundary at 0.05, corresponding to an error range of 10% and 5% respectively. For example, in Figure 11, the relative MOE drops below a threshold of 0.1 at a sample size of approximately 50 for grid 668 in July 2008. At this sample size, the mean chl-a value for that specific year–month–grid combination has a 10% error range (= 0.1 on the Y-axis) with 95% confidence. Likewise, the relative MOE for grid 453 in June drops below 0.05 at a sample size of about 200, so the error range for this year–month–grid combination is 5% at that sample size.

Despite small differences between year–month–grid combinations, the relative MOE curves in Figure 11 show a comparable pattern. Therefore, general sample size thresholds for each year–month–grid combination in the entire SNS area were assessed. We randomly sampled 100 grids for each year–month combination with a sample size greater than 100. For each grid, chl-a measurements were sampled up to a size of 800, with replacement, and the relative MOE was calculated (as shown in Figure 11). An exponentially decaying function was then fitted to each year–month–grid combination:

relative MOE \(= a \times \text{sample size}^b\)

with the nonlinear least squares (nls) function in R, where the a and b parameters were estimated by the model. Based on the fitted values, the sample size corresponding to the 0.1 and 0.05 MOE thresholds was determined for each year–month–grid. The means of these sample sizes serve as the overall confidence class boundaries for the SNS area. A threshold of 0.1 corresponds to a sample size of approximately 50 (‘moderate confidence boundary’), and a threshold of 0.05 corresponds to a sample size of around 200 (‘high confidence boundary’).

Per year–month–grid confidence classes can thus be assigned based on sample size. For subsequent aggregation per year–grid or per year classes were assigned numerical values of 1 (‘low confidence’), 2 (‘moderate confidence’) or 3 (‘high confidence’), where the overall mean determines the confidence rating. A score smaller than 1.5 represents a low confidence rating, between 1.5 and 2.5 a moderate confidence rating and a score higher than 2.5 a high confidence rating, as shown on the map in Figure 12. Results aggregated per year are shown in Table 7, including numerical scores. Some edge effects are present, as discussed in chapter 4.1, but these do not significantly affect the confidence rating.

The same method can be used for the MPM and RHPM coastal areas, where we first determined the MOE thresholds based on chl-a measurements and sample sizes. As a result of the higher variation in chl-a concentrations closer to shore, the threshold of 0.1 corresponds to a sample size of approximately 60 (‘moderate boundary’), and the threshold of 0.05 corresponds to a sample size of around 260 (‘high boundary’). Consequently, confidence rating is on average lower than in the SNS area (Figure 13 and Table 7).

4.5 Conclusions new method

We propose to first aggregate chl-a data of different methods on a 10x10 km grid scale. This approach allows for adjustments based on variations in sampling effort, while its impact on the overall yearly means is expected to be minimal. Additionally, we recommend a new, more objective method for assigning confidence ratings, where boundaries can be set based on acceptable error range.

The assessment of chl-a levels is essential for understanding eutrophication in the OSPAR maritime area, and the current method of combining in-situ and Earth Observation (EO) data provides a valuable, albeit imperfect, approach to estimating chl-a concentrations. One key issue identified in this analysis is the disparity between the temporal and spatial distributions of the in-situ and EO datasets, with in-situ samples being far less representative for the entire assessment area. Despite this, the current methodology often assigns significant weight to in-situ data, even when its confidence rating is low. This overemphasis may lead to inflated growing season means in years where in-situ measurements are coincidentally biased towards high chl a values, such as algal blooms.

Our findings suggest that the current weighting approach, especially the 50:50 weighting used in most years, can introduce bias into the overall chl-a means. In years with high in-situ chl-a concentrations, particularly during the early monitoring period, this bias leads to a moderate downward trend in chl-a levels from 1998 to 2020. However, when applying the alternative method presented in this report—where in-situ and EO data are aggregated on a year–month–grid basis—the trend is stable over time. This indicates that too much weight is given to the sporadic high in-situ values in the current method, and that the alternative method may provide a more accurate reflection of overall eutrophication trends. It is noteworthy that this new approach significantly diminishes the influence of in-situ measurements in monitoring eutrophication, bringing it close to negligible levels.

Furthermore, the introduction of confidence ratings based on the relative margin of error (MOE) offers a robust framework for assessing the reliability of chl-a estimates. By linking confidence ratings to sample sizes, we have provided a method for determining ‘low,’ ‘moderate,’ and ‘high’ confidence levels across different OSPAR assessment areas. The confidence rating thresholds for these classes can be determined per OSPAR area. For instance, the coastal the Meuse and Rhine plumes have higher sample size thresholds than the Southern North Sea due to increased variability in chl-a concentrations closer to shore. However, differences in confidence rating thresholds between these assessment areas were relatively small, meaning these thresholds may be determined more generally for the entire OSPAR maritime area. They could be established conservatively, including more grids with higher chl-a variation, but this would require data from more assessment areas. As an example the MOE threshold for a high rating was set to 0.05, and the threshold for a moderate rating to 0.1 in this report, however these thresholds can be adjusted to the acceptable error range. Moreover, data was aggregated on a 10x10 km scale, as it provided the sufficient sample size for reliable chl-a estimates, while minimizing edge effects and resolution losses. Optimal grid size should also be determined for other assessment areas.

In conclusion, the revised approach reduces the bias introduced by low resolution in-situ measurements and provides a clearer picture of eutrophication trends, supporting more informed decision-making for environmental management in the OSPAR maritime area.