Remote sensed SST and surface wind are used to verify the model accuracy. Wind observations are from the Advanced Scatterometer (ASCAT) ascending path dataset, on an approximately 12 km grid (Verhoef et al., 2012). The analysis of SST uses the Multi-scale Ultra-high Resolution (MUR) dataset, that is made of foundation night SST values, with a daily periodicity, on a 0.01° grid. It comprises MODIS, AVHRR, microwave and in-situ retrievals. A comprehensive description of MUR can be found in Chin et al. (2017). Both wind and SST are compared on the model grid, by interpolating the 12 km satellite wind product and averaging the 0.01° SST satellite product.

Wind and SST results are presented, side by side, in Figure 2 . Qualitatively, the mean WRF wind taken at 00h UTC ( Figure 2B ) is very similar to the ASCAT-12 km field ( Figure 2A ), with strong signatures of the coastal jet at Cape Ghir, and of tip jets downstream of Madeira and of the main Canary Islands. WRF wind speed bias ( Figure 2C ) is below ±0.5 m.s^-1 away from the coasts, slightly negative in the wakes and the upstream blocked flow of the islands, slightly positive in the tip jets and on the African coast. These differences mean that WRF perceives a slightly higher impact from topography on the wind field than observed by ASCAT, which is consistent with its higher resolution and with known ASCAT sampling issues near the coast. It is important to mention that a gap free satellite product that is available on a 0.25° grid (Bentamy et al., 2021) fails to represent the tip jets in some islands, smoothing out almost all jet-like features (not shown). To further assess the realism of the simulated surface wind, we compared it with data from two meteorological stations of the Portuguese Institute for Sea and Atmosphere (IPMA) located near the East and West flanks of Madeira, which revealed a bias of -0.79 and -0.21 m.s^-1 and a time correlation of 0.8 and 0.7, for the East and West meteorological stations, respectively. Moreover, the in-situ observed wind was also compared with satellite values and showed a similar or a higher bias than the one obtained with simulated data. For instance, in the West flank meteorological station the satellite bias is almost 3 m.s^-1 higher than that of the modeling wind (see Figure S1 ). This seems to suggest that some of the difference between model and satellite based data, as seen in Figure 2C , could result from the lack of resolution of satellite data.

For SST, model results are taken at 00h UTC. Observations ( Figure 2D ) indicate a strong signature of coastal upwelling, with a persistent filament offshore Cape Ghir and a consistent zonal SST gradient perturbed by the presence of the islands. Coastal features in ROMS ( Figure 2E ) are qualitatively similar to observations ( Figure 2D ), but more intense. There is a small warm bias ( Figure 2F ) offshore, only slightly exceeding 1°C outside the Canary current. Near the coast there is, however, a stronger cold bias, exceeding -3°C, in places, that is consistent with the positive coastal wind bias, which result in an intensified upwelling. The coastal cold bias is possible also partially explained by the gap filling in the MUR dataset in the presence of low clouds, that tend to intensify in conditions of coastal upwelling (Samelson et al., 2021). Overall, the mean computed SST bias is +0.28°C. For SST validation, we also used in-situ data retrieved from oceanographic buoys maintained by Puertos del Estado Institute. The buoys are located near La Palma, Gran Canaria and Tenerife (see islands locations in Figure 1B ), and revealed low cold biases of -0.30°C, -0.31°C and -0.54°C, respectively.

The atmospheric flow drives the regional ocean circulation by wind stress, and the patterns of wind stress curl force regions of upwelling and downwelling in the coastal ocean. The simulated wind stress curl, represented in Figure 3 , has a lot of small-scale noisy features, but is mostly characterized by a slight negative curl in regions away from coastal effects, a strong positive curl on the African coast, and strong dipolar features downstream of the islands. The background negative curl is the signature of the Azores subtropical anticyclone. The strong positive curl on the African coast is due to the cross-coast wind drop-off, as also found in other major upwelling systems, such as California (e.g., Renault et al., 2016a). A distinctive feature of the CANUS is the presence of islands with significant height, generating the strong curl dipoles, therefore confirming a rule of thumb in the upwelling near the islands by a -0.30 time-correlation between wind stress curl and SST nearby Gran Canaria west flank.

The offshore extension of the mean coastal positive curl anomaly is highly variable along the African coast, with potential impact in the regional ocean. Following Renault et al. (2016a), we define the wind drop-off extension as the distance from the coast of a reference contour line of the wind stress curl, here chosen as 3×10^-7 N.m^-3 (cf. Figure 3 ). Different processes may control that extension: coastal topography, coastline orientation in relation to the prevailing wind, and maybe transient processes occurring along the coast (Dorman, 1985; Beardsley et al., 1987; Dorman, 1987) or larger scale subtropical processes unconfined to the coast (Taylor et al., 2008). The maximum drop-off extension is found south of Cape Ghir, coinciding both, with the wake of the southern edge of the Atlas Mountains, and a significant perturbation of the coastline geometry. The latter has a direct impact on the wind direction relative to the coast, a parameter that will control the efficiency of Ekman pumping to generate upwelling. As shown in Figures 2A, B, C, E , the regions of intensified coastal upwelling coincide with sections of the coast which are almost parallel to the prevailing wind.

The spatial distributions of the summer mean latent and sensible heat fluxes ( Figure 4 ) show the combined effect of coastal upwelling and low-level jets. Both fluxes are positive in the offshore region, where the ocean is warming and moistening the atmosphere. The latent heat flux presents maxima near most islands (except the two closest to the coast) coinciding with the tip jets and decreases to very low values close to the continental coast. The sensible heat flux is intensified in the Madeira eastern tip jet but is, however, very small near the Canaries and negative near the African coast, where the atmosphere is being cooled, indicating that the SST cooling by upwelling is more important than the increased transfer coefficients at the low-level jets.

The strong signature of both, the African coast and of the different islands, in the atmospheric circulation, with quasi-permanent maxima of low-level jets at preferred locations, even with time varying oscillations in their intensity, suggests the need to look at the spatial structure of the atmosphere and of the ocean boundary layers. The atmosphere-ocean interactions lead to a strongly perturbed mean spatial distribution of the boundary layer in the atmosphere ( Figure 5A ) and ocean ( Figure 5B ). The atmospheric PBL height varies widely from near 200 m in the more active upwelling spots on the African coast to near 1400 m in the NW of the domain, with PBL subsidence downstream of the islands and main capes. The ocean Mixed Layer Depth (MLD) is roughly a mirror image of the former, with deeper MLD in regions of stronger wind, where PBL height is lower, an indication of the role of wind stress in driving the vertical mixing in the upper ocean.

All major coastal upwelling systems develop low-level jets, usually observed at about 400m height (Lima et al., 2022), that eventually may interact with the islands jets. The focus here is in the intra-seasonal variability of those jets, directly associated with the corresponding variability of the underlying ocean. By intra-seasonal variability, we mean all variability with timescales less than 3 months. Miranda et al. (2021) found that, in the case of the Madeira tip jets, the variability is controlled by oscillations in the Azores Anticyclone, through a variability of the Planetary Boundary Layer (PBL) height at a time scale of weeks. Figures 6 and 7 show that such process extends throughout the Canary system, affecting both the tip jets of the different islands and the African coastal jet. While the jets are quasi-permanent features of the summer circulation, their intensity varies at week to multiweek time scales by a large factor.

The intra-seasonal variability of the jets is larger at Madeira, the island less affected by the African coastal jet, but all jets seem to oscillate at comparable timescales, and they all occasionally relax to low intensity ( Figure 6 ). Considering the extension of CANUS, one would not expect a perfect synchronization between the jets, neither a simple connection to the PBL height at a single representative upstream location as found for Madeira (Miranda et al., 2021). Instead, Figure 7 proceeds to a more regionalized view by looking at composite means of the low-level wind speed for the lower tercile of the mean PBL height along the 31N parallel (Cape Ghir, cf. Figure 1B ) and the corresponding upper tercile. With a lower PBL height ( Figure 7A ), the coastal jet widens throughout the CANUS almost to Madeira longitude, with strong tip jets in all islands; with higher PBL height ( Figure 7B ) all jets attenuate, and the continental coastal jet is more evident close to Cape Ghir and south of the Canaries. The intense dependence of surface wind on PBL height is also evident in the terciles differences to the mean wind speed ( Figures 7C, D ). In these panels the highest differences (about +- 2 m.s^-1) are observed near the islands, particularly the Canaries, revealing the importance of island orography on the PBL height and consequently on the surface wind. Moreover, it also shows the importance of a high-resolution simulation to assess the islands effect, otherwise the islands orography and its impact would be too smooth.

The wide oscillation that takes place at the multiweek time scale in CANUS, is even more evident in Figure 8 , showing vertical wind speed composite means for the lower and upper terciles of PBL height at 31N and at 00h UTC. In these panels, the PBL height varies by a factor of 2 in the whole zonal section. The jet maximum always remains close to the top of the PBL. The abrupt lowering of the PBL near the coast fits the conceptual model proposed by Beardsley et al. (1987) from observations, at a time without the support of global reanalyses.

A more detailed view of the along coast wind ( Figure 9A ) reveals that most of the time the wind flows southwards along the African coast, with a clear diurnal cycle, apart from the bay south of Cape Ghir (30N) where a cyclonic circulation imposes a northward flow. Besides the diurnal cycle, visible in the power spectral density ( Figure 10 ) by a peak at 1.16 × 10^–5 Hz (24 h period) and of the two main sub-daily harmonics at 2.31 × 10^–5 Hz (12 hours), at 3.52 × 10^–5Hz (8 hours), it is also noted a peak at 9.26 × 10^–7 Hz, corresponding to a multi-weekly frequency. This latter peak is in agreement with the week and multi-week wind variability noted before in the surface wind time series ( Figure 6 ). The cape Ghir nearby high orography and the regional strong coast veer lead to a high variability in the along coast wind at 31N, noted by modes of oscillation clearly more intense than at other latitudes (green line in Figure 10 ), and by computed standard deviation and variance local maximums (not shown). The computed narrow 90% confidence band of each line (not shown), demonstrate the peaks statistical accuracy.

The evolution of the along coast wind concerning the full JAS period of 2019 ( Figure 9A ) also shows at least 3 northward propagating wind reversals, from south of Cape Ghir to the northern limit of the domain, lasting a few days. From 1979 until 2018, the northwest Africa, summer (JAS) Hovmöllers of ERA5 surface wind along the coast, show 2-3 northward wind reversals per summer season (not shown) despite the 30 km data spatial resolution, which shows the recurrence of this phenomena in the studied region.

However, the wind reversal events are not simultaneously observed along the coast. Indeed, a detailed view of the July event ( Figures 11A, B ) puts in evidence a clear northward propagation of the wind reversal accompanied by a less clear, but still visible, corresponding propagation of a positive anomaly in the PBL height. The wind speed perturbation propagates relatively to the ground at a mean speed of 3.8 m.s^-1, as shown by the star’s positions, in panel 9A, that show the maximum wind northward propagation.

Dorman (1985; 1987) and Beardsley et al. (1987) identified in the California current upwelling system the occurrence of similar short-lived episodes of near coast wind reversal, accompanied by changes in the atmospheric PBL height, propagating northward. These were described as episodes of upwelling relaxation and were attributed to coastal trapped atmospheric events.

Following the conceptual 2-layer model discussed in Dorman (1985), after analysis of potential temperature profiles in Cape Ghir zone, we computed a mean value of 293.5 K for the potential temperature of the lower layer (marine layer), that goes from ocean surface until a Brunt-Väisälä frequency square local maximum, and, a mean value of 302.5 K for the upper layer, that goes from the upper limit of the lower layer until an upper Brunt-Väisälä frequency square local maximum. These potential temperature mean layer values and the local mean 250 m PBL height, observed near cape Ghir (black rectangle in Figure 5A ) permitted to compute the Kelvin wave phase speed, following eq 4 of Dorman (1985) (eq S1 in the supplementary material ) with a value of 8.6 m.s^-1. But because the perturbation travels against a mean upstream flow of about 4.6 m.s^-1, between 30N and 33.2N, during the 10 days preceding the analyzed event (black rectangle in panel 9A), the perturbation ground-relative speed is near 4 m.s^-1, a value close to the simulated value (3.8 m.s^-1) of propagation to the north of maximum southerly wind.

Further evidence of a Kelvin wave is the geometry of the northward propagation event, clearly anchored at the coastline ( Figure S2 ) excluding the hypothesis of a synoptic disturbance (like a front) as the cause of the wind anomaly. Moreover, the offshore extension of the coastal anomaly matches reasonably well the theoretical offshore extension of a Kelvin wave (eq.5 of Dorman, 1985 and eq S2 in supplementary material ) that is about 53 km at 31N. Finally, the wind reversal coincides with transient positive anomalies of surface pressure ( Figure S2 ) as would be expected for a Kelvin wave (Dorman, 1985).

The control of the PBL height in an extended region ( Figure 5A ) is certainly unrelated to coastal trapped perturbations, which, as mentioned before, decay exponentially in the cross-shore direction. Indeed, the larger scale environment driving the variability of the CANUS is predominantly dominated by the Azores anticyclone. As shown in Figure 12 , built from 1979-2018 July-August-September ERA5 data, the lowering of the PBL height offshore Cape Ghir occurs when the anticyclone core is less intense, smoother and extends in ridge through the Bay of Biscay towards NE France and the English Channel, as noted by the displacement of 1020 hPa isobar, with an intensification of the geostrophic wind south-west of Iberia, leading to a stronger and wider coastal jet near Cape Ghir, extending almost to Madeira, and to intensified tip jets in both Madeira and the Canaries, as also clearly noted before in Figure 7 . In contrast, higher PBL heights are associated with an anticyclone more intense in its core and more constrained in its zonal extension. These results are very similar to the 40-year mean shown in Miranda et al. (2021), but with a different compositing approach.

Although not so intense as for surface wind, the coastal ocean surface temperature also shows variability, notably the one associated with wind reversals ( Figure 9B ). During the wind reversal event, of 11-13 July the coastal SST attains a relative maximum. Similarly, two other SST maxima in the same region, in the first and second fortnights of September are also coincident with the two other wind reversals events, with about one-two days of delay due to ocean inertia. In short, coastal trapped Kelvin waves are associated with periods of upwelling relaxation, as noted before in the California coast (Beardsley et al., 1987).

The cross-correlations of the SST at Cape Ghir with the surface wind at Cape Ghir, Gran Canaria, Madeira and La Palma ( Figure 13A ) show the co-variability of these two variables at four distinct locations where intense jets are usually observed. The maximum computed correlations, noted with a lag of one or two days, range from -0.37 at Madeira to -0.57 at La Palma. The lowest correlations are observed at Madeira jet, the one more distant from the Africa coast and consequently the less synchronized to the others, as noted before in surface wind time series ( Figure 6 ). These values show that a significant part of the SST variability near the Africa coast is due to the jets variability. We hypothesized that some of the unexplained variability may be due to larger scale ocean and atmosphere dynamics, a research topic that is beyond the scope of the present study due to the constrains imposed by the geographic limits of the computational domain. Nevertheless, the impact on ocean surface of the different wind regimes imposed by the jets is clearly noted in Figure 13B where is shown the SST difference between the wind speed upper tercile days and the lower tercile days at Cape Ghir. Similar figures for Madeira, Gran Canaria and La Palma jets show a similar pattern (not shown). As expected, near the coast a cooler ocean surface is observed when the jets are more intense, with values, as low as, -2°C. It is interesting to note that the cooler ocean surface is still verified near the west limit of the computational domain, a possible sign of the large-scale impact of the jets. The main exception to a cooler ocean surface during windy days occurs in the lee of the islands, due to a lower cloud cover that implies an intensified incident radiation, as verified in the simulated data (not shown).

The top-level ocean circulation ( Figure 14 ) varies on longer timescales than its atmospheric counterpart, but is characterized by an intense mesoscale structure, with counterrotating eddies and offshore filaments, as noted in previous studies (e.g., Barton et al., 2004). Two main features, whose origin strongly depends on the jets, characterize that circulation: the southward Canary current along the African coast, and the eddy train propagating westwards in the subtropics north of 30N. The latter was identified by Sangrà et al. (2009) as the Madeira corridor. The Canary Islands, and to a lesser extent Madeira, mark the origin of eddy trains propagating south-westwards into the tropics. These eddies are a response to land-atmosphere-ocean interactions, which lead to atmosphere and ocean vortex shedding, in the most dynamic regions of the tip jets. It is well known that ocean eddies have an important role in the zonal transport of the properties of the upwelling zone to the oligotrophic open ocean. For instance, Sangrà et al. (2009) estimated that the total primary production related to westward eddies may be as high as the total primary production of the Canary upwelling system. Due to the estimated long life (>3 months) of these mesoscale eddies, we will no further analyse them here, but this is certainly a relevant issue to examine in a longer high-resolution ocean-atmosphere coupled simulation.