3.1 Random forest regressor algorithm

As the aim here is to predict a continuous numerical value for sea-surface iodide, a regression approach is taken. As discussed in the introduction, previous approaches have been made to parameterise sea-surface iodide, and the most commonly used relationships employ sea-surface temperature as the predictor variable. Here we take a different multivariate and non-parametric approach, using the computationally cheap and interpretable random forest regressor (RFR) algorithm (Breiman, 2001; Pedregosa et al., 2011).

Random forest regression is based on finding a number of decisions trees, which predict the dependent variable. All of the trees contribute to the prediction, and they are collectively referred to as a “forest”. These trees can be explained as a record of the way the algorithm has linearly traversed a subset of the training data, splitting the data into two parts at each decision point or “node” in a way that minimised the internal differences of the parts. The best split is chosen between the available variables based on an error metric (e.g. mean square error), and this process is continued until a criterion of purity is reached or a minimum number of data points are left from a split. This is essentially a classification problem. The prediction of the forest is the mean value of the prediction of all of the different decision trees, which attempts to make the results more of a regression problem. More details of this approach can be found in Friedman et al. (2009).

This approach differs from previous approaches which have individually tested proposed relationships and selected the best-performing model(s) as a parameterisation (e.g. Table A1). Here, an algorithm uses the data it is provided to build a model that gives a prediction, and therefore it is the data themselves that define the model that is used to predict new values. A key difference of this approach is also that only a subset, the training set, is used to build the model, and the rest (or withheld set) is then used to test the performance of the model. Here we use 80  % of the data for the training set and use the remaining 20 % as the withheld set (also commonly referred to as the “testing set”).

To ensure that the models built are generalisable and mitigate overfitting, the random forest approach used here artificially increases the randomness within the forest (Pedregosa et al., 2011). This is done by randomly combining single decision trees by an approach referred to as “bootstrap aggregation” or “bagging” (Breiman, 2001; Tong et al., 2003). This additional bagging approach randomly samples observations within the training dataset and so mitigates overfitting of the trees to the dataset (Friedman et al., 2009). Furthermore, a stratified sampling approach was taken. Specifically, the overall dataset was split into quartiles according to iodide concentration value, and training data were randomly selected from each quartile. This approach was used to maintain the same statistical distribution in the training/testing data as the overall dataset.

Machine learning algorithms can generally be tuned to increase performance using settings called hyperparameters. However, random forests are known to generally perform well without tuning. The default hyperparameters were therefore used here (Pedregosa et al., 2011), except for increasing the number of trees (“n_estimators”) from 10 to 500. Mean square error (MSE) was used as the criterion for evaluating each split (also referred to as a “node”). The maximum number of “features” (the ancillary variables provided to the algorithm, such as temperature or nitrate concentration) considered when looking for the best split is set to the number provided to the algorithm. The number of splits a tree is allowed to make (“max_depth”) is not restricted, and further nodes are made until leaves contain less than two samples (“min_samples_split”) and a minimum of one (“min_samples_leaf”). All the random forest models are built using bootstrapping.

3.2 Error and uncertainty calculations

Understanding the errors and uncertainties in the global iodide distribution is important due to any sensitivities to this value within the modelled Earth system. We consider three sources of error in our predictions: the “dataset selection” error due to the splitting of the dataset into training and withheld parts, the “model selection error” due to the choice of dependent variables, and the “observational error” on the iodide measurements.

To quantify the range of the dataset selection error, we construct models from 20 pseudo-random splits of the dataset into training and withheld parts. The hyperparameters and input ancillary variables are kept the same for the generation of the 20 models, so that the only difference between the models is the training dataset. These 20 models are then used to predict the withheld data. Performance metrics (e.g. root mean square error (RMSE) and average absolute prediction) can then be calculated for each model. This gives a range of 20 values, which can then be converted to a percentage range as the error. This is done by dividing the largest range in predicted values for a model by the minimum predicted value, to give a maximum value for the range of error, and then taking the smallest range in the predicted values and dividing this by the maximum value, to give a minimum value for the error range. Significant differences between the model's performance metrics would suggest important sensitivity to the training/withheld dataset splits.

We define the “model selection” error as the uncertainty resulting from the choice of input ancillary variables. A number of combinations of input variables are possible in generating the models, and each will generate a different prediction. We quantify this error as the difference in performance against the withheld dataset and prediction value (e.g. average global value). Similarly to our calculation of the dataset selection error, this can be converted to percentage error by considering the range in these values and dividing them by minimum and maximum values.

For the observational error we refer to Chance et al. (2019a), who provide individual error estimates for each of the iodide observations in the data compilation. Over half (51 %) of the data points have an error of 5 % or less, and a further ∼25 % have an uncertainty in the range of 5 %–10 %. We therefore use a value of 10 % as a conservative estimate of the observational error.

3.3 Outlier identification and removal

Our dataset consists of values for ancillary variables and iodide concentration for all of the 1342 measurement locations in the observational dataset (Sect. 2). As discussed in Sect. 3.1, we split this dataset into two parts: (i) a training set for use in building and optimising models, and (ii) a withheld set to evaluate the models built. Particular care was taken to ensure the withheld and training datasets were representative of the entire dataset in the way the models are built, therefore improving performance and generalisability to unseen data (see Sect. 3.1).

We take a random forest regressor model built with variables that were intuitively assumed to give a reasonable ability to differentiate the observations (using depth, temperature, and salinity as the independent variables – abbreviated to “RFR(DEPTH+TEMP+SAL)” following Table 1). The RFR(DEPTH+TEMP+SAL) model was then used to explore the variation of error in the predictions using the dataset selection error approach described in Sect. 3.2. This builds multiple versions of the same model with different splits of training and test data and yields a distribution of root mean square error in the predicted iodide for withheld data as summarised in the final column of Table 2 and shown graphically in Appendix Fig. A1.

We define outliers here as values greater than the third quartile plus 1.5 times the interquartile range (Frigge et al., 1989). Removing these 49 values categorised as outliers (>309.5 nM) leads to a vast improvement in the RMSE error in the ensemble prediction from 95.1 to 37.6 nM (Table 2). This is shown graphically in Appendix Fig. A1, with the other subsets of the data explored (Table 2). This demonstrates that the high values are not well enough represented by the dataset to be able to be captured by the RFR approach. The removal of these high values from the dataset can also be justified as the driver for these concentrations is not yet well understood (Chance et al., 2014, 2019c; Cutter et al., 2018).

Removing these outliers reduces RMSE in the prediction with the 20 independent model builds from 48.2 to 2.3 nM (third quartile–first quartile). Once these outliers are excluded, more modest changes in average RMSE are then seen if models are built only using coastal or non-coastal data. Figure A1 also shows this is seen when removing lower salinity data (“Salinity ≥30 PSU and no outliers”), which is indicative of estuarine water. This highlights the strength in this approach's ability to predict iodide in different biogeochemical regions (i.e. not just coastal or non-coastal locations).

An additional removal of a single dataset of 19 observations from the Skagerrak strait (Truesdale et al., 2003) was made due to it exerting a disproportionate influence on iodide prediction in high northern latitudes ($\ge \mathrm{65}{}^{\circ }$ N), an area that is almost entirely unconstrained by local observations. We note that the Skagerrak is relatively unusual oceanographically, being an estuarine location with high ship traffic, and is considered unlikely to be an analogue for iodine speciation in the Arctic. This is decision is discussed further in the Appendix (Sect. A1), and the predictions made including this dataset are also provided in the shared output (Sect. 5).

From here, only the 1293 observational points excluding outliers and the data from the Skagerrak strait (Truesdale et al., 2003) are used.

3.4 Selection of ancillary variables and building an ensemble prediction

To decide which ancillary variables (temperature, salinity, etc.; see Table 1 and Sect. 2) should be used to predict sea-surface iodide concentration, RFR models were built and evaluated with different combinations of variables. Thirty eight combinations were considered (see first column of Appendix Table A2).

The top 20 performing models, based on their root mean square error against the withheld data, are plotted in Fig. 2, alongside existing parameterisations. The standard deviation for all predicted values is also shown to illustrate variation in the predictions. A complete list of the performance and of all models built here and their performance is given in the appendix (Table A2).

Figure 2Random forest regression (RFR) model performance (root mean square error (RMSE), blue) against the withheld data for the top 20 models on the left-hand y axis, along with values from the parameterisations from Chance et al. (2014) and MacDonald et al. (2014). The right-hand y axis is standard deviation of the prediction for the withheld data (orange). The top 10 performing models and the two exiting parameterisations considered here (Chance et al., 2014; MacDonald et al., 2014) are shown in bold. Parameterisations are ordered by their RMSE. Abbreviations are given in Table 1.

Download

The RMSE values in Fig. 2 show the increased skill present in the new predictions compared to the existing parameterisations. The RMSE improves from the 75.3 and 50.2 nM found for the Chance et al. (2014) and MacDonald et al. (2014) parameterisations, respectively, to 33.2–37.4 nM for the top 10 models created here. Only modest gains are seen in RMSE between models with three variables or more.

The best-performing model in the list is only marginally better than the 10th best-performing one; therefore there is not an obvious “best” performing set of ancillary variables. Thus going forward we use an ensemble prediction approach based on the mean value from an ensemble of the 10 top-performing models.

3.5 Model descriptions

Unlike many machine learning approaches, the random forest regressor algorithm is interpretable. The decision trees can be visualised to explain the main features driving the splits. Figure 3 shows schematically the whole regression approach taken here. Panel (a) shows single trees, of which 500 are built with the same input variables and then combined into forest (b). Then this forest is combined with the nine other top-performing models (made from different combinations of ancillary variables) to make an ensemble (c). The 10 predictions of (c) are then arithmetically averaged into a single prediction, which thus includes the predictions of 5000 trees with 10 different combinations of input variables. In Fig. 3a, the colour of a limb or “branch” following a node is given by the variable driving that split within the training dataset. For Fig. 3b and c it shows the percent of times that a variable drives that node within the forest. The value of the ancillary variable that sets the split is shown inside the circle (a, b, c). The thickness of the branch scales to the throughput of training dataset samples contained within that split. The trees are shown to a depth of five nodes for aesthetic reasons and due to increased divergence of the trees within a forest the deeper you go. However the trees themselves are unlimited in the depth they can reach.

Figure 3Schematic illustration of how (a) multiple decision trees are combined into (b) a forest and then combined into (c) a further 10-member ensemble. (a) shows individual trees in a forest. (b) represent a forest of 500 trees as a single figurative tree. (c) shows the 10 forests of 500 trees combined into a single prediction. The branches in plots (a)–(c) are coloured by the percentage of the decisions at a given node that are driven by a given variable. That value within the circle gives the value of the main ancillary variable driving a split. The thickness of branches gives the throughput of the dataset through a given node for single trees (a) or the average for plots of forests (b, c). The 10 forests shown as thumbnails in panel (c) are also shown in larger form in Appendix Fig. A5. Variable names are coloured as per the following coloured text: temperature (blue, ^∘C), depth (orange, metres), chlorophyll a (green, mg m⁻³), salinity (pink, PSU), nitrate (brown, µg m⁻³), mixed layer depth (MLD; purple, metres), phosphate (red, µg m⁻³), and shortwave radiation (grey, W m⁻²).

Download

The first and larger splits in the data at decision nodes in the models can be simply read, which can provide understanding of the main variables driving the initial and largest splits in the prediction. For all models in the ensemble, the initial split is driven by temperature, with a split occurring at around 21.1 ^∘C (with a standard deviation of 1.2 ^∘C). The data are then split by two further nodes from this, a left- and right-hand split (e.g. Fig. 3b). If depth or temperature is present as a variable, then they drive the majority of the next splits. If depth is not present as a variable, then either nitrate or mixed layer depth (MLD) is the most common variable to dictate the split in the data at the next node in the tree. Thus a qualitative way of interpreting the initial splits of the dataset would be to say that the model is primarily differentiating between warmer and shallower locations.

4.1 Prediction of iodide at observational locations

Figure 4 shows a point-by-point comparison between parameterised and observed iodide for the entire dataset, the withheld dataset, and the withheld coastal dataset and withheld non-coastal dataset. Predictions are shown for the ensemble random forest regressor approach described here and for both the Chance et al. (2014) and MacDonald et al. (2014) parameterisations. The root mean square errors of observed and predicted values are given in Fig. 4 and in Table 3.

Figure 4Regression plots showing comparisons between predicted values and observations in the entire (blue, N=1293) and withheld data (orange, N=259), withheld data classed as coastal (green, N=157), and the withheld data classed as non-coastal (pink, N=102). Solid lines give the orthogonal distance regression line of best fit. The dashed grey line gives the 1:1 line. Root mean square error for each line is annotated by subplot in nanomolar (nM).

Download

Table 3Statistics for observations and predictions by the ensemble prediction (RFR(ensemble)) and existing parameterisations against the entire dataset of observations. The root mean square error is shown against the entire and the withheld data.

Download Print Version | Download XLSX

The new ensemble prediction is the best performing, with a lower RMSE (35 nM) compared to the existing parameterisations (75 and 50 nM for the MacDonald et al., 2014, and Chance et al., 2014, respectively) for the withheld data. Both the new parameterisation and Chance et al. (2014) parameterisation are relatively unbiased (both have best-fit-line slopes of 0.84, against the withheld data), but the new parameterisation shows less noise than Chance et al. (2014). MacDonald et al. (2014) shows a significant low bias and significant noise. The improved skill from the RFR ensemble is consistent for both coastal and non-coastal observations.

Figure 5 shows comparisons between the probability distribution functions (PDFs) of the observed iodide and the predictions, together with the PDFs of the biases for the entire, coastal, and non-coastal withheld datasets. The PDF of the new parameterisation shows the greatest similarity to the observations. The PDF from Chance et al. (2014) shows a similar range to the observations and structure to the observations, whereas the PDF from MacDonald et al. (2014) shows again a significant underestimate. The bias plots show the new predictions are generally clustered around zero with a relatively narrow peak. Chance et al. (2014) is again roughly clustered around zero but shows a wider peak. The largest biases are found from MacDonald et al. (2014), which systematically underestimates observed iodide concentrations.

Figure 5Probability density function (bars) and Gaussian kernel density (lines) estimate of observations and predicted concentrations (left), and bias (right, model minus observations) in entire withheld dataset (upper, N=259), the withheld coastal dataset (middle, N=157), and the withheld non-coastal dataset (lower, N=102).

Download

The dataset selection error range, which shows the influence of the choice of how the dataset is split into training and withheld data on model prediction, is described in Sect. 3.2. Within the 20-member ensemble of different testing/withdrawn choices, the average variation in RMSE was 8.4 nM (5.9–11.02 nM), and in the range of average predicted values it was 6.1 nM (5.4–6.6 nM). This translates to a percentage error range of 13.9 %–36.3 % on the RMSE and 5.4 %–7.3 % on the average predicted value.

The model selection error, which is the influence of the different independent variables used, is described in Sect. 3.2. The difference in the average prediction of the 10 members of the ensemble is 1.8 nM (with a range of average prediction from 96.0 to 97.8 nM), and the range of the difference in model performance is 3.9 nM (33.2–37.2 nM). As a percentage this model selection translates to a percentage uncertainty on the RMSE of 10.6 %–11.9 % and on the average of 1.8 %–1.9 %.

The dataset selection and model selection compare to an error on the observations of ∼10 %. Uncertainty from dataset selection has a far greater effect on the prediction error than model selection. This is expected due to the small dataset size. The combined error in the prediction (dataset selection + model selection error) is either comparable to (7.2 %–9.2 % in terms of average prediction) or greater (24 %–48 % in terms of RMSE) than the observational error.

From this analysis we have shown that the new ensemble RFR model performs significantly better than those currently in the literature. We now turn to explore the predicted global distribution of sea-surface iodide using our ensemble model.

4.2 Global sea-surface iodide distribution

From the ensemble prediction system we calculate monthly global grids ($\mathrm{0.125}{}^{\circ }\times \mathrm{0.125}{}^{\circ }$, $\sim \mathrm{12.5}\mathrm{km}\times \mathrm{12.5}\mathrm{km}$) of sea-surface iodide using the gridded ancillary data (Sect. 2). The annual average spatial predictions are shown in Fig. 6 with the observations overlaid in circles. Similar to previous work, annual average maximum concentrations of 220 nM are found in tropical and coastal regions (e.g. Oceania and in the Caribbean/Gulf of Mexico), with the lowest concentrations in mid-latitude waters (22.4 nM). Seasonal variability is also seen within the monthly prediction (Appendix Fig. A3). However, this spatial and temporal variability is not well constrained by observations. For example, some of the highest concentrations are predicted for the South China Sea, a region without any observations (Fig. 6). Some features visible in the concentration field appear to be associated with deep bathymetric features (e.g. the higher concentrations over the Mid-Atlantic Ridge – Fig. 6), even though a physical explanation for such a link seems unlikely.

Figure 6Annual average predicted sea-surface iodide by the 10-member ensemble of models (RFR(Ensemble)), overlaid with iodide observations from Chance et al. (2019a) without outliers. Outliers are defined here as values greater than the third quartile plus 1.5 times the interquartile range (Frigge et al., 1989). Only locations that are entirely water are included in the spatial average. Earth raster and vector map data used are freely available from Natural Earth (http://www.naturalearthdata.com/, last access: 23 July 2019).

Summary statistics on the global predictions are shown in Table 4. These show that, as for comparisons at the observed locations (Sect. 4.1), the ensemble prediction is broadly in between the two existing parameters. The new ensemble model predicts a mean value of 106 nM (with members ranging from 102.3 to 108.8 nM), with predicted values from existing parameterisations ranging from 58.9 (MacDonald et al., 2014) to 128.1 nM (Chance et al., 2014).

Table 4Statistics on predicted global annual sea-surface values from new and existing parameters at a horizontal resolution of $\mathrm{0.125}{}^{\circ }\times \mathrm{0.125}{}^{\circ }$. Existing parameterisations (Chance et al., 2014; MacDonald et al., 2014) and the ensemble prediction are shown in bold.

Download Print Version | Download XLSX

The annual latitudinal average of these fields, together with predictions from Chance et al. (2014) and MacDonald et al. (2014), and the observations are shown in Fig. 7. Far greater structure is seen compared to the two existing parameterisations (Fig. 7) due to the multivariate and non-parametric ensemble approach used here. All parameterisations capture the broad observed feature of decreasing iodide from lower to higher latitude. The new predicted values lie between Chance et al. (2014) and MacDonald et al. (2014) in the tropics; however, within the polar regions, the new prediction is significantly higher than both of the previous parameterisations. The lower concentrations in the predicted values from MacDonald et al. (2014) for most of the global sea surface are clear.

Figure 7Predicted annual average sea-surface iodide plotted against latitude (lines), overlaid with observed concentrations (diamonds). Solid lines give mean values and shaded regions give (±) the average standard deviation. The standard deviation is the monthly standard deviation across a latitude between all 10 ensemble members (RFR(Ensemble)) or within a single prediction for existing parameterisations (Chance et al., 2014; MacDonald et al., 2014). Filled diamonds show non-coastal observations and non-filled ones show coastal values. Extent of x axis is shown for grid boxes that are entirely water.

Download

The range of dataset error is found for the 20 models with different training data splits, as described in Sect. 3.2. This gives an uncertainty in the form of a average range in predicted global mean surface iodide for all of the multiple builds of ensemble members of 4.0 nM (2.8–5.0) compared to a annual mean prediction of 106 nM. This maximum and minimum of this range in predicted values can then be divided by the minimum and maximum predicted global mean surface iodide values (98 and 109.3 nM, respectively) to give percent range of 2.5 % to 5.1 %. This is lower than that calculated for the individual locations of observations (Sect. 4.1) due to large global areas being similar in chemical and physical regimes compared to the subset of sampled locations within the observations.

The model selection error due to variability within ensembles' 10 members, generated with different independent variables, gives a global average surface concentration between 102.3 and 108.8 nM. This range in prediction gives a model selection error of 6.45 nM, which equates to 6.0 %–6.3 %. Like with the global uncertainty from dataset selection, the global value would be expected to be lower than the uncertainty at the specific locations of the observations (Sect. 4.1) due to the more homogeneous nature of the predicted areas. However, a greater variation is seen from different model predictions than within predictions for the observation locations. This highlights the importance of the different ancillary variables considered here and also therefore the strength gained from the ensemble approach taken here.

Within members of the ensemble, variation is modest except for two ensemble members which diverge north of $\ge \mathrm{65}{}^{\circ }$ N (Appendix Fig. A2). As noted earlier (Sect. 2), the values in this region are very poorly constrained by the observational dataset (Fig. 6).

In addition to the three errors described above, we also attempt to gain an understanding of the spatial uncertainty in the 10-model ensemble prediction. We do this by calculating the differences in the predicted spatial fields from the 10 ensemble members. Figure 8 shows the monthly average of the standard deviation of the 10 model ensemble as a percentage of the annual mean of the ensemble prediction. This is also shown in absolute terms in the Appendix (Fig. A4). Relative uncertainties are largest at the poles where predicted concentrations are lowest and where (at least in the Northern Hemisphere) very few observations are available to constrain the system. The southern oceans show an distinct pattern, where values close to coastal Antarctica appear well constrained but values further north appear poorly constrained.

Figure 8Annual average spatial percent uncertainty in predicted sea-surface iodide for the ensemble of models. Percent spatial uncertainty was calculated as the standard deviation in monthly average values for all models, divided by the annual mean. Values are limited to 25 % for contrast, but the maximum plotted value is 77 % in the northern high latitudes. Only locations that are entirely water are included in the spatial average. Earth raster and vector map data used are freely available from Natural Earth (http://www.naturalearthdata.com/, last access: 23 July 2019).