The present study is based on the analysis of two numerical simulations: an atmospheric simulation performed by the WRF model forced by ERA5 reanalysis (Hersbach et al., 2020), and an ocean simulation, one-way forced by low level WRF fields, with the ROMS model. The 40-year WRF simulation is described in Miranda et al. (2021), although here the analysis is limited to the 35-year period (1983–2017) due to the time window of the available ocean reanalysis, later described, and the need for an extended spin-up of the ocean model. The one-way coupling from WRF to ROMS is performed by providing the ocean model with three-hourly low-level fields from WRF for wind, relative humidity, air temperature, net shortwave and longwave radiation, pressure, and precipitation, the surface fluxes being computed by ROMS bulk formula. Both simulations (WRF and ROMS) used colocated nested grids, with an inner grid at 3 km and an outer grid at 9 km resolutions, with two-way interaction. The inner domain encloses Madeira and two nearby small islands located to the NW and SW of the main island, Porto Santo and Desertas, respectively. The limits of both computational grids are depicted in Figure 1A.

The initial and boundary conditions of the ROMS ocean model were obtained from the Simple Ocean Data Assimilation (SODA v3.4.2) (Carton and Giese, 2008), a global reanalysis with a 0.25° spatial resolution (approximately 24 km, near Madeira), with data every 5 days. That reanalysis assimilated observed data from the World Ocean Database (WOD) and the International Comprehensive Ocean–Atmosphere Dataset (ICOADS). It is forced by the atmospheric ERA-Interim reanalysis (Dee et al., 2011), using the Coupled Ocean–Atmosphere Response Experiment (COARE4) bulk formula. It was recently shown that SODA-v3.4.2. is an accurate reanalysis, particularly near Madeira (Carton et al., 2019), and while an improved product is not available with the new higher-resolution ERA5 dataset, it constitutes the best compromise for the present study.

At the outer ROMS computational grid boundaries, radiation and nudging conditions were imposed for tracers (temperature and salinity) and the baroclinic current (Marchesiello et al., 2001), a Shecheptikin condition for the barotropic current (Mason et al., 2010), and a Chapman condition for free-surface height (Chapman, 1985). Vertical turbulent mixing was parameterized with the Generic Length Scale (GLS) scheme (Umlauf and Burchard, 2003), with a k–ε closure, where k is the turbulent kinetic energy and ε is the dissipation. ROMS used 40 sigma levels, with an increasing vertical resolution near the ocean surface. The sea floor bathymetric grid is based in a 30 arc-second grid spacing data, from the General Bathymetric Chart of the Oceans (GEBCO; Becker et al., 2009). Several passes of a smoothing filter guaranteed a bathymetry r-factor (r = ∇h/h) below 0.2 (Haidvogel and Beckmann, 1999) to avoid pressure gradient errors associated with sigma coordinate calculations. A thorough description of the ROMS model is presented in Shchepetkin and McWilliams (2005, 2009).

The accuracy of 10-m wind can be assessed from reprocessed Advanced SCATterometer (ASCAT) data, available at a 12.5-km swath grid (Verhoef et al., 2012). Here, we used data from the ascending passage (descending passage gives similar results); that in Madeira occurs between 21 and 23 h UTC and WRF data at 21 h UTC, for the longest available period coincident with the numerical simulation, from 2007 until 2017. After remapping the simulated variables to the satellite grid, the remote-detected values were used to compute the following statistical parameters: bias, time correlation, and normalized standard deviation. The statistical parameters show the high accuracy of the simulated surface wind (Table 1), at four locations near Madeira: points 1 and 2 are located near the maximum tip-jets mean wind speed, point 3 in an upstream region not affected by the island influence on surface wind, and point 4 in the wake (see Figure 1B for exact position of points). The computed biases vary between −0.19 and −0.90, indicating a very slight underestimation of wind by the model. The time-correlation values are high, varying between 0.90 and 0.76, with the worst value in the wake. The simulated wind speed variability is close to observations.

The simulated SST was compared with the dataset created by the European Space Agency (ESA) Climate Change Initiative/Copernicus Climate Change Service (CCI/C3S), with a 0.05° resolution, combining SST retrievals from the Advanced Very High-Resolution Radiometer (AVHRR) and the Along Track Scanning Radiometer (ATSR) series of satellite sensors. For a full description of the dataset, the reader is referred to Merchant et al. (2019). As for surface wind, the statistical parameters show the robustness of the simulated values (Table 2), with an absolute bias lower than 0.23 °C in the four points, a very high synchronization between model and satellite data revealed by time correlations between 0.96 and 0.97, and a realistic simulated variability noted by a ratio of standard deviations of 1.04 or 1.05. The robustness of the simulated SST is also revealed by the mean spatial SST bias and root mean square error (Figure 2) with values always lower than 0.3 and 0.8 °C around Madeira, respectively. In both parameters, the highest values are located near the west flank, coinciding with high SST variability and gradients (shown latter) that are not accurately captured by daily 0.05° resolution satellite SST.

The presence of Madeira Island is especially clear in the perturbation it imposes to the wind. The seasonal mean wind around Madeira (Figure 3) displays a rather constant pattern, with northeasterly flow impinging the island throughout the year, leading to the establishment of two localized wind speed maxima near the east and west island tips. The upstream mean wind speed is stronger in winter (DJF) and weaker in summer (JJA). However, the tip-jets are most prominent in summer.

The establishment of the summer tip-jets was discussed in some detail by Miranda et al. (2021), and both its intraseasonal and decadal variability were found to be controlled by the regional PBL height, defined as the height of maximum stability in the low troposphere. Figure 4 shows the seasonal mean potential temperature gradient in the north upstream sampling point (Figure 1B), indicating a much more stable and lower height PBL in summer, lowering from about 1,700 m in winter (DJF) to about 1,000 m in summer (JJA). The lower summer value is below the height of the main island plateau, implying the need for flow deflection around the island with streamline squeezing near the island tips, a simple Bernoulli mechanism to accelerate the flow in the jets.

On the other hand, the seasonal mean wind directional constancy (Bromwich, 1989), computed as the ratio of the magnitude of mean wind vector to the mean wind speed:

represented in Figure 5, also shows a strong seasonality, with a much stronger constancy in summer (JJA) contributing to more frequent northeasterly flow and favoring the establishment of the tip-jets in their preferred locations.

Summer is also the season with the stronger wind speed variability in the jets, associated with a multiweek oscillation of the jet intensity, driven by a corresponding oscillation in the regional PBL height (Miranda et al., 2021). The seasonal mean intraseasonal variability is shown in Figure 6, with the standard deviation of wind speed exceeding 3.5 m/s in the east jet, slightly smaller values in the west jet, and much smaller values in the upstream flow. In the following section, we will concentrate on summer results.

Figure 6 also displays the intraseasonal summer variability of SST. While SST variability is driven by wind, there is a clear mismatch between wind variability, with two maxima in the two jets, and SST variability, with a single maximum under the west jet. This asymmetry between atmospheric forcing and the ocean response is explained by the asymmetry in the upper ocean circulation, later discussed, implying a stronger thermal response to wind variability near the west tip-jet.

The perturbed surface wind in the vicinity of Madeira imposes a very heterogeneous forcing on the upper ocean layers, through energy and momentum fluxes, which is reflected in the ocean surface temperature and 3D currents, particularly during summer, with higher wind speed in the tip-jets and lower wind speed in the wake. This wind environment favors the onset of cyclonic current anomalies near the west flank and anticyclonic anomalies near the east flank (Figure 7), which are the source of the frequent shedding of ocean eddies usually observed in the lee of Madeira (Caldeira et al., 2014). As also shown in Figure 7, this positive (negative) vorticity anomalies imply colocated upward (downward) currents near the west (east) flank. It is interesting to note that the small nearby islands produce qualitatively similar dipolar anomalies of vorticity, of a much smaller horizontal scale, as previously suggested by Caldeira and Sangrà (2012) (numerical) study.

The asymmetry in the vertical current is also observed in the chlorophyll concentration, assessed by remote sensing data, particularly in the 2016 summer (a summer with intense tip-jets) characterized by high concentration anomaly near the west flank (Figure 8). We hypothesize that this may be due to the local upwelling intensification observed in the simulation (not shown).

The effect of the tip-jets is clearly visible in the summer mean SST and mixed layer depth (MLD). The latter was computed considering a threshold of | Δρ| < ± 0.3 kg/m³ in the upper density vertical profiles (chosen considering the high stratified upper ocean layers near Madeira, during summer). In Figure 9A, the ocean surface is locally cooler at both jets, as one would expect from enhanced evaporative cooling, but cooling is much stronger near the west jet because extra cooling is provided there by local upwelling (in a stably stratified ocean), as shown in Figure 7, whereas the local downwelling (associated with the east anticyclonic eddies) produces no relevant signature. An asymmetric pattern is also observed in the MLD (Figure 9B) with shallower mixed layer in the west upwelling region and a deeper one in the east downwelling sector.

As shown in the MLD spatial distribution and ocean vertical currents, the tip-jets effect on the ocean goes well beyond the ocean surface. To analyze the undersurface effect in the upper ocean stability, Figure 10 shows vertical cross-sections of ocean static stability, given by squared Brunt–Väisälä frequency in the zonal and meridional directions, following the lines shown in Figure 9A. During summer, near the east jet (Figure 10A), a locally decreased vertical stability is shown by a deepening of the pycnocline. The opposite is observed near the west jet (Figure 10B). At both locations, the perturbations in the top 30 m are compensated by opposite perturbations below, with respect to surrounding water. Below 120 m, the stability field is almost flat (not shown). A detailed look shows that the maximum impact in the upper ocean (black triangles) is located near the maximum wind speed variability and not in the maximum mean wind (red triangles) (cf. Figure 9A), since the maximum wind variability coincides with maxima in ocean vorticity and vertical current, showing the importance of the tip-jets and wake for the upper ocean imposed perturbations. This is in line with the Couvelard et al. (2012) study, which suggested that there was a transfer of momentum (wind stress curl) from the atmosphere to the ocean, which contributed to the generation of the oceanic eddies at the island flanks. A rather different pattern is shown in the corresponding winter cross-sections (Figures 10C,D) with a much better mixed upper ocean and lower pycnocline (dashed line), probably due to the lower incident radiation during winter and the better ocean capacity to transport the heat to its interior. As a consequence, during this season, the tip-jets impact in the upper ocean is attenuated.

The shallow MLD observed in summer means a lower ocean heat capacity and consequently a higher SST sensitivity to atmospheric warming. However, there are, at least, two mechanisms that may lead to an attenuation of the ocean surface warming. One is an increased evaporative cooling that would be produced by tip-jet intensification (Miranda et al., 2021); the other is the intensification of the upwelling current associated to cyclonic vorticity near the west tip-jet (not shown), which may lead to a further attenuation of SST warming or even to local cooling. Both processes are controlled by the interdecadal variability of the tip-jets.

Figure 11 shows trends of summer 10-m wind, computed by the nonparametric Theil–Sen estimator (Sen, 1968, Theil, 1992), after filtering time series at each grid point by removing the mean annual cycle. This technique computes a central estimate of the trend and upper and lower estimates at the (chosen) 95% confidence level. In Figure 11, the central estimate is only shown when upper and lower estimates have the same sign. That is the case in both jets, which are found to have been intensifying slightly more than 0.7 m/s in 35 years.

Results using the same method (Theil–Sen) applied to summer SST are shown in Figure 12. In this case, all estimates are presented, as the underlying positive trend would mostly avoid negative estimates even from the lower estimate. The mean warming trend is found to be around +0.5°C/35 years but with a lot of spatial heterogeneity. Trends in the west jet are, however, much smaller with a minimum approaching zero, indicating that enhanced upwelling at those locations has been able to compensate most of global warming.

The ocean temperature trend at 20 m also shows high spatial heterogeneity, with zones where the warming is more intense than at surface and zones where the inverse is observed, particularly near the west flank upwelling zone. Interestingly, the mean ocean temperature warming trend at 20 m (Figure 12E) is slightly higher (about +0.06°C per 35 years) than on the surface, suggesting a decoupling between the surface processes expose to daily wind and thermal variations and the interior oceanic mixed layer, as suggested by Azevedo et al. (2021)