*3.1. Derivation of the Optimal Currents*

Following PIT09 and RS18, we derive the surface circulation in the Mediterranean Sea combining the CMEMS geostrophic velocities (GV hereinafter) and the SST fields. The Optimal reconstruction is based on the inversion of the ocean heat conservation equation in the mixed layer. The equation, indicated below, is inverted for inferring information on the sea surface currents:

$$\frac{\partial \text{SST}}{\partial \mathbf{t}} + \mathbf{u} \frac{\partial \text{SST}}{\partial \mathbf{x}} + \mathbf{v} \frac{\partial \text{SST}}{\partial \mathbf{y}} = \mathbf{F} \tag{1}$$

In Equation (1), (u,v) are respectively the zonal and meridional component of the ocean surface flow, (x,y) are the zonal and the meridional directions and F is the forcing term, containing the source and sink terms for the heat conservation equation. The resulting expressions for the optimally reconstructed velocity (uopt, vopt) are the following:

$$\mathbf{u\_{opt}} = \mathbf{u\_{bck}} + \mathbf{u\_0}\sin\phi + \mathbf{v\_0}\cos\phi = \mathbf{u\_{bck}} + \mathbf{u\_{COR}}$$

$$\mathbf{v\_{opt}} = \mathbf{v\_{bck}} - \mathbf{u\_0}\cos\phi + \mathbf{v\_0}\sin\phi = \mathbf{v\_{bck}} + \mathbf{v\_{COR}}\tag{2}$$

The Optimal currents are obtained adding the correction factors uCOR and vCOR to the background ("bck") geostrophic velocities. In practice, Equation (2) is the formal expression of the synergy between the altimeter-derived SSH and the SST observations for the sea surface currents reconstruction. We provide a brief illustration of the terms appearing in the correction factors of Equation (2), referring the reader to PIT09 for more details. The correction velocities u0 and v0 have the following expressions:

$$\begin{aligned} \mathbf{u}\_0 &= \frac{\mathbf{f}(\min(\boldsymbol{\beta}, \mathbf{q})) - \mathbf{f}(\max(\mathbf{a}, -\mathbf{q}))}{\mathbf{g}(\min(\boldsymbol{\beta}, \mathbf{q})) - \mathbf{g}(\max(\mathbf{a}, -\mathbf{q}))} \\ \mathbf{v}\_0 &= \mathbf{p} \mathbf{u}\_0 \end{aligned} \tag{3}$$

In addition, *φ* is a function of the SST spatial derivatives, as expressed by Equation (4):

$$\phi = \operatorname{atan}(\mathbf{A}/\mathbf{B})\tag{4}$$

where A = *∂x*SST, B = *∂y*SST and (x,y) are the generic zonal and meridional coordinates. Moreover, in Equation (3) the functions f, g (expressed as functions of a generic variable *γ*) and the coefficients p, q, *α* and *β* are given by the set of Equation (5):

$$\begin{aligned} \mathbf{f}(\gamma) &= -2(\mathbf{q}^2 - \gamma^2)^{3/2}/3; \\ \mathbf{g}(\gamma) &= \mathbf{x}(\mathbf{q}^2 - \gamma^2) + \mathbf{q}^2 \mathbf{a} \mathbf{s} (\gamma/\mathbf{q}); \\ \mathbf{p} &= \sin \phi \cos \phi (\sigma^2 \mathbf{v} - \sigma^2 \mathbf{u}) \mathbf{q}^{-2}; \\ \mathbf{q} &= \sqrt{\sigma^2 \mathbf{u} \sin^2 \phi + \sigma^2 \mathbf{v} \cos^2 \phi}; \\ \mathbf{a} &= (\mathbf{A} \mathbf{u}\_{\text{bck}} + \mathbf{B} \mathbf{v}\_{\text{bck}} + \mathbf{E} - \mathbf{h}) / \sqrt{\mathbf{A}^2 + \mathbf{B}^2}; \\ \beta &= (\mathbf{A} \mathbf{u}\_{\text{bck}} + \mathbf{B} \mathbf{v}\_{\text{bck}} + \mathbf{E} + \mathbf{h}) / \sqrt{\mathbf{A}^2 + \mathbf{B}^2}; \end{aligned} \tag{5}$$

where h is the error associated to the forcing term F; *σ*u, *σ*<sup>v</sup> are respectively the errors for the zonal and meridional background velocities, and E depends on the temporal derivatives of SST, i.e., E = *∂*tSST-F. The determination of the input parameters included in Equations (1)–(5) is illustrated in Sections 3.2 and 3.3.

#### *3.2. The Forcing Term F and the Computation of the Associated Error h*

In RS18, the forcing term *F* was approximated as the satellite-derived, low-pass filtered SST temporal derivatives, only retaining spatial scales larger than 500 km. This choice is due to the following two facts:


For the case of the Mediterranean basin, F was identified as the satellite SST temporal derivatives, *∂*tSST, low-pass filtered at 400 km. This was achieved filtering the *∂*tSST at several spatial scales and choosing the one that minimized the in-situ forcing estimates obtained with the OGS drifter-derived observations.

Using 30 years of drifting buoys-derived SST in the Mediterranean basin (data doi = 10.6092/7a8499bc-c5ee-472c-b8b5-03523d1e73e9), we computed an in-situ derived forcing term (Fdri), whose formal expression is given by Equation (6).

$$\mathbf{d}\_{\rm t} \mathbf{S} \mathbf{S} \mathbf{T}\_{\rm ins} = \frac{\partial \mathbf{S} \mathbf{S} \mathbf{T}\_{\rm ins}}{\partial \mathbf{t}} + \mathbf{u} \frac{\partial \mathbf{S} \mathbf{S} \mathbf{T}\_{\rm ins}}{\partial \mathbf{x}} + \mathbf{v} \frac{\partial \mathbf{S} \mathbf{S} \mathbf{T}\_{\rm ins}}{\partial \mathbf{y}} + \mathbf{w} \frac{\partial \mathbf{S} \mathbf{S} \mathbf{T}\_{\rm ins}}{\partial \mathbf{z}} = \mathbf{F}\_{\rm ins} \tag{6}$$

where the subscripts "ins" and "t" respectively stand for in-situ measured and derivative with respect to time.

For every drogued buoy trajectory, we evaluated the in-situ daily dtsst as

$$\mathbf{d}\_{\rm t} \mathbf{S} \mathbf{T}\_{\rm inv} = \left[ \mathbf{S} \mathbf{S} \mathbf{T}\_{\rm inv} (\mathbf{t} + \mathbf{d} \mathbf{T} / 2) - \mathbf{S} \mathbf{S} \mathbf{T}\_{\rm ins} (\mathbf{t} - \mathbf{d} \mathbf{T} / 2) \right] / \mathbf{d} \mathbf{T}$$

where dT = 1 day. Moreover, taking advantage of the 1986 to 2017 satellite-derived SST (dataset labeled with "2" in Section 2) we could compute the daily remotely sensed *∂*tSSTsat, including the filtered maps at scales ranging from 50 to 600 km. This operation enabled to derive the quantity L-*∂*tSSTsat, with L ∈ [50 km , 600 km]. Afterwards, each of the filtered L-*∂*tSSTsat was spatially and temporally colocated with respect to the drifters estimates, i.e., our benchmark. This allowed to generate two-dimensional maps of Root Mean Square Errors (RMSE) between the satellite and the in-situ-derived SST temporal derivatives. In the end, computing the mean, basin-scale RMSE between each of the L-*∂*tSSTsat and the dtSSTdri, we found that for L = 400 km such RMSE reached a minimum (Figure 1) letting us identify it as the optimal filtering scale. In addition, we assumed that the RMSE error map between the 400 km-DtSST*sat* and the *∂*tSST*dri* is an estimate of the forcing term error h appearing in Equation (2). A map of h, binned in 1◦ × 1◦ boxes, is given in Figure 2.

**Figure 1.** RMSE between the satellite and the in-situ derived forcing term as a function of the D*t*SST*sat* filtering scale.

**Figure 2.** RMSE between F and F*dri*, binned on 1◦ × 1◦ boxes.

The order of magnitude of the RMSE is 10−<sup>6</sup> ◦C·s−<sup>1</sup> over most of the Mediterranean Sea and its spatial variability is mostly correlated with the major dynamical features in the Mediterranean Sea (e.g., the Algerian current, the Liguro-Provenzal, the Atlantic Ionian Stream and Western Adriatic current) where the values can also exceed 1.5 ◦C·s<sup>−</sup>1.
