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

: 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 e ﬀ ective 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 ﬁelds of ocean temperature and chlorophyll- a concentration, two variables of primary importance for many upper-ocean bio-physical processes. Artiﬁcial neural networks can e ﬃ ciently 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 proﬁles and several di ﬀ erent combinations of satellite-derived surface predictors are used to identify the optimal model conﬁguration, 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, signiﬁcantly improving with respect to Mediterranean climatological data (MEDATLAS). 3D ﬁelds 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 proﬁles and their spatial and temporal variability. It thus represents an e ﬀ ective 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.


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. 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, 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.

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.
Remote Sens. 2020, 12, 4123 4 of 27 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.

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).
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 30 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.

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). Figure 1. Spatial distribution of the stations used for the training set (red triangle, ~70% of the total stations) and test set (blue triangle, corresponding to ~30% of the total stations). The selection between training and test stations was random.

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

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)(1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(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.

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].

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(2π * (day o f the year/365)) (1) day2 = sin(2π * (day o f the year/365)); Then, all inputs (except day1 and day2) were standardized by applying a z-score transformation as follow: where the µ and σ are respectively the mean and standard deviation for each co-predictor of the training dataset, while x i is a generic iinput 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 −3 up to values even considerably higher than 10 mgm −3 .
Obviously, once completed the MLP training, the outputs were converted into original values applying the z-score inverse formula based on reference µ and σ of training targets.

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].
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: with γ representing the performance ratio, a regularization parameter that defines the penalty imposed to the weights and biases, and = defining the mean of the sum of the squared weights and biases ( ). 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.  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: with γ representing the performance ratio, a regularization parameter that defines the penalty imposed to the weights and biases, and msw = 1 n n .

J=1
w 2 j defining the mean of the sum of the squared weights and biases w j . The regularization effect forces the network to adopt smoother weights. This modified Remote Sens. 2020, 12, 4123 8 of 27 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.

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 µ and σ), 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 µ and σ of targets of training dataset. The metrics used for this analysis are reported in the following Table 1: 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. 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 2 = 0.73; test, r 2 = 0.71), while the RMSE is relatively low with values of 0.20 mgm −3 for training and 0.21 mgm −3 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 −3 and test, MBE = 0.04 mgm −3 ) 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. 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 2 = 0.73; test, r 2 = 0.71), while the RMSE is relatively low with values of 0.20 mgm −3 for training and 0.21 mgm −3 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 −3 and test, MBE = 0.04 mgm −3 ) 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 2 = 0.71) is higher than r 2 = 0.69 of [48], and the relative percentage errors (MBE% = 17%; MeanAPE% = 48%; MedianAPE% = 30%) and RMSE = 0.21 mgm −3 are also lower (MBE% = 28%; MeanAPE%= 57%; MedianAPE% = 34%; RMSE = 0.26 mgm −3 ). 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 2 (0.71 against 0.63 of [48]) and RMSE (0.21 mgm −3 against 0.23 mgm −3 of [48]). Figure 4 shows RMSE profile computed on the test set as a function of depth. The higher errors (RMSE = 0.35 mgm −3 ) 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 −3 ) towards deeper layers, reflecting the lower chlorophyll variability and concentrations far from the surface [3].  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 2 = 0.71) is higher than r 2 = 0.69 of [48], and the relative percentage errors (MBE% = 17%; MeanAPE% = 48%; MedianAPE% = 30%) and RMSE = 0.21 mgm −3 are also lower (MBE% = 28%; MeanAPE%= 57%; MedianAPE% = 34%; RMSE = 0.26 mgm −3 ). 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 2 (0.71 against 0.63 of [48]) and RMSE (0.21 mgm −3 against 0.23 mgm −3 of [48]). Figure 4 shows RMSE profile computed on the test set as a function of depth. The higher errors (RMSE = 0.35 mgm −3 ) 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 −3 ) towards deeper layers, reflecting the lower chlorophyll variability and concentrations far from the surface [3].

Determination Coefficient
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). 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 −3 and 0.04-0.08 mgm −3 (Figure 5a,b) and other classes of higher surface concentrations as in 0.12-0.2 mgm −3 and 0.2-0.3 mgm −3 (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. 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 −3 and 0.04-0.08 mgm −3 (Figure 5a,b) and other classes of higher surface concentrations as in 0.12-0.2 mgm −3 and 0.2-0.3 mgm −3 (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. 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%.

Temperature Assessment
Both determination coefficients are high (training, r 2 = 0.94; test, r 2 = 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. 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 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 at higher surface temperatures (28-30 • C, Figure 8q-

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). Table 2. Different combination of co-predictors used in the sensitivity analysis of inputs in the MLP performance.

Depth
Chl SAT SST Day of the Year Lat Lon ADT Geostrophic U&V 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 −3 to 0.35 mgm −3 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]).

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 2 = 0.75 for MLP and r 2 = 0.32 for climatology, while the RMSE is 0.23 mgm −3 and 0.39 mgm −3 , 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 −3 ) notably lower than climatology, which may exceed 0.58 mgm −3 (Figure 10d). Similar results are found for the MAE (Figure 10e), whose magnitude

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 2 = 0.75 for MLP and r 2 = 0.32 for climatology, while the RMSE is 0.23 mgm −3 and 0.39 mgm −3 , 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.
[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. determination coefficient of r 2 = 0.95 and r 2 = 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 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 −3 ) notably lower than climatology, which may exceed 0.58 mgm −3 (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. determination coefficient of r 2 = 0.95 and r 2 = 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 −4 • C for the MLP.
Remote Sens. 2020, 12, x FOR PEER REVIEW 19 of 30 the MLP. The MBE displays the most evident difference, which is 0.3 °C for climatology against the −2 × 10 −4°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,d,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 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.

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).   [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 Figures 13b and 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.

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, biogeochemical cycle dynamics and ocean acidification processes.
As such, accurate estimations of surface and deeper ocean thermal properties become crucial for 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.

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 Figures 3 and 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.

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 −3 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 2 = Remote Sens.0.71 and MeanAPE = 48% against the r 2 = 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.  introducing, guiding and teaching her the machine learning techniques applied to the ecology and for always having supported her with his scientific knowledge. M. Sammartino thanks Marco Bellacicco for the fruitful review and comments that improved the final version of this paper. She also thanks the research group of the GOS of the Italian CNR, in particular Simone Colella, for his technical support and the helpful advices. We also thank Marco Bracaglia, Annalisa Di Cicco, Francesca Elisa Leonelli and Chiara Marullo for their useful discussions and suggestions.