Mapping the global distribution of lead and its isotopes in seawater with explainable machine learning

Lead (Pb) and its isotopes are powerful tools for studying the pathways of Pb pollution from land to sea and, simultaneously, investigating biogeochemical processes in the ocean. However, the scarcity and sparsity of in situ measurements of Pb concentrations and isotope compositions do not allow for a comprehensive understanding of Pb pollution pathways and biogeochemical cycling on a global scale. Here, we present three machine learning models developed to map seawater Pb concentrations and isotope compositions, leveraging the global GEOTRACES dataset as well as historical data. The models use climatologies of oceanographic and atmospheric variables as features from which to predict Pb concentrations, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb. Using SHapley Additive exPlanations (SHAP), we found that seawater temperature, atmospheric dust, atmospheric black carbon, and salinity are the most important features for predicting Pb concentrations. Dissolved oxygen concentration, salinity, temperature, and atmospheric dust are the most important features for predicting ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, atmospheric black carbon and dust, seawater temperature, and surface chlorophyll a for ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb. In line with observations, our model outputs show that (i) the surface Indian Ocean has the highest levels of pollution, (ii) pollution from previous decades is sinking in the North Atlantic and Pacific oceans, and (iii) waters characterised by highly anthropogenic Pb isotope fingerprints are spreading from the Southern Ocean throughout the Southern Hemisphere at intermediate depths. By analysing the uncertainty associated with our maps, we identified the Southern Ocean as the key area to prioritise in future sampling campaigns. Our datasets, models, and their outputs, in the forms of Pb concentrations, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb climatologies, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb climatologies, are made freely available to the community by Olivelli et al. (2024a; https://doi.org/10.5281/zenodo.14261154) and Olivelli (2025; https://doi.org/10.5281/zenodo.15355008).

1 Introduction

The low concentrations at which lead (Pb) is present in the ocean (on the order of parts per trillion) and the ubiquity of Pb contamination during sampling and sample processing did not allow for successful measurements of seawater Pb concentration and isotope composition until the late 1970s (Schaule and Patterson, 1981, 1983). Indeed, since the early 1900s, human activities have caused an increase in Pb concentrations in the environment, both on land and in the ocean, which has altered the background levels and the biogeochemical cycling of Pb (Boyle et al., 2014; Nriagu, 1979; Nriagu and Pacyna, 1988; Pacyna and Pacyna, 2001).

This increase in Pb concentrations has been accompanied by a shift in the Pb isotope compositions of seawater, reflecting the predominance of anthropogenically sourced Pb over naturally sourced Pb (Boyle et al., 2014). In fact, the relative abundances of stable Pb isotopes (²⁰⁴Pb, ²⁰⁶Pb, ²⁰⁷Pb, and ²⁰⁸Pb), expressed as ratios such as ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, can be used to identify and quantify the contributions of the different sources of Pb (Reuer and Weiss, 2002).

In the years and decades that followed the first successful measurements of Pb and its isotopes in seawater, sampling efforts mainly concentrated on the North Atlantic Ocean (Alleman et al., 1999; Boyle et al., 1986; Helmers et al., 1991; Pohl et al., 1993; Shen and Boyle, 1988; Sherrell et al., 1992; Véron et al., 1994, 1998, 1999). This basin, surrounded by the early developed economies of North America and western Europe, initially saw a sharp increase in pollution followed by a decrease thanks to environmental policies that phased out and banned the use of leaded petrol starting from the 1970s (Bridgestock et al., 2016; Weiss et al., 2003).

A major breakthrough in the understanding of the marine biogeochemistry and pollution of Pb on a global scale came with the GEOTRACES international marine geochemistry programme, which has run since 2006 and includes Pb concentrations and isotope compositions as key parameters of interest for measurement on all of its cruise sections (GEOTRACES Planning Group, 2006). However, despite these great efforts, the majority of the world's oceans remain unsampled for Pb concentrations and isotope compositions, especially in the Southern Hemisphere (Fig. 1).

Figure 1 Geographical distribution of all samples in the Pb concentration (a), ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb (c), and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb (e) datasets. Panels (a), (c), and (e) include all samples from the top 100 m. Panels (b), (d), and (f) represent the frequency distribution of sampling depths for panels (a), (c), and (e), respectively.

In this context, modelling studies represent a powerful tool for expanding our knowledge of the distribution and biogeochemical cycle of Pb, complementary to in situ observations. While deterministic models require a good understanding of the processes at play in the biogeochemical cycle of Pb, machine learning (ML) models do not require such prior knowledge as they are data-driven (Glover et al., 2011). Therefore, the latter can be regarded as a computationally efficient way of investigating the distribution of Pb and its isotopes in the global ocean without the need to explicitly parameterise causal relationships derived from fundamental knowledge of the biogeochemical cycle of Pb.

In recent years, ML algorithms have been employed to investigate several topics within chemical oceanography, including the distribution and cycling of tracers, such as iron (Huang et al., 2022), zinc (Roshan et al., 2018), copper (Roshan et al., 2020), barium (Mete et al., 2023), iodide (Sherwen et al., 2019), and nitrogen isotopes (Rafter et al., 2019). A thorough overview of the use and potential of ML in the ocean sciences is provided in Sonnewald et al. (2021).

In this study, we leverage the growing GEOTRACES dataset, in combination with climatologies of several oceanographic and atmospheric variables, to map the distribution of Pb and its isotopes in the global ocean using the ML XGBoost (eXtreme Gradient Boosting) regression algorithm. With this approach, we aim (i) to produce three-dimensional global maps of Pb concentrations and isotope compositions (${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb) for use as a benchmark for future levels of pollution and (ii) to identify areas where to concentrate sampling efforts in the upcoming years and decades by analysing the uncertainty of our model outputs. We make the resulting data products available to the community in the form of mean global climatologies and envisage updates to the models' outputs as new observations become available.

2 Materials and methods

Regression ML models require training on existing observations of a target variable (to be predicted by the model) alongside a number of features (i.e. ancillary variables) in order to make accurate predictions. Therefore, three different datasets were compiled to develop three separate ML models to map the global distribution of Pb concentrations, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, respectively. Each dataset includes information on 14 geographical, oceanographic, and atmospheric features (Table 1) and a target variable (either Pb concentration, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, or ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb). The features included in the datasets were chosen due to their proven predictive power in other ML studies focusing on ocean chemistry (Huang et al., 2022; Mete et al., 2023; Roshan et al., 2018, 2020; Sherwen et al., 2019) and their likely connection to the biogeochemical cycle of Pb.

The three datasets and models were built following the same procedure explained in Sect. 2.1 and 2.2. All data preparation, modelling, and analyses were done in Python 3.11.5 in a Linux operating system.

Table 1 List of the features used for dataset compilation and model development for Pb concentrations, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb.

^a Copernicus Marine Environment Monitoring Service. ^b Aerosol optical depth. ^c Copernicus Atmosphere Monitoring Service. ^d Variables used for dataset compilation and initial versions of the models but dropped for the final models (see more in Sect. S1 in the Supplement).

2.1 Data sources

2.1.1 Pb concentrations and isotope compositions

Observations of seawater dissolved and total dissolvable Pb concentrations, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb included in the GEOTRACES Intermediate Data Product 2021v2 (IDP2021v2; GEOTRACES Intermediate Data Product Group, 2023) form the basis of this study. Only GEOTRACES IDP2021v2 data with a quality control flag equal to 1 (equivalent to good data) were included in our Pb concentration, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb datasets. Of these good data points, eight samples from the southern sector of the Indian Ocean and the Southern Ocean were not included in the Pb concentration dataset as they had unusually high values compared to those collected at the surrounding locations. Additionally, one sample from the Arctic Ocean was not included in our dataset as it had a negative Pb concentration (−0.05 pmol kg⁻¹).

Additional data were included in our dataset to increase the number of observations that were available; these encompass published (Chen et al., 2023b; Chien et al., 2017; Lanning et al., 2023; Olivelli et al., 2023, 2024b) and unpublished studies (from the MAGIC laboratories at Imperial College London). As previous studies have shown that dissolved and total dissolvable Pb concentrations and isotope compositions are comparable in the open ocean (Bridgestock et al., 2018; Olivelli et al., 2023, 2024b; Schlosser et al., 2019), data for dissolved and total dissolvable Pb concentrations, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb were pooled together for the purpose of this work.

2.1.2 Oceanographic and atmospheric features

Global gridded climatologies of seawater temperature, salinity, density anomalies (σ₀), nitrate, phosphate, silicate, oxygen, apparent oxygen utilisation (AOU), and mixed layer depth (MLD) were obtained from the World Ocean Atlas 2018 product (WOA18; Garcia et al., 2019). The WOA18 has a 1° × 1° spatial resolution and 102 depth levels (spaced every 5 m between 0 and 100 m, every 25 m between 100 and 500 m, every 50 m between 500 and 2000 m, and every 100 m between 2000 and 5500 m). Chlorophyll-a concentrations in the surface ocean were obtained from the Copernicus Marine Service Global Ocean Biogeochemistry Hindcast product (0.25° × 0.25° spatial resolution) as monthly averages for the period between January 2005 and December 2017 and averaged over the entire period to create a climatology.

Dust aerosol optical depth (AOD) and black carbon AOD were included as most Pb enters the ocean via atmospheric transport. Dust was found to be an important source of natural Pb in the surface ocean in several regions, including the Atlantic Ocean (Bridgestock et al., 2016; Olivelli et al., 2023) and the Red Sea (Benaltabet et al., 2020). Black carbon, on the other hand, was chosen as an indicator of Pb emissions from industrial and high-temperature activities, which are known sources of anthropogenic Pb in the environment (Nriagu and Pacyna, 1988; Pacyna and Pacyna, 2001). Data for the years between January 2005 and December 2017 were obtained from CAMS global reanalysis (ECMWF Atmospheric Composition Reanalysis 4; EAC4) monthly averaged fields with a 0.75° × 0.75° spatial resolution, and they were averaged over the whole period to obtain a climatology.

2.1.3 Data handling and dataset preparation

The WOA18 grid was chosen as a reference for all feature and target variables in this study. Therefore, Pb concentration, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, chlorophyll a, dust AOD, and black carbon AOD values, which were either obtained as point observations or gridded at higher resolution, were re-gridded according to it. Lead data were binned to the 1° × 1° grid cell the samples were collected in and assigned to the depth level closest to the depth at which the samples were collected. All samples collected at depths lower than 5500 m were discarded, as 5500 m is the deepest depth level considered in the WOA18. The distance from the bottom depth was calculated as the difference between the bottom depth value reported for the sampling station and the WOA18 depth level to which the observation was assigned. Chlorophyll a, dust AOD, and black carbon AOD were re-gridded by assigning each 1° × 1° grid cell the feature value of the respective product cell located closest to the centre of the WOA18 gridded cell.

All samples in each 1° × 1° column were given the same dust AOD, black carbon AOD, surface chlorophyll a, and mixed layer depth values, regardless of their depths. This was done to assess whether characteristics of the surface ocean have an important impact throughout the water column. These features represent atmospheric inputs (dust AOD and black carbon AOD), surface primary production as a proxy for sinking particulate matter in the water column (surface chlorophyll a), and the thickness of the ocean layer that is directly influenced by the atmosphere (mixed layer depth).

In total, the Pb concentration dataset includes 9920 observations, the ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb dataset 2014, and the ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb dataset 2010 (Fig. 1). In the WOA18 grid, 9920 observations cover ∼ 0.3 % of the total gridded ocean volume, while 2014 and 2010 cover ∼ 0.06 %.

In the initial phase of model development, latitude and longitude were included as features and deconvoluted in a three-dimensional coordinate system to allow for continuity (Gregor et al., 2017; Huang et al., 2022). However, when including coordinates as features, the global maps of Pb concentrations and isotope compositions showed unrealistic spatial artefacts (see Sect. S1). For this reason, coordinates were excluded as features for further model development and are not included in the results presented here.

2.2 Model development and validation

The XGBoost non-linear regression algorithm (Chen and Guestrin, 2016) was chosen to develop three separate models for Pb concentrations, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb. Indeed, evidence shows that tree-based models consistently outperform deep-learning models in tabular-style datasets (Grinsztajn et al., 2022). Moreover, in recent years, XGBoost has been proven to perform well in a variety of applications, from finance to medicine and including the Earth sciences (Biass et al., 2022; Xie et al., 2024; Ye et al., 2023). XGBoost is an ensemble algorithm that consists of individual decision trees built in a sequential manner. In brief, each successive tree aims to reduce the errors perpetrated by its predecessor. The final prediction made by the model equals the weighted sum of the predictions made by all trees in the ensemble. XGBoost was preferred to the random forest algorithm because the sequential building of trees in XGBoost allows the algorithm to focus on poorly understood fields and because random forest is more prone to overfitting than XGBoost. Moreover, the random forest algorithm was found to perform worse than XGBoost for the Pb concentration and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb models and comparably for the ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb model (Sect. S2).

In addition to performing well in a variety of applications, tree-based models offer a higher level of interpretability compared to deep-learning and linear models (Lundberg et al., 2020). Of the different approaches that have been used to explain predictions of XGBoost models, SHapley Additive exPlanations (SHAP) has stood out in recent years as a unified approach to explain local predictions and gain a global understanding of a model's structure (Lundberg et al., 2020; Lundberg and Lee, 2017; Molnar, 2022).

The SHAP approach is based on game theory (i.e. features are considered to be players that contribute to a prediction) and allows us to interpret the importance of each feature of a given prediction by calculating the marginal contribution of each feature across all possible permutations. For each sample, a positive (negative) SHAP value indicates that a specific feature contributes towards increasing (decreasing) the final predicted value. Additionally, features with a larger mean absolute SHAP value are considered to contribute more to a given prediction compared to features with smaller mean absolute SHAP values.

In recent years, several studies in the Earth sciences have used SHAP values to interpret predictions made by tree-based and deep-learning models, providing insightful findings on atmospheric pollution (Hou et al., 2022; Qin et al., 2022; Stirnberg et al., 2021), ocean dynamics (Clare et al., 2022), and vegetation vulnerability (Biass et al., 2022). Here, we used the TreeExplainer method from the shap library in Python (Lundberg et al., 2020; Lundberg and Lee, 2017) to compute SHAP values and evaluate the contribution of each feature to predicted Pb concentrations, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb. Interactions between features, which are defined as additional combined feature effects after accounting for individual feature effects, were calculated using SHAP interaction values.

The first step of model development consisted of withholding a cruise transect in the Atlantic Ocean from the Pb concentration, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb datasets. This group of withheld samples, identified as the “geographic test set” (Fig. S1), included 216 samples for the Pb concentration model and 26 samples for the Pb isotope ratio models and was used to mimic areas of the ocean where in situ samples have never been collected and assess the ability of the model to generalise to these. The remainder of each of the three datasets was split into a training test consisting of 80 % of the data and a randomly selected test set consisting of 20 % of the data.

To identify the best model architectures, hyperparameters were tuned using 5-fold cross-validation on the training set using the grid search method. Cross-validation is a technique that allows us to evaluate the performance of a machine learning model during the hyperparameter tuning phase. It consists of splitting the training dataset into a specified number of subsets (k, in this case equal to 5) and using k − 1 subsets for training and 1 subset for testing. The training and testing process is repeated k times (i.e. k-fold cross-validation), and the performance of each hyperparameter configuration is calculated by averaging the scores obtained for each fold. Hyperparameters are a set of parameters that control the learning process of a ML model and are tuned to obtain the most optimal model performance. The grid search method consists of specifying a set of possible values for each hyperparameter and subsequently training and testing the models built with each unique combination of hyperparameter values. The hyperparameters tuned for our Pb concentrations, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb include the number of decision trees built and boosted (“n_estimators”), the rate at which the model learns information (“learning_rate”), the maximum number of split nodes in a tree (“max_depth”), the minimum number of samples a node must represent in order to be split further (“min_child_weight”), the percentage of features used to construct each tree (“colsample_bytree”), and the L1 and L2 regularisation terms (“reg_alpha” and “reg_lambda”, respectively). The two regularisation terms, L1 and L2, were only tuned for the Pb isotope ratio models in order to reduce the computational costs of hyperparameter tuning for the Pb concentration model, which was built using a much larger dataset than the ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb models. As the target variables of the three models (Pb concentrations, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb) were not uniformly distributed, a least-squares loss function was chosen as the preferred method for model optimisation, as squaring the error penalises the model more when the size of the error increases.

Different performance metrics were used to assess model performance during cross-validation and testing on the geographic and random test sets. These include the root mean square error (RMSE), mean average percentage error (MAPE), and coefficient of determination (R²) and are calculated as follows:

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n(y_{\mathrm{o}}-y_{\mathrm{p}}{)}^{\mathrm{2}}},\text{(2)}& {\displaystyle }\mathrm{MAPE}=\mathrm{100}{\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n\left|{\displaystyle \frac{y_{\mathrm{o}}-y_{\mathrm{p}}}{y_{\mathrm{o}}}}\right|,\text{(3)}& {\displaystyle }{R}^{\mathrm{2}}=\mathrm{1}-{\displaystyle \frac{\sum _{i=\mathrm{1}}^n(y_{\mathrm{o}}-y_{\mathrm{p}}{)}^{\mathrm{2}}}{\sum _{i=\mathrm{1}}^n(y_{\mathrm{o}}-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}}},\end{array}$$

where n is the number of samples; y_o and y_p are the observed and predicted values, respectively; and $\stackrel{\mathrm{‾}}{y}$ is the mean of the observed values.

Root mean square error was chosen as the primary evaluation metric for the two isotope ratio models. Indeed, RMSE is reported in the same unit as the target variables (in this case unitless), which makes it more intuitive to interpret the results. On the other hand, MAPE was chosen as the preferred metric with which to evaluate the performance of the Pb concentration models. This was done as Pb concentrations in the dataset vary between 0.10 and 95.00 pmol kg⁻¹, with an average value of 17.33 ± 13.05 (1 SD, n= 9920). While MAPE might be harder to interpret, as it is reported as a percentage rather than an absolute value with the same unit of measurement as the target variable, minimising the percentage error ensures that basins and locations where Pb concentrations are low are not overlooked by the model.

3 Results and discussion

3.1 Model development and performance

During hyperparameter tuning, we assessed model configurations by changing hyperparameter grid spacing and increasing granularity to get as close as possible to the best combination achievable. In the final phase of the tuning, 1280 different configurations were assessed for Pb concentrations, with 11 520 for both ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb. The hyperparameter space explored and the best hyperparameters identified for each of the three models are reported in Table 2. The variations in MAPE and RMSE values for different combinations of hyperparameters are visible in Fig. S2.

Table 2 Hyperparameter space explored for the XGBoost regression models. Bold values identify the combination of hyperparameters that returned the best model performance.

The mean MAPE values of the Pb concentration models calculated on the cross-validation set, random test set, and geographic test set were 22.7 ± 2.2 %, 20.7 ± 1.5 %, and 19.5 ± 4.1 % (n = 1280, 2 SD), respectively. The 10 best hyperparameter configurations for the Pb concentration model shared the same learning rate (0.01) and number of columns subsampled by each tree (0.7). The maximum differences in MAPE values between these 10 models were 0.05 % in the cross-validation set, 0.17 % in the random test set, and 0.82 % in the geographic test set.

The ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb models had mean RMSE values of 0.008 ± 0.004 (cross-validation), 0.007 ± 0.003 (random test set), and 0.006 ± 0.002 (geographic test set; n= 11 520, 2 SD). The 20 best-performing hyperparameter configurations shared the same learning rate (0.01), L1 regularisation (0), and number of columns subsampled by each tree (0.7). These models have maximum differences in RMSE values of 0.000024 in the cross-validation set, 0.000194 in the random test set, and 0.000298 in the geographic test set. Lastly, the ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb models had mean RMSE values of 0.007 ± 0.003 (cross-validation), 0.007 ± 0.003 (random test set), and 0.008 ± 0.003 (geographic test set; n= 11 520, 2 SD). The five best hyperparameter configurations have the same learning rate (0.01), number of columns subsampled by each tree (0.5), and L1 regularisation (0). The maximum differences in RMSE values between these models were 0.000003 in the cross-validation set, 0.000046 in the random test set, and 0.000069 in the geographic test set.

Overall, given the very comparable performances of the best hyperparameter configurations in the cross-validation set, random test set, and geographic test set for each of the three models described above, we built the final Pb concentration, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb models using the hyperparameter configuration that performed best in cross-validation.

The predictive performance of the final models was assessed on both the random and geographic test sets and is reported in Table 3. For all three models, the agreement between observed and predicted values (R² score) is higher for the random test set than for the geographic test set. However, MAPE and RMSE are not always better for the random test set than for the geographic one, indicating that the model is able to generalise well without drastically reducing its performance. Indeed, the Pb concentration model has a lower MAPE for the geographic test set (18.3 %) than for the random one (20.2 %), and the ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb model has a lower RMSE for the geographic set (0.005) than for the random test set (0.006).

Table 3 Performance of the best Pb concentration, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb models.

Model performance on the random test set was also assessed by splitting the data between ocean basins to identify which areas are better reproduced by the models (Fig. 2). The Pb concentration model achieves the highest R² in the Pacific Ocean (0.93) and the lowest MAPE in the Atlantic Ocean (14 %). The model fitting of Indian Ocean observations is overall very good (R²= 0.84). However, the MAPE is relatively high at 33 %, mostly due to a poor performance on data between 25 and 65 pmol kg⁻¹. Lastly, the two high-latitude basins have the worst R² values (Arctic Ocean: 0.49; Southern Ocean: 0.58) and high MAPE values (Arctic Ocean: 37 %; Southern Ocean: 29 %), due to the presence of outliers as well as a high density of low-concentration values (i.e. MAPE is higher for a given difference between modelled and actual values when observation values are closer to zero).

Figure 2 Basin-wise model performances assessed on all samples from each basin in the random test set. Top row: Pb concentration; middle row: ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb; bottom row: ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb.

The Pb isotope ratio models perform worst in the Arctic Ocean (${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb: R²= 0.52, RMSE = 0.010; ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb: R²= 0.55, RMSE = 0.011) and show discrepant performances in the Atlantic, Pacific, Indian, and Southern oceans. The ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb model performs best in the two basins, with the lowest number of samples in the random test, i.e. the Southern Ocean (R² = 0.93, RMSE = 0.004) and Indian Ocean (R²= 0.87, RMSE = 0.004). Samples from the Pacific Ocean are reproduced well (R²= 0.80, RMSE = 0.004), while the performance in the Atlantic Ocean is poorer (R²= 0.59, RMSE = 0.006). Lastly, for the ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb model, the Indian Ocean is the basin with the best performance (R²= 0.89, RMSE = 0.003), followed by the Atlantic Ocean (R²= 0.74, RMSE = 0.004) and the Pacific Ocean (R²= 0.57, RMSE 0.005).

3.2 Model explanation

SHAP values were calculated for all samples in the Pb concentration, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb datasets to maximise the interpretability of the models. In the next subsections, the four most important features of each model, as well as their interaction terms, are explained.

3.2.1 Pb concentrations

Temperature, dust AOD, black carbon AOD, and salinity are the most important features of the Pb concentration model (Table 4, Fig. 3). For temperature values above 5.8 °C, SHAP values are consistently positive, with warmer waters having a larger impact on predicted Pb concentrations (Fig. 3a). Black carbon AOD has the largest interaction with temperature, showing positive SHAP values in areas where both features are high. Taken together, these results suggest that surface and intermediate waters at low and middle latitudes (where black carbon AOD is highest) are associated with higher predictions of Pb concentrations. An exception to this trend is observed for samples with low temperatures and extremely low black carbon AOD, which have unexpectedly positive SHAP values for temperature (Fig. 3a). Analogously, SHAP values for black carbon AOD are positive for the lowest black carbon AOD values (< 0.001; Fig. 3c). The positive impact of extremely low black carbon concentrations in the atmosphere cannot be explained by a physicochemical process, as one would expect lower values of black carbon AOD to lead to lower Pb concentrations in the ocean. However, the positive SHAP values for temperature and black carbon AOD can be explained by the presence of samples with high Pb concentrations in the Southern Ocean (Fig. 1a), where black carbon AOD is lowest and seawater temperatures are below 5.0 °C (Figs. S3 and S4). Further evidence supporting the positive impact of the lowest black carbon AOD values on Pb concentrations in the Southern Ocean can be found in the interaction between black carbon AOD and salinity that is observed for salinity SHAP values (Fig. 3d). Indeed, this interaction indicates positive SHAP values for extremely low black carbon AOD and salinities around 34.0 ± 0.2 (Fig. S5), which are characteristics of the higher latitudes of the Southern Hemisphere.

Table 4 SHAP importance values for all features included in the Pb concentration and ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb models. Bold values indicate the four most important features for each model.

SHAP values are consistently positive for dust AOD values above 0.003 (Fig. 3b), and they show an interesting interaction with distance from bottom depth, which is especially visible for larger values of dust AOD (> 0.05). Indeed, as all grid cells at the same location were assigned the same dust AOD value throughout the water column, the interaction between dust AOD and distance from the bottom depth shows that samples more distant from the seafloor (i.e. closer to the surface) have higher SHAP values than their deeper counterparts (Fig. 3b). This finding agrees with observations of higher Pb concentrations at the surface ocean due to atmospheric deposition (Duce et al., 1991; Patterson and Settle, 1987), thereby suggesting that the model has learned from the data the key process of atmospheric Pb sourcing to the ocean.

Figure 3 SHAP values for temperature (a), dust AOD (b), black carbon AOD (c), and salinity (d) for the Pb concentration model. The colour of the dots represents the value of the feature that has the highest interaction (reported on the right-hand side of each plot). The y-axis values differ between the panels.

3.2.2 Pb isotope compositions

Seawater oxygen concentration, salinity, temperature, and dust AOD are the four most important features for the ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb model (Table 4; Fig. 4). Oxygen concentration values below 220 µmol kg⁻¹ have SHAP values between −0.004 and 0, and a sharp transition to positive SHAP values can be observed at oxygen concentrations of ∼ 240 µmol kg⁻¹ (Fig. 4a). The strong interaction between oxygen and nitrate concentrations suggests that areas of high dissolved oxygen and low nitrate concentrations, such as the Arctic and North Atlantic oceans (Figs. S6 and S7), have a large positive impact on the predicted ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb.

Figure 4 SHAP values for oxygen (a), salinity (b), temperature (c), and dust AOD (d – on a logarithmic scale) for the ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb model. The colour of the dots represents the value of the feature that has the highest interaction (reported on the right-hand-side of each plot). The y axis values differ between the panels.

SHAP values for salinity show an increasing trend as the latter varies between 32.5 and 37.0 (Fig. 4b). This is especially visible for surface and intermediate waters with densities (σ₀) between 25 and 35 kg m⁻³, indicating that the model learned well from the observations, which show higher Pb concentrations at shallower depths and a general decrease in Pb concentration values from tropical and subtropical to polar latitudes (Fig. 1a). SHAP values for temperature show a decreasing trend as waters get warmer (Fig. 4c). Interestingly, and expectedly from a theoretical perspective, this trend is reversed for density, which is the fifth-most important feature of the model and shows increasing SHAP values as waters become denser (Fig. S8). Both temperature and density also show a clear interaction with salinity, which can be explained by the fact that temperature and salinity are the two conservative parameters that contribute to seawater density. Overall, the SHAP values of and interactions between these three variables agree with the general view of Pb isotope ratios as tracers of water mass movements and ocean ventilation (Bridgestock et al., 2018; Frank, 2002; Henderson and Maier-Reimer, 2002).

Dust AOD shows a more variable distribution of SHAP values, which is better visualised on a logarithmic scale (Fig. 4d). Strongly positive SHAP values (up to 0.009) are visible for dust AODs lower than 0.004, followed by a rapid decrease to strongly negative SHAP values (up to −0.007) for dust AODs between 0.004 and 0.007. As dust AOD increases further, SHAP values first show a decreasing trend for dust AODs up to 0.12, whilst they increase afterwards (Fig. 4d).

Perhaps surprisingly, given the generally observed correlation between the ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb isotope ratios (Boyle et al., 2012; GEOTRACES Intermediate Data Product Group, 2023), the most important features for ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb (Fig. 5) do not closely match those of ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb (Fig. 4). In fact, the four most important features for the ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb model include black carbon AOD, dust AOD, seawater temperature, and surface chlorophyll-a concentration (Table 4; Fig. 5). Black carbon AOD has mostly positive SHAP values for feature values above 0.0045 and is negative otherwise (Fig. 5a). It has the strongest interaction with surface chlorophyll a and vice versa. Moreover, SHAP values for chlorophyll a are mostly negative for the lowest chlorophyll-a concentrations and range between 0.003 and −0.002 for chlorophyll-a concentrations above 0.13 mg m⁻³ (Fig. 5d). The interaction between these two features, however, does not show a clear pattern and cannot be interpreted intuitively, as both high and low chlorophyll-a concentrations have high and low SHAP values for black carbon AOD and vice versa (Fig. 5a and d). SHAP values of dust AOD (Fig. 5b) show a similar pattern to the ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb model, with lower dust AOD values having both the highest and lowest SHAP values. The interaction between dust AOD and black carbon AOD suggests that, when both features are low, dust AOD has a strong positive impact on the predicted ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb. By contrast, when dust AOD is low and black carbon AOD is higher than 0.004, dust AOD contributes to reducing predicted ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb values (Fig. 5b). Lastly, SHAP values for temperature also show a similar decreasing trend to those of the ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb model and have a strong interaction with dissolved nitrate concentrations (Fig. 5c). This interaction indicates that warmer waters with low nitrate content, such as those in tropical and subtropical areas (Figs. S4 and S7), are associated with a reduction in predicted ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb values.

Figure 5 SHAP values for black carbon AOD (a), dust AOD (b – on a logarithmic scale), temperature (c), and surface chlorophyll a (d) for the ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb model. The colour of the dots represents the value of the feature that has the highest interaction (reported on the right-hand side of each plot). The y-axis values differ between the panels.

3.2.3 SHAP model explanations and current understanding of marine Pb cycling

SHAP explanations for the Pb concentration, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb models align well with the current understanding of marine biogeochemical cycling of Pb with respect to the oceanographic predictors identified as the most important features of the three models. This indicates that the models are able to identify the key oceanographic patterns and processes that drive the distribution of Pb and its isotopes on a global scale. However, some regional explanations of dust and black carbon AODs for the Pb concentration and isotope composition models appear to contradict the theoretical intuition that dust and black carbon are sources of Pb for the ocean, with Pb isotope compositions broadly reflecting those of natural and anthropogenic land sources, respectively. In particular, the high SHAP values obtained for extremely low black carbon AODs in the Pb concentration model disagree with the hypothesis that low black carbon concentrations in the atmosphere should result in lower Pb concentrations in the ocean. Analogously, the highly positive SHAP values obtained for extremely low values of dust AOD for the ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb models and the consistently positive SHAP values obtained for relatively high black carbon AODs for the ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb model do not match the intuition that dust should contribute to shifting the Pb isotope composition of seawater toward more natural values and black carbon toward more anthropogenic ones.

Possible explanations for the partial mismatch between expected and obtained local SHAP values are provided in the following, but these will need validation through new campaigns at sea and further studies of the Pb isotope compositions of dust and black carbon sources in the atmosphere, such as the one by Nizam et al. (2020). A first explanation is that dust and black carbon are not the only sources of Pb in the global ocean, and this might have important implications at the regional scale. An example is the Southern Ocean, where ice melting and runoff could contribute significantly to the marine Pb budget of the area. Another possibility is that global climatologies of dust AOD and black carbon AOD do not capture the regional and temporal variability of the Pb isotope compositions of natural and anthropogenic Pb sources in the environment. Moreover, black carbon is emitted not only by fossil fuel consumption but also by forest fires and coal combustion. In detail, radiocarbon measurements show that the relative contributions from biomass and fossil fuel combustion to the composition of Asian black carbon can vary on a regional scale, where fossil fuel is the predominant source of black carbon on the Tibetan Plateau but biomass burning predominates for western India and the Maldives (Gustafsson et al., 2009; Li et al., 2016). As such, regional variations in the relative contributions of different atmospheric black carbon sources are expected to be associated with variations in the Pb isotope signatures associated with these black carbon emissions. Therefore, while the global climatologies of dust AOD and black carbon AOD used as model predictors are currently the best available proxies for atmospheric sources of Pb, they might not fully encompass the complexity of Pb sources in the environment, resulting in locally counter-intuitive SHAP explanations for the three models.

3.3 Global maps of Pb concentrations and isotope compositions in seawater

Global reconstructions of Pb concentrations, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb were made for each grid cell in the WOA18 using the same structure and ancillary features used for model development (Sect. 2.1).

Figure 6 shows the global distribution of Pb concentrations, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb at four depth levels (10, 1000, 2500, and 4000 m), arbitrarily chosen to represent surface, intermediate, deep, and bottom waters, respectively. These maps provide the first global-scale estimates of Pb concentrations and isotopes and can be used to further understand the dynamics and mechanisms that govern the large-scale distribution of Pb.

Figure 6 Global maps of reconstructed Pb concentrations (a), ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb (b), and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb (c). The four panels represent different depth levels (10, 1000, 2500, and 4000 m), with the white patches corresponding to the seafloor. The filled circles with white edges represent true data values for comparison with the modelled values in the background.

3.3.1 Pb concentrations

Lead concentrations vary between 2.05 and 77.65 pmol kg⁻¹. Generally, the concentrations are higher at the surface and show a decreasing trend with depth (Fig. 6a). The main exception to this trend is visible at intermediate and deep depths in the North Atlantic Ocean, where mapped Pb concentrations at 1000 and 2500 m are higher than at 10 m (Fig. 6a). However, subsurface maxima in Pb concentrations are also visible between 70 and 500 m in the North Pacific, in line with recent observations (Chan et al., 2024; Jiang et al., 2025; Wu et al., 2010; Zheng et al., 2019; Zurbrick et al., 2017). Both the general trend and the exceptions of the North Atlantic and, to a lesser extent, the North Pacific are in line with the expected distribution of Pb concentrations. Indeed, the majority of the Pb is sourced to the ocean via atmospheric wet deposition and then sinks towards the bottom, aided by particulate matter (Duce et al., 1991; Patterson and Settle, 1987). Therefore, high Pb concentrations are expected for the surface ocean, with the values decreasing with depth. However, extensive historical pollution from leaded petrol in North America and Europe (Kelly et al., 2009) led to high Pb concentrations in North Atlantic Ocean surface waters between the 1950s and 1980s. Lead concentrations eventually decreased due to environmental policies that phased out and banned leaded petrol, and Pb pollution in surface waters has decreased steadily since 1975 (Boyle et al., 2014; Bridgestock et al., 2016; Kelly et al., 2009). Therefore, the relatively high Pb concentrations in intermediate and deep waters reflect the sinking of elemental Pb pollution from earlier decades (Fig. 6a), in agreement with observations (Fig. 1a; Noble et al., 2015).

Mapped surface water Pb concentrations are highest in the Indian Ocean (average [Pb] = 32.0 ± 14.2 pmol kg⁻¹, 1 SD) and in the North Pacific (33.6 ± 11.3 pmol kg⁻¹, 1 SD), in line with recent observations from the two basins (Chen et al., 2023b; Echegoyen et al., 2014; Lanning et al., 2023). The North Pacific also shows high mapped Pb concentrations at depths of up to 1000 m (Fig. 6a), with subsurface maxima of [Pb] of up to 70 pmol kg⁻¹. These features can be reconciled with pollution from previous decades (Lanning et al., 2023; Wu et al., 2010), like in the North Atlantic, and subduction of mode and intermediate waters from their areas of formation in the north-western Pacific (Jiang et al., 2021).

Relatively high Pb concentrations are furthermore found at the surface level near South and Central America, western Africa, and Southeast Asia. Except for samples collected close to the coast in the South Atlantic Ocean, these areas have never been sampled before, and the high concentrations predicted by the model cannot be verified directly. However, the distribution of the predictor features shows that these are all areas with high seawater temperatures and relatively high black carbon and dust AODs (Figs. S4, S3, and S9). Therefore, following the model's interpretation with SHAP values (Sect. 3.2.1), these features must contribute greatly to the high predictions made by the model.

Lastly, the polar regions are characterised by the lowest Pb concentrations throughout the water column, with average Pb concentrations of 9.0 ± 3.1 pmol kg⁻¹ (1 SD) for the Arctic Ocean and 11.9 ± 3.6 pmol kg⁻¹ (1 SD) for the Southern Ocean. Additionally, both polar oceans have the least variable Pb concentrations with depth (Fig. S11), which denotes them as the areas of the global ocean that are least affected by anthropogenic pollution.

3.3.2 Pb isotope compositions

Global maps of ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb ratios show similar patterns at the four different depth levels that are analysed (Fig. 6b, c). Globally, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb ratios vary between 1.144 and 1.200, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb ratios vary between 2.420 and 2.479. Generally, lower ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb values are associated with Pb sourced from anthropogenic activities, while higher values reflect natural sources of Pb (Bollhöfer and Rosman, 2000, 2001). However, the exact isotope composition of natural and anthropogenic Pb varies between sources and basins, and there is no single reference value that is valid on a global scale due to the short residence time of Pb in seawater and the heterogeneity of global Pb sources. Recent work has, furthermore, shown that reversible scavenging can affect and modify the advected Pb isotope signature in areas of high suspended particulate matter (Lanning et al., 2023; Olivelli et al., 2024b). However, historically, Pb isotopes have been regarded as tracers of water mass movements and ventilation in the North Atlantic Ocean (e.g. Véron et al., 1998, 1999). We therefore additionally assessed the distribution of ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb across different ranges of the potential density anomaly (σ_θ; Figs. 6 and S12).

The median ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb value is highest for the Arctic Ocean (1.186) and lowest for the Indian Ocean (1.165; Fig. 7). The median values for the North and South Pacific, South Atlantic, and Southern oceans are very similar, varying between 1.168 and 1.171, while the North Atlantic Ocean has a higher median ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb of 1.180. Additionally, the Arctic Ocean shows the smallest difference between the first and third quartiles, indicating that Pb isotope compositions are rather constant throughout the geographical domain and the water column, while the Southern Ocean has the most variable distribution (Fig. 6).

Figure 7 Violin plots of the ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb distributions in the different ocean basins (top left panel) and within each basin at different potential density ranges. The white dashed line in each violin plot represents the median value, while the dotted lines represent the lower and upper quartiles (Q1 and Q3, respectively).

In the surface layer (10 m), the lowest ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb values are observed in the northern Indian Ocean, while the Arctic Ocean has the highest ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb values and coastal areas around North America and the Arctic Ocean and eastern Pacific Ocean have the highest ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb values (Fig. 6b, c). Additionally, the distributions of both isotope ratios show the presence of low values (${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb: 1.148–1.160; ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb: 2.420–2.437) at a depth of 1000 m with a core at the location of the Subantarctic Front (∼ 48° S), expanding horizontally to all basins in the Southern Hemisphere and vertically to the deep layer (2500 m). Low ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb values at intermediate depths have been observed in the Indian (Lee et al., 2015), South Atlantic (Olivelli et al., 2024b), and Southern oceans off the coast of Tasmania (GEOTRACES cruise GS01; unpublished data from the MAGIC group at Imperial College London). Historically, anthropogenic emissions from Australia, New Zealand, Chile, and South Africa have recorded the lowest Pb isotope ratios and are mostly associated with the use of broken-hill-type leaded petrol (with ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb as low as 1.060 and 2.328, respectively; Bollhöfer and Rosman, 2000). Our maps, therefore, suggest that Antarctic Intermediate Water (AAIW), which sinks to depths at latitudes between 50 and 60° S across the Southern Ocean, in the eastern South Pacific Ocean, and near the Drake Passage, has the greatest anthropogenic signature of all of the non-surface water masses considered (Fig. 6b, c). Further evidence for the highly anthropogenic signature of AAIW can be found by analysing the distribution of Pb isotope compositions across density layers. Indeed, AAIW, which has a potential density anomaly between 27.0 and 27.4 kg m⁻³, corresponds to the layer with the lowest median ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb in the South Pacific, South Atlantic, and Indian Ocean (Fig. 7). As leaded petrol was phased out in the abovementioned countries and their neighbours, before the time of sampling on GEOTRACES cruises, the widespread anthropogenic signature observed at 1000 m in the Southern Hemisphere either reflects persisting pollution from previous decades that was deposited in the Southern Ocean and is transported northward or the remobilisation of previously deposited Pb on land and subsequent atmospheric deposition in the formation areas of AAIW.

In the North Pacific, North Atlantic, and Southern oceans, the median ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb values increase with depth and potential density ranges (Figs. 6b, c and 7). By contrast, in the South Atlantic, Indian, and, to a lesser extent, South Pacific oceans, a decrease in ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb values can be observed between surface and intermediate water masses, generally followed by an increase as water masses become denser (Fig. 7). The mapped distributions of Pb isotope compositions, taken together with those of Pb concentrations, therefore indicate that the deepest and densest water masses are least affected by anthropogenic pollution on a global scale, as they have the lowest Pb concentrations and highest ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb ratios. On the other hand, surface and intermediate waters are most affected by pollution, which is discernible (i) in the Northern Hemisphere from Pb concentration data and to some extent from Pb isotope compositions, although natural and anthropogenic sources affecting those basins have a relatively small difference in Pb isotope fingerprints, and (ii) in the Southern Hemisphere and northern Indian Ocean from the low ratios of ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb.

3.4 Uncertainty quantification

Three different sources of uncertainty were identified for our model outputs and are described and, where possible, quantified in this section. The first source of uncertainty relates to the observed values of Pb concentrations and isotope compositions, the second to the ancillary variables used as features in the models, and the third to the random split of the data into a training set and a test set.

All observations of Pb concentrations, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb ratios, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb ratios are characterised by a fundamental level of uncertainty associated with measuring Pb and its isotopes in seawater. While the GEOTRACES programme ensures data quality of the highest standard through a thorough process of intercalibration, it cannot eliminate uncertainty that is intrinsic to sample processing and analysis. In the case of Pb concentrations, individual sample uncertainties arising from replicate measurements generally range between 0.1 % and 80 %, and on average they are around 10 % to 25 % (relative 1 SD; e.g. Boyle et al., 2020; Echegoyen et al., 2014). Inter-laboratory uncertainties arising from the use of different sample processing methodologies and instrumentations and assessed on intercalibration samples, however, typically show better reproducibility. In fact, the relative standard deviation of the North Atlantic GEOTRACES and SAFe reference samples is between 3.5 % and 9.4 % (1 SD; GEOTRACES, 2009, 2013). For Pb isotope compositions, individual sample uncertainties, assessed either as replicate measurements of the sample itself or extrapolated from replicate measurements of the NIST 981 standard reference material, are generally around 0.2–0.5 ‰ for both ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb (e.g. Boyle et al., 2020; Bridgestock et al., 2016). Inter-laboratory calibration exercises of Pb isotope composition measurements of seawater have returned uncertainties ranging between 0.1 ‰ and 3.5 ‰ for both isotope ratios considered (Boyle et al., 2012). So, if we consider the MAPE to be a reflection of the uncertainty associated with the models' predictions, predicted values of Pb concentration have an uncertainty (20.2 % random test set and 18.3 % geographic test set; Table 3) that falls within the average values for replicate measurements (10 %–25 %). Similarly, the uncertainty associated with predicted ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb values (MAPE: 1.0 ‰–3.0 ‰; Table 3) falls within the range observed for inter-laboratory calibration.

Pooling together dissolved and total dissolvable Pb concentration and isotope composition data creates an additional layer of uncertainty. Our choice is backed by the evidence that, in the open ocean, dissolved and total dissolvable Pb data are comparable (Bridgestock et al., 2018; Olivelli et al., 2023, 2024b; Schlosser et al., 2019). However, this might not be the case in estuaries and coastal areas, where boundary exchange and isotopic re-equilibration are believed to play an important role in determining the observed Pb concentrations and isotope compositions (Chen et al., 2016, 2023a). Therefore, new observations, particularly from coastal areas and other sparsely sampled locations, would provide valuable data for the validation of our modelling results. The uncertainty associated with the models' features cannot be quantified directly and contributes to the models' MAPE values, but it is nevertheless expected to be small and not create systematic bias. Indeed, the data products used in this study, i.e. WOA18 for oceanographic variables, CMEMS's Global Ocean Biogeochemistry Hindcast product for chlorophyll a, and CAMS's ECMWF Atmospheric Composition Reanalysis 4 for black carbon AOD and dust AOD, are all subjected to stringent quality control.

Figure 8 Coefficients of variation for reconstructed Pb concentration (a), ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb (b), and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb (c) values. As in Fig. 4, the different panels represent different depth levels (10, 1000, 2500, and 4000 m) and the white patches correspond to seafloor locations. Please note that the colour scales differ for the three panels.

Lastly, to quantify the uncertainty arising from the random splitting of the data into a training set and a test set, the best Pb concentration, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb model architectures were initiated with 100 different random train–test splits while keeping the hyperparameters constant. The means and standard deviations of the model ensemble predictions for the global maps of Pb concentrations and isotope compositions were used to calculate the coefficients of variation for the mapped Pb distributions (Fig. 8). The coefficient of variation was calculated as the ratio between the average value of each cell across the 100 different ensemble members and the associated standard deviation. For both Pb concentrations and isotope compositions, the coefficients of variation are highest in the Arctic and Southern oceans, especially in the surface and intermediate layers (Fig. 8). Both the Arctic and Southern oceans are characterised by the lowest Pb concentrations observed, and therefore any variability in the outputs of the ensemble would cause a larger coefficient of variation as the mean concentrations are low. However, the spatial variability in the observed Pb concentrations in the Southern Ocean is much larger than in the Arctic Ocean. Indeed, high concentrations are found at the surface and in the water column near the Antarctic Peninsula, likely caused by glacial meltwater inputs, while much lower values are visible in the Indian and Australian sectors. Therefore, the different splits between the training and test sets lead to much more variable model outputs for the Southern Ocean than for its Arctic counterpart. The high latitudes are also the areas with the largest coefficients of variation for ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, arguably because of the significant presence of both anthropogenic and natural Pb. Based on the variability of the outputs of the three model ensembles, it can be argued that at-sea sampling efforts should focus on the Southern Ocean to obtain a clearer understanding of the distribution and cycling of Pb and its isotopes in that area and to reduce the uncertainty associated with our modelled outputs.

4 Model and data applications

More than 20 years ago, Henderson and Maier-Reimer (2002) used a general circulation model to investigate the natural marine cycle of Pb and its isotopes. As these authors did not consider inputs from anthropogenic activities, here we present the first study that provides global maps of the distribution of Pb concentrations and ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb ratios.

In addition to providing a detailed demonstration of the applicability of machine learning to marine (isotope) geochemistry, the new models and their outputs can be used to (i) benchmark future levels of pollution, (ii) inform sampling strategies and campaigns, and (iii) compare outputs from different models, whether process-based or data-driven. Lead concentration and isotope composition data from our 12-year climatologies provide useful information on the levels of marine pollution in terms of both sources and spatial extent. This insight can be used to assess the impact of human activities and any potential environmental policies implemented in the future. Both the mapped distributions and their associated uncertainties can be used by the community to support any process-oriented or large-scale sampling campaigns. In particular, the great uncertainty associated with Pb concentrations in the Southern Ocean and the widespread distribution of low ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb values at intermediate depths throughout the Southern Hemisphere suggest that these areas should be prioritised in the future. Lastly, the 1° × 1° data products developed will provide a useful and robust comparison for any global and regional modelling studies that attempt to parameterise and reproduce the processes driving the biogeochemical cycle of Pb, which still requires future work.

5 Code and data availability

All scripts generated to build the Pb concentration, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb datasets and the corresponding models are available in Olivelli (2025; https://doi.org/10.5281/zenodo.15355008). The gridded climatologies of Pb concentrations, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb, as well as the results of the model ensemble for uncertainty calculations, are available in Olivelli et al. (2024a, https://doi.org/10.5281/zenodo.14261154).

6 Conclusion

We successfully generated global climatologies of Pb concentrations, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb ratios, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb ratios using the non-linear regression algorithm XGBoost, trained and tested on high-quality Pb data collected as part of the GEOTRACES programme. We found that Pb concentrations are highest in the surface layer of the Indian and North Pacific oceans and at intermediate depths in the North Atlantic Ocean, and they are lowest in the Arctic Ocean. The ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb ratios are lowest in the surface Indian Ocean and at intermediate depths in the Southern Hemisphere, and they are highest throughout the water column in the Arctic Ocean. Taken together, the distributions of Pb concentrations, ${}^{\mathrm{206}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb ratios, and ${}^{\mathrm{208}}\mathrm{Pb}{/}^{\mathrm{207}}$Pb ratios indicate that (i) the surface Indian Ocean has the highest levels of pollution, (ii) pollution from previous decades is sinking in the North Atlantic and Pacific oceans, and (iii) a highly anthropogenic fingerprint originating in the Southern Ocean is spreading at intermediate depths throughout the Southern Hemisphere with AAIW. The lower model performance and higher uncertainty associated with reconstructions of the Southern Ocean suggest that this is a key area when prioritising future sampling campaigns to expand our understanding of biogeochemical processes that drive the distribution and cycling of Pb and its isotopes in the ocean. More broadly, our approach demonstrates the applicability of machine learning to marine (isotope) geochemistry, even when data are scarce and very sparse, as with Pb isotope compositions. Lastly, our model output will provide a useful benchmark for future levels of pollution and a valuable comparison for process-based models that aim to improve understanding of the biogeochemical processes governing the marine cycling of Pb.