Rainfall Enhancement Downwind of Hills Due to Stationary Waves on the Melting Level and the Extreme Rainfall of December 2015 in the Lake District of Northwest England

Abstract

This paper investigates how stationary gravity waves generated by flow over orography enhance rainfall, with particular attention to the role of induced waves in the melting level. The findings reveal a new mechanism by which gravity wave flow focuses precipitation, amplifying rainfall intensity downwind of hills. This mechanism, which depends on the differential velocities of rain and snow, offers fresh insights into how orographic effects can intensify rainfall. A two-dimensional diagnostic model based on linear gravity wave theory is used to investigate the record-breaking rainfall of December 2015 in the Lake District of northwest England. The pattern of ascent is shown to have a qualitatively good fit to that of the Met Office’s operational high-resolution UKV model averaged over 24 h, suggesting that orographically excited stationary waves were the principal cause of the rain. Precipitation trajectories imply that a persistent downstream elevated wave caused by the Isle of Man supported a spray of seeding ice particles directed towards the Lake District, and that these grew whilst suspended in strong upslope flow before being focused by the undulating melting-level into intense shafts of rain.

1. Introduction

It has long been observed that orography has a profound influence on rainfall enhancement and distribution, with ranges of hills, especially those exposed to the prevailing wind, often being wetter than lower-lying or more sheltered areas. Given the importance of the intensity and spatial distribution of rainfall, e.g., to river flow modelling and flood warning systems, mechanisms of enhancement have been much studied conceptually, observationally, analytically, and with numerical models, for the UK and more widely [1,2,3]. One of the first mechanisms to be described was the seeder–feeder effect [4], whereby pre-existing, seeding raindrops sweep out cloud droplets formed in a capping, or feeder, cloud, which is continuously replenished by the ascent of moist air over a hill. Since then, many authors have dealt with the mechanism [5,6]. The seeder–feeder mechanism has been identified as the most important process by which precipitation is enhanced over the moderate (<1000 m) hills of England and Wales [5,7,8]. Descriptions generally invoke some non-orographic origin for the independently existing seeder precipitation, such as frontal ascent, although Browning et al. [9] discuss self-seeding by hills via the release of potential instability, whilst Robichaud and Austin [10] point out that hills of long wavelength (half width more than about 20 km, depending on wind speed) may allow a long enough residence time for cloud droplets to grow to precipitation size within the capping cloud itself, leading to what is effectively a separate category of orographic enhancement. Whilst the precipitation maximum in orographically enhanced rain tends to occur some way up the slope on the windward side of the hill, the effects of wind drift can cause the largest amounts to fall beyond the summit of a hill or range of hills [5,11].

Other authors have pointed to the importance of stationary gravity waves in determining the distribution of precipitation in areas of complex terrain [12]. Colle [13] used two-dimensional simulations to show how the upstream tilt of a gravity wave from a hill crest can displace precipitation distribution in this direction, whilst Stout et al. [14,15] found that differential advection in gravity wave motion preferentially concentrates particle deposition in distinct zones.

Although there has been much progress in the field of orographic precipitation in the sixty years since Bergeron’s exposition of the seeder–feeder effect [4], notably through advances in high-resolution modelling in the last couple of decades, there remain challenges in understanding important processes and improving our conceptual models of them.

Whilst Minder et al. [16] examined influences on the melting level in semi-idealised situations on windward mountain slopes over which there was orographically enhanced precipitation, the role of waves in the melting level in causing precipitation enhancement seems to have received no attention. The study detailed in this paper was conducted to investigate this, using the case of extreme rainfall in December 2015, when the UK 24 h rainfall record was broken at Honister Pass in the Lake District of northwest England. The interest was not in microphysical enhancement mechanisms, but rather in what role gravity wave flow, and in particular wave-induced undulations in the melting level, may have played. A secondary goal was to examine the extent to which a stationary elevated gravity wave may have contributed to the seeder precipitation.

2. Materials and Methods

To examine the extreme rainfall event, 24 h mean forecast fields were taken from the operational 09:00 UTC run of the Met Office’s UKV model on 4 December 2015. This is a 1.5-km-grid-length formulation of the non-hydrostatic primitive equation Unified Model [17]. Output was extracted principally in the form of cross sections of mean values for the 24 h to 18:00 UTC 5 December. These were aligned along the mean wind above the boundary layer through Honister and consisted of horizontal and vertical wind, along with temperature, at full horizontal resolution and on 14 vertical levels. Data from this model relating to moist variables were not available to the author, so to fill this gap, and to provide some verification for the UKV to justify its use as ‘truth’, use was made of the Copernicus European Regional Reanalysis (CERRA) dataset (model and pressure level data available at the Climate Data Store, catalogue entries DOI: 10.24381/cds.7c27fd20 and DOI: 10.24381/cds.a39ff99f). This is produced by a 5.5 km formulation of the Harmonie–Aladin assimilation NWP system with boundary conditions from the ECMWF ERA 5 reanalysis system [18]. Like the UKV data, the CERRA output was used in 24 h averaged form. Being of lower resolution, this reanalysis dataset lacked much of the critical short-wave orographically induced detail seen in the UKV temperature fields and was therefore considered insufficient on its own. It has been stated [3] that when the length scale of valleys and hills is around 10 km, a grid spacing of order 1 km is needed for proper simulation.

A simple gravity wave model, described below, was formulated for the purposes of examining the event. The objective was to provide evidence for the extent to which stationary waves may have initiated and enhanced the rain, by comparing its vertical velocity field with that of the UKV. Experiments with the model also allow attribution of the precipitation genesis and enhancement mechanisms to different topographic features. UKV vertical and horizontal winds are used to generate plausible representative precipitation trajectories from the likely source of the precipitation that fell in the Honister area to investigate kinematic non-microphysical enhancement mechanisms arising from trajectory focusing. For this purpose, a simple Lagrangian model is used, whereby initially regularly vertically spaced precipitation particles are advected with UKV horizontal and vertical winds, their assigned fall speeds relative to the surrounding air changing with their phase transition across a melting level. At each eight-second time step, UKV winds and temperatures are interpolated onto the computed position of the particle being followed, the winds used together with fall speed to compute its position eight seconds later. The extent to which these trajectories diverge or converge in different areas is used to infer depletion or concentration of rain.

A mathematical treatment, largely given in Appendix A, is used to derive equations from first principles that quantify how stationary gravity waves, along with their adiabatically mirrored counterparts on the melting level, can focus rainfall. The kinematic Lagrangian model is applied, but this time with idealised wind and melting-level structure defined, in order both to visualise the enhancement mechanisms for simple cases and to provide quantitative validating results against which to check the derived equations. For the latter purpose, in order to maximise accuracy for comparison with the equations, its time step is shortened progressively as precipitation particles approach the critical surfaces of the melting level and ground, reaching a minimum of 0.01 s.

2.1. Gravity Wave Model

The two-dimensional gravity wave model is based on the analysis of stationary waves, which are characterised by a phase velocity in the upstream direction that is equal and opposite to the horizontal wind (U), allowing them to maintain a fixed position relative to the orographic forcing. By using the linearised momentum, continuity, and thermodynamic equations, and considering wave solutions with phase velocities that counterbalance the flow within which they exist, and with length scales small enough that Earth’s rotational forces are negligible compared to buoyancy forces, the following equations [19] for stationary waves are derived:

$$w=W\exp \left(i\left(kx+mz\right)\right),$$

(1)

$$m^2=N^2/U^2-k^2,$$

(2)

where k and m are horizontal and vertical wavenumbers and N is the buoyancy frequency given by

$$N=\sqrt{{\displaystyle \frac{g}{T}}\left({\displaystyle \frac{\partial T}{\partial z}}+\mathsf{\Gamma }\right),}$$

(3)

where Γ is the adiabatic lapse rate, g is acceleration due to gravity, and T temperature in K.

In the model, orography is decomposed into sinusoidal Fourier components, each exciting stationary waves. For each, the boundary value of W, the amplitude of the vertical velocity forced at the surface, is derived from the product of the horizontal wind and the maximum slope of the terrain component. The stationary wave’s vertical wavenumber m is calculated from Equation (2) and fed into Equation (1), along with the k value appropriate to the Fourier component, to extend W upwards from the surface. Amplitude is increased by a factor of $\sqrt{{\rho }_s/\rho }$, where ρ is density and ρ_s density at the surface. Positive m² leads to proportionality with cos mz through Euler’s formula, and a wave tilt from the vertical given by m/k; but if m² is negative, it results in a term in exp (−|m|z) and, thus, vertical decay (evanescence). Contributions from each Fourier component are combined to arrive at a field of vertical velocity in the x–z plane. The group velocity effects fall out via constructive interference of the resulting waves.

The model was run with 20 vertical levels at a spacing of 500 m, and at a maximum wavenumber of 60 along its 300 km length, giving a smallest resolved wavelength of 5 km. Winds, temperatures, and buoyancy frequency N², as well as being averaged over 24 h, are also spatially averaged on vertical levels for the entire 300 km cross section and, therefore, vary only with height. Only the wind component in the plane of the cross section is used. There is a difficulty in specifying a suitable value of N, static stability being a key parameter in determining the structure of waves. The standard expression assumes an unsaturated atmosphere, in which case the dry adiabatic lapse rate is used for Γ in Equation (3). However, with such large values of rainfall in the case study, it was considered that latent heat release should not be ignored. Whilst it is possible to specify the equivalent buoyancy frequency for a saturated atmosphere, unsaturated layers will occur, especially in zones of descent, and even in largely saturated regions, inhomogeneity of humidity at smaller spatial and temporal scales will likely result in some intermediate value being appropriate to a treatment of time-averaged fields. A pragmatic choice was made to use a value of N in which Γ in Equation (3) takes a value part way between the saturated and unsaturated value, experiments showing that a 55–45% weighted mean of moist and dry values produced a good result. For full accuracy, virtual temperatures should be used in the calculation of both dry and moist Γ, and terms in the vertical variation in total water mixing ratio should be included in moist Γ [20]. However, it is considered that at the temperatures in the case study, these terms would be very small. The limitations of application of such a model to a real case include its two-dimensionality, the linearity of the equations, and the fact that they are strictly valid in vertically uniform conditions. These are discussed in more detail in Section 4.3.

It is recognised that there are formulations with fewer compromises and that can deal better with vertical non-uniformity and non-linearity [21]. However, as demonstrated in the following section, in the case study, the output aligns reasonably well with that of the primitive equation UKV model. Utilising a simple model provides significant benefits. Although fully featured models like the UKV produce more precise simulations, their complexity can obscure causal relationships and complicate the attribution to specific processes. Conversely, simple models can offer valuable insights into the phenomenon and allow attribution not only to atmospheric processes but also to individual topographic features.

2.2. Case Study

In the 24 h period up to 18:00 UTC on 5 December 2015, 341.4 mm fell at Honister in the Lake District of Cumbria, northwest England, setting a new UK 24 h rainfall record, return period calculations suggesting this to be more unusual than a one-in-a-thousand-year event [22]. This event and exceptionally high rainfall totals recorded by other nearby gauges in the catchment of the river Derwent resulted in severe flooding in Cockermouth and Keswick. Sixteen kilometres further east, Brotherswater, in the catchment of the river Eden, recorded 293 mm in 24 h, this river registering its highest level on record, causing flooding in Carlisle [22].

The synoptic situation is shown in Figure 1a,b. A slow-moving pattern maintained a warm, moist airstream over the area through the period, originating at low latitudes in the North Atlantic. Lying on the warm side of the polar front, thermal wind was relatively weak, so winds through much of the troposphere were unidirectional and strong. An exceptionally strong zonal index was present at 18:00 UTC on 4 December, evidenced by a low south of Iceland of 943 hPa and a high in the Bay of Biscay of 1037 hPa, the flow strength and direction showing little change over England through the period. Secondary inner frontal zones analysed crossing the area of interest may well have contributed deep layers of moisture. Synoptic-scale dynamical forcing for ascent does not appear to have been high, but the long, narrow corridor of moisture extending from the Atlantic over the UK can be identified with a warm conveyor belt/atmospheric-river-type structure [23]. Hourly rainfall totals at Honister for the wettest 24 h (Figure 1d) showed a remarkable uniformity for most of the period, suggesting a predominantly steady-state mechanism rather than the passage of transient forcings, and the prolonged nature of the event seems to be a key reason for its severity. Figure 1c gives reported rainfall totals more widely. A 48 h period is used, since some gauges only report every 24 h, and the event spanned two standard rain days. It is likely that a very large proportion of the values given (90% in the case of Honister) fell within 24 h.

Figure 2 shows a more detailed topographic map of the Lake District, a roughly circular region of hills around 50 km across, rising to over 950 m and incised by glacial valleys radiating from a central region, within which Honister is located.

Figure 3 shows the orientation of the baseline of the gravity wave model and the cross-section plots. The closest sounding station to it is around 60 km away at Albemarle (A in Figure 2). The two soundings for 00 and 12 UTC on 5 December are interpolated onto levels at 10 hPa intervals and averaged to be representative of the middle period of the 24 h (Figure 4). Strong, unidirectional flow with a stable layer above hill-top level favour marked gravity wave activity. Note also the transition to unstable, near-adiabatic-lapse-rate conditions above 6000 m/500 hPa, which suggest diminishing wave activity in the upper troposphere.

2.2.1. Interpretation of Gravity Wave Structure from Model Cross Sections

The vertical velocity structures apparent in the case study are examined using output from the UKV and GWM and interpreted in the framework of stationary gravity wave theory. As a consequence of Equations (1) and (2), there is a critical wavelength L = 2π U/N below which waves are vertically aligned and decay in amplitude with height (i.e., exhibit evanescence). Stationary wave components of greater wavelength have an upstream tilt of wave phase with height and downwind transmission of energy along the group velocity vector, which slopes upwards at m/k from the horizontal. Values of U and N in the lower layers give L~22 km. The Isle of Man, inducing a 30 km wavelength stationary wave, transmits wave energy downstream at a slope of m/k~0.8. However, the reduction in stability above 6000 m means this energy becomes trapped, leading to a maximum amplitude in the downwind portion at x = 85 km, z = 6000 m, where both UKV and GWM produce over 0.8 m s⁻¹ (Figure 5).

In the context of frontal ascent, this value would be considered large, e.g., Heymsfield [24] reported values of 0.01–0.1 m s⁻¹ in warm frontal settings, 0.2 m s⁻¹ in warm occluded fronts, and 0.5 m s⁻¹ in convective regions of fronts. The vertical velocity pattern is characteristic of the response to an isolated hill [25,26], and can result in a stream of cold-topped cloud downwind of a hill, often with a sharp upwind edge. When seen in satellite imagery, it is a familiar sign to operational meteorologists of a stationary, or standing, wave. Elsewhere across the domain, an elevated zone of ascent of this sort is not in evidence, presumably because of interference by adjacent features in the more complex orography of the Lake District and Pennines, where the evanescent short wavelengths dominate.

Trapping leads to a low-amplitude oscillation in GWM extending further downwind of the Isle of Man at 6000 m, manifested in a weak minimum at x = 95 km and a weak maximum x = 110 km, but with no obvious counterparts in the UKV.

The low-amplitude surface-based undulation apparent either side of the Isle of Man in GWM is at the scale of the shortest resolved wavelength (5 km) and is considered to be grid-scale noise, an artefact of the numerical treatment. Being of short wavelength, it is evanescent.

Fourier analysis of the orography for the whole of the domain shows the maximum amplitude at a wavelength of 77 km, representing the coarse components of the Isle of Man, Lake District, and Pennines. Being of long wavelength, it generates a rearward-tilting zone of ascent extending back over the sea towards the Isle of Man, shown more clearly in the coarse resolution GWM experiment in Figure 6. The change to vertically aligned phase seen in this figure is due to the transition to lower static stability and strengthening winds aloft associated with the trapping of wave energy previously referred to.

Experiments were conducted with the GWM in which orographic features were selectively excluded in order to attribute the individual, at-a-distance influence of the Isle of Man, Lake District, and Pennines. These confirm that the Lake District and, to a smaller degree, the Pennines contributed to a sloping zone of ascent tilting back across the Irish Sea towards the Isle of Man, and that the Isle of Man caused the elevated column of ascent at x = 85 km.

In general, the UKV tends to confine areas of descent to lower altitudes, such as those over the Isle of Man and Lake District. This may be linked to the transition to lower stability above 3000 m and the result of an inadequacy in the GWM’s equations, which are valid for slowly varying U and N. Some other detail differences are seen between the two models, especially over the eastern side of the Lake District, where the orography varies strongly perpendicular to the baseline of cross section, amplifying the effect of the small normal wind component.

However, the generally good correspondence between the pattern of ascent in this output and that of the UKV, despite the limitations of the GWM, strongly suggests that gravity waves were the main driver of vertical velocity over the period of extreme rainfall and, therefore, of the rain itself. It is acknowledged that the time averaging of the UKV data may mask diabatic or other effects that occurred randomly through the domain, or that dipoles of ascent and descent could pass through the region, leaving little trace when averaged out. However, the unusual steadiness of rain rate at Honister (Figure 1d) through most of the period argues against the importance of any such effects. The column of elevated ascent above 2000 m at x = 85 km, caused by the Isle of Man, is thought to be influential for rainfall over the Lake District and will be examined in the next section.

2.2.2. Precipitation Trajectories and Enhancement

Rather than attempting to model precipitation growth mechanisms, this study concentrates on establishing plausible representative precipitation trajectories to establish the likely origin of the precipitation that fell in the Honister area, as well as non-microphysical enhancement mechanisms based on trajectory convergence, i.e., focusing of rainfall. The term ‘snow’ will be used for all types of frozen precipitation, e.g., in the term ‘snow fall speed’.

Microphysical precipitation growth mechanisms that lead to heavy rain require the ascent of saturated air. The sustained column of rapidly ascending air caused by the Isle of Man, coinciding with moist conditions, could potentially generate significant cloud and precipitation. An extended upstream cross-section analysis suggests this column is the only region of ascent that could account for Lake District rainfall. While strong ascent immediately over the Lake District’s windward hills is crucial for final growth stages, it occurs too close to allow sufficient time for initial precipitation formation.

Although lee descent over the Isle of Man may have dried the air somewhat, evaporation of pre-existing cloud water or ice could mitigate this effect. To assess the precipitation-producing potential of saturated ascent, the following equation is used, adapted from Equation (2) by Grant et al. [27]:

$$C\approx w{\displaystyle \frac{\partial r_s}{\partial z}},$$

where C is the condensation/deposition rate, and r_s is saturation specific humidity. This assumes that vertical motion is much more significant than horizontal for generating condensation and that the rate of change of supersaturation following the motion is much smaller than C, implying vapour converts to ice or liquid upon exceeding saturation (if it is argued that supersaturation should not be ignored, C can be taken to mean the rate at which the total water substance, liquid, solid, and vapour, is made available for precipitation growth processes). A field of C was computed from UKV vertical velocity multiplied by the vertical gradient of volumetric saturation specific humidity (Figure 7a). The distribution pattern, rather than absolute quantitative accuracy, is key here. Due to the temperature- and pressure-dependent reduction in water vapour carrying capacity, the maximum occurs at lower levels than that of the vertical velocity field, centred around 4000 m/−6 °C. The CERRA reanalysis profile of specific snow content at the location of the maximum (Figure 7b) is consistent with this finding. It shows a peak near 4000 m, decreasing slowly above and more rapidly below, in line with the calculated condensation rate (specific snow water content refers to the mass of aggregated ice crystals that can fall as snow).

A range of precipitation particles with varying sizes, forms, and fall speeds is likely to have been present, but nominal particles were traced on the assumption of steady, average fall speeds. Given that the greatest sensitivity in trajectory was to specified snow fall speed, values of 1.0, 1.5, and 2.0 m s⁻¹ were used and precipitation back trajectories calculated from the nearest rain landfall to Honister. These met the column x = 85 km, as shown in Figure 7a at 3200 m (1 m s⁻¹), 3850 m (1.5 m s⁻¹), and 4300 m (2 m s⁻¹). Temperatures at these locations are −2.5 °C, −5.5 °C, and −8.0 °C, respectively. The latter two are close to the peak potential condensation/deposition rates and snow content. For the main trajectories considered later, a 1.5 m s⁻¹ fall speed was specified for snow. This is based on an average between 1 m s⁻¹ for smaller ice precipitation and aggregation snowflakes and 2 m s⁻¹ for graupel, as observed in Doppler radar studies [28,29].

Evaporation of ice particles is particularly efficient [30]. Therefore, it is important for their survival that they remain in an environment saturated with respect to ice once formed. It can be seen in Figure 5a that the whole zone between the Isle of Man and the Lake District along the track of the ice particles is characterised by mean ascent of around 0.2 m s⁻¹, which would likely not only maintain saturation, but facilitate continued growth in a mixed-phase cloud. As mentioned earlier, the rearward-sloping influence from the long wave component forced by the Isle of Man, Lake District, and Pennines contributed to this ascent. The CERRA reanalysis data show cloud to be present along the whole trajectory path until a few hundred metres above the surface, with greatest cloud liquid and ice content just below 1500 m (Figure 8).

A rain fall speed of 6 m s⁻¹ was specified. This is based on values from Best [31], who found a mean drop size of 2 mm for rain whose rate is 10–15 mm h⁻¹, a size associated with a terminal velocity of around 6–7 m s⁻¹ [32]. No allowance was made in computing trajectories for aerodynamic forces on raindrops, although this is expected to make very little difference to the final trajectory, since the rain stage lasts about 4 min compared to the overall time of 36 min from source to landfall around Honister. Yuter et al. [33], using a disdrometer, found 0.5 °C marked a sharp transition from snow to rain, and this temperature is used to define the melting level in the precipitation trajectories.

2.2.3. Near-Isothermal Layer

A stable, near-isothermal layer is evident in the lower troposphere in the UKV and CERRA fields, as shown in the profile in Figure 9, close to the upwind end of the cross section. It is shown later that this layer is an important feature for rainfall focusing, so some attention is given here to its structure and position. Figure 9 also serves to provide some verification for the UKV forecast, which shows a close correspondence to the CERRA reanalysis data. The layer, with a lapse of less than 1 °C km⁻¹ several hundred metres in depth and lying below the melting level over the Irish Sea and Isle of Man, can be seen in Figure 10, trajectories showing it to move with the flow, intersecting the melting level over the Lake District and Pennines, where it weaves around the melting level. It seems that melting of snow from upstream precipitation could have played a part in forming this layer, a process discussed in Section 4.1. Note that its position relative to the melting level varies with vertical motion, and that ascent causes the melting level to sink and the stable layer to rise, with the opposite tendencies with descent. Both the Albemarle sounding data in Figure 4 and the model profiles in Figure 9 show the stable layer to lie at different distances below the melting layer, a result of their position at lower elevations, whereas over much of the Lake District, the stable layer extends either side of the melting level.

Low stability within the boundary layer is likely accentuated by vigorous mixing through mechanical turbulence forced by strong winds, but is also a common feature of profiles beneath isothermal snowmelt layers [6,34]. Regardless of origin, the near-isothermal layer seems likely to have had an important effect on the heavy rain at Honister, as will be advanced in Section 4.1.

2.2.4. Rainfall Focusing

Precipitation trajectories given in Figure 11 were calculated from release points with a vertical spacing of 80 m and arranged in a column spanning the zone of peak precipitation-producing potential as shown in Figure 7. The most intense areas of ascent and descent are on a small enough scale that they show some evanescence, although still with marked effect on the melting level, here defined by the 0.5 °C isotherm, which closely mirrors the topography, suggesting that adiabatic temperature changes are the strongest determinant of its spatial variation, warming induced by descent causing it to rise, cooling on ascent causing it fall. A large rain-free zone is suggested downwind of where strong ascent over the rising ground causes both a dip in the melting level and a levelling out of snow trajectories. Areas of descent then cause the melting level to rise and snow to fall more quickly to meet it, suggesting striking enhancement where trajectories converge. These trajectories are, of course, merely representative of particles with a spectrum of fall speeds for both rain and snow from the values used here, of 6 m s⁻¹ and 1.5 m s⁻¹, respectively. Enhancement, calculated on the basis that rain rate is inversely proportional to horizontal distance apart from neighbouring trajectories (Equation (4)), gives a factor of 4.0 compared with precipitation rates over the sea, but at a distance of some 6 km downwind of Honister (at 154 km), which itself lies within the largely rain-free zone. It is of note that this maximum is very close to another very wet site, Thirlmere, which concurrently set a new UK record of 405 mm over two consecutive rain-days.

Although the UKV has a horizontal grid length of 1.5 km for its dynamics, filtering coarsens the resolution of orography [35], and unsurprisingly, there are discrepancies between the model orography on the line of the cross section and its real-world counterpart, as given in Figure 12. This shows Honister to lie on a slope that drops from an area of high ground about 500 m above sea level. A few kilometres upstream of this is a range of mountains, including Pillar, at 892 m, which is less than 1 km to the north of the cross section at P. Between these two areas of high ground lies the steep-sided valley of the river Liza (L), with an abrupt drop of some 500 m to valley bottom within around 1000 m horizontal distance. The width of the valley at this point is only about 3 km. Therefore, it can only be expected to be poorly resolved at best and, indeed, it is not apparent in the cross section. Figure 12 also shows 48 rainfall totals reported by stations close to the line, the large majority of which fell in 24 h.

The valley of the river Liza (L) widens northwestwards, and to find a UKV cross section that captures it, another (cross section II) was extracted in the same orientation, but 3 km to the northwest of cross section I. Precipitation particle trajectories were then calculated as before and are shown in Figure 13. Here, the drop into the valley at around 146 km is better resolved, giving a marked upward step in the melting level directly above. The resulting rain falls on the projected position of Honister (150 km) with a similar, although slightly lower, level of enhancement (a factor of 3.8 at 152 km). Thirlmere, at 158 km, has a factor of 3.1.

Of course, augmentation of rainfall by wash-out of droplets from the feeder cloud, essentially a cloud physics problem, is not included in this analysis, and neither is the growth by riming whilst ice particles are in near suspension over the windward slopes as supercooled cloud droplets rise around them. It is of note that the trajectories allow very little rain to fall above 400 m on the windward slope of the first range of hills due both to the lifting of snow and the depression of the melting level in strong ascent, where standard seeder–feeder enhancement might normally be expected to be maximised. This supports the notion that the feeder cloud would be less depleted by wash-out on the windward side than in a situation with a higher melting level, allowing it to persist onto the lee side of the first ridge. Not only would this offset any tendency for processes to become dry adiabatic on descent, it would allow more wash-out enhancement here.

Both cross sections show marked enhancement towards the eastern end of the Lake District between 170 and 180 km, due to descent from the ridge (F in Figure 12), which includes such high peaks as Fairfield (873 m) and Helvellyn (950 m). Whilst it might be expected that there should be less scope for low-level enhancement of rain through washout given high upstream rainfall, Brotherswater, which lies 5 km south of the line of cross section I at 174 km, recorded 293 mm in 24 h and 372 mm in 48 h.

The convergence of precipitation trajectories emanating at different heights is indicative of precipitation in the given depth of column between trajectory starting points being spread over a smaller surface area, increasing the rain rate at the expense of reduction elsewhere, but leaving unaltered the total area integrated precipitation. In the next section, the process is idealised and examined analytically.

2.3. Theory

Consider snow with fall speed w_s in an environment of uniform horizontal velocity U that turns to rain on encountering the melting level. This is an idealised situation in which there is no evaporation or sublimation of precipitation, any modification of the environmental temperature by melting snow is neglected, and snow all melts on the same melting surface, ignoring the fact that frozen hydrometeors of different sizes and compositions can persist for different distances below the 0 °C isotherm [16]. If the melting level is fixed and slopes upwards along the flow, the intensity of precipitation is increased, since the resulting rain originating in column Δz is concentrated in a smaller horizontal area Δx_r, as can be seen in Figure 14.

With a horizontal melting level as in Figure 15a, the intensity of precipitation is the same above and below the melting level, again assuming no evaporation or sublimation. Whilst the sloping melting level causes increased water mass flux within the area, the total precipitation integrated over the whole region remains the same. This is in contrast with microphysical enhancement mechanisms, which sweep out cloud water and thereby redistribute extra water to the ground. For this simplified scenario, the enhancement factor, E, can be given by the ratio of the two cross-sectional areas over which precipitation generated in the column of height Δz falls:

$$E={\displaystyle \frac{\mathsf{\Delta }{x}_s}{\mathsf{\Delta }{x}_r}}$$

To solve for E, use is made of the identity

$${\displaystyle \frac{\mathsf{\Delta }{x}_s}{\mathsf{\Delta }{x}_r}}={\displaystyle \frac{\mathsf{\Delta }z}{\mathsf{\Delta }{x}_r}}/{\displaystyle \frac{\mathsf{\Delta }z}{\mathsf{\Delta }{x}_s}}$$

and the fact that in a base situation of no vertical wind and horizontal melting level (as in Figure 15a), snow particles fall with a slope given by

$${\displaystyle \frac{\mathsf{\Delta }z}{\mathsf{\Delta }{x}_s}}={\displaystyle \frac{{-w}_s}{U}}$$

Note that fall speeds are treated as negative because they are in the direction of decreasing z, so that −w_s in the equation above is a positive quantity. Rearranging and taking limits, as Δ → 0, the enhancement of rain rate over the base situation of no vertical wind and horizontal melting level can be given as

$$E=-{\displaystyle \frac{\partial z}{\partial x_r}}{\displaystyle \frac{U}{w_s}}$$

(4)

In Figure 15 are four simplified scenarios based on a wind blowing falling precipitation from the left. The melting level remains invariant with time and topography, where it is not flat, it is shown in black. Precipitation particle trajectories shown were calculated with a simple Lagrangian model in which an abrupt change in fall speed occurs at the geometrically determined melting level. As before, there is a single well defined melting surface and no evaporation or sublimation. In order to visualise the separate mechanisms that contribute to the adiabatic case, effects such as the tendency for the melting level to sink along the flow due to extraction of latent heat by melting snow are excluded, although they will be introduced in the subsequent analytical treatment.

In Figure 15a, which is considered as the base state, there is no vertical wind and a horizontal melting level. Horizontal distances between streams of falling rain and snow are the same and the precipitation rate remains unenhanced and uniform at the base rate.

Figure 15b represents a range of hills that give rise to stationary gravity waves and modulate the melting level, lowering it through adiabatic cooling in areas of ascent, raising it in areas of descent (in this case, the environmental lapse rate is half the adiabatic lapse rate. It is shown later that this gives the melting-level variation the same amplitude as that of the wave motion). Precipitation is concentrated into discrete zones of heavier rain separating dry zones. The relatively low melting level depicted here causes rainfall on the ground to be most concentrated on the lee side of the hill, although this is not an intrinsic characteristic of this scenario, since a higher melting level could cause it to fall on the windward side of downstream hills. Rainfall modulation can be seen as the aggregated effect of the two scenarios in Figure 15c,d.

In Figure 15c, there is no vertical wind, but sinusoidally varying diabatic heating and cooling imposed by some external agent causes the melting level to undulate as in Figure 15b, giving rise to dry and varyingly wet zones, as snow is preferentially intercepted by the melting level, where it slopes upwards along the flow as in Figure 14. This will be referred to as the melting-level slope effect.

In Figure 15d, there is gravity wave motion, but also sinusoidally varying diabatic heating and cooling, which exactly counteracts adiabatic temperature change, leaving a horizontal melting level; a very unrealistic scenario, but one that is contrived to eliminate the melting-level slope effect and isolate a separate modulating mechanism, whereby zones of ascent delay the descent of snow, causing it to become more concentrated in the descending phases. This will be referred to in this article as the bunching effect and is similar to the mechanism investigated by Stout et al. [14,15], although theirs lacked a melting transition.

To examine these scenarios, the simplest wave case is taken, in which the flow pattern is constant with height, as would be the case when the orographic wavelength is equal to the critical wavelength (explained in Section 3.1). As in studies such as those of Hobbs et al. [11], it is assumed that the precipitation particle moves relative to the airstream within which it is embedded only downwards and at its terminal velocity. In doing so, any unbalanced aerodynamic forces are neglected. Use is made of a stability parameter γ, defined as the environmental lapse rate expressed as a proportion of the adiabatic lapse rate, i.e.,

$$\gamma ={\displaystyle \frac{-\partial T/\partial z}{\mathsf{\Gamma }}},$$

where Γ is the adiabatic lapse rate, saturated or dry according to the situation.

Since there is a discontinuity at the melting level, the approach taken is to apply Equation (4) in two stages—firstly, to calculate enhancement along the melting level, and then to extend it from there down to the ground. Derivations, which include a diabatic term to represent the cumulative effect of melting snow on the melting level, are given in Appendix A, with the results summarised here. It is shown that enhancement at the ground over base rate is given by

$$E_g=\left({\displaystyle \frac{w_r+w_g}{w_s}}\right)\left({\displaystyle \frac{{w_m+\gamma (w}_s-D)}{{w_m+\gamma (w}_r-D)}}\right),$$

(5)

where w_r and w_s are the fall speeds of rain and snow, respectively, w_g is the vertical wind at the ground due to the slope of the terrain in the direction of the horizontal wind, and w_m is the vertical wind experienced by a precipitation particle at the melting level. D is the instantaneous diabatically induced local melting-level-height tendency, where the particle crosses the melting level. In other words, it is the local rate of change of the melting-level height; measured in metres per second, it is generally negative and of much smaller magnitude than w_r and w_s. The terms in the inner brackets are effectively fall speeds relative to the melting level, which can be approximated by the absolute fall speeds.

As the vertical temperature profile around the melting layer approaches isothermal, γ → 0, the diabatic terms disappear and the vertical wind at the melting level cancels to give

$$E_g\to {\displaystyle \frac{w_r+w_g}{w_s}}$$

Since w_g, the vertical wind at the ground is normally small compared with w_r, and persistent melting snow inclines the vertical temperature profile towards isothermal, the first-order relationship can be stated as

$$E_g\approx {\displaystyle \frac{w_r}{w_s}}$$

(6)

Taking the example of snow fall speed of −1 m s⁻¹, a rain fall speed of −5 m s⁻¹, enhancement at the ground tends towards a factor of 5, given a near-isothermal environment around the melting level. If there is +0.2 m s⁻¹ (ascent) at the ground due to motion up the slope of a downstream hill, this reduces enhancement to a factor of 4.8, whilst −0.2 m s⁻¹ gives 5.2. Of course, these enhancement values are purely kinematically derived and separate from those due to cloud physics processes, such as feeder cloud washout.

It is shown in Appendix A that the melting-level slope and bunching effects are separately equal when γ = 0.5 (the situation depicted in Figure 15b), and that at higher static stability, the melting-level slope effect dominates, equalling the full effect (diabatic and adiabatic) at γ = 0 (isothermal).

3. Results

3.1. Validation of Analytical Derivation

Equation (5) was tested against the simple time-integrated Lagrangian finite difference model used to create Figure 15 and Figure A1, with D allowed to take on values other than zero. It computes trajectories of falling and melting precipitation, with the vertical wind geometrically defined according to the simplest vertically uniform gravity wave, a uniform horizontal wind, and a response in the melting-level consistent with γ. The method amounts to a finite difference solution in which the trajectory lines are determined by following the motion of precipitation particles and iteratively calculating their position one time step in advance according to their fall speeds and the horizontal and vertical winds at their location. Figure 16 shows a scatter plot of the model results against Equation (5) for a series of different experiments in which varying values of w_r, w_s, D, U, γ, hill height, and length were specified. For these experiments, D varied between 0 and −1.5 m s⁻¹, and γ between 0.8 and 0.2. The very good fit to the analytical values (within 0.3%) gives high confidence in the mathematical derivation of the enhancement equation and in the argument that it gives correct results for idealised situations.

3.2. Applcation of Equations to Case Study

Equations (5) and (6) are also tested on the case study by seeing to what extent they explain the convergence of the precipitation trajectory lines shown in cross sections I and II (Figure 11 and Figure 13). The input values and results are summarised in

Table 1. Input values and enhancement at the ground (E_g), according to Equations (5) and (6), compared with those measured from trajectories for cross sections I and II. See text for comments regarding the estimated value of D.

Parameter	Cross Section I	Cross Section II
wm	−1.2	−0.7
wg	0.0	−0.2
wr	−6.0	−6.0
ws	−1.5	−1.5
D	−0.15	−0.15
γ	0.075	0.075
Eg Equation (5)	3.2	2.9
Eg Equation (6)	4.0	4.0
Eg from trajectories	4.0	3.8

. The minimum resolved temperature lapse in the stable layer (0.6 °C km⁻¹) was used to calculate γ. Consistent with treatment of stability in the GWM, a 55–45% mix of saturated and dry adiabatic lapse rate was used, giving γ = 0.075. The diabatic contribution from snow melt, D, cannot be estimated directly, but as will be discussed in Section 4.2, it would generally be expected to be much less than the snow fall speed w_s, and a value of an order of magnitude less was chosen, although it is shown that the result has a very weak dependence on this.

Feeding the tabulated values into Equation (5) gives values of E_g = 3.2 and 2.9 for cross sections I and II, which are 80% and 76%, respectively, of those suggested by the trajectories (E = 4.0 and 3.8). Rounding values to two decimal places, the same result is given if diabatic effects are ignored (D = 0 m s⁻¹), and giving D a value of −1.5 m s⁻¹, the maximum possible value, only changes the computed values slightly, to 3.1 and 2.8, respectively.

There are a few possible explanations for the mismatch between Equation (5) and the measured values. One is that that the trajectory offset measurements give finite-difference estimates that are subject to errors associated with numerical approximation. Another is that the implicit assumption that the model vertical velocity arises only from movement within the two-dimensional plane leads to errors. More significantly, the theoretical treatment that leads to Equation (5) assumes that waves and horizontal wind are invariant with height. As previously noted, the wave forced by the orographic component of highest amplitude (Figure 6) tilts backwards due to its long wavelength, which could have introduced an extra factor, as suggested by Stout et al. [14].

As shown in Figure 12, the rain gauge closest (3 km) to the cross section on the coastal upstream side recorded 24 mm in 48 h, while Honister registered 382 mm during the same period. A back-trajectory from the coastal site meets the column of ascent at x = 85, z = 3450 m, where the condensation rate and snow amounts are not far off their peak values. Making a first-order simplification that the difference between the two is solely due to orographic modulation, we could characterise the Honister reading as being orographically enhanced by a factor of 15.9 (382/24). Given that the focusing processes dealt with seem to have contributed a factor of about four, this suggests that microphysical processes, such as seeder–feeder washout of cloud droplets, may also have accounted for a factor of about four, implying a comparable magnitude of contribution from each.

4. Discussion and Summary

4.1. Interpretation of Enhancement Equations

Irregularities on the bed of a stream give rise to stationary waves on the water surface [36]. The resulting differential refraction of light leads to stationary patterns of light and dark on the substrate. This effect has some similarities to the melting-level focusing mechanism embodied in equation (5) and can be used as an analogy, albeit while recognising that the physical mechanisms at play are fundamentally different.

The enhancement Equations (5) and (6) relate to situations in which the thermal pattern is fixed, with air flowing through the isotherms as in stationary gravity waves. They are not applicable to structures such as frontal surfaces, where the thermal gradients move largely with the wind and are therefore crossed neither by air parcels nor, horizontally, by precipitation particles. Indeed, with the frontal case, at least in the frontogenetic phase, the transverse thermally direct circulation [37] means that the slope of the melting level has an opposite relationship to the vertical velocity, with a higher melting level and warmer air associated with ascent.

If γ takes on a negative value, the inversional profile indicates that rather than a situation of snow turning to rain, rain turning to snow, or freezing rain, is to be expected, and the established relationships break down. With γ ≥ 1, the atmosphere is statically neutral or unstable and gravity waves do not form. Consequently, the applicability range is defined by 0 ≤ γ < 1.

The enhancement equations derived relate strictly only to an environment with constant horizontal wind speed and with flow patterns that do not vary in the vertical. In the context of wave motion, this is the simplest type of wave, whose phase remains constant in the vertical and that exhibits no change in amplitude with height, consistent with the critical wavelength (defined in Section 2.2.1). Stout et al. [14] did not deal with a melting transition, but found in their treatment of the concentrating effect of the advection field that intensity was increased by the upstream tilt of a wave with height and decreased by evanescence, i.e., a reduction in wave amplitude with height. If only considering the enhancement across the melting level, Equation (5) is not invalidated by evanescence, since for the melting-level slope effect, which dominates at a smaller γ, it is only the amplitude at the melting level that matters, although the upstream tilt found with larger wave components could further increase enhancement.

Because peak vertical velocity is proportional to wavenumber (k) through the relationship w_max = ±kAU, where A is the wave amplitude and U horizontal wind, the largest vertical velocities tend to be driven by the shorter waves. However, these waves, being generally smaller than the critical wavelength, lose amplitude with height through exponential evanescence. Therefore, the lower the melting level, the greater the enhancement potential through larger values of w_m in Equation (5), although this effect becomes less important as static stability decreases towards isothermal. A lower melting level also results in less drift, keeping rain tied more to the lee slopes and, therefore, bringing a positive contribution to enhancement from w_g.

Since air cannot flow into the ground, w_g is only non-zero where the ground slopes along the flow, this component contributing to enhancement on lee slopes and reduction on windward slopes. This treatment has dealt with rain that falls onto ground whose contours follow the shape of the airflow, but in the case of trapped lee waves, the pattern may persist well downwind of the forcing topography, and with it, the enhancement.

4.2. Diabatic Effects

The diabatic effect, whereby melting snow extracts from the environment the heat necessary to make the phase change from ice to water, over time results in a cooling tendency through the melting layer towards 0 °C, which can frequently lead to an isothermal or near-isothermal layer several hundred metres deep [6,34]. Even modest warming from forced descent, e.g., down the lee side of a hill, in these circumstances can produce a steep downwind rise in the melting level. The resulting melting-level slope effect will cause rapid interception of snow by the melting level, leading to the enhancement of rain at the surface, which can be approximated by the rain fall speed divided by the snow fall speed. Even in situations in which the cumulative effect of chilling in the melting layer is not enough to turn the profile isothermal, it will act to decrease γ and, therefore, increase the enhancement according to Equation (5). Such a layer was in evidence in the case study (Figure 4, Figure 9, and Figure 10), and its vertical position relative to the melting layer changed as it flowed over the topography.

The study of Minder et al. [16] on the modelled influences on the depression of the melting level on windward mountain slopes showed that, in addition to the adiabatic effects, orographically enhanced precipitation resulted in a lowering of the melting level through two diabatic effects: the latent heat extraction process mentioned above, and changes in the microphysical melting distance, whereby greater downward penetration of frozen hydrometeors below the 0 °C isotherm is expected where those hydrometeors become larger. All three effects were shown to vary with mountain size, static stability, surface temperature, humidity, and the microphysics scheme employed. However, in Minder et al.’s study, the melting level intersects the ground, which constitutes a significant difference with that presented in this paper. Here, the elevated melting level means that strong ascent, rather than increasing precipitation at the surface, reduces it as snow is lifted. Conversely, more precipitation crosses the melting level in areas of downstream descent, the steepness of the adiabatically induced slope of the melting level being somewhat offset by the diabatic tendency to lower it along the flow. In this way, the diabatic effects tend to operate in opposition to, rather than reinforcing, the adiabatic changes.

The quantification of the diabatically induced rate of change in the height of the melting level, represented by D in Equation (5), is not straightforward and varies depending on a number of factors, including the vertical temperature profile [6,38] and the extent to which air is blocked by topography [39]. Atlas et al. [38] advanced a formula relating the depth to which the 0 °C isotherm extends downwards in unblocked flow due to melting snow, but as Stoelinga et al. point out [6], this formula fails to take account of the countering effect of warm air mixing up from below; this, in turn, is encouraged by instability generated in the layer below by the latent cooling in evidence in the case study (Figure 9).

Because it comes about through the operation of snow as it overtakes the melting level, D might generally be expected to have a magnitude that is significantly smaller than the snow fall speed w_s. Wexler [40] cites a well-marked case in which an elevated melting level in weak, unblocked flow was lowered by 4000 feet (~1200 m) in about 2 h by heavy precipitation, equating to a value of D of around −0.17 m s⁻¹, about an order of magnitude less than w_s. (w_s − D) and (w_r − D) are, in effect, the precipitation fall speeds relative to the melting level, both lessening in magnitude as D increases, proportionally more in (w_s − D) than in (w_r − D), since the magnitude of w_r exceeds that of w_s; therefore, the overall effect of diabatic processes is to reduce enhancement in areas of descent. However, given that D << w_s, like many atmospheric processes, the adiabatic case (D = 0) can be considered as a good approximation; and in any case, the reduction is strongly modulated by static stability via the γ parameter, and as static stability increases towards isothermal, the importance of D diminishes further, vanishing when γ = 0, i.e., with an isothermal profile.

4.3. Gravity Wave Model Limitations

There are several theoretical limitations that should be kept in mind when applying the gravity wave model, whose output is shown in Figure 5 and Figure 6:

• Linearity: The linearised governing equations are theoretically most valid for low hills, but still form a reasonable approximation when the Froude number exceeds one [41]. In this case study, the Froude number for the lower layers, given by U/NH (where H is hill height), is approximately three. However, the linear treatment of trapped waves can still significantly underestimate the amplitude of shortwave components [21], whilst non-linear features associated with trapped waves, such as rotors, can lead to unsteadiness.
• Vertical uniformity: The equations assume that U and N are invariant with height. According to Shutts [42], the equations may be applied to an atmosphere with varying U and N if these parameters change sufficiently gradually relative to the dominant gravity wave wavelength. For the case study, the winds show little variation up to 8000 m, but the stability (N) does vary, which may compromise accuracy.
• Two-dimensionality: The model treats hills as infinite ridges regardless of their extent perpendicular to the cross section. It fails to account for effects from orography outside the cross-section, such as lateral wave deflection [41,43]. Additionally, wind components normal to the cross section may induce lower boundary vertical velocities that differ from the two-dimensional model, especially where terrain varies most in this direction.

Despite these points, the model seems to have performed well enough in the case study for it to be useful, e.g., in the confident attribution to the Isle of Man of the zone of persistent, elevated ascent upstream of the Lake District.

4.4. Summary

The record-breaking rainfall of December 2015 in the Lake District has been examined. A simple gravity wave model gave a similar vertical velocity pattern to that of the operational high-resolution model time-averaged over 24 h, suggesting that rainfall was principally driven by gravity wave motion. A zone of strong elevated ascent caused by the Isle of Man was correctly determined to have been a steady source of seeding ice particles to generate rain over the Lake District. A distinct stable layer is in evidence, intercepting and amplifying undulations in the melting layer that mirror the orography and have a strong focusing effect on computed precipitation trajectories over central and eastern parts of the Lake District, including Honister.

The analytical treatment prompted by this observation has drawn attention to the importance of gravity-wave-induced adiabatic variations in the height of the melting level for the modulation of rainfall, which operates via Equation (5) and, in the limit of an isothermal vertical profile, enhances rain by a factor of rain fall speed divided by snow fall speed (Equation (6)). The resulting enhancement differs from conventional mechanisms in concentrating rain, preferentially in areas of descent. The factors favouring the mechanism include the following:

• A small temperature lapse around the melting level, ideally isothermal;
• A large ratio of rain fall speed to snow fall speed;
• Strong descent at (and, therefore, strong upward slope of) the melting level;
• Strong downslope wind at the surface;
• Proximity of the melting level to the ground.

Strong descent at the melting level and on the ground are, in turn, encouraged by strong horizontal wind speed. Attention has been drawn to the importance of a model’s resolution of narrow valleys, since these can be key features for generating intense shafts of rain from the melting level, which hit the ground downstream.

Equation (5) also embodies the modulating effect of the bunching together of precipitation trajectories through differential advection in wave motion, whilst diabatic effects, such as the extraction of latent heat by melting snow, act to lower the melting level downwind. Both these effects tend to zero, as the vertical temperature profile approaches isothermal.

Equation (5) has been checked against an idealised Lagrangian model and found to be accurate. When applied to the case study, it was found to give about 75–80% of the enhancement value suggested by the computed precipitation trajectories. The additional concentrating effect of the upstream tilt of the larger wave components in evidence in the gravity wave model output may have been responsible for the rest of the enhancement.

In situations that favour strong lee enhancement, precipitation on windward slopes is lessened, as snow is prevented from reaching the melting level by strong ascent. The effect of this is to reduce wash-out of the feeder cloud upwind of the peak, increasing the potential for such seeder–feeder enhancement downwind of the peak. The ratio of the rainfall total at Honister to that at an upstream coastal site indicates, for the case study, that the contributions from the mechanisms embodied in Equation (5) and cloud physical processes are of a similar order of magnitude. The conditions that favour the mechanism enumerated above are not uncommon, especially in the winter half of the year, but further studies would be needed to establish how representative these relative proportions are of the generality of orographic precipitation events in this region and elsewhere. The mechanism also fits with the fact that the largest rainfall totals are often well to the lee of the first range of hills. This is counter to the simplistic expectation encouraged by standard conceptual models that the wettest places should be on the upwind slopes of the most exposed hills. In the case study, the prolonged alignment of the Isle of Man with Honister along the strong flow seems to have been the key element that maximised the seeder precipitation and made the event so unusual.

The effects treated here are, of course, implicit in a primitive equation model’s representation of the atmosphere and the hydrometeors within it, and in that sense, they do not need to be accounted for separately. However, they provide evidence of an enhancement mechanism that seems to have been overlooked, that of the redistribution of precipitation by the kinematics of stationary waves and, more especially, by the response in the melting level.