The 19–21 June 2013 Alberta flood was the second costliest ($6 billion CAD) natural disaster in Canadian history, trailing only the 2016 Fort McMurray, Alberta, Canada, wildfires. One of the primary drivers was an extreme rainfall event that resulted in 75–150 mm of precipitation in the foothills west of Calgary, Canada. Here, the mesoscale dynamics and thermodynamics that contributed to the extreme rainfall event are elucidated through high-resolution numerical model simulations. In addition, terrain reduction model sensitivity experiments using Gaussian smoothing techniques quantify the importance of orography in producing the extreme rainfall event. It is suggested that the extreme rainfall event was initially characterized by the formation of a surface cyclone on the eastern side of the Canadian Rockies due to quasigeostrophic (QG) mechanisms. Orographic processes and diabatic heating feedbacks maintained the surface cyclone throughout the event, extending the duration of both easterly upslope flow and QG forcing for ascent in the flood region. The long-duration ascent and associated condensational heat release in the flood region vertically redistributed potential vorticity, anchoring and further extending the duration of the surface cyclone, upslope flow, and the rainfall. Although the magnitudes of ascent and precipitation were smaller in 10% and 25% reduced terrain simulations, only a terrain reduction of greater than 25% drastically altered the location and magnitude of the heaviest precipitation and the associated physical mechanisms.
The June 2013 Alberta flood was the second most expensive natural disaster in Canadian history, resulting in more than $6 billion CAD in losses (Milrad et al. 2015; Pomeroy et al. 2016; Teufel et al. 2017). Only the 2016 Fort McMurray, Alberta, Canada, wildfires were more costly. The flood was caused by multiple factors (Pomeroy et al. 2016; Teufel et al. 2017), including an extreme rainfall event during 19–21 June (Milrad et al. 2015). An unusual synoptic-scale evolution characterized the extreme rainfall event (Milrad et al. 2015), which resulted in a Rex block (Rex 1950) over the northern United States and Canadian Rockies. The cutoff low portion of the Rex block helped to spur surface cyclogenesis on the eastern slopes of the Rockies and associated anomalously strong long-lived easterly upslope flow into the southern Alberta foothills region during the extreme precipitation event (Milrad et al. 2015).
As detailed in Milrad et al. (2015) and Pomeroy et al. (2016), the heaviest precipitation occurred in the foothills west of Calgary (CYYC), near Banff (CWZG) and Canmore, Alberta, Canada. For the remainder of this paper, we refer to the foothills west of Calgary as “the flood region,” outlined with a triangle between the three cities (stars) in Figs. 1 and 2. The long-lived heavy rainfall resulted in increased streamflow downstream toward Calgary, forcing an evacuation of downtown on 21 June (Milrad et al. 2015; Pomeroy et al. 2016; Teufel et al. 2017). As in other infamous flash floods in the Front Range of the Rocky Mountains such as the 2013 Great Colorado Flood (Lavers and Villarini 2013; Gochis et al. 2015; Morales et al. 2015) and the 1976 Big Thompson Flood (Maddox et al. 1978; Caracena et al. 1979), differential cyclonic vorticity advection (CVA), lower-tropospheric warm-air advection (WAA), and easterly upslope flow all played important roles in providing forcing for ascent (Milrad et al. 2015). These ascent forcing mechanisms acted to release instability (Maddox et al. 1979; Schumacher and Johnson 2005, 2006, 2008, 2009; Milrad et al. 2015).
In this study, we present high-resolution numerical model simulations and sensitivity experiments. Our specific objectives are to do the following:
Further elucidate the mesoscale processes associated with the extreme rainfall event, using an ingredients-based methodology for flash flooding (Doswell et al. 1996; Gyakum 2008). Specifically, these ingredients are as follows: lift, moisture, and instability. In addition, we investigate the importance of rainfall duration to the event magnitude.
Quantify the importance of the local orography to ascent forcing and cyclone maintenance, using numerical model sensitivity experiments that employ terrain smoothing techniques.
The remainder of the manuscript is organized as follows: section 2 details the data used and model setup, section 3 summarizes the numerical simulation results, while section 4 presents insights gained from the terrain smoothing model sensitivity experiments. Finally, a discussion of the most pertinent findings and conclusions is presented in section 5.
2. Data and model configuration
We compared our model precipitation to the Environment and Climate Change Canada Canadian Precipitation Analysis (CaPA), which has a 15-km grid spacing and 6-h temporal resolution (Mahfouf et al. 2007). Given that it was our primary dataset in Milrad et al. (2015), we used the National Centers for Environmental Prediction (NCEP) Climate Forecast System Reanalysis (CFSR) to verify that the model accurately reproduced the synoptic-scale structures associated with the Alberta flood. The CFSR is run on a T382 spectral resolution (~38-km grid spacing) and obtained on a 0.5° global grid, with a 6-h temporal frequency (Saha et al. 2010). After 2010, the CFSR has been produced in real-time as the CFSv2 (Saha et al. 2014). Finally, all plots except Fig. 1 were created using the General Meteorological Package (GEMPAK), version 7.2.0, updated from the original package devised by Koch et al. (1983).
b. Terrain sensitivity numerical experiments
We performed a series of terrain smoothing experiments in order to quantify the impact of the Alberta topography on the extreme rainfall event using version 3.6 of the Advanced Research version of the Weather Research and Forecasting (WRF) Model (WRF-ARW; Skamarock et al. 2008). All simulations employed a two-way nested 27/3-km domain with 48 vertical levels centered over CYYC (Fig. 1). We initialized the simulations at 1200 UTC 19 June and ran them for 48 h, with NCEP Global Forecast System (GFS) 0.5° analyses [NOAA National Centers for Environmental Information (NCEI) 2016] used as initial and lateral boundary conditions. In 1- and 4-km grid spacing WRF simulations of the 2013 Great Colorado Flood, Schwartz (2014) also used the GFS 0.5° analyses as initial and lateral boundary conditions. Similar to Flesch and Reuter (2012), we initialized the WRF simulation 12–24 h before the onset of the heaviest precipitation, at 1200 UTC 19 June. All of the precipitation fell within the 48-h simulation period, and approximately 85% fell after 0000 UTC 20 June (Milrad et al. 2015; Pomeroy et al. 2016).
We used the Kain–Fritsch cumulus parameterization scheme (Kain and Fritsch 1990; Kain 2004) for the 27-km domain, while the 3-km domain was convection allowing. We parameterized microphysical processes with the Morrison scheme (Morrison et al. 2009) and the Mellor–Yamada–Nakanishi–Niino Level 2.5 (MYNN2.5; Nakanishi 2001; Nakanishi and Niino 2004, 2006) planetary boundary layer scheme was used in conjunction with the MYNN surface layer scheme and the Unified Noah land surface model (Chen et al. 1996, 1997; Koren et al. 1999; Chen and Dudhia 2001; Ek et al. 2003).
We based the magnitudes of the smoothing on the elevation of the two tallest peaks (2386 and 2395 m) of the Canadian Rockies located to the immediate west of the flood region. Subsequently, we applied a Gaussian smoother to successively reduce the height of these peaks by 10%, 25%, and 40% within both the 27- and 3-km domains. Although the Gaussian smoother can introduce noise at the boundaries of the inner domain, we made the 3-km nest intentionally large to limit the impact of this noise on the flood region, which is located near the center of the domain (Fig. 1).
While a variety of other techniques have been employed to perform terrain reduction sensitivity experiments (e.g., Flesch and Reuter 2012; Schumacher et al. 2015), the application of a Gaussian smoother eliminates the ringing effect within the spatial domain, which can occur when applying techniques such as Fourier transforms. Specifically, ringing effects are the result of the Fourier transforms’ decaying oscillations, which lack high frequencies, unlike the Gaussian smoother that smoothly tapers to zero. Signal processing with a hard cutoff can result in ringing artifacts in the filtered signal. In addition, applying a Gaussian smoother allowed for the examination a spectrum of reduced terrain regimes, providing more in-depth analysis than studies that apply, for example, a uniform 75% terrain height reduction above a certain threshold (Flesch and Reuter 2012). In regions where the reduced terrain height was below the height of the original surface, the numerical grids were extrapolated vertically to the surface by the WRF.
3. Numerical simulation results
a. Large-scale characteristics
Figure 2 compares WRF 3-km simulated precipitation to the CaPA. For a comparison of both the WRF and CaPA precipitation amounts to station gauge data, we refer the reader to Milrad et al. (2015) and Pomeroy et al. (2016). With respect to the CaPA, the WRF 3-km total precipitation is 10–20 mm too small in the flood region near Banff and Canmore and 40–50 mm too large north of the flood region (Figs. 2a,b). Most of the total precipitation differences in the flood region between the WRF and the CaPA occurred on day 2 (1200 UTC 20 June–1200 UTC 21 June), when the WRF precipitation maximum was located too far to the north (Figs. 2g–i).
Although the WRF and CaPA are not an exact match, primarily with respect to the location of the precipitation maximum, we have confidence that the WRF is accurately simulating the dynamic and thermodynamic mechanisms responsible for the event (sections 3 and 4). In previous work on convection and rainfall simulations (e.g., Lombardo and Colle 2013), we compared WRF surface temperature and moisture, vertical profiles, and synoptic and mesoscale dynamics and thermodynamics to observations to ensure that the simulations were acceptable. Though the WRF simulation of the QPF structure is important, it is also crucial for the WRF to simulate these metrics for us to have confidence in the simulations. In the current case, the simulated synoptic- and mesoscale mechanisms and structures were checked extensively with observations and CFSR data, and were found to generally be a good match.
To demonstrate that the WRF accurately reproduces the synoptic-scale characteristics of the extreme rainfall event, Figs. 3 and 4 compare the WRF 27 km to the CFSR for 500-hPa geopotential height and mean sea level pressure, respectively. Figure 3 examines quasigeostrophic (QG) forcing for ascent using the Q-vector form of the inviscid, adiabatic QG omega equation:
where is the constant Coriolis parameter, σ is the static stability parameter, ω is the vertical velocity, and the sense of the forcing for vertical motion is related to the divergence of the Q-vector, Q. Areas of Q-vector convergence are associated with QG forcing for ascent.
The Q-vector form of the QG omega equation is useful in that there is no contradiction between QG forcing terms, which can frequently occur with the full QG omega equation. For example, from the full QG omega equation, differential CVA and the horizontal Laplacian of WAA are both associated with QG forcing for ascent. However, collocated differential CVA and the horizontal Laplacian of cold-air advection (CAA) result in ambiguous QG forcing for vertical motion. We calculated Q-vector divergence in the 600–400-hPa layer such that it was representative of QG forcing for midtropospheric ascent, but entirely above the highest terrain peaks. Other layers (e.g., 500–300 and 700–400 hPa) were tested and found to be similar. In addition, all Q-vector diagnostics shown in this paper were filtered through a 12-point Gaussian-weighted smoothing function in GEMPAK, because the WRF (27 km) and CFSR (0.5°) grid spacing were too fine for unfiltered Q-vector diagnostics to be of reasonable magnitudes (Barnes et al. 1996).
The WRF 27 km and CFSR both show a cutoff 500-hPa cyclone over the northern U.S. Rockies (Fig. 3), which forms the southern portion of the Rex block (Milrad et al. 2015). The WRF 27 km shows that over the flood region, Q-vector convergence is observed throughout the event (Figs. 3b,d,f). The CFSR shows Q-vector convergence in the flood region at 1200 UTC 20 June and 0000 UTC 21 June (Figs. 3c,e), although the magnitudes are smaller than in the WRF 27 km. Particularly at 1200 UTC 20 June (Figs. 3c,d), there is larger Q-vector divergence on the western slopes of the Rockies in the WRF 27 km than in the CFSR, although both clearly indicate QG forcing for descent.
In composite analyses of heavy mesoscale precipitation events, Junker et al. (1999) found that regions of heavy precipitation were collocated with areas of both high equivalent potential temperature (θe) advection and geostrophic WAA. However, high-θe advection was of a larger magnitude, suggesting positive moisture advection in addition to WAA. Figure 4 shows high-θe advection in the flood region throughout the duration of the extreme rainfall event. In addition, the surface cyclone on the eastern slopes of the Rockies (Fig. 4) is not what is typically referred to as lee cyclogenesis. Lee cyclogenesis on the eastern slopes of the Rockies occurs with downslope-oriented westerly flow aloft, while in the Alberta flood the upper-tropospheric flow is easterly (Fig. 3). The cyclone observed on the eastern slopes of the Rockies from 0000 to 1200 UTC 20 June in Montana, Alberta, and Saskatchewan (Figs. 4a–d), Canada, forms in response to differential CVA (Milrad et al. 2015) associated with the cutoff upper-tropospheric cyclone located over Idaho and Montana during the same time period (Figs. 3a–d).
b. Dynamics and thermodynamics
Figure 5 shows that Q-vector convergence is present in southern Alberta from 0000 UTC 20 June (Fig. 5a) to 0000 UTC 21 June (Fig. 5e), with the exception of 1200 UTC 20 June (Fig. 5c), when Q-vector divergence is observed in the flood region. Milrad et al. (2015) found similar results using the CFSR and showed that during 0600–1200 UTC 20 June, only relatively light precipitation fell. Junker et al. (1999) and Maddox et al. (1978) found that in flash flood events, geostrophic WAA is typically a larger contributor to QG forcing for ascent than differential CVA. To that end, geostrophic WAA is present in the flood region throughout the event (Fig. 5), with the largest magnitudes in the foothills at 1800 UTC 20 June and 0000 UTC 21 June (Figs. 5d,e), collocated with Q-vector convergence.
Figure 6 shows high-θe advection in the flood region throughout the majority of the extreme rainfall event. At 1200 UTC 20 June, high-θe advection is coming from the east-northeast, on the north side of the surface cyclone (Fig. 6c), and collocated with the heaviest precipitation (Fig. 2). The high-θe advection maximum observed north of the flood region at 1200 UTC 20 June (Figs. 6c,d) may explain why the WRF 3 km oversimulated precipitation to the north on day 2 (Fig. 2). Buzzi et al. (1998) and Lin et al. (2001) also observed high-θe advection within a narrow upslope low-level jet (LLJ) in case studies of heavy orographic precipitation.
Figure 5e shows that starting at 1200 UTC 20 June, Q-vector convergence is present on the northern side of the surface cyclone in southern Alberta. The cyclone on the eastern slopes is also persistent throughout the entirety of the extreme rainfall event (Figs. 5 and 6). In full-terrain WRF simulations of two Alberta extreme rainfall events in 2005, Flesch and Reuter (2012) similarly found a long-lived surface cyclone on the eastern slopes of the Rockies. In terms of cyclone formation, the traditional form of the QG omega equation (Bluestein 1992, p. 328) dictates that differential CVA is associated with QG forcing for ascent. This in turn results in surface convergence, leading to an increase in cyclonic vorticity via the QG vorticity equation (Bluestein 1992, p. 327), and thus a stronger surface cyclone. The long duration of the cyclone in the Alberta flood contributed to the magnitude of the extreme rainfall event (section 4).
Figure 7 presents southwest–northeast-oriented cross sections along the blue line in Fig. 6a. At the start of the extreme rainfall event (0000 UTC 20 June), the flood region is characterized by a convectively neutral lower troposphere (Fig. 7a). At 0600 UTC 20 June (Fig. 7b), the lower troposphere in the flood region is convectively stable, with descent overhead and ascent located to the east. Convective instability is first observed over the flood region and collocated with ascent at and above 700 hPa at 1200 UTC 20 June (Fig. 7c). After 1200 UTC 20 June, convective instability becomes more pronounced across the flood region for the remainder of the rainfall event (Figs. 7d–f).
4. Role of orography
a. Terrain smoothing and precipitation
Figure 8 shows WRF 3-km elevation maps of the four smoothing thresholds employed for the sensitivity experiments. The spectrum of smoothing thresholds includes unsmoothed (hereafter control, Fig. 8a), 10% smoothed (Fig. 8b), 25% smoothed (Fig. 8c), and 40% smoothed (Fig. 8d). In the flood region, the maximum elevation is reduced from approximately 2400 to 1800 m in the 25% smoothed simulation (Fig. 8c) and 1400 m in the 40% terrain smoothed simulation (Fig. 8d).
Figure 9 shows that the 48-h total precipitation in the 10% smoothed simulation is similar to the control (Figs. 9a,b). In the 25% smoothed simulation, maximum precipitation accumulations in the foothills are reduced by 30–40 mm, while the precipitation maximum on the western slopes of the Rockies in southeastern British Columbia increases slightly in areal extent (Fig. 9c). The most substantial precipitation differences occur when the terrain is smoothed to 40% of its original value; Fig. 9d shows a precipitation maximum in southeastern British Columbia, which corresponds to the elevation maximum (Fig. 8d). A secondary precipitation maximum is seen in Alberta, but far to the northeast of the flood region. The 40% smoothed terrain resembles a gently sloping plateau, which “peaks” to the west of the actual peaks of the Canadian Rockies (Fig. 8d). Consequently, the 40% smoothed simulation exhibits a more horizontally disperse precipitation shield (Fig. 9d), which Buzzi et al. (1998) also found in simulations of heavy orographic precipitation events.
b. Dynamics and thermodynamics
For the sake of brevity and to demonstrate the largest differences in dynamics between smoothing thresholds, we primarily focus on three smoothing thresholds in this section: control, 25% smoothed, and 40% smoothed. Considering 12-h accumulated precipitation totals, the control and 25% smoothed simulations are similar at 0000 UTC and 1200 UTC 20 June (Figs. 10a,b,d,e). By 0000 UTC 21 June, there is incrementally less precipitation over the flood region in the 25% and 40% smoothed simulations (Figs. 10f,i) than in the control (Fig. 10c). In addition, the magnitude of the precipitation maximum to the north of the flood region is also reduced, and in the case of the 40% smoothed simulation, located farther to the northeast of the flood region. These observations are representative of smaller day-2 precipitation totals in the flood region at larger smoothing thresholds relative to the control, and suggest the importance of orographic processes during the latter part of the event. In the 40% smoothed simulation, a cyclone extends zonally across the Alberta–British Columbia border (Figs. 10h,i), likely explaining smaller precipitation totals in the flood region than in the control and 25% smoothed simulations. Meanwhile, at 0000 UTC 21 June (middle of day 2), the control and 25% smoothed simulations feature a surface cyclone on the eastern slopes of the Rockies (Figs. 10b–f), which acts to focus east-northeasterly upslope surface flow, and results in more precipitation in the foothills flood region. These findings further suggest that orographic processes are important to the maintenance and anchoring of the surface cyclone east of the highest terrain.
Figure 11 shows high-θe advection over the flood region for all smoothing thresholds at 0000 UTC 20 June (Figs. 11a,c,e), as the surface cyclone is located over Montana. During day 2, low-θe advection is present in the flood region in the 25% and 40% smoothed simulations (Figs. 11d,f), while neutral-θe advection is present in the control, although high-θe advection is observed immediately to the west. More importantly, in the 40% smoothed simulation the cyclone stretches across the Rockies, indicating that the strongest high-θe advection is in southeastern British Columbia (Fig. 11f), collocated with the precipitation maximum (Fig. 9d).
Mesoscale circulations over the flood region, including frontogenesis and its contributions to the 3D circulation, also contributed to the extreme nature of the event. Petterssen (1936) and Keyser et al. (1988) defined 2D frontogenesis as the rate of change over time of the horizontal potential temperature gradient:
Frontogenesis is calculated from the WRF 27-km control and 40% smoothed simulations (Figs. 12 and 13) only. We do not use the WRF 3 km for this comparison because deconvolving the flow at 3-km grid spacing is untenable, since geostrophic and hydrostatic balance conditions are not met. Figure 12 presents a plane view of the 800–700-hPa layer-averaged 2D frontogenesis, potential temperature, winds, and the most unstable lifted index. The most unstable lifted index is calculated by taking the parcel in the lowest 2 km with the largest value of θe and lifting it to 500 hPa. Two key times in the frontal dynamics, 0600 and 2300 UTC 20 June, are instructive in elucidating the impacts of terrain.
At 0600 UTC 20 June, the control and 40% smoothed simulations (Figs. 12a,c) both show a westward-oriented ridge of higher potential temperatures extending along and just north of the flood region. In addition, both simulations reveal a region of frontogenesis associated with confluence in the flow on the northern periphery of the cyclonic circulation. Although the thermal structure of the two simulations looks similar, the easterly “jet” in the control shows a decrease in speed as well as a small deflection toward lower pressure as the flow approaches the terrain barrier, resulting in enhanced frontogenesis along the eastern slopes of the terrain. Presumably, this decrease in speed is a function of a loss in horizontal kinetic energy as the flow ascends the barrier in a stably stratified atmosphere, and the deflection is consistent with the disruption of geostrophic flow.
While the flow is slightly deflected, it does not result in barrier jet formation or significant cold air damming (e.g., Bell and Bosart 1988). We examined the WRF 3-km control 10-m wind and actual surface observations from the flood region (not shown) to further verify that there was no evidence of a barrier jet. We suggest the lack of a barrier jet is due to air within the upslope flow being too unstable. This is supported by the negative lifted index values located just east and northeast of the flood region (Fig. 12). Although convective available potential energy (CAPE) in the flood region is relatively small (section 4b), air ingested into the system is characterized by conditional instability.
An investigation of a possible barrier jet can be formalized by looking at the ratio of the kinetic energy to the stability, the Froude number, as seen in Eq. (4):
where U is the mean wind speed of the flow, H is the depth of the barrier, and N is the Brunt–Väisälä frequency given by Eq. (5):
where g is the gravitational acceleration and θ is potential temperature. In general, “dammed” flows are characterized by values of Fr < 1 (e.g., Colle and Mass 2000). Here, the mean 10-m wind speed in the upslope jet is approximately 25 m s−1, the barrier height is approximately 1500 m, and the Brunt–Väisälä frequency squared is roughly 15 × 10−5 s−2, which yields a Froude number of approximately 1.4. This Froude number value is marginal for barrier jet formation (e.g., Bell and Bosart 1988).
Later in the event, the differences become more pronounced; at 2300 UTC 20 June, frontogenesis has all but disappeared in the 40% smoothed simulation (Fig. 12d), while an area of frontogenesis is still situated along the eastern slopes of the Rockies in the control (Fig. 12d). This difference in frontogenesis can be explained by both the differences in the flow field as well as the potential temperature structure. In the 40% smoothed simulation (Fig. 12d), the cyclonic circulation has shifted westward and is elongated in an east–west direction relative to the control (Fig. 12b), which shifts the region of confluent and convergent flow to the southwest into British Columbia. In the control (Fig. 12b), the cyclonic circulation is confined to Alberta, as is the confluence zone in the southwest quadrant of the cyclone. Furthermore, the potential temperature gradient remains anchored along the eastern slopes, as the warm air has failed to cross the terrain due to adiabatic cooling of the air in the upslope region and diabatic impacts of precipitation situated along the eastern slopes.
The frontogenesis in the flood region in the control is dominated by the advection of potential temperature by the ageostrophic wind (not shown), suggesting that semigeostrophic theory (e.g., Martin 2006, p. 203) is more suitable than QG theory when considering the implications of frontogenesis and the secondary circulations. Regions of frontogenesis are characterized by thermally direct circulations with air rising on the warm side of the thermal gradient and sinking on the cold side. Beyond frontogenesis, the magnitude and shape of the secondary circulations is dictated by both the static stability as well as the geostrophic shear profile. For a full discussion of semigeostrophic flow and the Sawyer–Eliassen equation, see Martin (2006, p. 204) and Holton and Hakim (2013, p.285).
The cross sections in Fig. 13 show the relationship of the secondary circulations to the frontogenesis. Plotted are 2D frontogenesis, the tangential circulation, , for which dashed lines represent small negative (i.e., weakly stable) values, and the condensational heating rate (CHR). The CHR is calculated as a finite approximation of Eq. (6):
where L is the latent heat of condensation; pb and pt are the pressures at the surface and 100 hPa, respectively; ω is vertical motion; and rs is the saturation mixing ratio. The expression in any given layer is automatically set to zero if ω is positive or if the relative humidity is <90%.
In Fig. 13, the orientation of each cross section is shown by the corresponding thick solid blue line in Fig. 12, and the cross sections are roughly oriented along the thermal gradient in the control. While the circulations should be illustrative of the impacts of the frontogenesis, the circulations themselves are not strictly a balanced response to the frontogenesis and include the full model ascent.
Figures 13a and 13c compare the control and 40% smoothed circulations for 0600 UTC 20 June, when QG forcing for ascent was dominant (Milrad et al. 2015). In both runs, there is strong frontogenesis located near the eastern slopes of the terrain (Figs. 13a,c). The frontogenesis in both cases results in a region of ascent on the warm side of the potential temperature gradient. While there is evidence of the descending branch of the circulation in both runs, it is more pronounced in the control (Fig. 13a). Upslope ascent is seen along the eastern slopes of the terrain (Fig. 13a), in addition to the broadening of the frontogenesis in the lower troposphere. The broadening of the lower-tropospheric frontogenesis subsequently results in another region of ascent to the north of the branch anchored just north of the highest peaks. While the updrafts are slightly stronger in the 40% smoothed simulation, the control has a broader region of ascent, with the strongest updrafts and more intense condensational heating concentrated in the 800–600-hPa layer.
The differences between the two runs are more pronounced at 2300 UTC 20 June (Figs. 13b,d). The frontogenesis in the 40% smoothed simulation (Fig. 13d) has a small vertical extent in the flood region and the resulting vertical circulation is nearly zero, as shown by arrows only slightly deviated from horizontal. Meanwhile, the control (Fig. 13b) shows a region of frontogenesis tilted southwestward with height, along with an ascending vertical branch of the circulation and condensational heating. The control also shows larger static stability over and to the west of the highest terrain, relative to the 40% smoothed simulation. This is likely associated with both adiabatic cooling and the diabatic redistribution of potential temperature associated with the diabatic heating above 825 hPa (Fig. 13b), further discussed below.
An importance difference between the 40% smoothed and lesser thresholds on day 2 is the lower-tropospheric winds. In the control and 25% smoothed simulations (Figs. 14b,d), the winds from the surface to 600 hPa are distinctly more northeasterly and stronger than they are in the 40% smoothed simulation (Fig. 14f). Because the Canadian Rockies are oriented from southeast to northwest (Fig. 1), northeasterly winds are perpendicular to the terrain, suggesting there is more orographic ascent at smaller smoothing thresholds. Finally, we compared the WRF winds in Fig. 14 to the nearest radiosonde soundings at Stony Plain, Alberta, and Glasgow, Montana, and surface observations (not shown). Although both radiosonde stations are more than 200 km away from the flood region, they were enveloped in the easterly upslope flow. We concluded that the WRF wind directions are reasonable, but the wind speeds are likely slightly too strong in the control, particularly near the surface. For example, the surface wind at Stony Plain was 30 kt (1 kt = 0.5144 m s−1), and observed wind at CYYC was 20 kt.
Interpolated values of surface CAPE at Canmore are plotted on each panel in Fig. 14. In the control (Figs. 14a,b), the CAPE values are similar to the CFSR values found by Milrad et al. (2015) and confirm that the extreme precipitation event occurred in a low-to-moderate CAPE environment (e.g., Maddox et al. 1979). Thus, given enough saturation and lift to the level of free convection, buoyant ascent is likely occurring at 0000 UTC 20 June (Fig. 14a), but not at 0000 UTC 21 June (Fig. 14b), where the midtroposphere is quite dry. Overall, the WRF 3-km soundings demonstrate that both moisture and conditional instability were present during the first part of the event, particularly in the lower troposphere. Finally, some saturation is still observed in the control and 25% smoothed simulations at 0000 UTC 21 June.
Figure 15 shows cross sections as in Fig. 7, but for each smoothing threshold, with the respective terrain plotted at the bottom. The 25% smoothed simulation shows smaller values of ascent in the flood region (Figs. 15c,d) compared to the control simulation (Figs. 15a,b), as also found by Galewsky and Sobel (2005) in model sensitivity experiments of heavy orographic precipitation in California. In addition, at 0000 UTC 20 June (Figs. 15a,c), the control simulation has a larger descent maximum on the downslope (western) side of the Rockies, suggestive of the impact of the higher terrain peaks on vertical motion.
Both the control and 25% smoothed simulations show a convectively neutral lower troposphere within the ascent in the flood region at 0000 UTC 20 June (Figs. 15a,c), transitioning to strong convective instability at 0000 UTC 21 June (Figs. 15b,d). However, the magnitudes of both lower-tropospheric θe (328 vs 326 K at 0000 UTC 21 June) and convective instability are larger in the control simulation. This suggests a more favorable environment for heavy precipitation since warmer, moister air (larger θe) is present at the surface and the atmosphere is more convectively unstable. In the 40% smoothed simulation, the strongest ascent is located largely to the west of the flood region (Figs. 15e,f). At both times, the magnitude of ascent in the flood region decreases with increased terrain smoothing (Fig. 15).
In Fig. 16, we investigate the relationship of the condensational heating to the development of the low-level circulation centered south of Calgary by 0000 UTC 21 June (Fig. 11b) using the WRF 27 km. The WRF 27 km was chosen simply for ease of interpretation relative to the noise in the WRF 3-km fields even after smoothing, although the results are qualitatively similar. Figure 16 also includes a diabatic tendency of potential vorticity (PV), which is calculated by taking advantage of the conservation of Ertel’s PV on isentropic surfaces (Holton and Hakim 2013, p. 121), as defined in Eq. (7):
where g is the gravitational acceleration, is relative vorticity on an isentropic surface, f is the Coriolis parameter, and is the change in potential temperature with pressure. Since for frictionless and adiabatic flow, the local time tendency of PV is equal to the horizontal advection of PV in isentropic coordinates [i.e., ]. Any differences with time can then be interpreted as a diabatic PV tendency. This is then compared with the column total CHR as well as the instantaneous PV field to facilitate an understanding of the relationship between the near-surface PV field (and consequent circulation) and the precipitation field. Figure 16 focuses on the 4 h prior to and including 0000 UTC 21 June, as this is the period where the impact of the orography is most apparent in the simulations. As the circulation is relatively shallow, an isentropic surface close to the ground (306 K) is utilized, resulting in some of the PV calculations being missing near the highest peaks in the control runs, although this does not impede the interpretation of the results.
At 2000 UTC 20 June, a large band of condensational heating in the control run is evident in a zonally oriented band north of the flood region, along with a local PV maximum of 1.5 PVU (1 PVU = 10−6 K kg−1 m2 s−1) (Fig. 16a). The heating maximum is in situated in the PV gradient north of the maximum in the region of the strongest implied easterly flow. Although it is not shown, the origin of this PV anomaly lies in a condensational heating plume that originated just northeast of Calgary about 6 h earlier before traveling westward to its current position. However, of crucial importance is the southward extension of condensational heating along the border of Alberta to a local maximum in the extreme southwest corner of Alberta. Collocated with this local maximum is a diabatic PV tendency of approximately 0.15 PVU h−1. At this point, a coherent PV maximum in this region is not evident, although we do see a localized “trough” in the PV field. By 2200 UTC 20 June (Fig. 16c), the region enclosed by the 1-PVU contour in southern Alberta has expanded and we continue to see a diabatic PV tendency of 0.15 PVU h−1 situated underneath a region of condensational heating.
It is after this point, however, that we see an apparent feedback between the evolution of the near-surface PV field and the precipitation, as suggested by the CHR. As the local PV maximum south of Calgary strengthens to over 1.5 PVU by 2300 UTC (not shown), the associated low-level easterlies north of the feature strengthen, resulting in stronger upslope flow and larger heating rates, which in turn produces a larger diabatic PV tendency, with values exceeding 0.4 PVU h−1. This area of heating and subsequent diabatic PV tendency is maintained through 0000 UTC 21 June in the control simulation, as the PV maximum increases to greater than 1.75 PVU and the heating (and implied precipitation) is situated along the immediate eastern slopes of the highest peaks (Fig. 16e).
The above evolution is in stark contrast to the evolution in the 40% smoothed simulation during the same time period. Before 2100 UTC 20 June, the evolutions are somewhat similar (not shown) leading to a setup where a zonally extended PV maximum is situated to the north of the stations of interest (Fig. 16b). However, unlike in the control run, condensational heating is occurring to the west of the Alberta border in a region of WAA associated with the synoptically forced PV maximum to the south. Accordingly, in Fig. 12 we showed that the westward progression of warm air was impeded by the terrain in the control run, but not in the 40% smoothed simulation. Between 2100 UTC 20 June and 0000 UTC 21 June (Figs. 16b,d,f), the PV field is distorted in response to the condensational heating and the subsequent PV tendency, but the region is more diffuse and situated well west of the flood region. The diffuseness of the pattern is important in that the resultant PV maximum fails to achieve the intensity as measured by the gradient in PV, which is crucial to modulating the strength of the wind field.
Physically, the feature evolution described above is perhaps best understood in the context of PV thinking (e.g., Hoskins et al. 1985; Morgan and Nielsen-Gammon 1998). The condensational heating associated with ascent in the foothills flood region acts to redistribute isentropes and PV in the vertical. The tendency in PV due to diabatic processes depends on the orientation of the three-dimensional absolute vorticity vector as it relates to the three-dimensional gradient of diabatic heating (e.g., Martin 2006). Considering that the dominant impact arises from the vertical gradients of condensational heating, the PV increases below the level of maximum heating, and decreases above. The strengthened PV below the condensational heating maximum acts to induce/enhance the cyclonic vorticity (Hoskins et al. 1985), thus contributing to the maintenance of the surface cyclone on the eastern slopes of the Rockies. From an isobaric, QG perspective, the condensational heating maximum forces ascent, and the associated near-surface convergence produces a surface cyclonic vorticity increase and enhanced cyclogenesis.
In summary, a positive feedback occurs between synoptic-scale and mesoscale processes: the initial cyclogenesis is likely a result of QG processes (i.e., differential CVA), causing increased upslope flow in the flood region on the north side of the surface cyclone. Subsequently, the condensational heating associated with the ascent acts to maintain the cyclone through the vertical redistribution of PV. As a result, the durations of QG, frontogenetical, and orographic ascent increase within the upslope flow in the region of heaviest precipitation.
Feedbacks among ascent, static stability, and diabatic heating were also observed by Colle and Mass (2000) and Galewsky and Sobel (2005) in studies of heavy orographic precipitation events. Using model sensitivity experiments, Morales et al. (2015) found that condensational heating intensified a preexisting cyclonic mesovortex during the 2013 Great Colorado Flood that had formed due to large-scale dynamic forcing in a region where vortex formation on the eastern slopes of the Rockies is common, often referred to as the Denver Cyclone. In other words, low-level relative vorticity increased due to the diabatic feedback mechanism described here for the Alberta flood. The results here in the region of strongest condensational heating (Fig. 16) are strongly suggestive of a similar phenomenon.
c. Quantifying orographic ascent
Vertical motion produced purely by the orography (ωoro) is calculated by
where ρ is density, g is the gravitational acceleration, and is representative of vertical velocity in height coordinates (i.e., w) due to orography, calculated using the 10-m vector wind. Written in Bluestein (1992, p. 301) and used by Lin et al. (2001) and Milrad et al. (2015), Eq. (8) allows us to quantify the contributions of orographically forced vertical motion to total vertical motion. Although Eq. (8) diagnoses vertical motion produced directly by the terrain, orography can influence vertical motion above summit elevation, depending on mountain shape and atmospheric stability.
Figures 17 and 18 quantify the contributions of orographic ascent () to total near-surface ascent. At 0000 UTC 20 and 21 June, larger values of orographic ascent are located over the flood region in the control simulation than in the 25% and 40% smoothed simulations. Circulation vectors in the 25% smoothed simulation show upslope flow in the foothills (Figs. 17c,d), but are generally of a lesser magnitude than in the control simulation (Figs. 17a,b); this is supported by larger values of orographic ascent in the control. In the 40% smoothed simulation, the circulation is considerably weaker with a smaller vertical component in the flood region, leading to weaker orographic ascent, particularly after 1200 UTC 20 June (Figs. 17e,f).
Figure 18 presents the fraction of total near-surface vertical motion accounted for by orographic vertical motion. Our focus here is primarily on the western part of the flood region (i.e., Banff and Canmore), in which the duration and magnitude of ascent and precipitation were larger (Milrad et al. 2015; Pomeroy et al. 2016). In the control simulation, accounts for more than 80% of at Banff and Canmore at both times shown (Figs. 18a,b). Progressively smaller fractions of comprise at larger smoothing thresholds at 0000 UTC 20 June (Figs. 18c,e), while at 0000 UTC 21 June, the 25% and 40% smoothed are similar to each other, but less than the control. At 0000 UTC 20 June, and throughout the majority of the event (not shown), accounts for 60%–80% of in the 25% smoothed simulation (Fig. 18c), while in the 40% simulation, it is only 40%–60% (Fig. 18e). These results suggest that the role of the orography in producing ascent is reduced when terrain heights are reduced by at least 25%.
In summary, WRF terrain sensitivity experiments suggest that the Alberta terrain had two primary impacts on the extreme rainfall event that caused the Alberta flood:
5. Discussion and conclusions
The June 2013 Alberta flood was the second costliest natural disaster in Canadian history. An extreme rainfall event during 19–21 June primarily drove the event, which caused 75–150 mm of precipitation in the foothills west of Calgary. The primary objectives of this study were to use 1) high-resolution WRF simulations to further investigate the mesoscale and orographic mechanisms that caused the extreme rainfall event, and 2) terrain smoothing model sensitivity experiments to quantify the importance of the topography.
As found in the CFSR-based study of Milrad et al. (2015), an upper-tropospheric Rex block over western North America resulted in surface cyclogenesis in Montana and southeastern Alberta. The surface cyclone initiated and sustained upslope easterly and northeasterly lower-tropospheric flow into the flood region. Using the ingredients-based approach of Doswell et al. (1996) for flash flooding, we found ascent (i.e., QG, frontogenesis, orographic, and buoyant), moisture (i.e., high-θe advection and a saturated lower troposphere), and convective instability in the flood region, with conditionally unstable air to the east ingested by upslope flow. We therefore suggest that the heavy precipitation was triggered by ascent in the presence of convective (potential) instability, which was primarily sustained by orographic ascent in the presence of an abundance of water vapor being transported into the region by east-northeasterly winds.
Our model sensitivity experiments included four terrain-smoothing thresholds: unsmoothed (control), 10%, 25%, and 40% smoothed. Results show that the magnitudes of accumulated precipitation and ascent decreased gradually from the control to the 25% smoothed simulation. However, the ingredients for heavy precipitation were still similar in the 25% smoothed simulation, albeit with smaller magnitudes. In contrast, the 40% smoothed simulation featured a surface cyclone that stretched across the true location of the Canadian Rockies, resulting in reduced orographic ascent, an unsaturated lower troposphere, and convective instability displaced westward from ascent, particularly later in the rainfall event. In fact, the precipitation maximum in the 40% smoothed simulation occurred in British Columbia on the true western slopes of the Rockies, collocated with the 40% smoothed simulation elevation maximum. Finally, we showed in our control simulation that the long duration of the surface cyclone in southeastern Alberta acted to extend ascent in the flood region throughout day 2 of the rainfall event. Diabatic redistribution of PV on the eastern slopes of the Rockies was the primary mechanism for sustaining the surface cyclone.
In conclusion, the regional orography had both direct (i.e., upslope ascent) and indirect (i.e., increased duration of the surface cyclone on the eastern slopes) impacts on the magnitude of the extreme rainfall event. The spectrum of smoothing thresholds used for our sensitivity experiments suggest that slightly reduced terrain can still result in orographic ascent and increased precipitation totals. For the Alberta flood, terrain reductions >25% result in a near-complete removal of orographic processes, including the diabatic PV feedback mechanism discussed in section 4b.
This research has been supported by grants from the Natural Sciences and Engineering Research Council of Canada (Discovery, and Climate Change and Atmospheric Research). Thanks to NOAA and NCAR for access to reanalysis data, Environment and Climate Change Canada for station precipitation data, and for providing the CaPA dataset. Special thanks to Dr. Kelly Mahoney of NOAA for her invaluable assistance with the postprocessing of the WRF simulations and to Dr. Kevin Brown for his insight and assistance of terrain smoothing procedures. Finally, thanks to the editors Dr. Russ Schumacher and Dr. David Schultz, three anonymous reviewers, and Dr. Lynn McMurdie of the University of Washington, all of whom provided helpful comments and critiques of earlier versions of this manuscript.