An Artificial Neural Network to Infer the Mediterranean 3D Chlorophyll-a and Temperature Fields from Remote Sensing Observations

Abstract

Remote sensing data provide a huge number of sea surface observations, but cannot give direct information on deeper ocean layers, which can only be provided by sparse in situ data. The combination of measurements collected by satellite and in situ sensors represents one of the most effective strategies to improve our knowledge of the interior structure of the ocean ecosystems. In this work, we describe a Multi-Layer-Perceptron (MLP) network designed to reconstruct the 3D fields of ocean temperature and chlorophyll-a concentration, two variables of primary importance for many upper-ocean bio-physical processes. Artificial neural networks can efficiently model eventual non-linear relationships among input variables, and the choice of the predictors is thus crucial to build an accurate model. Here, concurrent temperature and chlorophyll-a in situ profiles and several different combinations of satellite-derived surface predictors are used to identify the optimal model configuration, focusing on the Mediterranean Sea. The lowest errors are obtained when taking in input surface chlorophyll-a, temperature, and altimeter-derived absolute dynamic topography and surface geostrophic velocity components. Network training and test validations give comparable results, significantly improving with respect to Mediterranean climatological data (MEDATLAS). 3D fields are then also reconstructed from full basin 2D satellite monthly climatologies (1998–2015) and resulting 3D seasonal patterns are analyzed. The method accurately infers the vertical shape of temperature and chlorophyll-a profiles and their spatial and temporal variability. It thus represents an effective tool to overcome the in-situ data sparseness and the limits of satellite observations, also potentially suitable for the initialization and validation of bio-geophysical models.

1. Introduction

In recent years, the assessment and monitoring of the marine environmental status has received ever-growing attention, due to the potentially critical impact of ongoing natural and human-induced changes on related ecosystem functioning and services. A deeper understanding of marine ecosystem evolution and dynamics, however, would require 3D observations of key essential ocean variables at different temporal and spatial scales and much wider and regular coverage than presently achievable.

Indeed, most of the existing data, especially those providing biological/bio-optical/bio-chemical parameters, are collected by in situ sampling either through coastal monitoring programs and time-limited oceanographic cruises, or by fixed platform, as moored buoys, and by autonomous instruments, such as Biogeochemical Argo (BGC-Argo) [1,2,3,4,5,6]. These data are eventually able to provide accurate descriptions of local conditions along the water column, but they are clearly not sufficient to describe the processes occurring over the wide range of temporal and spatial scales involved/affected by undergoing changes [7].

On the other hand, the space-time coverage and accuracy of surface observations from satellite platforms has been continuously improving [8,9,10]. Remotely-sensed measurements have significantly changed our way to investigate marine ecosystem dynamics: the retrieval of chlorophyll-a concentration (hereafter chlorophyll) from ocean color measurements, used as a proxy of biomass, revolutionized our understanding of the biological-physical interactions and opened new possibilities for quantitative assessments of basin scale variations [11,12,13]. Nevertheless, one important limit of the remote sensing data is that they sense only the surface layer and are not able to observe properties at depth. In fact, the detection capability of satellite observations depends on the band and measurement principle considered. For ocean color data, the measured signals result from the interaction of light with the water constituents, decreasing exponentially with depth. The attenuation coefficient, modulated by the presence of phytoplankton pigments, non-algal particles (NAP), and colored dissolved organic matter (CDOM) [14], determines the depth of the euphotic zone (Z_eu) (defined as the depth where 99% of the surface light is attenuated), limiting the amount of energy available for the photosynthesis. Satellite observations can only reveal the signals integrated between the surface and the first optical depth, thus sensing even shallower layers [15].

As such, most of the variability of marine bio-physical variables remains unobserved. Noteworthy, the synergic use of satellite and in situ data already proved as an effective way to exploit the advantages of both, using satellite data to complement the accuracy of punctual in situ observation for a description of ocean’s vertical structure at higher spatial and frequency resolution, either at regional and global scales. The quantitative extrapolation of essential ocean variables from surface to deeper layers can be pursued with several techniques, which span from purely empirical models to more recent machine learning approaches. For physical variables, as temperature and salinity, different empirical methodologies have been proposed. Some of them include (multi)linear regression models [16,17,18,19,20,21,22,23], as, e.g., in Jeong et al. [23], where the subsurface temperature is inferred by exploiting its high correlation with surface fields, Coupled Pattern Reconstruction [20], Empirical Orthogonal Function reconstruction, [21,24,25] and Gravest Empirical Modes [26], where empirical mode decompositions are used. As a counterpart, for biological variables (e.g., chlorophyll concentration), Morel and Berthon [27] and Uitz et al. [28] proposed empirical parametrizations for the vertical biomass retrieval in different surface trophic regimes.

Nevertheless, these types of inversion models can often come across erroneous estimates, due to the complexity and non-linearity of the system, especially in those areas characterized by non-trivial vertical stratification and small-scale phenomena, that can easily determine different vertical profile shapes corresponding to the same surface observation of a single variable.

More recently, investigations on the potential use of artificial intelligence (AI) tools to address this challenge have shown promising results. For ocean physical applications (at local and global scales), they include: the combination of a pre-clustering process and the use of a neural network to determine e.g., the Subsurface Temperature Anomaly (STA) using Sea Surface Temperature (SST), Sea Surface Height (SSH) and surface wind observation [29]; the use of self-organizing map (SOM) applied on SST, SSH and Sea Surface Salinity (SSS) data from Argo gridded monthly anomaly datasets in order to infer the STA [30]; the application of support vector machine on satellite SST, SSH and SSS anomalies for the estimation of the STA in the Indian Ocean [31], at global scale, the use of an extreme gradient boosting (XGBoost) for the STA retrieval from surface observations [32]; the use of a multi-layer network or multilayer convolutional long- and short-term memory (M-convLSTM) to predict directly the subsurface temperature field from surface data [33,34], the use of SST, (SSS) and Absolute Dynamic Topography (ADT) to predict hydrographic profiles and related uncertainties through a deep LSTM combined with a dropout strategy [35].

Chlorophyll vertical reconstruction with machine learning was first carried out by Richardson [36], who tested an unsupervised neural network to characterize the vertical shape of chlorophyll from satellite data and the estimation of primary production in the coastal Benguela upwelling system. A similar unsupervised network, combined with Hidden Markov Models, is proposed in [37], while a cascade neural network, trained with a quasi-Newton method, is presented in [38]. Conversely, supervised neural networks were applied in [39,40,41,42] to provide primary production estimates. Similar supervised approaches, based on Multi-Layer-Perceptron (MLP), were implemented in two works [43,44] aimed to reconstruct the vertical bio-optical and bio-geochemical variables from Bio-Argo data combined with satellite measurements, still requiring concurrent vertical observations of physical variables. Sauzède et al. [45] proposed a network for the retrieval of water-column nutrients from satellite and BCG-Argo measurements. Similarly, a Bayesian approach has been proposed to augment in situ observations by inferring vertical profiles of variables that are not directly measured [46].

All these works exploit two key points: the high frequency and coverage of satellite observations and the power of AI to reveal non-linear relations [47]. In fact, considering the complexity of marine ecosystems and the reduced number of observed environmental and biological variables, the capability of machine learning techniques to recognize the implicit relationship between specific variables, without imposing a priori models, allows to overcome some of these limitations, which would result impossible to solve through poorly constrained deterministic models [47].

In this work, we develop and test a new model to directly project chlorophyll and temperature surface data at depth, based on a deep feedforward network. The proposed model, focused on the Mediterranean Sea, extends previous developments described in [48]. As detailed in the following sections, it represents a major evolution with respect to our previous work, because:

(1)

it introduces a different network architecture (different input vector definition, different network depth, different optimization algorithm);

(2)

it introduces a new set of input variables (with dedicated sensitivity tests): Absolute Dynamic Topography and U & V geostrophic velocities;

(3)

both training and prediction are effectively based on satellite data and not on in situ surface values;

(4)

for the first time, the technique is applied to retrieve and analyze full 3D monthly climatological fields in the Mediterranean Sea.

Specifically, we propose a regional model to jointly infer the vertical chlorophyll and temperature fields at once. A two-hidden-layer feedforward neural network is employed and trained directly on satellite data and not anymore on in situ surface values. Here, we focus on the Mediterranean Sea, which is characterized by very complex bio-geochemical processes and physical dynamics (and interactions among them), that make the prediction of the ocean interior structure from surface observations a real challenge. Results of the reconstruction are validated against in situ data over fully independent matchup points, significantly improving with respect to climatologies. The network is successively applied to reconstruct and analyze the Mediterranean seasonal climatologies of 3D chlorophyll fields.

2. Materials and Methods

Artificial neural networks require a training and test dataset to verify not only model performances but also to assess their ability to generalize the unseen data. Here, a Multi-Layer-Perceptron (MLP) is used for the concurrent retrieval of the subsurface chlorophyll and temperature profiles from satellite observations plus some ancillary information. Specifically, the network training and test were based on a matchup database where each record includes surface (satellite plus latitude, longitude and date) variables and in situ vertical profiles when simultaneously available.

In the following sections, we provide: a description of the in situ and satellite datasets used; details on MLP, the pre- and post-processing of the input/output data; a description of the network implementation and its training.

2.1. In Situ Database and Quality Control

The in situ database includes all concurrent profiles of temperature and chlorophyll collected during 20 oceanographic cruises carried out in the Mediterranean Sea from 1998 to 2015 (same as in [48]). Most of these data were acquired and processed by the GOS group (Global Ocean Satellite monitoring and marine ecosystem studies) of the Institute of Marine Sciences (ISMAR) of Rome of the Italian National Research Council (CNR).

More in detail, the chlorophyll profiles were acquired with a fluorometer and then calibrated with bottle samples that were concurrently sampled in each cruise.

All data passed the quality controls based on SeaDataNet standards (see quality control manual at https://www.seadatanet.org/Standards/Data-Quality-Control) and further requirements stated for the specific goal of this work. The same procedure adopted in [48] was applied on either in situ chlorophyll and temperature profiles and included the following points:

• Visual check (to guarantee the consistency of the database);
• Check of the first acquisition depth (accepted if ranging between 3 and 4 m to avoid any noise on the first acquired measure);
• Check of missing points along the profiles: a specific-case evaluation was done observing each sample. The profile showing too many missing points were discarded, otherwise they were linearly interpolated.
• Check maximum acquisition depth: only those profiles with a depth greater or equal to 150 m were included.

The quality control procedure selected 1213 over 2710 profiles suitable for the subsequent MLP development. From this data-pool, about 30% was randomly extracted and used for testing the network performance, as an independent dataset, once implemented the network (Figure 1).

2.2. Matchup Satellite Database

The network, here proposed, takes as inputs different surface remote sensing variables.

All Mediterranean satellite observations used in this work were downloaded from Copernicus Marine Environment Monitoring Service (CMEMS) (https://marine.copernicus.eu/) and cover the time window from 1998 to 2015.

Satellite daily chlorophyll (Chl_SAT) maps are regional products interpolated (L4) at 1 km of resolution from the merging of multi sensor data (SeaWiFS, MODIS-Aqua, NPP-VIIRS). The chlorophyll field is computed by applying two algorithms regionalized for the Mediterranean Sea and diversified for water types: the MedOC4 [49] applied for open waters (Case1) and AD4 [50] for the optically complex waters (Case2). The membership of each pixel to the different water types is determined by comparing its satellite spectrum to the average water type spectral signature obtained from the MedOC4 [49] and the CoASTS [51] in situ databases, respectively. The processing of intermediate pixels between Case1 and Case2 follows the [52] methods. Finally, the interpolation procedure is based on the DINEOF algorithm described in [53].

The Mediterranean Sea Surface Temperature (SST) products are produced by the CNR starting from the reprocessed Pathfinder V5.3 (PFV53) AVHRR data covering the 1982–2018 period and combined them with a bias-corrected version of the CMEMS NRT L4 data up to 2017 to provide a consistent time series of interpolated (L4) field at a resolution of 0.0417° × 0.0417°. The interpolation is pursued through an Optimal Interpolation algorithm (for more details about this product see the Quality Information Document; http://marine.copernicus.eu/documents/QUID/CMEMS-SST-QUID-010-021-022.pdf and [54,55,56]).

The Absolute Dynamic Topography (ADT) and absolute geostrophic currents (zonal and meridian components, U & V respectively) are additional variables included in the altimeter satellite gridded data of Sea Level Anomalies (SLA) available in the CMEMS web-portal, produced by the production unit of SL-CLS of Toulouse in France. The ADT is the instantaneous height above the Geoid and it is obtained from the sum of the SLA_N and MDT_N (Mean Dynamic Topography, [57], where N, here, is the mean reference period (1993–2012). The ADT is a delayed time optimally interpolated product at a resolution of 0.125° × 0.125° processed by the DUACS multimission altimeter data processing system and coming from the merging of different altimeter missions (for more details about the processing see the Quality Information Document at https://resources.marine.copernicus.eu/documents/QUID/CMEMS-SL-QUID-008-032-062.pdf or http://duacs.cls.fr). The absolute geostrophic currents are obtained by the SLA and ADT product, above mentioned.

The value for each satellite matchups was obtained by computing the median value of the cloud free pixel within a 3 × 3 box, for a total number of 1213 matchups.

2.3. Multi-Layer-Perceptron (MLP)

In this work, a Multi-Layer-Perceptron (MLP), [58,59] is developed starting from the nntool of MATLAB^® version 9.7 (R2019b) [60]. The MLP consists of a series of layers connected to each other from the input to the final output, through one or more hidden layers. Each layer contains the computing elements (neurons or nodes) characterized by a non-linear activation function and connected to the subsequent layer with weighted edges.

More specifically, the network used in this work is a Feedforward Neural Networks which can be employed for any kind of input/output mapping [61]. During the training, the weights are iteratively adjusted following the minimization of a cost function as the quadratic error between the observed and predicted outputs [44]. This iterative process continues until a minimum is reached, through the technique of error backpropagation [59,62]. In our network, the training is performed by applying the Scaled Conjugate Gradient backpropagation (SCG). This is an optimization method based on conjugate directions but with no line search at each iteration. By using a step size scaling mechanism, this method resulted to be faster than other second order algorithms [63].

2.4. Data Pre-and Post-Processing

The candidate input variables selected for the concurrent prediction of chlorophyll and temperature at different depths (here represented by the pressure levels) are the following: the depth, latitude, longitude, day of the year, the satellite SST, Chl_SAT, ADT and U & V components of geostrophic currents retrieved from satellite topography. Differing from [48], here, three new satellite-derived variables have been included and different configurations have been tested.

The choice of the co-predictors is related to the nature of the examined problems and to the limited number of variables available to characterize the ocean state. In our model, the time variable (day of the year) and geographical coordinates (latitude and longitude) are expected to account for the seasonal evolution of the examined parameters and for the spatial variability represented by the different bio-geographies of the Mediterranean Sea [64]. Moreover, ADT, U & V components provide information of the surface ocean dynamics, and might thus help to better retrieve vertical stratification and related environmental conditions. Zonal and meridional surface currents can have an impact on surface chlorophyll distribution in terms of horizontal advective processes and correlated nutrient transport, while the ADT, coupled to the SST, provides information on the steric variations that can be related to modifications of the internal density structure (e.g., [65,66,67,68]). All these predictors were tested together with Chl_SAT to optimize the MLP performance in the Mediterranean Sea.

The entire dataset was randomly divided in two subsets: one for the training and internal validation (70% of total) and the remaining (about 30% of the total) for the independent test of the MLP performance.

Before the training, the Chl_SAT and in situ vertical chlorophyll were log-10 transformed. Additionally, to take advantage of the annual cycle, the sampling date was projected on circular coordinates as follows:

$$day1=cos\left(2\pi \ast \left(day of the year/365\right)\right)$$

(1)

$$day2=sin\left(2\pi \ast \left(day of the year/365\right)\right);$$

(2)

Then, all inputs (except day1 and day2) were standardized by applying a z-score transformation as follow:

$$z\left(x_i\right)=\left(x_i-\mu \right)∕\sigma$$

(3)

where the $\mu$ and $\sigma$ are respectively the mean and standard deviation for each co-predictor of the training dataset, while $x_i$ is a generic i- input value.

The z-score transformation is a linear normalization of the inputs, widely employed in the MLP development, arranging the inputs and the target outputs within similar ranges of values [69,70]. At the same time, the log-10 transformation adopted for the chlorophyll concentrations is strictly related to the natural distribution of this variable in the Mediterranean Sea, which can vary from very small concentration as e.g., 0.01 mgm⁻³ up to values even considerably higher than 10 mgm⁻³.

Obviously, once completed the MLP training, the outputs were converted into original values applying the z-score inverse formula based on reference $\mu$ and $\sigma$ of training targets.

2.5. MLP Training

The structure of the network was defined on the base of different tests using one or two hidden layers, with a number of neurons varying between 5 and 36. To choose the final configuration, other hyper-parameters were varied: the number of valid fails before stopping the training (early stopping in the following), the application or not of regularization strategies, the use of a linear or hyperbolic tangent sigmoid as transfer function of the output layer and different number of epochs (e.g., varying from 1000 to 8000). After several trials, the architecture with minimum error on the test set and minimum number of neurons was selected as optimal. The structure of the network was thus set to (10-12-8-2) (Figure 2), which means one input layer with 10 inputs, two hidden layers with 12 and 8 nodes, respectively, and one output layer including the two output variables (chlorophyll and temperature). The transfer function chosen for both hidden layers was a hyperbolic tangent sigmoid, while a linear function was applied for the output layer. During the training phase the targets were scaled between a [−1, 1] interval. This choice ensures an equal treatment of the predictions and training targets in the cost function, without favoring the targets with the largest value ranges. Finally, to improve the model generalization, the entire training dataset was replicated to ensure a more balanced selection of those profiles less represented in the dataset [71].

In order to improve the network generalization, the training was carried out following the so called early stopping criterion: during the training phase, the dataset is randomly divided in two additional subsets: the first is properly the training set (80%) which is used for computing the gradient and updating the weights and biases, while the second is called validation set (20%) and it is used for an internal validation of the network and leading its learning. In fact, when the error computed on the entire validation set starts to increase and deviates from the mean-square error (MSE) computed on the whole training set, that’s still decreasing, the learning phase is stopped and the weights at the minimum of the validation error are saved; in our training the number of validation fails was set to 1. This means that the training was stopped at the first divergence between the validation and training MSE, respectively. This technique is widely used to avoid as much as possible the overfitting [47,58].

To further improve the generalization of the network, an additional method was adopted, the so-called regularization [70]. This is based on the modification of the performance function (usually represented by the MSE computed on the training set), by adding a term that consists of the mean of the sum of squares of the network weights and biases:

$$msereg=\gamma \ast msw+\left(1-\gamma \right)\ast mse$$

(4)

with γ representing the performance ratio, a regularization parameter that defines the penalty imposed to the weights and biases, and$msw=\frac{1}{n}{\displaystyle \sum }_{\dot{J}=1}^n{w}_j^2$defining the mean of the sum of the squared weights and biases $\left(w_j\right)$. The regularization effect forces the network to adopt smoother weights. This modified MSE, in addition to the normalization of the outputs within the interval [−1, 1], enables to have smaller weights and biases and thus prevent the overfit [60]. In our case, the performance ratio was set to 0.01, after several tests with values in the range between 0.001 and 2.

3. Results

3.1. MLP Performance Evaluation and Input Sensitivity Analysis

The MLP performance was assessed on the test set by comparing the predicted chlorophyll and temperature deep values with the in situ ones. As described before, the test set represents 30% of the total database (364 stations accounting for 53,872 single depth individual data values). The inputs were standardized and processed as done during the training (according to the reference values of $\mu$ and $\sigma$), and then passed to the network. The outputs were then transformed by applying the inverse formula of the z-score scaling, according to the reference $\mu$ and $\sigma$ of targets of training dataset. The metrics used for this analysis are reported in the following

Table 1. Reference statistical parameters employed for the assessment of the Multi-Layer-Perceptron (MLP) performance with respect to the observed data. N is the number of observations and x is the observed value.

Determination Coefficient	$r^2={\left(\frac{{\sum }_i\left(x_i-\overline{x}\right)\left(Mode{l}_i-\overline{Model}\right)}{\sqrt{{\sum }_i{\left(x_i-\overline{x}\right)}^2{\sum }_i{\left(Mode{l}_i-\overline{Model}\right)}^2}}\right)}^2$
Root Mean Squared Error	$RMSE=\sqrt{\frac{1}{N}\sum _{i=1}^N{\left(x_i-Mode{l}_i\right)}^2}$
Mean Squared Error	$MSE=\frac{1}{N}\sum _{i=1}^N{\left(x_i-Mode{l}_i\right)}^2$
Mean Bias Error	$MBE=\frac{1}{N}\sum _{i=1}^N\left(x_i-Mode{l}_i\right)$
Mean Absolute Percentage Error	$MeanAPE=\frac{1}{N}\sum _{i=1}^N\left|\left(\frac{x_i-Mode{l}_i}{x_i}\right)\right|\times 100$
Median Absolute Percentage Error	$MedianAPE=median\left|\left(\frac{x_i-Mode{l}_i}{x_i}\right)\right|\times 100$
Percent Mean Bias Error	$MBE\%=\frac{1}{N}\sum _{i=1}^N\left(\frac{x_i-Mode{l}_i}{x_i}\right)\times 100$
Mean Absolute Error	$MAE=\frac{1}{N}\sum _{i=1}^N\left|x_i-Mode{l}_i\right|$

:

3.1.1. Chlorophyll Assessment

Figure 3 shows the density scatterplots and statistics of the comparison between the observed and predicted chlorophyll values on training (Figure 3a) and test (Figure 3b) dataset, respectively.

The MLP predicts the target chlorophyll values with a good accuracy. The statistic on each panel is comparable for both datasets, suggesting the network is no-overfitting during the learning phase. The determination coefficient is high in both cases (training, r² = 0.73; test, r² = 0.71), while the RMSE is relatively low with values of 0.20 mgm⁻³ for training and 0.21 mgm⁻³ for the test, indicating the good fit of the model on the observed data and the closeness of the predictions to the real data. The biases (training, MBE = 0.03 mgm⁻³ and test, MBE = 0.04 mgm⁻³) are also low, while the negative values of the MBE% (training, MBE% = −20% and test, MBE% = −17%) suggests an over-all overestimation of the MLP with respect to the targets. This is also observable in the scatterplots of Figure 3 where, for lowest chlorophyll concentrations, the MLP tends to overestimate the observed data, showing a reduced prediction variability with respect to the targets.

Figure 3 shows that most of the data fall along the 1:1 line, with a reduced scatter. The performance of this network improves with respect to the shallower network set-up in [48], that was trained on in situ data alone. In fact, the test validation statistics show a clear advantage of the new network architecture: the determination coefficient (r² = 0.71) is higher than r² = 0.69 of [48], and the relative percentage errors (MBE% = 17%; MeanAPE% = 48%; MedianAPE% = 30%) and RMSE = 0.21 mgm⁻³ are also lower (MBE% = 28%; MeanAPE%= 57%; MedianAPE% = 34%; RMSE = 0.26 mgm⁻³). Furthermore, even if the test set used is not exactly the same in the two models, the comparison between both reconstructions based on satellite surface data and the in situ profiles also highlights an improvement in terms of r² (0.71 against 0.63 of [48]) and RMSE (0.21 mgm⁻³ against 0.23 mgm⁻³ of [48]).

Figure 4 shows RMSE profile computed on the test set as a function of depth. The higher errors (RMSE = 0.35 mgm⁻³) are observed between 20 m and 60 m, which corresponds to the portion of the water-column where major chlorophyll variations occur, mostly driven by the evolution of the Deep Chlorophyll Maximum (DCM). DCM occurrence is determined by the complex dynamics of light, temperature and nutrients availability within the euphotic layer [72]. This increment of the errors at the DCM depths is reduced when observing the MBE% profile (not shown), where, being relative, the errors are impacted not only by the DCM variability but also by higher chlorophyll concentrations. On the contrary, the errors noticeably decrease close to zero (RMSE = 0.1 mgm⁻³) towards deeper layers, reflecting the lower chlorophyll variability and concentrations far from the surface [3].

In order to summarize the capability to predict subsurface chlorophyll shape from surface data, all in situ and their corresponding MLP profiles were grouped according to different satellite surface biomass ranges (following [28,43]) and then compared. As done in [48], all the profiles falling within a specific surface biomass class were averaged, and the mean profile and standard deviation of observations (black line and shaded grey area in Figure 5) were compared with those of simulated data (blue line and blue shaded area in Figure 5).

There is a notable agreement between modelled profiles and observations, remarkably evident for some specific surface categories, as e.g., in low ranges of 0.001–0.04 mgm⁻³ and 0.04–0.08 mgm⁻³ (Figure 5a,b) and other classes of higher surface concentrations as in 0.12–0.2 mgm⁻³ and 0.2–0.3 mgm⁻³ (Figure 5d,e). The predicted profiles with their standard deviations are overlapped on real data in most of the cases, especially at shallower and deeper depths than those of the DCM.

The MLP performance is reduced for highest surface concentrations where the predicted profile appears smoother than the observed one (Figure 5h). However, the MLP forecast correctly recognizes the biomass peak position along the water-column for all classes, though its magnitude can, in some cases, be underestimated.

3.1.2. Temperature Assessment

Figure 6 shows the density scatter plots and statistics of the comparison between the observed and predicted temperature values on training and test dataset.

As for chlorophyll, the MLP performances on temperature are promising. The MLP displays an MBE of 0.01 °C for training and −0.01 °C for the test and the MedianAPE = 1.9%, which, being less sensitive to the outliers, results lower than the MeanAPE = 2.7%.

Both determination coefficients are high (training, r² = 0.94; test, r² = 0.93), while the RMSE of 0.67 °C for training and of 0.71 °C for the test, partly suffers from the scatter observable in Figure 6. The data cloud fit closely the 1:1 line for the entire temperature ranges, with an increase of the scatter for intermediate temperature intervals (Figure 6, 20–26 °C). Also, for temperature, the training and test errors are comparable, once more suggesting that the network is not overfitting.

Figure 7 provides the RMSE profile computed on the test set at different depths. At surface, the RMSE is 0.7 °C. The error peaks are observed at $~$19–20 m attaining values of RMSE = 1.13 °C. Then, the error decreases, quite constantly, down to 80 m, reaching very low values (RMSE = 0.3 °C) at bottom layers (about 150 m). The curve shows the maximum variability between 20 m and 80 m, where usually the mixed-layer depth (MLD) has significant variations, clearly impacting also the DCM position.

Comparing the statistical results obtained from the MLP temperature validation with those shown in other works is not trivial, as most of them refer to global reconstructions or different specific regional oceans (e.g., North Atlantic, Southern Ocean), only few of which are focused on the Mediterranean basin, and also in those cases limited to specific sub-basins and fixed locations (e.g., [20,21]). Overall, most reconstructions show a RMSE close to our (0.7 °C) with values spanning in the range between 0.45 and 1.2 °C [35,73,74,75]. Along the vertical, our MLP errors show similar shape and values of the same order of magnitude of those obtained for different seasons in [75].

As for chlorophyll, in situ and simulated profiles were divided in different classes and averaged to compare the mean profile shape of observations (black line in Figure 8) with those predicted by the network (blue line Figure 8). The classification was based on an incremental step of 1 °C, resulting in 18 different surface categories.

The comparison of MLP outputs with observations is very good in almost all surface temperature ranges. The best MLP performance is observed in cases of low temperatures (e.g., from 12 °C to 18 °C, Figure 8a–f), usually associated with a homogenous and mixed water column. An overlap of the two profiles is observed for all classes down 100 m depth, where the field becomes more stable. There is a slight discrepancy in the first meters of the class 12–13 °C (Figure 8a) which is probably due to the number of profiles available in the training set for this type of surface temperature interval. The same divergence is observed for higher temperature classes (23–25 °C and 28–30 °C, Figure 8n,o,s,t) as a result of two possible causes: the wide variability of the field in these temperature ranges, also suggested by the large standard deviation (represented by the shaded area in Figure 8o,p,t), and the relatively low number of observations available in that range. For higher surface temperatures, for which the presence of a thermocline is more evident, the MLP is able to recognize the depth level where the temperature discontinuity occurs. The magnitude of the thermocline is accurately predicted for most of the classes, with the exception of those referred to 23–25 °C (Figure 8n,o) in correspondence of higher standard deviation. As expected, profiles with a more evident thermocline correspond to higher surface temperatures (28–30 °C, Figure 8q–t), crossing a wide temperature range from 26 to 27 °C in surface to 14–15 °C at 150 m with an inflection point, usually, at about 20 m of depth. For most of the classes the predicted surface value and the observed one are very close, with some exceptions for classes that are less represented in our dataset, such as, e.g., 29–30 °C (Figure 8t).

3.1.3. Sensitivity Analysis of Inputs in the MLP Performance

In order to evaluate the impact of each MLP predictor and define the optimal network architecture, several tests were carried out on different combinations of input surface data (see

Table 2. Different combination of co-predictors used in the sensitivity analysis of inputs in the MLP performance.

	Depth	ChlSAT	SST	Day of the Year	Lat	Lon	ADT	Geostrophic U&V
TEST 1	●	●	●	●	●	●	●	●
TEST 2	●	●	●	●	●	●	●
TEST 3	●	●	●	●	●	●
TEST 4	●	●	●	●	●
TEST 5	●	●	●	●		●
TEST 6	●	●	●		●	●

).

Figure 9 shows the chlorophyll and temperature RMSE profiles with different colors identifying each input combination. For both variables, the curve with lowest RMSE is that corresponding to the configuration of the network, which uses all co-predictors. It is followed by the blue line where the two components (U & V) of geostrophic currents are neglected. However, a slight divergence of dark and blue line is observable in the case of the chlorophyll predictions (Figure 9a), which is not so evident in the temperature retrieval (Figure 9b), suggesting a higher sensitivity of the vertical biomass structure to the geostrophic currents. For both variables, the addition of the ADT determines a strong impact and reduction of the errors, with a RMSE that decreases from 0.44 mgm⁻³ to 0.35 mgm⁻³ and 1.4 °C to 1.13 °C for chlorophyll and temperature, respectively (see red line in Figure 9a,b). The second most influencing variable is the day of the year (green line, in Figure 9a,b), which better captures the signals characterized by a clear annual periodicity [65]. As expected, its impact is evident for the temperature and less for chlorophyll. In fact, Mediterranean Sea surface temperature variability is dominated by a marked annual cycle (e.g., [56]). This signal, though in conjunction with other physical forcing, clearly modulates the vertical stratification and mixing of the water column, directly impacting also the interior temperature structure. Finally, latitude and longitude inputs seem to have the same impact on both predicted variables, with some exceptions for latitude exclusion in temperature modelling that determines a higher RMSE either at surface and bottom layers (yellow line, Figure 9b). As expected, the model sensitivity to the use of different inputs is high from surface to 60 m for both variables and then decreases for chlorophyll, as represented by the closeness of the colored curves in Figure 9a. This is linked to the nature of the two predicted variables: the chlorophyll usually shows larger variations in the first 60 m of the water column, all over the basin, accounting for the limited light-penetration depth and nutrient availability; while, it decreases in concentrations and variability toward the bottom [3,48]. On contrary, the temperature can be strongly variable down to the bottom, reflecting many other physical processes beyond large scale surface heat exchanges, that clearly affect the different areas of the basin (e.g., main currents’ instabilities, (sub)mesoscale features, deep water formation events [76]).

4. Comparison with MEDATLAS Climatology

In this section, an assessment of the MLP performance is carried out comparing predicted 3D fields to the MEDATLAS seasonal and monthly climatologies [77,78]. As a reference database of Mediterranean and Black Sea regions, the MEDATLAS climatology is widely used for the approximation and definition of the vertical profiles of several bio-physical ocean variables at monthly mean and seasonal scales. This database originated from the cooperation of Mediterranean bordering countries within the context of the MEDAR/MEDATLAS-II project, and it is based on the collection and aggregation of multi-disciplinary, in situ hydrographic and bio-chemical measurements ([77], http://modb.oce.ulg.ac.be/backup/medar/contribution.html).

Differing from [48], where an annual climatology was used as reference for the assessment of the vertical chlorophyll estimation, here, seasonal and monthly means were used for the MLP chlorophyll and temperature comparison, respectively. The assessment was carried out on the test set.

The results of chlorophyll comparison are given in Figure 10. The determination coefficient is r² = 0.75 for MLP and r² = 0.32 for climatology, while the RMSE is 0.23 mgm⁻³ and 0.39 mgm⁻³, respectively. The higher accuracy against observed values of the MLP than climatology is also noticeable from the scatterplot in Figure 10a,b where the scatter of the data cloud is reduced for the MLP predictions with respect to climatological estimates.

Generally, all the statistical parameters result in a better accordance of the MLP with the observations with respect to the climatology. This is particularly true for the MeanAPE% which registered a value (MeanAPE% = 49%) one order of magnitude lower than that obtained for climatology (MeanAPE% = 111%). Same behavior can be observed also in terms of vertical variation (Figure 10c,d,f,) with some exception over 80 m for the MBE curve, where the climatological error is closer to zero with respect to MLP.

The RMSE curve shows an MLP error (0.35 mgm⁻³) notably lower than climatology, which may exceed 0.58 mgm⁻³ (Figure 10d). Similar results are found for the MAE (Figure 10e), whose magnitude and shape are very close to those described in [48]. The other MLP statistics and those obtained in [48] are very similar too, even if a direct comparison is hindered by the random selection of stations in the test set and the use in this work of seasonal and monthly means, instead of an annual climatology. As expected, for both datasets, the higher errors occur between surface and 50–60 m of depth, reflecting the strong chlorophyll variability and the large variation of the DCM magnitude and position in this portion of water-column. For RMSE and MAE (Figure 10d,e), the two lines converge to lower errors at the bottom layer.

The results of MLP temperature and climatology comparison are given in Figure 11. A determination coefficient of r² = 0.95 and r² = 0.92 is found for MLP and climatology, respectively. Nevertheless, observing all the other statistical parameters, the climatology clearly shows less accurate results. The RMSE and MeanAPE are 1.13 °C and 4.1% against those of 0.75 °C and 2.8% of the MLP. The MBE displays the most evident difference, which is 0.3 °C for climatology against the −2 × 10⁻⁴ °C for the MLP.

In Figure 11c, the MBE line highlights a significant improvement of the MLP with respect to climatology between 0 and 75 m, attaining surface values of MBE = 0.22 °C for MLP vs. MBE = 1.03 °C for climatology.

For the RMSE (Figure 11d), the climatology shows values over 1.5 °C against the 1.2 °C of the MLP, with maximum errors at the surface (RMSE = 1.6 °C), where the temperature pattern can vary more rapidly. The divergence among the errors of MLP and climatology continues up to 80 m of depth and then decreases towards deeper layers (Figure 11c–e).

The comparison with MEDATLAS demonstrates that the climatology fails to resolve the finer spatial and temporal scales as those retrieved with the MLP. In fact, our MLP matches the real observations with significantly higher accuracy, representing an effective alternative to climatology to overcome the discontinuity and sparseness of the in-situ observations of temperature and chlorophyll.

5. MLP Application on Satellite Data at Basin Scale

In this section, an example of the MLP reconstruction obtained from a satellite climatology is given. Figure 12 shows the MLP basin prediction of the chlorophyll and temperature fields at 3 m of depth. More in detail, the left side represents the monthly longitudinal mean (Figure 12a,c), while the right side shows the corresponding monthly latitudinal mean (Figure 12b,d).

The analysis of the images highlights a good consistency of the surface MLP prediction with the well-known pattern of both variables in the Mediterranean Sea. From a latitudinal point of view (Figure 12a), maxima chlorophyll values are observed at high latitudes (40–44°N) in March and April, when usually large phytoplankton blooms occur and become visible from space [3,64]. This chlorophyll pattern is consistent with that of temperature, which, according to the winter deepening of the mixed layer depth [76], shows the lowest values at the same higher latitudes (40–44°N, Figure 12c). As a mirror, the temperature pattern shows a gradual increase from April to November, along all latitudes (Figure 12c). In fact, according to the warming of surface layers towards the Eastern basin, maxima values are observed from August to September at low latitude (30–36°N, Figure 12c) and high longitudes (30–35°E, Figure 12d) [56,79]. At surface and latitude over than 44.5°N (Figure 12a), the chlorophyll concentrations remain high for all the year, maybe reflecting the influence in the longitudinal average of the more productive coastal waters, as e.g., North Adriatic Sea, that usually shows intraseasonal high biomass levels, observable also from satellite [80]. At Alboràn Sea and Algerian Basin longitudes (from −6°E to 10°E, Figure 12b), the chlorophyll shows high values in winter-early April (from January until April), a decrease in summer (from June to September) and again a slight increase in autumn. This behavior is in accordance with other works [3,81] and is strictly related to the particular dynamic of this area. Nevertheless, the high chlorophyll values registered in summer at very low longitudes (lower than −5°E, Figure 12b) suggest a bias in the MLP prediction in this basin, maybe due to the absence of in situ data useful for the training and the complex hydrodynamics that occur in the Alboràn Sea.

Figure 13 gives the MLP reconstruction of temperature and chlorophyll at 50 m of depth. For chlorophyll, the most evident pattern is a decreasing concentration gradient from North to South, with minima at latitudes lower than 36°N, especially in the months from May to October (Figure 13a). As a counterpart, the temperature shows low values in winter early-spring in those areas of vertical dynamic as that of the North-Western Mediterranean Sea (40–44.5°N, Figure 12c; Figure 13c), where a seasonal thermocline is absent in this period, allowing the uplift of nutrients from the bottom layers. At basin scale, at this depth (50 m), the chlorophyll resulted unbalanced, with higher concentration in the western basin than in the eastern. In fact, as described also in [3], the vertical position of the DCM varies toward the Eastern side, reaching depths higher than $~$70 m, and therefore, becoming less observable at 50 m. This pattern is furthermore highlighted by the chlorophyll discontinuity at 20°E in Figure 13b and Figure 12b and relates to the nature of Levantine basin that shows a deeper DCM along the year with respect to the more productive areas as that of South-Western and North-Western [3,82].

Last study case is reported in Figure 14 that shows the MLP chlorophyll and temperature reconstruction at 100 m of depth. Here, both temperature and chlorophyll patterns are very smooth with slight increase of chlorophyll concentrations in summer (Figure 14a). As described in [3,83], from spring to summer, the DCM remain constrained in a narrow band at $~$100 m along the longitudinal gradient, with a small increase in the Eastern basin, as highlighted by the higher values at 20°E onward observable in panel b of Figure 14.

The results of this reconstruction figure out a promising capability of the network to recreate the principal patterns of both variables. The analysis has demonstrated how these models can embed the seasonal variation at basin and sub-basin scales, just using only surface information.

6. Discussion

The evolution of chlorophyll concentration and temperature 3D fields are expected to be strictly linked, both reflecting the changes driven by ongoing global warming. Phytoplankton is a key player in the marine food webs and monitoring its dynamics is crucial for ecosystem functioning and services assessment [84]. Phytoplankton photosynthetic activity occurs via intracellular pigment as chlorophyll and the dynamics of photosynthetic cells respond to fluctuations of temperature, upper ocean turbulence and nutrient availability [85]. Indeed, climate changes acting on these environmental factors can lead to modifications in the phytoplankton community structure and dynamics, both directly, i.e., through physiology and species interactions, and indirectly, i.e., by modifying water column stratification and nutrients and light availability. The increase of temperature in the water column may affect the organism physiology but also modifies the stratification of the euphotic zone. More stratified waters reduce the mixing in the column limiting the nutrient supply at surface and, thus modulating the algal biomass content, commonly approximated by the chlorophyll concentration [86]. Therefore, retrieving 3D chlorophyll and temperature distribution is key to better monitor phytoplankton blooms, representing a hot spot for primary production, of high relevance for the understanding of ecosystem functioning, fisheries, bio-geochemical cycle dynamics and ocean acidification processes.

As such, accurate estimations of surface and deeper ocean thermal properties become crucial for the understanding of the mechanisms of cause/effect that relate the biomass with the physical processes (vertical motion, horizontal transport, and mixing). Moreover, the modelling of subsurface and deeper ocean temperature allows to better describe the important role which the temperature and other related physical variables (salinity and water density) have in coastal processes, heat exchange, air-sea interaction and climate change [87]. The synergic analysis of different marine environmental variables can lead to an improvement in the understanding of physical and biological interactions in the upper marine layers and their evolution over time and space. In recent years, different techniques have been proposed to exploit the relationship between the sea surface and subsurface parameters to reconstruct the interior structure of the oceans, based on empirical or advanced statistics. Machine learning methods are also emerging as fast and suitable techniques to find a solution to this kind of problems [87].

In fact, even within a complex system as the Mediterranean Sea, neural network models as the one presented here proved able to map the latent relationships between vertical chlorophyll profiles and proper combinations of surface variables (Chl_SAT, SST, ADT, U & V). These results further highlight the utility of such approaches, especially for regional cases, where the ecological dynamics can be very complex.

The error statistics of our MLP predictions are better for the 3D temperature field than for chlorophyll (see MBE%; MeanAPE% and MedianAPE% in Figure 3 and Figure 6). This is easily explained by the complex ecological and biological processes that modulate phytoplankton abundance throughout the water column (e.g., grazing), which cannot be predicted by (relatively) simple environmental data. Nevertheless, a further improvement of our work can be expected by including new satellite-derived data (such as sea surface salinity) as additional predictors in the MLP network, as soon as sufficiently accurate products will become available at regional scale.

As a future work, it would be also very interesting to carry out an intercomparison between the performance of our model and that obtained following other machine learning approaches in the prediction of 3D Mediterranean chlorophyll and temperature structure from regional surface observations.

7. Conclusions

Accurate monitoring of the 3D ocean state is a crucial requirement for marine environment preservation and climate change impact assessment. The use of satellite sensors has notably improved the frequency and resolution of surface ocean observations. However, these measurements are only representative of the surface layers, and cannot provide direct information on the variability of several ocean variables along the water-column. On the contrary, field activities provide accurate vertical estimates, but with low time and spatial coverage. For this reason, the combination of satellite and in situ observations to extrapolate information on the deeper layers through artificial intelligence tools/deep learning approaches is expected to play a key role in the next future.

In this work, we propose a machine learning approach to reconstruct the 3D fields of ocean temperature and chlorophyll concentration in the Mediterranean Sea, starting from surface measurements only. A novel deep Multi-Layer-Perceptron (MLP) network is applied to satellite data to concurrently infer the subsurface structure of these two variables. The model was validated against independent in situ test dataset. The relative errors demonstrated significant improvements in the temperature reconstruction (MBE% = −0.26%, and MeanAPE = 2.72%), but also better predictions for chlorophyll (MBE% = −17 and MeanAPE = 48%). Along the water column, a maximum RMSE of 1.2 °C for temperature and 0.4 mgm⁻³ for chlorophyll was found within the range of about 10–60 m depth, namely where the highest variability is observed in both variables. The results of the temperature model validation agreed with most of similar approaches applied on other ocean/seas [20,75]. The statistics and general fit of chlorophyll validation on test set data demonstrated an improvement with respect to that already obtained in [48], with an r² = 0.71 and MeanAPE = 48% against the r² = 0.69 and MeanAPE = 57% of the network trained on in situ surface values and showed in [48].

Aiming to understand the influence of each co-predictor in the 3D vertical reconstruction and optimize the network architecture, a careful sensitivity analysis was carried out. The configuration of the network that includes all the input variables (day of the year, Lat, Lon, Chl_SAT, SST, ADT and geostrophic U & V) showed the lowest RMSE, both for chlorophyll and temperature. ADT had the strongest impact on the prediction of both variables, as an indirect proxy for the thermocline dynamics, followed by the day of the year, which, as expected, resulted of particular importance for temperature vertical shape.

MLP clearly outperformed MEDATLAS climatology when compared to fully independent test data. This suggests that, for an analysis of processes occurring at finer spatial and temporal scales, the use of our model would provide much better predictions, embedding the wide basin and sub-basin variability of the temperature and chlorophyll, at surface and at deeper layers. This was also supported by the analysis of the MLP application to satellite monthly climatology over 18 years, reliably reconstructing the main patterns of both variables in the basin, both at the surface and along the water column. The predicted 3D fields of chlorophyll showed higher chlorophyll concentrations in those areas characterized by an intense vertical mixing e.g., Gulf of Lion and coastal ones, with an Eastward decrease in surface and a consequent deepening of the DCM. Conversely, the modelled temperature showed lower temperature in winter in the Western basin, with a homogenous vertical pattern and a more pronounced thermocline in summer season in the Eastern basin. The capability of our network to recognize the DCM and thermocline position for the several surface trophic and thermal regimes in different seasons, highlighted its suitability to integrate the high variability of the temperature and chlorophyll vertical fields and its responsiveness to the complex subsurface dynamic of the Mediterranean basin.

Machine learning techniques applied to satellite estimates demonstrate huge and still only minimally exploited potentialities, both as predictive models and to better initialize and validate numerical bio-geophysical models. Indeed, even if a continuous effort is still required for the collection of accurate in situ data at even higher spatial and temporal scales, artificial neural network represent a valid alternative to provide an estimated profile where in situ sampling are missing or difficult.

As a future perspective, it will be of a great interest to test additional predictors, for instance, directly ingesting radiances instead of satellite chlorophyll, or including input from sea surface salinity satellite products. Future works could also extend this approach to the modelling of other ocean bio-chemical and physical variables (e.g., particulate organic carbon, salinity etc.). This will be possible once improved/new satellite products will be developed at regional scale, allowing to develop a cost-effective method to fill the gap of information between the 2D and 3D ocean systems.