Extreme Rainfall over Complex Terrain: An Application of the Linear Model of Orographic Precipitation to a Case Study in the Italian Pre-Alps

: Intense meteorological events are the primary cause of geohazard phenomena in mountain areas. In this paper, we present a study of the intense rainfall event that occurred in the provinces of Lecco and Sondrio from 11 to 12 June 2019. The aim of our work is to understand the effect of local topography on the spatial distribution of rainfall and to attempt the reconstruction of a realistic rainfall ﬁeld relative to that extreme event. This task represents a challenge in the context of complex orography. Classical rain-gauge interpolation techniques, such as Kriging, may be too approximate, while meteorological models can be complex and often unable to accurately predict rainfall extremes. For these reasons, we tested the linear upslope model (LUM) designed for estimating rainfall records in orographic precipitation. This model explicitly addresses the dependence of rainfall intensiﬁcation caused by the terrain elevation. In our case study, the available radio sounding data identiﬁed the convective nature of the event with a sustained and moist southern ﬂow directed northward across the Pre-Alps, resulting in an orographic uplift. The simulation was conducted along a smoothed elevation proﬁle of the local orography. The result was a reliable reconstruction of the rainfall ﬁeld, validated with the ground-based rain gauge data. The error analysis revealed a good performance of the LUM with a realistic description of the interaction between the airﬂow and local orography. The areas subjected to rainfall extremes were correctly identiﬁed, conﬁrming the determinant role of complex terrain in precipitation intensiﬁcation.


Introduction
In alpine environments, hydrogeological hazards are mainly triggered by heavy rainfall events that can occur across the seasons in different ways [1][2][3]. The autumn is the most critical season, where extratropical cyclone structures can strike the alpine range frequently, affecting its southern flank slope [4][5][6]. These episodes can cause extensive damage, injuries, and deaths due to the possible triggering of critical events, such as floods and landslides [7,8]. In the literature, several approaches have been adopted to better understand the cause-effect relation between heavy rainfall and hydrogeological issues.
One of the most important is the elaboration of the intensity-duration curves where precipitation features are intended as predictors for the event triggering. These curves correlate the intensity and the duration of a precipitation event, defining a critical rainfall event with respect to a threshold above which landslide failures are highly expected [9][10][11][12][13]. For floods and flash floods, the event intensity is assessed through the evaluation of the correspondent return period of the rainfall, under the hypothesis of iso-frequency among floods and precipitation [1,[4][5][6]14,15]. As a result, the information acquired about rainfall characteristics is fundamental for the correct assessment of these hydrogeological events.
All these approaches have different levels of approximation but require high quality input data. These data are retrieved from rain gauges time series. Unfortunately, rain gauges networks are not always uniformly spatially distributed, and they only provide single-point information. To overcome this limitation, compromises are adopted to deal with ungauged areas. In some cases, the nearest rain gauge is taken into account directly; however, more generally, a rainfall field is reconstructed from a station network, applying automatic interpolation methods, such as the Inverse Distance Weighted (IDW) or Kriging methods [16][17][18][19][20].
Interpolation techniques work well across flat areas where the terrain elevation is nearly constant or very smooth [18]. The simplest ones, such as IDW, are geometry-based where interpolated values depend only on the rain gauge distance. In this case, the availability of a sufficiently dense network is necessary to obtain a reasonable interpolation; however, this is sometimes not possible in mountain environments. The dependence of the rainfall field on elevation is not taken into account. In this light, geostatistical techniques, such as Kriging External Driven (KED) may be a viable solution: the classical Kriging could be conditioned (driven) considering an external drift, such as the elevation [18,21,22].
KED yields good results when a strong relationship between a meteorological variable and elevation is present. The spatial interpolation of temperature is an example of a successful application of KED [23,24]; however, this cannot be extended to precipitation. Precipitation follows neither an additive error model nor a Gaussian distribution (typical of temperatures), which are both prerequisites for a rigorous application of most geostatistical methods [18,25]. Instead, precipitation follows a multiplicative error model, and its distribution is positively skewed [25,26]. Therefore, the Kriging performances are low, especially in the reconstruction of daily and sub-daily rainfall fields.
To avoid these problems, several authors have proposed extensions of these interpolation methods. One is the PRISM model (Parameter Elevation Regressions on Independent slopes Models) [27][28][29]. The former is based on a weighted climate-elevation regression function that acknowledges the dominant influence of elevation on precipitation. To operate, for each meteorological station, there are assigned weights that account for other physiographical factors in addition to elevation, such as the topographic exposure or coastal proximity, which affect the climate at a variety of scales. The PRISM model was widely tested all over the world for other meteorological parameters, obtaining good responses by the scientific community for climatic data interpretation [28][29][30].
However, the PRISM model was not successfully applied to interpret localized events, such as the orographic thunderstorm triggered in watershed basins [31]. In these cases, Limited Area Models (LAMs) for local weather forecasting represent an important data source [5,14,[31][32][33]. Their affordability is still increasing worldwide particularly in nowcasting meteorology, and they are continuously updated and trained in the prediction of atmosphere dynamics [31,32]. Historically, one of the main drawbacks of these models has been the complexity of their parameterization regarding the cloud microphysics, which is the "engine" responsible for the formation of precipitation [5,31,32].
Difficulties have risen across mountain areas where the morphology could dramatically influence the rainfall processes, even at a local level [34]. Result variability is also dependent on the microphysical scheme implemented inside the model [5]. Currently, according to [35], the lack of predictability after a few hours represents a current challenge aspect for small-scale precipitation events. In particular, for scales smaller than~60 km, predictability is lost rapidly within the first 6 h of the forecast, which is the typical scale length of the intense precipitation phenomena. For this reason, the design and the use of a convection-permitting ensemble appears promising [35,36]. However, also in these cases, meteorological outputs should be always compared and validated through appropriate error analysis with rain gauges records on the ground [14,25,26].
Bearing in mind that a unique technique to determine rainfall amounts for ungauged regions does not exist, especially for remote places in mountain territories, we focused on a model able to reconstruct a reliable rainfall field. The aim of our study was to attempt to overcome the poor accuracy of classical interpolation techniques, while avoiding the complexities of running a LAM. For these reasons, we tested the application of the "linear upslope model" (LUM) of orographic precipitation proposed by [37][38][39][40]. Three aspects drove our choice: the rigorous but rather simple mathematical formulation that is inspired by the complex physics of meteorological models; the flexible applicability that required limited spatial data, such as the local morphology, i.e., elevation, and available radiosonde data [31,40]; and the extreme event of orographic rainfall that was recorded in the case study area.
The idea behind the linear upslope model is to reproduce, with sufficient accuracy, the mechanism of rainfall over complex morphology without considering the cloud microphysics. In our case study, we experimented with a typical situation of orographic precipitation where the territory of the provinces of Lecco and Sondrio ( Figure 1) experienced an average rainfall depth of 110 mm in 13 h [41], with an extreme that reached 210 mm. That episode suggested the application of LUM. The aims were twofold: to better understand the role of orography in rainfall intensification and to reproduce a realistic rainfall field for the considered area.
Bearing in mind that a unique technique to determine rainfall amounts for ungauged regions does not exist, especially for remote places in mountain territories, we focused on a model able to reconstruct a reliable rainfall field. The aim of our study was to attempt to overcome the poor accuracy of classical interpolation techniques, while avoiding the complexities of running a LAM. For these reasons, we tested the application of the "linear upslope model" (LUM) of orographic precipitation proposed by [37][38][39][40]. Three aspects drove our choice:  the rigorous but rather simple mathematical formulation that is inspired by the complex physics of meteorological models;  the flexible applicability that required limited spatial data, such as the local morphology, i.e., elevation, and available radiosonde data [31,40]; and  the extreme event of orographic rainfall that was recorded in the case study area.
The idea behind the linear upslope model is to reproduce, with sufficient accuracy, the mechanism of rainfall over complex morphology without considering the cloud microphysics. In our case study, we experimented with a typical situation of orographic precipitation where the territory of the provinces of Lecco and Sondrio ( Figure 1) experienced an average rainfall depth of 110 mm in 13 h [41], with an extreme that reached 210 mm. That episode suggested the application of LUM. The aims were twofold: to better understand the role of orography in rainfall intensification and to reproduce a realistic rainfall field for the considered area. The model is described in detail in the Materials and Methods section together with a presentation of the case study. In the Results section, we first present an elaboration of the radiosonde parameters used for model initialization, followed by the outcomes of the rainfall field reconstruction and an error analysis regarding the rain gauge measurements. Finally, our Discussion and Conclusion of the back-analysis study are reported with our comments. The model is described in detail in the Materials and Methods section together with a presentation of the case study. In the Results section, we first present an elaboration of the radiosonde parameters used for model initialization, followed by the outcomes of the rainfall field reconstruction and an error analysis regarding the rain gauge measurements. Finally, our Discussion and Conclusion of the back-analysis study are reported with our comments.

The Linear Upslope Model (LUM)
LUM evaluates the spatial evolution of a precipitation P (mm) triggered by the local orography h (m). The idea behind the model formulation is the following: when a humid airflow rises along a slope, it starts to condense, cooling adiabatically and triggering rainfall on the up-slope flank of the mountain range. Beyond the mountain peak, the airflow begins to descend and dries out due to adiabatic warming, causing a decrease in rainfall. This simple conceptualization neglects all the cloud microphysics and airflow dynamics that characterize the complex atmosphere interaction with the local orography. For example, the perturbation of airflow trajectories after crossing a mountain peak are not explicitly simulated by the model; however, the water vapor depletion due to generated rainfall is taken into account [37,40].
P is expressed as a function of only two terms: the Water Vapour Flux and the local terrain slope. Water Vapour Flux (WVF) is an indicator of the airmass moisture and represents the integrated flux of water substance across vertical section of the atmosphere [39,40]. The WVF module is expressed in kg s −1 m −1 and is computed as: where ρ is the air density kg/m 3 , q v is the specific humidity %, and U(z) is the horizontal component of the wind vector m/s retrieved from balloon sounding where atmospheric parameters are measured with altitude. U(z) attributes the direction to the WVF vector. The airflow evolution across space is described by the atmosphere continuity equation (Equation (2)) where the horizontal divergence of Water Vapour Flux vector (WVF) is equal to the difference among two terms: the precipitation P, that depletes WVF, and the surface evaporation E (mm), that refills WVF.
Under the hypothesis of a steady state and a mono-dimensional airflow directed along positive x-axes, the WVF module can be expressed simply as a function of the spatial coordinate x. Equation (2) can be rewritten as Equation (3) making explicit the change of the WVF module along an x-path in two different locations, A and B [40]: The P and E components can be expressed as a function of the WVF vector. Following the assumption proposed by [37,39], an airflow can rise up a mountain flank following the same slope as the underlying terrain. The airflow exhibits a vertical wind velocity, W, as a function of the horizontal component U(z) and terrain slope. Bearing in mind the mono-dimensional approximation, the elevation h is represented by a profile h(x) along the positive x direction.
If the airflow rises following a moist-adiabatic lapse rate, i.e., a temperature gradient of around 6.5 • C/1000 m, and it is mostly saturated at all levels, the condensation rate, CR, kg s −1 m −1 can be evaluated as the vertical integral of the wind velocity W(x,z) m/s, the air density at saturation ρ sat kg/m 3 , and the saturation scale height H −1 sat m −1 .
The saturation scale height H −1 sat for water vapour ascending moist-adiabatically is approximated using the Clausius-Clapeyron theory, and can be expressed as a function of the latent heat of water, L = 2.5 × 10 6 J kg −1 , the gas constant for vapour R V = 461 J kg −1 K −1 , the temperature at the ground T 0 ∼ = 20 • C = 293 K, and the moist-adiabatic lapse rate, γ. Adjusting the previous Equation (5b), the precipitation rate, P, can be rewritten from the CR as in Equation (6) highlighting its dependence from the WVF vector and the terrain slope. Equation (6) states that, if the terrain elevation is increasing in the direction of the incoming upslope wind, the condensation of the air mass can generate precipitation in that area.
On the other hand, when the WVF vector is approaching the downslope, the contribution due to the gradient is the opposite and the air mass tends to evaporate, due to the adiabatic warming (Equation (7)).
LUM was implemented following these steps: 1.
An initial condition of WVF vector was defined to initialize the model at coordinate x = 0.

2.
Considering the local slope evaluated from the elevation profile and the WVF vector, P and E were estimated using Equations (6) and (7). 3.
The continuity equation (Equation (3)) was considered to retrieve the new value of the WVF vector at coordinate x > 0. 4.
The operation at points 2 and 3 was repeated until the end of the elevation profile.
The mechanism depicted by LUM is rather realistic even though simplified. The suitability has been tested on several occasions by authors who applied the methodology to different domains all over the world [37][38][39][40].

The Case Study of 11-12 June 2019
During the night between 11 and 12 June 2019, an extreme convective rainfall event occurred in the upper part of the Lake of Como, affecting the territory of the provinces of Lecco and Sondrio ( Figure 1). Rather persistent and auto-regenerating thunderstorms started during the evening of 11 June around 8 p.m. and did not dissipate completely until 9 a.m. the following day. A huge amount of rainfall fell with an average rainfall depth of 110 mm in 13 h [41]. The town of Premana experienced a total of 210 mm in 13 h that corresponds to a precipitation with a return period of 200 years. This datum has been estimated considering the local "Intensity-Duration-Frequency" (IDF) from [41]. All rainfall depths recorded in the study area and reported from local environment agency ARPA Lombardia [41] are listed in Table 1. Figure 1A,B depict the location of the study area in the Lecco and Sondrio Provinces (Lombardy Region, Italy) and the local rain gauge network from Arpa Lombardia. This event caused serious damage. As shown in Figure 1, thunderstorms followed a narrow corridor path along a south-north section, striking the area near the town of Premana. Due to the extreme precipitation, a sudden flash flood in the Varrone river hit Dervio, a town located downstream to Premana. The Pagnona dam risked a dam break due to the reservoir overtopping, and thus Dervio's inhabitants were temporarily sheltered. In Sondrio province, several shallow landslides and debris flows affected the slopes of Mount Legnone in the municipality of Delebio and Piantedo, and a road in the S. Giacomo Filippo municipality collapsed due to Liro river flood.
All these events happened in the same area and in a relatively short time. For these reasons, particular attention was posed to the analysis of meteorological maps of that episode [42]. Looking at a regional scale, the position of air masses in Figure 2A,B show a low pressure centred on western France. This configuration was responsible of the warm and humid air advection coming from Mediterranean Sea in the direction of the central sector of the Alps. The temporal persistency of heavy rainfalls was caused by the stationarity behaviour of the air masses that were not advected eastward faster as it is typically seen during a thunderstorm episode (Figure 2A,B shows a not appreciable difference over 12 h).
That humid airflow was also sustained by the presence of intense jet streams as reported in Figure 2C. The core velocities reached values around 180-190 km/h at 300 hPa geopotential height contributing to increasing the air upslope motions due to upper layer divergence [5,43]. The straight path followed by rainfall extreme events depicted in Figure 1B is, therefore, defined by the upstream currents that advected thunderstorm cells northward ( Figure 2C,D). Premana. Due to the extreme precipitation, a sudden flash flood in the Varrone river hit Dervio, a town located downstream to Premana. The Pagnona dam risked a dam break due to the reservoir overtopping, and thus Dervio's inhabitants were temporarily sheltered. In Sondrio province, several shallow landslides and debris flows affected the slopes of Mount Legnone in the municipality of Delebio and Piantedo, and a road in the S. Giacomo Filippo municipality collapsed due to Liro river flood. All these events happened in the same area and in a relatively short time. For these reasons, particular attention was posed to the analysis of meteorological maps of that episode [42]. Looking at a regional scale, the position of air masses in Figure 2A,B show a low pressure centred on western France. This configuration was responsible of the warm and humid air advection coming from Mediterranean Sea in the direction of the central sector of the Alps. The temporal persistency of heavy rainfalls was caused by the stationarity behaviour of the air masses that were not advected eastward faster as it is typically seen during a thunderstorm episode (Figure 2A,B shows a not appreciable difference over 12 h).
That humid airflow was also sustained by the presence of intense jet streams as reported in Figure 2C. The core velocities reached values around 180-190 km/h at 300 hPa geopotential height contributing to increasing the air upslope motions due to upper layer divergence [5,43]. The straight path followed by rainfall extreme events depicted in Figure  1B is, therefore, defined by the upstream currents that advected thunderstorm cells northward ( Figure 2C,D). The exceptionality of this event could be explained only considering the local effect of orography on airflow motion. According to [4,6,24,30,34], two more meteorological forcings were identified as responsible for that intensity. The first regards the local lowlevel convergence of moist southerly flows as reported in Figure 3, which shows wind The exceptionality of this event could be explained only considering the local effect of orography on airflow motion. According to [4,6,24,30,34], two more meteorological forcings were identified as responsible for that intensity. The first regards the local lowlevel convergence of moist southerly flows as reported in Figure 3, which shows wind trajectories and velocities at 2000 m. The second is the presence of the orographic barrier of Orobie Alps (Figure 1) that was hit by an incoming southerly flow triggering the upslope movement of airmasses. These two elements were simulated in more detail considering the linear upslope model presented. trajectories and velocities at 2000 m. The second is the presence of the orographic barrier of Orobie Alps ( Figure 1) that was hit by an incoming southerly flow triggering the upslope movement of airmasses. These two elements were simulated in more detail considering the linear upslope model presented.

Results
In this section, the application of the LUM to the case study is reported. In the first and second part, the initial conditions of LUM were retrieved from radiosonde data, and hypotheses on model applicability are discussed. In the third part, the model results are validated through error analysis comparing with rain gauge recordings.

Convective Nature of the Rainfall Event
Radiosondes are a valuable instrument for the exhaustive definition of the 3D structure of the atmosphere [5,31,32]. They are equipped with a small aerostatic balloon and launched every 12 h from selected sites around the globe. These balloons carry a small box containing meteorological sensors for the vertical profiling [46] of quantities, such as the pressure (hPa), elevation (m), temperature (°C), relative humidity (%), dew point temperature (°C), wind velocity (m/s), and wind direction (°).
In our case study, sounding data collected near the Linate Milan Airport (60 km distance from Lecco, shown in Figure 1) were provided by University of Wyoming [46]. Two radio soundings were examined: the 12 h UTC of 11 June and the 00 h UTC of 12 June (Figure 4). The first was intended as a precursor of the atmosphere state favourable to trigger the extreme rainfall. The second was indeed representative of the airflow dynamic that occurred during the rainfall. Therefore, both soundings were used for assessing the initial condition WVF0 for LUM.

Results
In this section, the application of the LUM to the case study is reported. In the first and second part, the initial conditions of LUM were retrieved from radiosonde data, and hypotheses on model applicability are discussed. In the third part, the model results are validated through error analysis comparing with rain gauge recordings.

Convective Nature of the Rainfall Event
Radiosondes are a valuable instrument for the exhaustive definition of the 3D structure of the atmosphere [5,31,32]. They are equipped with a small aerostatic balloon and launched every 12 h from selected sites around the globe. These balloons carry a small box containing meteorological sensors for the vertical profiling [46] of quantities, such as the pressure (hPa), elevation (m), temperature ( • C), relative humidity (%), dew point temperature ( • C), wind velocity (m/s), and wind direction ( • ).
In our case study, sounding data collected near the Linate Milan Airport (60 km distance from Lecco, shown in Figure 1) were provided by University of Wyoming [46]. Two radio soundings were examined: the 12 h UTC of 11 June and the 00 h UTC of 12 June (Figure 4). The first was intended as a precursor of the atmosphere state favourable to trigger the extreme rainfall. The second was indeed representative of the airflow dynamic that occurred during the rainfall. Therefore, both soundings were used for assessing the initial condition WVF 0 for LUM.
To facilitate the interpretation of the atmosphere state, convective parameters, such as CAPE, CIN, and LI ( [5,34,47]), and height parameters, such as LCL, LFC, and EL ( [5,32,34]) have been retrieved from radio soundings: CAPE is the acronym of Convective Available Potential Energy, and it is an indicator of atmospheric instability, a necessary condition severe weather hazard. CIN stands for Convective Inhibition and represents the amount of energy required to overcome the negatively buoyant energy that the environment exerts on an air parcel. The latter is not the opposite of CAPE, but the two indexes give complementary information: when CAPE is high, above 1000 J kg −1 and CIN is low, less −30 J kg −1 , the probability of convective triggering is high. LI (Lifted Index) is the temperature difference between the environment T e and air parcel lifted adiabatically T p at a given pressure height, usually 500 hPa.
LCL is the Level of Condensation and represents the height where the cloud formation is started. LFC (Level of Free Convection) is the height where the air parcel can rise-up without external forces and thunderstorms systems can form. EL (Equilibrium Layer) is the height where vertical motion is stopped and generally corresponds to the thunderstorm cap or anvil.  To facilitate the interpretation of the atmosphere state, convective parameters, such as CAPE, CIN, and LI ( [5,34,47]), and height parameters, such as LCL, LFC, and EL ( [5,32,34]) have been retrieved from radio soundings:  CAPE is the acronym of Convective Available Potential Energy, and it is an indicator of atmospheric instability, a necessary condition severe weather hazard.  CIN stands for Convective Inhibition and represents the amount of energy required to overcome the negatively buoyant energy that the environment exerts on an air parcel. The latter is not the opposite of CAPE, but the two indexes give complementary information: when CAPE is high, above 1000 J kg −1 and CIN is low, less −30 J kg −1 , the probability of convective triggering is high.  LI (Lifted Index) is the temperature difference between the environment Te and air parcel lifted adiabatically Tp at a given pressure height, usually 500 hPa.  LCL is the Level of Condensation and represents the height where the cloud The values reported in Table 2 confirm the instability conditions of the atmosphere column. This fact is highlighted primarily by Brunt-Vaisala frequency N 2 , which assumes negative values [5,32]. Looking at the convective parameters, the first sounding is the most critical. The former exhibits the highest CAPE and the lowest CIN, which correspond to conditions of "moderate" instability. LI is low, depicting a situation where severe storms with strong lifting mechanisms have a high probability to occur. The second sounding exhibits lower values of CAPE and LI and higher values of CIN since the thunderstorms have been already triggered and have started to dissipate the atmosphere energy.
Considering the height parameters, the thunderstorm-prone structure of the atmosphere can be well recognized. Comparing the two soundings, this configuration remained stationary for 12 h, confirming the previous hypothesis coming from air masses analysis. The LCL was settled slightly above the upper boundary layer limit, and the LFC was around 2000 m, at the top of the mountain peaks of the Orobie Alps. Therefore, convective cells could have been triggered easily above the mountain range and evolved freely up to the EL that was settled around 10 km.
In summary, the convective nature of the extreme rainfall event was clearly described by the radio sounding even if located slightly far away from the area studied. All the analysed parameters indicated the instability of the atmosphere but were not sufficient to explain the severity of the phenomena recorded over the Alps. The soundings data provide us with two essential parameters required by LUM: the initial conditions and the precipitation efficiency ratio.

The Initial Conditions and the Precipitation Efficiency Ratio
The continuity equation of LUM (Equations (2) and (3)) requires an initial condition. As specified in the Materials and Methods section, we described the model under two hypotheses: a steady state condition and a mono-dimensional geometry along the x coordinate. Therefore, the initial condition for the incoming airflow is defined by the vector WVF 0 calculated at location x = 0. Looking at Table 3, the WVF 0 vector was retrieved for each sounding described previously, plus a third one recorded on 12 June 2019 12 h UTC. As reported in Equation (4), the module of WVF vector is represented by the integral over an air column height. The integral limits were assessed looking at the sounding data: they were comprised between the LCL (lower) and EL (upper) levels. For each sounding, validity windows were defined as 6 h backward and 6 h forward from the central reference time. Then, considering that the rainfall event started at 8 p.m. of 11 June until 9 a.m. of 12 June, the two soundings of 12 June 2019 00 h UTC and 12 June 2019 12 h UTC were selected to compute the average module of WVF 0 at location x = 0.
Inside LUM, the microphysics of rainfall processes is not explicitly considered since an instantaneous transformation between condensate quantity in the cloud (CR) and rain, P is assumed [31,39,40]. Therefore, a correction term should be introduced to take into account the efficiency of these processes that can significantly reduce the effective rainfall amount [48]. The precipitation efficiency, E p , is expressed as the ratio of the P atm , rainfall computed by a generic atmospheric model and P gauge , the rainfall recorded by rain gauges on the ground: Several authors correlated E p with different parameters of radio sounding, and a rather useful relationship can be found considering the wind shear. The former represents the wind velocity difference between the ground and a reference height of 6000 m [5,32]. Using the relationship proposed by [48] in Figure 5, from the three-radio sounding, a reasonable precipitation efficiency value was assessed to be around 0.3 (Table 4).
Inside LUM, the microphysics of rainfall processes is not explicitly considered since an instantaneous transformation between condensate quantity in the cloud (CR) and rain, P is assumed [31,39,40]. Therefore, a correction term should be introduced to take into account the efficiency of these processes that can significantly reduce the effective rainfall amount [48]. The precipitation efficiency, Ep, is expressed as the ratio of the Patm, rainfall computed by a generic atmospheric model and Pgauge, the rainfall recorded by rain gauges on the ground: Several authors correlated Ep with different parameters of radio sounding, and a rather useful relationship can be found considering the wind shear. The former represents the wind velocity difference between the ground and a reference height of 6000 m [5,32]. Using the relationship proposed by [48] in Figure 5, from the three-radio sounding, a reasonable precipitation efficiency value was assessed to be around 0.3 (Table 4). Table 4. Wind shear and precipitation efficiency for the three radio soundings.

Linear Upslope Model Applied to the Case Study
The application of the LUM to our case study required further simplifications to make the problem more tractable. The first regards the mono-dimensional approximation of the domain geometry. The event analysis has depicted a clear south-north corridor followed by thunderstorms ( Figure 1). That rectilinear path can ideally define a bidimensional domain that slices the atmosphere along the horizontal and vertical coordinates x and z. Bearing in mind that the quantity WVF is vertically integrated, the problem is reduced to one dimension.

Linear Upslope Model Applied to the Case Study
The application of the LUM to our case study required further simplifications to make the problem more tractable. The first regards the mono-dimensional approximation of the domain geometry. The event analysis has depicted a clear south-north corridor followed by thunderstorms (Figure 1). That rectilinear path can ideally define a bidimensional domain that slices the atmosphere along the horizontal and vertical coordinates x and z. Bearing in mind that the quantity WVF is vertically integrated, the problem is reduced to one dimension.
For an accurate representation of the ground elevation h along the thunderstorm path depicted in Figure 1, a group of five parallel traces was extracted from a Digital Elevation Model (DEM). The DEM has a 5 m cell resolution and was freely downloaded from geoportal of Lombardy Region [49]. The obtained traces were then averaged in a unique profile to derive a representative value of h(x) along the south-north section. As reported by [31,37,40,50], the LUM results were not well correlated with very detailed DEMs. The airflow was generally not sensitive to a local abrupt slope change but was instead conditioned by the average height field. Consequently, we resampled the original DEM at a coarser resolution of 500 m to obtain a smoother representation of the local terrain orography.
Inside the LUM, the presence of the boundary layer (BL) near the ground was not explicitly taken into consideration. BL is a portion of airmass that tends to be at rest in contrast to the free atmosphere [5,32]. The height of BL is difficult to estimate with a closed formulation as it depends on several variables and has a typical diurnal cycle [5,32,34,51,52]. It tends to be completely depleted after intense rainy events due to the increased turbulence caused by downdrafts [53,54]. We included this parameter in our simulation, nevertheless. Using the soundings data and looking at temperature inversion in proximity of the ground, the BL height was estimated. As shown in Table 5, the diurnal fluctuation of the BL was in the order of 1000 m, and this represents the range considered for the study. The BL was included in the LUM formulation proposing a new elevation profile that can differ significantly from the real terrain. In particular, the areas settled under the BL height were automatically flattened because they were assumed to be excluded from rainfall formation processes. Figure 6 shows the comparison between the real terrain profile and the adjusted one considering BL. A reference value of boundary layer of 1000 m was adopted for the Val Padana side, and an indicative value of 600 m was used for the Valtellina floodplain. The latter was decreased considering the possibility that local valley winds are responsible for its rapid depletion [31,52,55]. For an accurate representation of the ground elevation h along the thunderstorm path depicted in Figure 1, a group of five parallel traces was extracted from a Digital Elevation Model (DEM). The DEM has a 5 m cell resolution and was freely downloaded from geoportal of Lombardy Region [49]. The obtained traces were then averaged in a unique profile to derive a representative value of h(x) along the south-north section. As reported by [31,37,40,50], the LUM results were not well correlated with very detailed DEMs. The airflow was generally not sensitive to a local abrupt slope change but was instead conditioned by the average height field. Consequently, we resampled the original DEM at a coarser resolution of 500 m to obtain a smoother representation of the local terrain orography.
Inside the LUM, the presence of the boundary layer (BL) near the ground was not explicitly taken into consideration. BL is a portion of airmass that tends to be at rest in contrast to the free atmosphere [5,32]. The height of BL is difficult to estimate with a closed formulation as it depends on several variables and has a typical diurnal cycle [5,32,34,51,52]. It tends to be completely depleted after intense rainy events due to the increased turbulence caused by downdrafts [53,54]. We included this parameter in our simulation, nevertheless. Using the soundings data and looking at temperature inversion in proximity of the ground, the BL height was estimated. As shown in Table 5, the diurnal fluctuation of the BL was in the order of 1000 m, and this represents the range considered for the study. The BL was included in the LUM formulation proposing a new elevation profile that can differ significantly from the real terrain. In particular, the areas settled under the BL height were automatically flattened because they were assumed to be excluded from rainfall formation processes. Figure 6 shows the comparison between the real terrain profile and the adjusted one considering BL. A reference value of boundary layer of 1000 m was adopted for the Val Padana side, and an indicative value of 600 m was used for the Valtellina floodplain. The latter was decreased considering the possibility that local valley winds are responsible for its rapid depletion [31,52,55]. After these considerations and adjustments, LUM was applied to the case study. The continuity equation was integrated explicitly to evaluate the variation of P(x) and E(x) along the northward direction x. The spatial evolution of the WVF vector exhibited a depletion of 55% after the transit over the mountain region. The WVF0 were around 540 After these considerations and adjustments, LUM was applied to the case study. The continuity equation was integrated explicitly to evaluate the variation of P(x) and E(x) along the northward direction x. The spatial evolution of the WVF vector exhibited a depletion of 55% after the transit over the mountain region. The WVF 0 were around 540 kg s −1 m −1 near Lecco, and its final value dropped down to 260 kg s −1 m −1 near the Spluga Pass, 75 km northward.
In Figure 7, the spatial evolution of P (blue) and E (red) cumulated over 13 h are reported. Comparing the graph with the elevation profile, it could be highlighted that LUM works in accordance with the physics of the orographic rainfall mechanism that led to precipitation formation in the windward portion of the mountain (upslope) and the evaporation in the lee side (downslope). This is particularly evident for Mount Legnone, which is the highest peak of the entire cross section (1500 m smoothed and 2650 m effective).
Geosciences 2021, 11, x FOR PEER REVIEW 12 of 18 kg s −1 m −1 near Lecco, and its final value dropped down to 260 kg s −1 m −1 near the Spluga Pass, 75 km northward. In Figure 7, the spatial evolution of P (blue) and E (red) cumulated over 13 h are reported. Comparing the graph with the elevation profile, it could be highlighted that LUM works in accordance with the physics of the orographic rainfall mechanism that led to precipitation formation in the windward portion of the mountain (upslope) and the evaporation in the lee side (downslope). This is particularly evident for Mount Legnone, which is the highest peak of the entire cross section (1500 m smoothed and 2650 m effective). The LUM results of P(x) and E(x) appear to be physically correct; however, we noticed that the rainfall amounts were under/overestimated at certain points and rather irregularly distributed. Therefore, to obtain a realistic representation of the rainfall pattern, a smoothing function was applied to P(x). We tested three different smoothing windows related to the meso-gamma scale of the studied phenomena: 5, 10, and 20 km ranges, as reported in Figure 8. We considered a smoothing window of 10 km that better fit the observed values. This also led us to believe that 10 km was the length scale of the thunderstorms that occurred in the region during the event [5,32,53]. The green continuous line in Figures 7 and 8 represents the smoothed function for the precipitation P(x): the values of the rainfall field are now more uniformly distributed with respect to the initial prediction ( Figure 7). The results have already considered the reduction ratio related to the precipitation efficiency (0.3). The LUM results of P(x) and E(x) appear to be physically correct; however, we noticed that the rainfall amounts were under/overestimated at certain points and rather irregularly distributed. Therefore, to obtain a realistic representation of the rainfall pattern, a smoothing function was applied to P(x). We tested three different smoothing windows related to the meso-gamma scale of the studied phenomena: 5, 10, and 20 km ranges, as reported in Figure 8. We considered a smoothing window of 10 km that better fit the observed values. This also led us to believe that 10 km was the length scale of the thunderstorms that occurred in the region during the event [5,32,53]. The green continuous line in Figures 7 and 8 represents the smoothed function for the precipitation P(x): the values of the rainfall field are now more uniformly distributed with respect to the initial prediction ( Figure 7). The results have already considered the reduction ratio related to the precipitation efficiency (0.3).
In Figure 9A,B, the rain-gauge amounts of the nearest weather station were plotted against the rainfall profile, showing a good accordance among the two datasets. The highest value recorded at the Premana rain gauge was well represented by LUM, indicating that in the surrounding area the upslope motion of the air likely conditioned the dynamic of the incoming WVF flux significantly. ranges, as reported in Figure 8. We considered a smoothing window of 10 km that better fit the observed values. This also led us to believe that 10 km was the length scale of the thunderstorms that occurred in the region during the event [5,32,53]. The green continuous line in Figures 7 and 8 represents the smoothed function for the precipitation P(x): the values of the rainfall field are now more uniformly distributed with respect to the initial prediction (Figure 7). The results have already considered the reduction ratio related to the precipitation efficiency (0.3). indicating that in the surrounding area the upslope motion of the air likely conditioned the dynamic of the incoming WVF flux significantly.
In Table 6, we report the error quantifications of the smoothed rainfall field with respect to the reference rain gauge stations located along the profile. The absolute error was generally confined around few millimetres of rain, with a maximum of 40 mm across the Valtellina flood plain. Looking at the relative errors, an average was assessed at nearly 20%, which is rather low considering the model adopted and its approximations. At the point of Premana, a small error (around 30 mm) with respect to the total rainfall amount was estimated.

Discussion
As we have seen from the results listed in Table 6, the errors were maintained under control with respect to the reference rain gauges, giving a realistic reconstruction of the rainfall field. In contrast to the classical interpolation techniques, we also gained a physical interpretation of the phenomena under twofold aspects. The first regards the analysis of the soundings data that can depict the atmosphere state before and during the event, highlighting potential critical situations, as in the present case. The second consists of a better understanding of orographic rainfall mechanism in relation to several parameters In Table 6, we report the error quantifications of the smoothed rainfall field with respect to the reference rain gauge stations located along the profile. The absolute error was generally confined around few millimetres of rain, with a maximum of 40 mm across the Valtellina flood plain. Looking at the relative errors, an average was assessed at nearly 20%, which is rather low considering the model adopted and its approximations. At the point of Premana, a small error (around 30 mm) with respect to the total rainfall amount was estimated.

Discussion
As we have seen from the results listed in Table 6, the errors were maintained under control with respect to the reference rain gauges, giving a realistic reconstruction of the rainfall field. In contrast to the classical interpolation techniques, we also gained a physical interpretation of the phenomena under twofold aspects. The first regards the analysis of the soundings data that can depict the atmosphere state before and during the event, highlighting potential critical situations, as in the present case. The second consists of a better understanding of orographic rainfall mechanism in relation to several parameters that can influence it. We can say that LUM can be more "didactic" compared with complex LAM models, and its simple implementation permits making the problem more tractable.
One of the strengths is the relatively small amount of input data required for the implementation. A lack of data is generally a problem in mountain catchments; however, as we have seen in our case study, the LUM was correctly initialized using only the sounding data and topographic elevation. Unfortunately, even if the terrain elevation is available worldwide with good resolution, sounding data are typically single-point data. In this particular case, the radio soundings were taken some 60 km away from the starting point of the simulation. This approximation was considered acceptable for defining the initial conditions of the model, i.e., the WVF 0 .
The estimation of WFV 0 , could be a non-trivial task at some locations or even not sufficiently precise to feed the linear model. It is more likely that the sounding is located too far away from the studied area, i.e., hundreds of kilometres. In that case, the atmosphere's vertical profile should be inferred through comparing different data sources, especially considering the vertical atmosphere profiles retrieved from reanalysis models or using satellites. The information detected by infrared sensors and LiDAR scanners from geostationary satellites are now providing encouraging results in terms of WVF estimation. These techniques are currently adopted for the study of atmospheric rivers [56,57]. As a matter of fact, increasing the data availability regarding atmospheric moisture flows may help in this way.
Another point of discussion concerns the hypothesis of linearity under which the model operates. In fact, trying to represent or simulate a chaotic system, such as the atmosphere, with a liner model could appear as nonsense. In our view, this is partially true. As described by the model authors, the model contains a strong simplification of the processes occurring in the atmosphere. Microphysics is completely neglected, and rainfall formation is reduced to a simple water mass balance of the WVF flow. This is clearly not true in reality, and the process is simulated by LAM. The complexity of LAMs sometimes does not permit effective control of the solution especially over complex terrain and when the orography abruptly influences the rainfall generation. LUM permitted us to highlight the critical parameters that affect the final outcome of the model, neglecting second-order influences. In the case study, the test of LUM moved on in this direction, accepting some simplification but sill presenting a reasonable rainfall field reconstruction. These simplifications are discussed here in detail: The mono-dimensional domain is a strong idealization of the event that occurred. The sounding data compared to a local LAM depicted an event that developed northward following a narrow cone extension. The thunderstorm corridor had a clear starting point, an average width around 10-15 km but to an extent of 100 km. Therefore, also considering the low-level wind convergence (Figure 3), a mono-dimensional reduction was sufficiently realistic. The resampling of orography is an operation that is generally adopted in atmospheric models. For the Alps, a rectangular shape range with a 2000 m average slope was considered in the past as a sufficient representation of the morphology in regional atmospheric-dynamic models [6,32,34,50]. Currently, LAMs can assimilate orography at higher resolutions, but the smoothing operations are still necessary [58,59]. In our case, the high sensitivity of LUM to the terrain profile required that the local morphology be smoothed to obtain realistic simulations. The topographic influence on incoming airflow can generate gravity waves and turbulences that can also perturb the airflow dynamic along the vertical [5,34,60]. These secondary effects likely played a significant role in the spatial redistribution of the rainfall, especially behind the peak of Mount Legnone where the results showed the highest errors with underestimation of around 40 mm. Using the linear model, the airmass uplift triggered by orography was the predominantly simulated process, and the others were confined to the second order. For these reasons, the downslope dynamics are poorly described due to high non-linearities that may occur in these processes. The estimation of the boundary layer height cannot be computed explicitly, and this represents a significant uncertainty that should be treated carefully. However, in our opinion, this quantity should be considered when LUM is adopted. As we determined for this case study, BL was essential to determine the portion of the atmosphere that can contribute to the effective rainfall generation. The motivation is the following: due to surface friction, low atmosphere layers are maintained at rest and do not experience any upslope motion until the BL is completely eroded. In our case study, the evidence was confirmed looking at three sounding data where, for layers comprised in BL, the wind velocities were sensibly reduced, and their directions were not aligned with the WFV airflow. If these layers are included in the computation of WFV 0 , they could lead to a sensible overestimation of the initial conditions due to their high concentration in water vapour. For these reasons, we tested this BL adjustment, and the results demonstrated good improvements in the LUM simulation.
Bearing in mind these necessary simplifications, the final result was rather encouraging. The model highlighted how the orography played a determinant role in intensifying the precipitation rate locally. As was observed in analysing the reconstructed precipitation field, the critical value recorded at Premana rain gauge was matched correctly. In the upper part of Valchiavenna, the model outcome was in good accordance with the rain gauge measurements. Likely, hidden errors were also contained in the reference rain gauges, particularly due to their relative distance from the cross section considered, such as in the case of Barzio station. For this reason, an extension of LUM for working with 2D domains could be a further improvement of the analysis to fulfil these uncertainties. We concluded that the linear upslope model produced a reliable reproduction of this exceptional rainfall event.

Conclusions
The aim of this paper was to reconstruct a realistic rainfall field through the application of a simple linear upslope model. Over complex terrain, this task is required in many fields, in particular for hydrogeological risk assessment. Interpolation techniques that consider only the data recorded by rain gauge networks may fail as they do not consider orography as critical factor for precipitation intensification. On the other hand, LAMs offer a good solution to investigate rainfalls; however, their complexity is highly data-demanding, and, especially in mountain environments, solutions could often be too approximate.
For these reasons, we attempted to apply the linear upslope model proposed by [37], which represents an intermediate alternative. To test the model application, a case study was chosen due to its severity, which was intensified by two main causes: the favourable meteorological conditions of the atmosphere and the interaction with the local orography. Considering the known limitations and the simplification adopted, LUM was determined to be suitable for the purpose.
The rainfall field reconstruction was able to fit the data obtained from the reference rain gauges, particularly around the critical areas where the highest amount of rainfall was recorded. In addition, this allowed us to increase our understanding of this type of complex meteorological phenomena, discovering the importance of BL for the assessment of the initial conditions related to the WVF vector. Consequently, the role of airflow interactions with orographic terrain was more clearly depicted in relation to the extreme localized rainfall, which was able to trigger several hydrogeological hazards. In this light, LUM represents a valuable instrument for a rapid but comprehensive analysis, Water Vapour Flux as demonstrated in this critical case study.