The signature of the tropospheric gravity wave background in observed mesoscale motion

How convection couples to mesoscale vertical motion and what determines these motions is poorly understood. This study diagnoses profiles of area-averaged mesoscale divergence from measurements of horizontal winds collected by an extensive upper-air sounding network of a recent campaign over the western tropical North Atlantic, the Elucidating the Role of Clouds-Circulation Coupling in Climate (EUREC⁴A) campaign. Observed area-averaged divergence amplitudes scale approximately inversely with area-equivalent radius. This functional dependence is also confirmed in reanalysis data and a global, freely evolving simulation run at 2.5 km horizontal resolution. Based on the numerical data it is demonstrated that the energy spectra of inertia gravity waves can explain the scaling of divergence amplitudes with area. At individual times, however, few waves can dominate the region. Nearly monochromatic tropospheric waves are diagnosed in the soundings by means of an optimized hodograph analysis. For one day, results suggest that an individual wave directly modulated the satellite-observed cloud pattern. However, because such immediate wave impacts are rare, the systematic modulation of vertical motion due to inertia–gravity waves may be more relevant as a convection-modulating factor. The analytic relationship between energy spectra and divergence amplitudes proposed in this article, if confirmed by future studies, could be used to design better external forcing methods for regional models.

1 Introduction

The effects of large-spatial-scale and long-timescale processes on convection are well understood, such as shallow convection prevailing in geographic locations of free-tropospheric subsidence associated with the Hadley cell (e.g., Albrecht et al., 1995; Norris, 1998; Wood and Hartmann, 2006). These links between convective regimes and circulation averaged over seasonal timescales are often the basis for prescribing an external forcing in limited-area models (Albrecht et al., 1979; Sobel and Bretherton, 2000; van der Dussen et al., 2016). Cloudiness in conventional climate models is highly sensitive to the seasonal and synoptic-scale variations in the atmospheric state, but clouds do not vary adequately on timescales shorter than a day (e.g., Nuijens et al., 2015). Whereas models struggle to capture sub-daily cloud variability, which is pronounced on scales of a few tens to several hundreds of kilometers, several observational studies have established direct links between sub-daily variability of cloudiness and mesoscale – i.e., 𝒪(100) km – vertical motion. For instance, Vogel et al. (2020) found that the shallow convective mass flux is regulated primarily by mesoscale vertical motion at cloud base. George et al. (2020) showed that lower-tropospheric mesoscale vertical velocity explains a larger fraction of variations in cloud base cloud fraction than local humidity or stability. Mauger and Norris (2010) reported variations in near-surface divergence to be the most important driver of cloud cover in the subtropical eastern North Atlantic on timescales shorter than 12 h, whereas boundary layer stratification was relevant for timescales of 12–48 h. Their study relied on satellite observations as well as on reanalysis data to complement the meteorological fields. The still relatively poor understanding of how clouds couple to circulation on scales of 20–200 km was the main motivation for recent efforts to design observations to measure the vertical profiles of area-averaged divergence on sub-daily timescales (e.g., Bony et al., 2017).

By deriving vertical profiles of divergence from sondes dropped along circular flight paths of 200 km diameter, Bony and Stevens (2019) showed that wave-like structures with vertical wavelengths of several kilometers, autocorrelations of multiple hours, and amplitudes on the order of $\mathrm{2}\times {\mathrm{10}}^{-\mathrm{5}}$ s⁻¹ dominated the divergence structure in the free troposphere over the tropical Atlantic. Stephan et al. (2020a) confirmed such features in soundings from the 2006 Tropical Warm Pool–International Cloud Experiment (TWP-ICE; May et al., 2008), which took place near Darwin, Australia. Furthermore, they tested if the characteristics of divergence variability might be consistent with the gravity wave dispersion relation and concluded that gravity waves may serve to explain the temporal, vertical, and horizontal variability observed in area-averaged horizontal divergence. Their study, however, is limited in two ways. First, they only diagnosed the average properties of gravity waves passing the domain, not how variable the wave characteristics were. Second, the TWP-ICE campaign, as well as the circle flights of Bony and Stevens (2019), used constant-sized areas with a diameter of about 200 km and thus could not be used to quantify the dependence of divergence characteristics on area size.

In this study we analyze the extensive radiosounding network (Stephan et al., 2021) of the Elucidating the Role of Clouds-Circulation Coupling in Climate (EUREC⁴A) campaign (Stevens et al., 2021). A network composed of a station at Barbados, four ships, and two aircraft collected more than 2700 vertical profiles of winds, temperature, and relative humidity over the course of 33 d. Since the ships were moving, this network sampled area-averaged divergence over different-sized regions. We derive area-averaged divergence using a three-dimensional variational analysis (described in Sect. 2.2). Independently of this, we diagnose the properties of gravity waves by applying a hodograph analysis to the horizontal winds measured by radiosondes (described in Sect. 2.3).

One aim of the present study is to interpret the properties of mesoscale divergence variability in the context of the tropospheric background wave spectrum, which we diagnose from the ERA5 reanalysis (Hersbach et al., 2020) and a global 2.5 km horizontal-resolution simulation initialized on the starting date of the campaign. It is well known that the global kinetic energy spectrum follows a k⁻³ power law at synoptic scales of 1000–4000 km (e.g., Boer and Shepherd, 1983), which is related to slow balanced dynamics, and a $k^{-\mathrm{5}/\mathrm{3}}$ power law at mesoscales (e.g., Nastrom and Gage, 1985), which is mostly associated with fast unbalanced dynamics and divergence. Hence, one may expect a link between the spectral characteristics of the global energy spectrum and the characteristics of regional area-averaged divergence. A related question is to what extent the tropospheric wave background at given times and locations is saturated versus being dominated by nearly monochromatic waves. Several case studies have shown that individual gravity waves can organize convection (e.g., Shige and Satomura, 2001; Stechmann and Majda, 2009; Lane and Zhang, 2011). Therefore, a second aim of this study is to extract the properties of individual gravity waves from the soundings to quantify their variability, and to test if there are any direct wave imprints visible in patterns of convection.

In Sect. 2 we describe the observations, numerical data, and methods of our analyses. The questions of if and how global energy spectra are linked to mesoscale divergence are discussed in Sect. 3.1, whereas Sect. 3.2 focuses on local waves at specific times. Results are summarized and discussed in Sect. 4.

2 Data and methods

2.1 Observational data

To discuss potential links between vertical profiles of mesoscale divergence and convection, we show images from the Geostationary Operational Environmental Satellite-16 (GOES-16), which is centered on 75.2^∘ W. Infrared (10.3 µm, Band 13) images of the GOES-16 Advanced Baseline Imager (ABI; Schmit et al., 2017) are provided by NASA's Worldview application (https://worldview.earthdata.nasa.gov, last access: 24 November 2020). Information on the cloud bottom and top levels derived from the GOES-16 ABI (NOAA, 2018) is distributed by the Langley SATCORPS group of NASA.

Diagnostics of mesoscale divergence and gravity wave characteristics are based on measurements of surface-launched radiosondes, and dropsondes released from aircraft, taken during the EUREC⁴A campaign in January and February 2020. The experimental design of EUREC⁴A involved 85 dropsonde circles from aircraft flights and regular releases of surface-launched radiosondes from five platforms. A combined data set with 895 dropsonde soundings from the German Aerospace Center's HALO aircraft and 322 dropsonde soundings from NOAA's WP-3D will be made available as part of the Joint dropsonde-Observations of the Atmosphere in tropical North atlaNtic large-scale Environments (JOANNE) data set described elsewhere. We use the quality-controlled Level 2 data with a vertical resolution of 10 m. HALO's circle typically consisted of 12 sondes and is centered at $\sim \mathrm{13.3}{}^{\circ }$ N, −57.7^∘ E, with a diameter of 200 km, as shown in Fig. 1b. The WP-3D circles, of which there are only 14, were mostly flown further east at varying locations within 11–16^∘ N, 59–50^∘ W. HALO data cover the surface up to the flight level of 9 km; WP-3D data cover the surface up to lower and variable flight levels with a maximum altitude of 7 km. We ingest all available dropsonde and surface-launched soundings into a three-dimensional variational analysis (3D-Var; described in Sect. 2.2) to estimate the area-averaged divergence. A separate 3D-Var analysis is performed just on the HALO circles to estimate divergence on this smaller scale.

Figure 1 (a) Equivalent radius of the areas sampled by the sounding network (blue) and the HALO circles (red) over the period of the campaign. All points fulfill $L/R<\mathrm{2}$ (see the text for details). (b) For 31 January 2020 at 16:00 UTC an example of a polygon spanned by radiosoundings (blue) and the location of dropsondes from the HALO circle (red).

In addition to using them in 3D-Var, we analyze Level 2 surface-launched soundings for gravity waves. These data, which cover the ground to the stratosphere (typically 25 km) in steps of 10 m, are described in Stephan et al. (2021) and are publicly available at https://doi.org/10.25326/62. Unlike the circle flights, the continuous data coverage allows us to filter them in time to a specific frequency interval, as we describe in Sect. 2.3. These radiosondes were launched regularly (usually 4 hourly) from a network consisting of the Barbados Cloud Observatory (BCO) and four ships within 6–16^∘ N, 60–51^∘ W: the R/V L'Atalante, R/V MS-Merian, R/V Meteor, and R/V RH-Brown. From 8 January to 19 February, the network launched a total of 848 radiosondes from the surface and also recorded the descent for 82 % of them. Since a typical ascent took about 90 min, this resulted in high-frequency coverage during normal operations. The Meteor remained nearly stationary at a longitude of 57^∘ W and moved north- and southward between 12.0–14.5^∘ N. The RH-Brown moved approximately orthogonal to Meteor's sampling line between the BCO (13.16^∘ N, 59.43^∘ W) and the Northwest Tropical Atlantic Station for air–sea flux measurements buoy (NTAS) at 14.82^∘ N, 51.02^∘ W. The MS-Merian and Atalante ventured further to the south to a minimum latitude of $\sim \mathrm{6.5}{}^{\circ }$ N and were often in close proximity, i.e., <200 km apart. The ships' routes are displayed in Fig. 2 of Stephan et al. (2021).

2.2 3D variational analysis

The computation of divergence and vorticity from unfiltered sounding data is performed with an analysis tool called 3D-Var, described in López Carrillo and Raymond (2011). It ingests spatially nonuniform profiles of temperature, horizontal wind, mixing ratio, and pressure and interpolates them onto a regular grid of the user's choice. Following a weighted least-squares algorithm, 3D-Var then minimizes the misfit between the computed fields and the gridded input measurements. This minimization problem is additionally constrained by the anelastic mass continuity equation. We use a $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ longitude–latitude grid for analyzing data from the full network of soundings and a $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ longitude–latitude grid for our separate analysis of data from HALO's circle flights. Since we are interested in variability on vertical scales >1 km, a vertical grid spacing of 200 m is chosen, covering the surface up to 16 km.

3D-Var is applied every hour between 18 January and 16 February (720 h). From the analysis of Bony and Stevens (2019) divergence autocorrelations for the HALO circles are expected to exceed 3 h and to increase with the area over which divergence is averaged. Therefore, to improve data density, we ingest, for each given time, data within $\pm (\mathrm{\Delta }t=\mathrm{90}$ min). This choice results in nearly continuous 3D-Var results between 18 January and 16 February, as can be seen in Fig. 1a, which shows the equivalent radius, computed as $R=\sqrt{A/\mathit{\pi }}$. The area A, in which the data assimilation is performed, is limited to the region with good data coverage, to avoid extrapolation. For each application of 3D-Var, this area is defined by the boundaries of a polygon that encloses all available soundings and whose corners are grid points with full vertical coverage in all variables. For 31 January, Fig. 1b shows an example of such a polygon. The blue polygon is the mask for the $\mathrm{0.5}{}^{\circ }\times \mathrm{0.5}{}^{\circ }$ all-soundings assimilation, and the red triangles mark the polygon bounds of the $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ HALO-circle assimilation. We discard all times at which the aspect ratio of a polygon's perimeter L and the equivalent radius R (as defined above) is $L/R>\mathrm{2}$, to limit the analysis to areas with little distortion. For the blue polygon in Fig. 1b, for instance, $L/R=\mathrm{1.3}$. The output of 3D-Var contains the 3D fields of divergence and vorticity within the polygon bounds. We average them horizontally to compute vertical profiles of divergence and vorticity.

2.3 Hodograph analysis

The hodographs of the horizontal wind perturbations associated with low-frequency gravity waves describe ellipses whose properties can be determined using the Stokes parameters, as in Eckermann and Vincent (1989). Properties that follow immediately are the ratio of the minor to the major axis, ϵ; the orientation of the major axis relative to the zonal direction, Θ; and the degree of polarization of the ellipse, d. Large d indicates that a sounding contains coherent waves rather than noise. Therefore, we iteratively search for time–altitude windows in the radiosonde measurements that give large values of d, as described in the following.

We analyze data from each platform separately, with one exception: over the studied period, the MS-Merian launched fewer radiosondes (118) than the others ships (Meteor: 203; RH-Brown: 169), which was partly due to her proximity to the Atalante (139 radiosondes). Thus, to achieve continuous data coverage, we include selected soundings collected by the MS-Merian with those of the Atalante; specifically, soundings of the MS-Merian that were taken within 100 km of the Atalante are added to the Atalante's data set. This choice is reasonable as the drift of a sonde can be on the order of 100 km. From now on, we will refer to this merged data set as the combined data set. The BCO is the only spatially fixed platform. They launched 182 sondes, whose horizontal wind measurements from 21 January at 02:00 UTC to 17 February at 12:00 UTC are shown in Fig. 2a and b.

Figure 2 (a) Zonal and (b) meridional wind measured by ascending and descending soundings at the BCO. (c) Zonal and (d) meridional wind perturbations after filtering the raw measurements to $\mathrm{500}\mathrm{m}<{\mathit{\lambda }}_z<\mathrm{10}\mathrm{km}$ and $\mathrm{3}\mathrm{h}<\mathit{\tau }<\mathrm{53}$ h.

Missing values in each data set are filled by interpolating in time to a 1-hourly resolution. A Lanczos filter is applied to horizontal winds and temperature profiles to retain vertical wavelengths of $\mathrm{500}\mathrm{m}<{\mathit{\lambda }}_z<\mathrm{10}$ km and periods of $\mathrm{3}\mathrm{h}<\mathit{\tau }<\mathrm{53}$ h, where 53 h corresponds to the inertial period at the BCO. Figure 2c and d show the resulting horizontal wind perturbations for the BCO. Coherent structures are clearly visible in both wind components. Oftentimes they exhibit a downward propagation with time, which is characteristic of upward-propagating gravity waves.

At each time we compute the quadrature spectrum between the zonal and meridional wind perturbations, u^′(z) and v^′(z). Then we search the resulting spectrum, which is a function of time and vertical wave number m, for local maxima in m that persist in time for at least 10 h. The vertical profiles of u^′(z) and v^′(z) are then filtered to include only the peak vertical wave number and the adjacent wave numbers. A first guess for the associated time window, [$t_{\mathrm{s}}^{\mathrm{0}},t_{\mathrm{e}}^{\mathrm{0}}$], is obtained from the positions of the adjacent minima. At this point, the vertical wavelength assigned to the event is determined from the power-weighted mean wave number of the selected time–wave number window. The initial altitude interval $[h_{\mathrm{s}}^{\mathrm{0}},h_{\mathrm{e}}^{\mathrm{0}}]$ is [0,18] km. The degree of polarization d is next computed for the initial window and all possible combinations of shorter windows contained in it. The window defined by [t_s,t_e], [h_s,h_e], where $t_{\mathrm{s}}^{\mathrm{0}}\le t_{\mathrm{s}}<t_{\mathrm{e}}\le t_{\mathrm{e}}^{\mathrm{0}}$ and $h_{\mathrm{s}}^{\mathrm{0}}\le h_{\mathrm{s}}<h_{\mathrm{e}}\le h_{\mathrm{e}}^{\mathrm{0}}$, with the greatest total number of data points and a polarization d≥0.8 is saved for subsequent analysis. If no window inside [$t_{\mathrm{s}}^{\mathrm{0}},t_{\mathrm{e}}^{\mathrm{0}}$], $[h_{\mathrm{s}}^{\mathrm{0}},h_{\mathrm{e}}^{\mathrm{0}}]$ with d≥0.8 can be found, then the peak is discarded. Otherwise ϵ and Θ are computed using the Stokes parameters.

The angle Θ of the major axis of the ellipse counterclockwise from the eastward direction is the direction of wave propagation, but with an ambiguity of 180^∘. To resolve this ambiguity, we also compute this angle as in Evan and Alexander (2008):

$$\begin{array}{}\text{(1)}& {\mathrm{\Theta }}_{T_{\text{+90}}^{\prime }}={\tan }^{-\mathrm{1}}\left({\displaystyle \frac{\stackrel{\mathrm{‾}}{u^{\prime }{T}_{\text{+90}}^{\prime }}}{\stackrel{\mathrm{‾}}{v^{\prime }{T}_{\text{+90}}^{\prime }}}}\right),\end{array}$$

where $T_{\text{+90}}^{\prime }$ is the value of the temperature perturbation after shifting the phase by $+\mathrm{90}{}^{\circ }$ using a Hilbert transform. Finally, the ratio ϵ of the minor to the major axis is related to the intrinsic frequency Ω of the gravity wave as

$$\begin{array}{}\text{(2)}& \mathit{ϵ}={\displaystyle \frac{f}{\mathrm{\Omega }}}-{\displaystyle \frac{k_{\mathrm{h}}m\stackrel{\mathrm{‾}}{V_z}}{[m^{\mathrm{2}}+\mathrm{1}/(\mathrm{4}{H}_{\mathrm{s}}^{\mathrm{2}}\left)\right]\mathrm{\Omega }}},\end{array}$$

where f is the Coriolis parameter, H_s the density scale height, k_h the magnitude of the horizontal wave number vector, and $\stackrel{\mathrm{‾}}{V_z}=\stackrel{\mathrm{‾}}{\left(\mathrm{d}V/\mathrm{d}z\right)}$ the mean shear of the horizontal wind in the direction perpendicular to Θ. To compute $\stackrel{\mathrm{‾}}{V_z}$, we first average V over [h_s,h_e] and [t_s,t_e].

The dispersion relation for inertia–gravity waves (e.g., Hankinson et al., 2014) is

$$\begin{array}{}\text{(3)}& {\mathrm{\Omega }}^{\mathrm{2}}=f^{\mathrm{2}}+{\displaystyle \frac{N^{\mathrm{2}}{k}_{\mathrm{h}}^{\mathrm{2}}}{m^{\mathrm{2}}+\mathrm{1}/\left(\mathrm{4}{H}_{\mathrm{s}}^{\mathrm{2}}\right)}}-{\displaystyle \frac{\mathrm{2}f{k}_{\mathrm{h}}m{\stackrel{\mathrm{‾}}{V}}_z}{m^{\mathrm{2}}+\mathrm{1}/\left(\mathrm{4}{H}_{\mathrm{s}}^{\mathrm{2}}\right)}},\end{array}$$

where N is the Brunt–Väisälä frequency. Combining Eq. (2) with Eq. (3) gives

$$\begin{array}{}\text{(4)}& {\mathrm{\Omega }}^{\mathrm{2}}\left({\mathit{ϵ}}^{\mathrm{2}}-a^{\mathrm{2}}\right)-\mathrm{2}\mathrm{\Omega }f\mathit{ϵ}\left(\mathrm{1}-a^{\mathrm{2}}\right)+f^{\mathrm{2}}\left(\mathrm{1}-a^{\mathrm{2}}\right)=\mathrm{0},\end{array}$$

where $a^{\mathrm{2}}=\stackrel{\mathrm{‾}}{V_z^{\mathrm{2}}}/\left(N^{\mathrm{2}}\right[\mathrm{1}+\left(\mathrm{4}{m}^{\mathrm{2}}{H}_{\mathrm{s}}^{\mathrm{2}}{)}^{-\mathrm{1}}\right])$. Solving for $\frac{\mathrm{\Omega }}{f}$, this results in

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\mathrm{\Omega }}{f}}={\displaystyle \frac{\mathit{ϵ}\left(\mathrm{1}-a^{\mathrm{2}}\right)\pm {\left[{\mathit{ϵ}}^{\mathrm{2}}{\left(\mathrm{1}-a^{\mathrm{2}}\right)}^{\mathrm{2}}-\left({\mathit{ϵ}}^{\mathrm{2}}-a^{\mathrm{2}}\right)\left(\mathrm{1}-a^{\mathrm{2}}\right)\right]}^{\mathrm{1}/\mathrm{2}}}{{\mathit{ϵ}}^{\mathrm{2}}-a^{\mathrm{2}}}}.\end{array}$$

After determining Ω, k_h follows from the rearranged dispersion relation (3):

$$\begin{array}{ll}\text{(6)}& {\displaystyle }{k}_{\mathrm{h}}=& {\displaystyle }{k}_{\text{disp}}={\displaystyle \frac{fm\stackrel{\mathrm{‾}}{V_z}}{N^{\mathrm{2}}}}{\displaystyle }& {\displaystyle }\pm {\displaystyle \frac{\mathrm{1}}{N}}{\left[{\displaystyle \frac{f^{\mathrm{2}}{m}^{\mathrm{2}}}{N^{\mathrm{2}}}}{\stackrel{\mathrm{‾}}{V_z}}^{\mathrm{2}}+\left({\mathrm{\Omega }}^{\mathrm{2}}-f^{\mathrm{2}}\right)\left(m^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}{H}_{\mathrm{s}}^{\mathrm{2}}}}\right)\right]}^{\mathrm{1}/\mathrm{2}},\end{array}$$

where by definition k_h is positive in the direction of wave propagation.

To obtain the ground-based frequency ω, we compute at each altitude the quadrature spectrum between u^′(t) and v^′(t), QS(t,ν), using the S transformation (Stockwell et al., 1996):

$$\begin{array}{}\text{(7)}& QS(z;t,\mathit{\nu })=\mathrm{\Im }{\left.\left[S{T}_u(t,\mathit{\nu })\cdot S{T}_v^{*}(t,\mathit{\nu })\right]\right|}_z,\end{array}$$

where ST_u(t,ν) is the S transformation in the time dimension of u^′(t) at a given z, and ν denotes frequency. Then, ω corresponds to the peak in $QS(z;t,\mathit{\nu })$ after integrating over [h_s,h_e] and [t_s,t_e].

In principle we could estimate ω from the background wind in the direction of wave propagation, U_h, using $\mathit{\omega }=\mathrm{\Omega }+k_{\mathrm{h}}{U}_{\mathrm{h}}$. However, due to the strong wind shear (see Fig. 2a and b), these results would be highly sensitive to the altitude range selected to compute U_h. As the Doppler shift equation assumes U_h(z)=const, it is not clear if U_h is best approximated as an average across an altitude range or the wind at a particular level. Moreover, as U_h is the projection of the horizontal wind along the wave propagation direction, an additional uncertainty is introduced by the uncertainty in Θ. For this reason we compute ω with the alternative method described above.

2.4 Normal-mode decomposition of numerical data

We derive global wave spectra as well as vertical profiles of regionally averaged divergence and vorticity from the ERA5 reanalysis and a numerical simulation. For both data sets we analyze 3-hourly horizontal winds between 20 January–29 February 2020, matching the EUREC⁴A period. The hodograph analysis is not applied to the numerical data, because we found the vertical resolution to be too coarse.

ERA5 is produced and made publicly available by the European Centre for Medium Range Weather Forecasts (ECMWF). The original data with a horizontal resolution of ∼30 km are stored on 137 hybrid sigma/pressure levels from the surface up to an altitude of 80 km (0.01 hPa).

The numerical simulation run at 2.5 km horizontal resolution is performed with the ICON model (Zängl et al., 2015). It was initialized at 00:00 UTC on 20 January 2020 from an operational analysis of the atmospheric state produced by the ECMWF and forced at the lower boundary with 7 d running-mean sea surface temperatures, which are also taken from an ECMWF analysis. Ninety vertical levels extend from the surface up to 75 km with a Klemp-type implicit Rayleigh damping (Klemp et al., 2008) of vertical winds starting at 44 km. The simulation is freely evolving and integrated until 00:00 UTC on 1 March 2020.

We re-grid the data from both ERA5 and ICON to a regular n256 Gaussian grid with 1024×512 points in longitude and latitude, respectively, corresponding to a resolution of 39 km at the Equator. In the vertical we interpolate to 68 hybrid sigma/pressure levels between the surface and 10 hPa (∼31 km).

The data at each time step are then subjected to a normal-mode function decomposition, using the freely available software package MODES, a detailed description of which is given in Žagar et al. (2015). The modal decomposition projects the three-dimensional fields of geopotential height and horizontal wind onto an orthogonal set of basis functions. These functions are solutions to the Rossby and inertia–gravity wave (IG) dispersion relationships, respectively (Kasahara, 1976). Since the former are vorticity-dominated, we label them as ROT. The orthogonality of the basis functions allows filtering out specific wave modes, including the option of a subsequent inversion back to physical space. This allows us to derive the local wind fields inside the EUREC⁴A domain, as well as corresponding divergence and vorticity fields, associated with specific modes.

3 Results

3.1 Effects of the global wave spectrum on mesoscale divergence

Figure 3a and b show how vertical profiles of area-averaged divergence vary in time. The strong vertical variability in the profiles derived from the array of all soundings and the circle, respectively, is consistent with previous observations by Bony and Stevens (2019) and Stephan et al. (2020a). Vorticity varies on longer scales in time–altitude space (Fig. 3c and d). For ERA5 and ICON we display the corresponding divergence and vorticity profiles associated with IG and ROT modes, respectively (Fig. 4). These profiles are computed for the same geographical region and an equivalent radius of ∼260 km. This radius corresponds to the smallest radius that allows us to compute a centered area average, i.e., using odd numbers of grid points in each direction on the n256 Gaussian grid, while at the same time exceeding the truncation scales applied in the normal-mode decomposition in both horizontal directions (truncation meridional wave number n=200; zonal wave number k=320). Divergence associated with ROT modes is negligible; ROT modes possess only 2.4 % (ERA5) and 3.3 % (ICON) of the variance contained in the IG modes of Fig. 4a and b. Similarly, IG modes possess only 3.8 % (ERA5) and 3.1 % (ICON) of the variance contained in the ROT modes of Fig. 4c and d. Figure 4 also shows strong variability at short scales in divergence and at long scales in vorticity, with similar magnitudes to those observed. Notably, 3D-Var, which ingests a 3 h running mean of measurements, and MODES, which is independently applied at single time steps, yield divergence as well as vorticity anomalies that exhibit temporal persistence with vertically coherent structures. From the comparison between observations and numerical data we conclude that by far the largest fraction of divergence variability inside the sounding network is due to IG waves.

Figure 3 (a, b) Divergence and (c, d) vorticity anomalies derived from 3D-Var applied to (a, c) the sounding network and (b, d) the HALO circles.

Figure 4 (a, b) Divergence and (c, d) vorticity anomalies associated with (a, b) IG and (c, d) ROT modes in (a, c) ERA5 and (b, d) ICON for an equivalent radius of 262 km. The area is centered on ∼12^∘ N, −57^∘ E, matching the approximate location of the EUREC⁴A sounding network (compare Fig. 1b).

Because the EUREC⁴A sounding network sampled areas of different sizes, we can for the first time address how observed magnitudes of area-averaged divergence variability change with area. Figure 5 shows scatterplots of divergence amplitudes multiplied by the corresponding area-equivalent radius versus the equivalent radius for three different tropospheric layers. We here define divergence amplitude as the difference between the maximum and minimum value found in a layer. The running means (thick blue dots) follow lines that are approximately horizontal, albeit with a slight upward slope at h=6–9 km. Thus, divergence amplitudes at equivalent radii of 100–300 km scale approximately inversely with radius.

Figure 5 Divergence amplitude anomaly (defined as maximum minus minimum at (a) 6–9 km, (b) 3–6 km, and (c) 1–3 km altitude) multiplied by equivalent radius, as a function of equivalent radius. Based on divergence from 3D-Var applied to the sounding network. Thick (thin) blue dots show the running mean (±1 standard deviation), using R intervals of 25 km.

ERA5 and ICON also show such an inverse scaling of divergence amplitudes (Fig. 6). As expected from the predominant IG contribution to divergence, this scaling is due to the scaling of the IG-driven divergence field. While the lower bound in R is limited by our choice of normal-mode expansion parameters, we can compute divergence for areas that exceed the sizes sampled during EUREC⁴A. Figure 6 shows that the approximate R⁻¹ scaling of divergence amplitudes remains valid across the entire mesoscale, up to at least the synoptic scale of about 1000 km (zonal wave number 38). ERA5 and ICON produce qualitatively similar results, but we notice in Fig. 6 that ERA5 has slightly larger divergence amplitudes than ICON. To shed light on this difference and, more importantly, to find an explanation for the scaling with radius, Fig. 7 examines the global energy spectra in ERA5 and ICON and their decomposition into IG and ROT modes.

Figure 6 Divergence amplitude anomaly (defined as maximum minus minimum at (a) 6–9 km, (b) 3–6 km, and (c) 1–3 km altitude) in (a, c, d) ERA5 and (b, d, f) ICON. The vertical lines mark ±1 standard deviation. The solid magenta lines in panels (b, d, f) are theoretical predictions based on the spectra shown in Fig. 7 (see the text for details). The dashed magenta line uses a spectral slope of $k^{-\mathrm{5}/\mathrm{3}}$ instead of k⁻² at k>100.

Figure 7 (a, b) Energy spectra as a function of zonal wave number from MODES applied to (a) ERA5 and (b) ICON. The (black) total energy spectrum is decomposed into (red) ROT modes and (blue) IG modes. Dashed black lines show spectral slopes of k⁻¹, $k^{-\mathrm{5}/\mathrm{3}}$, and k⁻³; magenta lines in panel (b) show $k^{-\mathrm{1}/\mathrm{2}}$, $k^{-\mathrm{5}/\mathrm{3}}$, and k⁻², where the transition between $k^{-\mathrm{5}/\mathrm{3}}$ and k⁻² is at k=100 and marked by a thin vertical magenta line.

In both data sets the total energy follows a k⁻¹ power law at synoptic scales that transitions to a steeper slope at k≈8. For k>8 the ROT spectrum in ERA5 scales as k⁻³, and the IG spectrum follows on average a $k^{-\mathrm{5}/\mathrm{3}}$ slope but with substantial deviations from a constant slope. Towards the largest k there is a further steepening in ERA5, which was also reported for ERA-Interim and points to a lack of energy at small scales in the reanalysis (Žagar et al., 2015). The spectral behavior in ICON is similar to ERA5, but the slopes for different regimes are more sharply defined. The IG curve follows a $k^{-\mathrm{1}/\mathrm{2}}$ slope for k<8 and then $k^{-\mathrm{5}/\mathrm{3}}$ up to k≈100. At k>100, there is a slight steepening in ICON as well, such that the slope between k≈100–150 is k⁻². The smallest scale shown in Fig. 6 is k=153. The three prominent slopes of $k^{-\mathrm{1}/\mathrm{2}}$, $k^{-\mathrm{5}/\mathrm{3}}$, and k⁻² are marked with magenta lines in Fig. 7b. Figure 8 shows the energy contained in IG modes in ERA5 relative to ICON. For k<100 the energy in ERA5 exceeds the energy in ICON, and vice versa for k>100. Most of the divergence variability considered here can be attributed to the more energetic large scale, i.e., k<100, explaining the larger divergence amplitudes in ERA5.

Figure 8 Fraction of the energy in IG modes in ERA5 relative to ICON as a function of zonal wave number (ERA5 divided by ICON).

As it turns out, the spectral slopes can be used to predict the magenta lines in Fig. 6b; i.e., they can explain how divergence magnitudes vary with equivalent radius. Consider a gravity wave with wave vector $(k,l,m)$, frequency ω, and phase ϕ. The associated wind perturbations can be expressed as

$$\begin{array}{}\text{(8)}& (u,v,w)=(A_x,A_y,A_z)\sin (kx+ly+mz-\mathit{\omega }t+\mathit{\varphi }),\end{array}$$

where the amplitudes A_i are functions of the wave vector. The associated horizontal divergence scales with k and l as follows:

$$\begin{array}{}\text{(9)}& \text{Div}\sim (k{A}_x+l{A}_y)\cos (kx+ly+mz-\mathit{\omega }t+\mathit{\varphi }).\end{array}$$

As we will demonstrate in Sect. 3.2, the IG waves inside the EUREC⁴A domain were mainly propagating zonally, so that we can assume divergence to scale as Div∼kA_x(k). Near-zonal propagation is also expected, because IG waves can become equatorially trapped (Wheeler and Kiladis, 1999), and it is these equatorially trapped east- or westward-propagating modes that we identify with the normal-mode decomposition. How A_x(k) depends on k is contained in the identified spectral slopes of IG modes in Fig. 7b: A_x scales as $k^{-\mathrm{1}/\mathrm{4}}$ for k<8, as $k^{-\mathrm{5}/\mathrm{6}}$ for k=8–100, and as k⁻¹ for k=100–150. To derive the relationship between these spectral slopes and the divergence amplitudes, we follow the ansatz that each averaging scale $R^{*}=\mathrm{2}\mathit{\pi }/k^{*}$ feels the divergence on all scales $R>R^{*}$, whereas divergence at scales $R<R^{*}$ has little or no effect. Under this assumption, the scaling of divergence amplitudes A_D(k^*) can be approximated by the sum over all k smaller or equal to $k^{*}=\mathrm{2}\mathit{\pi }/R^{*}$:

$$\begin{array}{}\text{(10)}& A_D\left(k^{*}\right)\sim \sum _{k=\mathrm{1}}^{k=k^{*}}k\cdot A_x\left(k\right)\sim \sum _{k=\mathrm{1}}^{k=k^{*}}k\cdot k^{\mathit{\sigma }\left(k\right)},\end{array}$$

with $\mathit{\sigma }\left(k\right)\in \mathit{\{}-\mathrm{1}/\mathrm{4},-\mathrm{5}/\mathrm{6},-\mathrm{1}\mathit{\}}$, as indicated above. Computing this sum and multiplying by a constant to match the IG point at R=1044 km in Fig. 6b, d, and f results in the solid magenta lines drawn in these panels. Evidently, at 3–6 and 6–9 km altitude, the match is excellent. At 1–3 km Fig. 6f misses the decrease in divergence magnitude at small R. This might be due to an introduction of additional energy at short horizontal scales by convection, which tends to have tops at 1–2 km during EUREC⁴A. We did not perform any calculation for ERA5, as ERA5 clearly lacks energy at small scales and the ICON simulation is more realistic in this regard. However, as noted above, even in ICON the slope of the IG spectrum steepens slightly at k=100. Since it is not clear if the change in spectral slope from $k^{-\mathrm{5}/\mathrm{6}}$ to k⁻² at k=100 is physical, we repeat the calculation using a slope of $k^{-\mathrm{5}/\mathrm{6}}$ for all k≥8. This results in the dashed magenta lines of Fig. 6b, d, and f. Now, for R<392 km there is no longer a drop in amplitude, but amplitudes remain approximately constant. This behavior better matches the observations (Fig. 5). The validity of Eq. (10) is further discussed in Sect. 4.

3.2 Effects of individual waves on mesoscale divergence

In the previous section we related the characteristics of the global IG spectrum to the statistical properties of variability in area-averaged mesoscale divergence. We can now explain divergence magnitudes in a statistical sense, but the longevity of structures identified in Figs. 3 and 4 suggests that individual waves may dominate the local divergence field at many times. It is important to note that the result obtained in the previous section does not require a saturated wave spectrum. Evidently, the local spatiotemporal predominance of one or a few waves would change the scaling law in the given period and region, but this is consistent with the relatively large standard deviations in Figs. 5 and 6. Therefore, we next investigate if the generic IG waves that appear to govern the statistics of mesoscale divergence can directly modulate convection through their individual impact on mesoscale vertical motion.

To study local wave characteristics, we use the iterative hodograph algorithm, which is designed specifically to detect nearly monochromatic waves in soundings (see Sect. 2.3). The algorithm extracted 34 wave events with polarization d≥0.8 and with an average polarization of 0.84. All events and extracted wave characteristics are listed in Table 1. Most detected waves have vertical wavelengths of 1.2–4 km with a mean of 2.6 km and highly variable horizontal wavelengths with a mean and standard deviation of 𝒪(450) km. Intrinsic periods range from 6–45 h with a mean of 21 h, and ground-based periods from 11–50 h with a mean of 28 h. The average propagation direction, computed with two different methods, is eastward.

Table 1 For all identified wave events, columns list the (1) time interval t₁–t₂ as “start date/start hour–end date/end hour” with hour in UTC, (2) data set S (R: RH-Brown; M: Meteor; B: BCO; C: combined), (3) altitude interval h₁–h₂ used for the computation of Stokes parameters, (4) degree of polarization d, (5) wave propagation direction Θ counterclockwise from the eastward direction computed with the Stokes parameters, (6) wave propagation direction ${\mathrm{\Theta }}_{T_{\mathrm{90}}}$ obtained with Eq. (1), (7) vertical wavelength λ_z, (8) horizontal wavelength λ_h, (9) intrinsic period τ^*, and (10) ground-based period τ obtained from the quadrature spectrum.

The identified wave characteristics do not show any evidence for systematic changes over the period of the campaign. This is consistent with a generally broad but unsaturated tropospheric IG background. If the wave background were saturated, then it would be impossible to identify waves with a high degree of polarization. Instead, at most times, we are able to find predominant waves (Table 1).

A day on which gravity waves may have had an effect on the cloud field is 31 January. Figure 9 shows hourly infrared satellite images of the campaign region on this day. A smooth modulation of the cloud field is clearly visible in the form of alternating southwest-to-northeast-oriented bands of dry and cloudy areas. The crests of the cloudy areas are topped by narrow stripes of bright clouds. One of the larger cloud clusters seen in the images is located inside the HALO circle. The top and bottom height of this cloud are ∼2.3 and ∼1 km, respectively. In the cloud layer the mean zonal wind is −4.3 m s⁻¹ and the mean meridional wind is −0.1 m s⁻¹. Thus, the motion of the cloud cluster across the circle (diameter of 200 km) matches the estimated advection of 100 km in 6.5 h.

Figure 9 GOES-16 images from 31 January 2020, on a regular longitude–latitude grid. The top edge of the images measures 1358 km, the bottom edge 1410 km, and the sides 889 km. Shown are brightness temperatures from the C13 channel. Green outlines show the Caribbean Islands, and red triangles the locations of HALO's dropsondes.

The 31st of January is also the day when the wave algorithm detected the most distinct wave event. Three out of the four data sets from the surface-based stations contain a wave signal with a high degree of polarization (0.86, 0.86, 0.87; Table 1). Notably, the three detected waves agree broadly in terms of the estimated horizontal wavelengths (480, 340, 440 km), which match well the 400 km spacing of the respective moist or dry regions in Fig. 9. In addition, the long ground-based periods (33, 27, 34 h) are consistent with the nearly stationary appearance of the cloud field in Fig. 9.

The propagation direction indicated by the algorithm, however, is more due eastward ($-\mathrm{14}{}^{\circ }$, $-\mathrm{28}{}^{\circ }$, $-\mathrm{11}{}^{\circ }$) than one would expect from the alignment of the cloud pattern. It could be that vertical propagation and refraction of the wave are partly responsible for this difference, because the ships were located further to the south, where a patterning of clouds is not visible (not shown). To obtain an additional estimate of wave parameters that allows for a more direct comparison with Fig. 9, we apply a hodograph analysis to the 66 HALO dropsonde profiles collected on 31 January over the course of multiple circle flights. Their locations are marked by red triangles in Fig. 9. We pre-select the altitude interval of 2–5 km in order to include the upper portion of the cloud layer while avoiding noise near the surface. Since the entire duration of the flight is less than 8 h, we compute wind perturbations by removing the average of all profiles instead of filtering in time. The resulting horizontal wavelength is 421 km with a propagation angle $\mathrm{\Theta }=-\mathrm{43.1}{}^{\circ }$, in agreement with Fig. 9.

The good match of the satellite image with wave properties found independently in data from four stations is a strong indication for a causal link. The 31st of January was the only day that showed a potential wave imprint on the cloud field that was detectable by eye. Thus, we conclude that visible effects of individual waves on convection occur rather seldom, at least during boreal winter in the studied region.

4 Conclusions

This study extended the analysis of recent novel measurements of mesoscale divergence variability (Bony and Stevens, 2019) by diagnosing vertical profiles of mesoscale divergence in the extensive sounding network of the 2020 EUREC⁴A field campaign. For the first time, we could measure how divergence magnitudes depend on the area over which averages are computed, which may be relevant to sub-daily variability of cloudiness (Mauger and Norris, 2010; Vogel et al., 2020; George et al., 2020). Observed divergence magnitudes during EUREC⁴A scale approximately inversely with the area-equivalent radius, a result that is also true for the ERA5 reanalysis and a 2.5 km global numerical simulation. By applying a normal-mode decomposition to the numerical data, we demonstrated that the major fraction of mesoscale divergence variability is due to inertia–gravity (IG) waves, confirming the hypothesis put forward by Stephan et al. (2020a).

Moreover, we derived an analytic relationship between the wave number dependence of global IG energy spectra and the scaling of divergence amplitudes with area. The essential ansatz for this Eq. (10) is that IG waves on all scales larger than a considered area contribute to the divergence averaged over this area. The data presented in this study follow Eq. (10) very closely. Yet the relationship is not trivial from an intellectual point of view, as one might expect some cancelation of the local divergence induced by some IG waves with the local convergence associated with other IG waves. Either this cancelation has no effect on the overall scaling of divergence amplitudes (note that Eq. 10 only predicts their scaling with area, not their absolute magnitude) or the cancelation is negligible. The second possibility is more likely to be true if the tropospheric IG wave background is not saturated but locally composed of few waves. Persistent features in filtered-wind and divergence profiles suggest that the latter possibility may apply.

To gain more information about the properties of local waves, we extracted nearly monochromatic wave events from soundings by means of an iterative hodograph analysis. For most days, clear IG signals could be found, and on one day the identified wave likely had a direct effect on the organization of clouds, as evidenced by comparison to satellite images. Thus, although seldom, direct wave imprints on convection do occur. This suggests that the IG wave background is indeed not saturated. However, the systematic modulation of vertical motion due to IG waves of different scales may be more relevant for convection. More research is required to quantify the relative importance of the stochastic wave forcing for shaping patterns of convection, compared to other well-known local influences, such as lower-tropospheric stability (Medeiros and Stevens, 2011), relative humidity (Slingo, 1987), and sea surface temperature (Qu et al., 2015).

Future studies should repeat the measurements analyzed in this study for different regions, especially different latitudes, as well as regions over land, to test the universal applicability of the proposed relationship between energy spectra and divergence amplitudes. If confirmed, then our results imply that all scales are relevant for mesoscale variability, which has consequences for regional models that introduce an artificial truncation scale, and could be used to design better external forcing methods for regional models.