Super Storm Desmond: a process-based assessment

In the subsequent sections, while all notation is introduced within the text, we include table 2 in the supplementary information (hereafter SI) for convenience when referring to terms featuring within equations.

We used winter half-year (October–March, 1979–2016) six-hourly data at 1° × 1° spatial resolution from the ERA-Interim (ERA-I) reanalysis dataset (Dee et al 2011) to contextualise the atmospheric circulation during Storm Desmond. ERA-I data have been used previously to identify and characterise ARs (Ramos et al 2016), and found to agree well with other reanalysis products when used for this purpose (Lavers et al 2012). To diagnose moisture transport, specific humidity (q), zonal (u), and meridional (v) wind components were extracted for the domain 100 °W–30 °E, 20 °N-80 °N at pressure levels 1000 to 300 hPa. Data from higher than 300 hPa in the atmosphere are neglected because little moisture is transported above this level (Ramos et al 2016).

We also obtained six-hourly q, u, v, air temperature (T) and vertical velocity on all model levels for a global domain for winter half year, 1979–2016. These data were extracted to trace the history of air masses arriving in the BI during AR events (see section 2.3). To this end, we additionally employed global boundary-layer height (blh) from the ERA-I surface fields, and surface pressure (ps) from the ERA-I model levels' archive. Sea-surface temperatures (SSTs) were similarly extracted for a global domain (at 1° × 1° spatial resolution) from the HadISST dataset (Rayner et al 2003) for the longer period 1870–2016. These data were used to interpret the role of SSTs in modulating the humidity of traced air parcels (sections 2.3 and 4).

To characterise and contextualise the AR associated with Desmond, all ARs affecting the BI during the winter half year were identified using the vertically-integrated vapour transport (IVT):

$$\begin{equation} {\rm{IVT}} = \frac{1}{g}\mathop \int \limits_{1000}^{300} q{({u^2} + {v^2})^{1/2}}{\rm d}p \end{equation} \tag{ 1 }$$

where g is the acceleration due to the Earth's gravity, p denotes air pressure, and the remaining terms retain their meaning from section 2.1.

Equation (1) was applied six-hourly to the ERA-I data. The criteria used to classify ARs were adapted from Lavers et al (2012), who similarly classified ARs in the vicinity of the BI; the specifics are briefly recounted here. First, grid cells with IVT greater than the 528.2 kg m⁻¹ s⁻¹ threshold employed by Lavers et al (2012) were identified as possible components of ARs. If this condition was satisfied along 10.5 °W between 50 and 60 °N (the approximate western fringe of the BI; figure 1(a)), a connected-neighbours algorithm was implemented to join all adjacent cells exceeding the IVT threshold (evaluated using eight cell connectedness). If these touching cells extended to at least 30.5 °W, then the BI was defined as experiencing an AR during the respective time step. However, we only retained persistent ARs with the highest potential impact. These were defined as events lasting at least 24 hours, and whose latitude of maximum IVT along 10.5 °W remained within 4.5° of the initial IVT maximum at this line of longitude. The 4.5° latitude tolerance means that the mid-point of ARs can deviate at most ~500 km from their starting position, implying sustained influence over a given region, which likely translates to greater hydrological impact. The minimum duration extent applied here (24 hours) is greater than employed by (Lavers et al 2012), but was preferred for identifying ARs that were potentially most hydrologically severe. Moreover, preliminary screening of winter 2015/2016 ERA-I data indicated that Desmond's AR satisfied this time criteria.

AR intensity for each event was defined via two metrics: (i) the sum of the maximum IVT along 10.5 °W (50 °N-60 °N) (hereafter 'TOTAL'); and (ii) the peak rolling 24 hour mean of the maximum IVT along the same arc of longitude (hereafter 'MAX24'). The former metric places weight on the total amount of moisture delivered over the lifetime of an AR, whereas the latter focuses on intensity over a 24 hour window. We also estimated the rarity of MAX24 for Desmond by calculating its return period using a Gumbel distribution, whose parameters were determined from our sample using maximum likelihood estimation. We selected the Gumbel distribution because initial screening indicated it had the minimum Akaike information criterion amongst a pool of candidate right-skewed distributions (see the SI for details of this provisional screening). The suitability of the fitted Gumbel distribution was evaluated using a Kolmogorov–Smirnov (KS) test, with p-values estimated using a Monte Carlo simulation (see Clauset et al (2009) for a description of this procedure). It is necessary to determine p-values via simulation because we estimate the Gumbel distribution's parameters using sample data, which has the effect of increasing the probability of a type II error, as critical values of the KS test statistic are biased high (Steinskog et al 2007).

The return period (R_p) for Desmond's IVT was estimated by evaluating:

$$\begin{equation} {R_p} = \frac{1}{{\zeta \left[ {1 - {\rm{CDF}}\left( {{\rm{IVT}}} \right)} \right]}} \end{equation} \tag{ 2 }$$

in which ζ is the mean number of ARs per year, and CDF is the Gumbel cumulative distribution function, given by:

$$\begin{equation} {\rm{CDF}}\left( {{\rm{IVT}}} \right) = {{\rm e}^{ - {{\rm e}^{ - \left( {\frac{{{\rm{IVT}} - \mu }}{\sigma }} \right)}}}} \end{equation} \tag{ 3 }$$

where the location and scale are denoted by µ and σ, respectively. Physically, the mean of the Gumbel distribution is a linearly increasing function of both µ and σ, whereas the variance is proportional to the square of σ. We fitted the Gumbel distribution over the 1981–2010 climate period, and calculated the uncertainty on R_p using the profile likelihood method to estimate the 95% confidence interval (Coles 2001).

To provide process insight into Desmond's AR, we separately computed the column-mean wind speed (WS), evaluated at the latitude of the maximum IVT, for the four time steps constituting MAX24:

$$\begin{equation} {\rm{WS}} = \frac{g}{{1000 - 300}}\mathop \int \limits_{1000}^{300} {({u^2} + {v^2})^{1/2}}{\rm d}p. \end{equation} \tag{ 4 }$$

We also computed the total column water (TCW) for the same period and latitude:

$$\begin{equation} {\rm{TCW}} = \frac{1}{g}\mathop \int \limits_{1000}^{300} q{\rm d}p \end{equation} \tag{ 5 }$$

which yields kg of water vapour per square meter and is equivalent to mm. The average of the four WS and TCW values were then used to characterise each AR.

Backwards air mass tracking enables identification of sources and sinks of moisture that influence atmospheric humidity and hence the intensity of ARs (Gimeno et al 2012, Nieto et al 2007, Ramos et al 2016, Winschall et al 2014). To this end, we used the Lagrangian analysis tool 'LAGRANTO' (Sprenger and Wernli 2015) to calculate 10 day back trajectories of air in BI ARs. We then applied the moisture accounting routine of (Sodemann et al 2008) to identify source locations of moisture entrained in the traced air parcels. Finally, the role of SST variability in modulating AR's TCW was explored by calculating uptake-weighted SSTs from the HadISST dataset at sites of moisture uptake. Further details of how these analyses were implemented, including the calculation of weighted SSTs, are provided in SI.

Note that wherever we report the significance of temporal trends, their magnitudes were computed using the Seil–Then slope estimation method (which is robust to outliers; Sen 1960, Theil 1950), and p-values were evaluated using a 1000-trial block bootstrap simulation to account for possible autocorrelation. For each trial, the Seil–Then slope was calculated for the series of randomly reordered blocks and p was defined as the fraction of times the absolute magnitude of these slopes exceeded the magnitude of the trend in the original series. Block lengths were determined for this procedure using the implicit equation given in Wilks (2011).

We detected 81 ARs during the 37 winter half years (figure 2). Desmond's AR is shown in figure 3 for each 6 hourly slice of passage (04/12/2015 12:00 to 06/12/2015 00:00). Maximum IVT along 10.5 °W was attained on 5 December at 12:00, and the 24 hour period of peak IVT intensity (corresponding to MAX24) occurred from midnight on 04/12/2015 to midnight on 05/12/2015.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The other 80 ARs provide context to Desmond. According to TOTAL, Desmond's AR was ranked third since 1979, surpassed only by events in 2011 and 1991. However, by MAX24, Desmond's AR was the most severe on record (table 1). The 1991 event is the only AR other than Desmond to feature in the top three of both metrics. This AR caused heavy orographic rainfall in the southern Pennines (Northern England), leading to river discharges with return periods estimated in excess of 100 years (Institute of Hydrology 1992). The algorithm also detected the well-known AR in November 2009 (ranked ninth according to both metrics), an event that was also associated with flooding in both Ireland and the UK (Adger et al 2013, Lavers et al 2011).

The unusualness of the 24 hour IVT associated with Desmond was explored further by estimating the return period of MAX24 according to a Gumbel distribution. Figure (4) demonstrates the appropriateness of the Gumbel distribution, with points clustered around the 1:1 line. Confidence in the appropriateness of this distribution is increased by the KS test results, which gave little support (p = 0.67) for rejecting the Gumbel as our sample's parent population. Application of equation (2) suggested a return period of 12–66 years (95% confidence interval), with a most likely estimate of 26 years.

Zoom In Zoom Out Reset image size

The representativeness of this return period rests on the MAX24 random variables being independent and identically distributed (iid). Independence seems likely, given that the minimum separation between ARs was 78 hours. It is more challenging to verify the identically distributed assumption. Global mean air temperatures have risen over the period of our AR sample (Hartmann et al 2013), and models indicate increasing AR intensity with climate warming (Lavers et al 2013, Ramos et al 2016), implying that the identically distributed assumption may be inappropriate. However, we found no significant temporal trends in means (p = 0.36) or standard deviations (p = 0.13) computed on five-year running samples of MAX24 (1979–2016). Thus, the iid assumption seems to be a reasonable approximation, and we interpret the return period as representative for the climate of the recent past.

Further insight into the unusualness of Desmond's MAX24 relative to the rest of the ARs was explored by assessing WS and TCW separately over the 24 hours of maximum intensity (table 2). Neither of these quantities was unprecedented, but both approached the upper end of their distributions, with empirical non-exceedance probabilities of 88.9% and 92.6% for TCW and WS respectively (figure 5). TCW for Desmond is perhaps more unusual when the timing is considered, as only two of the top ten events by TCW occured outside of October and November, which is physically consistent with the higher SSTs that prevail during these autumn months. By contrast, Desmond was less seasonally unusual in terms of WS, as all but one of the top ten events by this metric occurred in the core winter months (December to February), with six falling in December. The result of concurrent extremes in these quantities means that Desmond's AR is ranked top if the products of the 24 hour average WS and TCW are assessed (not shown). Whilst this result is consistent with the analysis of MAX24, it should not be regarded self-evident, as the latter depends on covariations between q and u/v (equation (1)), which are not captured by equations (4) and (5).

Zoom In Zoom Out Reset image size

Back trajectories from the region of maximum IVT are shown in figure 6. Across all ARs, incident air parcels evidently spend much of their time to the west of the BI. The majority of the air parcels' time is, however, spent over the North Atlantic and Desmond was no exception in this regard (figure 6(b)). For all tracks, we could attribute an average of 50% of moisture to source regions (see the SI for context regarding 'unattributed' moisture), somewhat lower than that reported by (Sodemann et al 2008). We explore reasons for this in section 4. Most moisture uptake occurs over the North Atlantic (figure 6(c)), with the location of the maxima following the path of the Gulf Stream extension near 40 °N.

Zoom In Zoom Out Reset image size

To explore SST influence on TCW we computed the mean SST at uptake locations for each AR during the period of MAX24. Figure 7 indicates a highly significant correlation (r = 0.62; p < 0.001) between mean SST and corresponding mean TCW, with a sensitivity (β) of 1.16± 0.17 mm °C⁻¹. We calculate the lowest mean SST experienced by air entrained in all ARs to be 13.1 °C, and the highest to be 22.6 °C, translating to a difference in TCW of around 11 mm (43% of the mean TCW across all ARs) using this regression equation. Air parcels entrained in Desmond's AR encountered a mean SST of 19.7 °C, higher than 75% of all other ARs assessed. Evaluating the regression function at this x-coordinate (assuming normally distributed residuals whose standard deviation is 3.15 mm; see figure 7 annotation), yields a conditional exceedance probability of TCW as high as experienced during Desmond (30 mm; hereafter DTCW) of 18%, which is 55% greater than the marginal probability (of 11.1%; given by 100%−88.9%; see figure 5).

Zoom In Zoom Out Reset image size

We have, therefore, demonstrated empirically that SSTs experienced by incident air parcels exert control on the TCW. This raises an interesting question about the possible role of climate change in increasing the likelihood of attaining DTCW. Variability in SSTs experienced by air entrained in ARs may result from (1) the positions of the boundary-layer tracks themselves (for example, tracks spending more time over the lower latitude North Atlantic are more likely to experience warmer SSTs); and (2) temporal variability in SSTs, resulting from internal variability at a range of time scales, and also from externally-forced climate warming. Under constant (1), long-term increases in SSTs would be expected to increase the likelihood of DTCW (see figure 7).

To explore the role of climate warming in contributing to Desmond's high TCW, we implemented an experiment that isolated (2). We took the observed air parcel trajectories (and uptake locations/amounts) from Storm Desmond and iteratively replaced the observed 2015 SST field with SSTs from all years in the HadISST dataset (1870–2015). The mean SST experienced by the air parcels for each of these alternative scenarios was then calculated from SI equation (3), with the resulting SSTs used to compute conditional exceedance probabilities (pr₂) for DTCW, using the regression equation (figure 7). These probabilities were compared to the reference value (pr₁) obtained for the 2015 SST field (18%—see above). An upward trend in the ratio pr₂/pr₁ (hereafter 'PR'—the probability ratio; Fischer and Knutti 2015) would result from increasing SSTs along Desmond's boundary layer air parcel trajectories. This experiment is warranted because the result is far from a forgone conclusion: global mean December SSTs have risen by almost 0.4 °C century⁻¹, increasing to more than 0.5 °C century⁻¹ if the trend is computed from the early 20th century (figure 8(a)). However, warming has been much suppressed over a large area of the North Atlantic, with cooling even apparent in the region of the sub-polar gyre (figures 8(a) and (b)). Accordingly, around 32% of Desmond's moisture uptake occurred from regions with SSTs below their 1951–1980 mean during December 2015 (figures 8 (c) and (d)). Our experiment therefore permits insight into how much warmer the track-averaged SSTs could have been in other years with less pronounced cool SST anomalies. We adopt the 1951–1980 climate normal period because it has been used widely as a reference, and the Atlantic Multidecadal Oscillation (AMO) (Kushnir 1994) was approximately neutral during this time.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

This experiment indicates that uptake-weighted mean SSTs along Desmond's track have risen since 1870, with a significant upward trend of 0.37 °C ± 0.12 °C century⁻¹ (figure 9). This rate is approximately equal to the global mean (figure 8 (a)) and slightly less than the trend in mean North Atlantic SSTs (0.40 °C ± 0.04 °C century⁻¹: calculated over the region −80 °E–5 °E, 0 °N–70 °N using weights based on grid cell surface area, accounting for unequal cell sizes). The increasing SSTs along Desmond's track translate into a significant rise of 17% century⁻¹ in PR (~25% since 1870). However, despite the increasing trend, pronounced high- and low-frequency variability in along-track SSTs means that, for identical air parcel trajectories, higher SSTs (and inferred TCW) should have been expected for Desmond's AR in 15 other years since 1870, and the rolling 30 year mean SST was actually at a peak during the middle of the 20th century (centred on 1947).