# Long-term foehn reconstruction combining unsupervised and supervised
learning

Reto Stauffer
Universität Innsbruck &Achim Zeileis
Universität Innsbruck &Georg J. Mayr
Universität Innsbruck  0000-0002-3798-5507 0000-0003-0918-3766
0000-0001-6661-9453<EMAIL_ADDRESS>
Reto Stauffer, Achim Zeileis, Georg J. Mayr Foehn winds, characterized by
abrupt temperature increases and wind speed changes, significantly impact
regions on the leeward side of mountain ranges, e.g., by spreading wildfires.
Understanding how foehn occurrences change under climate change is crucial.
Unfortunately, foehn cannot be measured directly but has to be inferred from
meteorological measurements employing suitable classification schemes. Hence,
this approach is typically limited to specific periods for which the necessary
data are available. We present a novel approach for reconstructing historical
foehn occurrences using a combination of unsupervised and supervised
probabilistic statistical learning methods. We utilize in-situ measurements
(available for recent decades) to train an unsupervised learner (finite
mixture model) for automatic foehn classification. These labeled data are then
linked to reanalysis data (covering longer periods) using a supervised learner
(lasso or boosting). This allows to reconstruct past foehn probabilities based
solely on reanalysis data. Applying this method to ERA5 reanalysis data for
six stations across Switzerland and Austria achieves accurate hourly
reconstructions of north and south foehn occurrence, respectively, dating back
to 1940. This paves the way for investigating how seasonal foehn patterns have
evolved over the past 83 years, providing valuable insights into climate
change impacts on these critical wind events. Foehn, reconstruction,
unsupervised, mixture model, supervised, climate, trend Reto Stauffer
Universität Innsbruck
Faculty of Economics and Statistics
_and_ Digital Science Center
Universitätsstr. 15
6020 Innsbruck, Austria
E-mail:

##  1 Introduction

Foehn winds are downslope winds on the leeward side of mountains and can be
found all around the world in areas with pronounced topographical features
that impede the airflow, such as the European Alps, the Southern Alps in New
Zealand, mountain ranges along the Mediterranean Sea, or the Rocky Mountains.
Depending on the region, these winds are given specific names such as Santa
Anna winds (Southern California), Chinook (Rocky Mountains), Mistral (Southern
France), Bora (Croatia), or Foehn (Central Europe and New Zealand). For
historical reasons ‘foehn’ has become a synonym for this type of terrain-
induced wind phenomena.

Often, foehn is characterized by a sharp increase in wind speed and sudden
changes in temperature and relative humidity, which can have a strong
influence on the local climate and the people living in the affected areas.
While foehn is often associated with a mild (and typically dry) climate,
strong foehn events can also cause extensive damage to vegetation and man-made
structures. In some areas it is not uncommon for strong foehn gusts to
overturn trucks or vans, or for airports and harbors to be closed due to
unsafe conditions. In addition, the strong and dry winds can kindle and spread
domestic fires and wild fires (Schoennagel _et al._ , 2004; Reinhard _et al._
, 2005; Zumbrunnen _et al._ , 2009) or affect the development of Antarctic ice
shelves (e.g., Cape _et al._ 2015; Elvidge _et al._ 2020).

Although the conceptual model of foehn is well established (Armi and Mayr,
2007; Mayr and Armi, 2008; Armi and Mayr, 2011; Richner and Hächler, 2013), it
cannot be measured directly. During the last decades several (semi)-automatic
methods have been developed which allow to differentiate ‘foehn’ and ‘no-
foehn’ events based on in-situ measurements from automated weather stations
(AWSs). Among the frequently used algorithms are Widmer’s föhn index (Widmer,
1966; Courvoisier and Gutermann, 1971) based on Fisher’s linear discriminant
analysis to distinguish between two or more distinct classes, and enhanced
versions of it (e.g., Jansing _et al._ 2022). Other studies use decision-based
or tree-based methods to classify foehn events (e.g., Dürr 2008; Speirs _et
al._ 2013; Cape _et al._ 2015; Turton _et al._ 2018; Datta _et al._ 2019;
Elvidge _et al._ 2020; Laffin _et al._ 2021; Francis _et al._ 2023; Laffin _et
al._ 2021; Francis _et al._ 2023), all of which are deterministic methods,
where the thresholds have often been selected manually. To overcome these
limitations, Plavcan _et al._ (2014) proposed a method based on finite mixture
models for automatic and fully parametric probabilistic foehn classification.
To perform the classification, all methods require AWS measurements with high
temporal resolution (ideally sub-hourly), which are typically only available
for recent decades. While this allows to classify and analyze foehn when the
AWS provides sufficient data, there is no information on foehn prior to the
installation of the AWS, nor when there have been outages or the AWS has been
decommissioned.

Additional information on the atmospheric conditions from numerical
(re-)analysis models are an excellent source to complement the in-situ
measurements. Reanalysis data sets are typically produced by physically based
numerical weather prediction (NWP) models and sophisticated data assimilation
schemes that use all available observations to estimate the ‘best known’
atmospheric state. An example is the global reanalysis data set ERA5 (Hersbach
_et al._ , 2023a, b) from the European Centre for Medium-Range Weather
Forecasts (ECMWF). ERA5 provides global hourly 4-dimensional atmospheric
conditions back to 1940. However, this comes with its own challenges. First,
foehn is not explicitly modeled by the NWP models. Second, due to technical
and computational limitations, the reanalysis can only approximate the real
world and cannot resolve small-scale atmospheric processes and topographic
features, which would be important for small-scale phenomena such as foehn.

One way to overcome these limitations is to combine AWS measurements and
reanalysis data using statistical or machine learning techniques. A
classification of foehn based on AWS data serves as the response
(target/outcome/labels) for a supervised model that uses reanalysis data as
explanatory variables (inputs/covariates). Once the relationship between the
two sets of data has been learned, the models can be used to predict the
expected state (‘foehn’ or ‘no foehn’) for periods for which reanalysis data
is available but AWS measurements are not. This is also known as (statistical)
downscaling or post-processing and has been used for foehn modeling in some
variations: For nowcasting foehn at a station in Switzerland (Altdorf),
Sprenger _et al._ (2017) use data from a local 7 km analysis data set
(COSMO-7) – but only to get consistent inputs for their adaptive boosting
algorithm (AdaBoost). Laffin _et al._ (2021) use both ERA5 and the Regional
Atmospheric Climate Model 2 (RACMO2; Wessen and Laffin 2022) combined with
tree-based gradient boosting (XGBoost; Chen and Guestrin 2016) to predict
foehn on the Antarctic Peninsula. XGBoost is also used by Mony _et al._ (2021)
to investigate future changes to foehn frequency in Switzerland.

This study proposes a novel probabilistic approach that combines unsupervised
and supervised machine learning methods to bridge the gap between in-situ
Automatic Weather Station (AWS) measurements and ERA5 reanalysis data to
diagnose foehn. This combined approach allows us to reconstruct long-term,
high-resolution foehn conditions over several decades. ERA5 data enables us to
reconstruct the probability of foehn occurrence with hourly temporal
resolution dating back to 1940, long before AWS were installed.

Our approach also allows to identify potential long-term changes in foehn
occurrence in the European Alps from this high-resolution reconstruction. The
validity of this combined approach is demonstrated by applying it to six
stations located both north and south of the main Alpine ridge to show the
method’s effectiveness for both north and south foehn wind situations.

Figure 1:  Flow chart for the new combined foehn reconstruction algorithm.

Figure 1 shows a schematic representation of the proposed approach. For all
six locations, data are available from an AWS at the target location as well
as from a nearby mountain station for the last 14–22 years (Sec. 2.1) on a 10
min temporal interval. A Gaussian mixture model (unsupervised learning; Sec.
3.1) is used for foehn classification. The result is then aggregated into
binary time series (‘foehn’/‘no foehn’) with an hourly temporal resolution to
match the resolution of the ERA5 data used in the next step. After combining
the different data sources, supervised learning (Sec. 3.2) is used to find the
relationship between a variety of interpolated and derived variables from ERA5
(Sec. 2.2) and the classified events. Once these statistical models have been
estimated, foehn can be reconstructed (Sec. 3.3) for the whole period from
1940–2022, allowing to investigate possible trends and/or seasonal changes
over the last decades (Sec. 3.4).

##  2 Data

Section 2.1 describes the measurement data utilized for foehn classification
along with the study area and the target stations. Section 2.2 explains the
reanalysis data set and its pre-processing for the supervised learning method,
along with the reconstruction process.

###  2.1 In-situ measurements

This study utilizes data from six AWSs situated across Switzerland and the
western part of Austria, all positioned at the bottom of valleys known to be
affected by foehn winds. Four of these stations are located north of the main
Alpine ridge, while two are located in the canton of Ticino (Switzerland)
south of the main Alpine ridge. Whilst the stations north of the Alps are
prone to south foehn, the stations south of the Alps are known for the
presence of north foehn.

An additional AWS upstream near the crest of the main Alpine range improves
the accuracy of the foehn classification (Plavcan _et al._ , 2014). The
stations Innsbruck and Ellbögen utilize data from Sattelberg, the remaining
four stations in Switzerland use observations from station Gütsch.

All stations provide data on mean wind speed, wind direction, air temperature,
and relative humidity at a 10 min temporal resolution. Table 1 displays the
locations and data availability of these stations, while Figure 2 depicts a
map illustrating their geographical position and the surrounding topography.

| Type | Location | Data availability
---|---|---|---
$\bigtriangleup$ Gütsch (Andermatt)1 | crest | 46.653N/8.616E 2286m | 2005-01-01–2023-12-31 (95.3%)
Altdorf1 | south | 46.890N/8.620E 438m | 2005-01-01–2022-12-31 (78.1%)
Montana1 | south | 46.290N/7.460E 1423m | 2005-01-01–2022-12-31 (77.1%)
Comprovasco1 | north | 46.460N/8.935E 576m | 2005-01-01–2022-12-31 (85.8%)
Lugano1 | north | 46.004N/8.960E 205m | 2005-01-01–2022-12-31 (90.2%)
$\bigtriangleup$ Sattelberg2 | crest | 47.011N/11.479E 2107m | 2006-01-01–2022-12-31 (75.0%)
Ellbögen2 | south | 47.200N/11.430E 1080m | 2006-01-01–2022-12-31 (92.0%)
(Universität) Innsbruck3 | south | 47.260N/11.385E 578m | 2009-06-21–2022-12-30 (99.1%)
Table 1: Station type, location and data availability; begin and end date plus
percent available within period. Observations provided by the Swiss national
weather service (1; MeteoSwiss), the University of Innsbruck (2) and the
Austrian national weather service (3; GeoSphere Austria). Four stations are
used to model south foehn, two to model north foehn and two serving
information at the mountain crest ($\bigtriangleup$; cf. Type).

###  2.2 ERA5 reanalysis

This study makes use of ERA5 reanalysis data, which is publicly accessible via
the Copernicus climate data store (Hersbach _et al._ , 2023a, b). ERA5 offers
four-dimensional gridded data with an hourly temporal resolution (starting
from 1940) on a spatial grid of $0.25^{\circ}\times 0.25^{\circ}$ ($\sim
28\,km\times 20\,km$ for Central Europe).

Figure 2:  Location of the six target stations (circles/C) and two crest
stations (triangles); Innsbruck and Ellbögen are shown combined due to their
close proximity ($<8\,km$). The solid black lines represents the main Alpine
ridge. On top, the neighborhood ‘star’ used for interpolation is shown,
exemplarily labeled on the most eastern station (details see App. A).

The data from 90 different fields (30 single level fields and 60 pressure
level fields; see App. A) are bilinearly interpolated to the geographical
location of the six stations of interest. Based on the 90 interpolated values,
a series of derived variables is calculated, including vertical temperature
gradients and level thickness, resulting in a total of 155 variables. These
155 variables are referred to as the ‘direct’ variable set as they solely rely
on information retrieved directly from the geographical location of the
corresponding target station.

Since foehn is driven by large-scale synoptic processes, it might be
insufficient to rely only on data at the target location. Therefore,
additional information from the surrounding is incorporated by extracting data
from ERA5 at a series of neighboring points in a ‘star’ formation around the
target location. In combination with the direct variables, a list of
additional derived variables is calculated, such as spatial temperature
differences and spatio-temporal pressure differences, among many others.
Finally, the first and second order harmonics of the day of the year are
included to capture seasonal variation. In total, this yields 497 variables:
four harmonics, 155 direct variables, 136 spatial variables, 120 temporal
variables and 82 spatio-temporal variables (e.g., temporal changes of spatial
differences). This expanded set is referred to as the ‘full’ variable set.

Figure 2 shows the target locations (`C`; center) and their neighboring points
used. These neighboring points are positioned geographically relative to the
target station upstream (`U`) and downstream (`D`) of the main foehn wind
direction as well as to the right (`R`) and left (`L`) of it. While the
interpolated information from the target location itself (`C`) is always used
as possible covariate for the statistical models, the values interpolated at
the neighboring points are not directly employed but are instead used for the
calculation of the derived/augmented variables. Further details regarding the
construction of the neighboring ‘star’ can be found in Appendix A.

##  3 Methodology

Section 3.1 introduces the unsupervised learning model used for foehn
classification, followed by the supervised learning models in Section 3.2. The
results from Section 3.2 are then used to reconstruct hourly foehn occurrence
over the past decades (Sec. 3.3) which is analyzed employing season-trend
decomposition in Section 3.4.

###  3.1 Unsupervised learning: Mixture model for foehn classification

As direct foehn measurement do not exist, the data from the AWSs are currently
unlabeled. Therefore, a mixture model is employed to distinguish between the
‘foehn’ and ‘no foehn’ events:

$\text{Pr}_{\text{obs}}(\text{foehn})=f_{\text{mix}}(\text{measurements}).$
(1)

$\text{Pr}_{\text{obs}}(\text{foehn})$ denotes the posterior probability for a
‘foehn’ observation at a specific time and station, which is modeled as a
function ($f_{\text{mix}}()$) of the 10 min in-situ measurements (Sec. 2.1).
We employ a two-component Gaussian mixture model (Grün and Leisch, 2008) with
concomitants, closely following the method proposed by Plavcan _et al._ (2014)
implemented in the _R_ package `foehnix` (Stauffer, 2023).

The prerequisite condition for an observation to be used for classification is
that the wind direction at the location of interest falls within the
prevailing foehn direction at the target location, and that wind from a
specific direction is also prevalent at the corresponding crest station at the
same time (details in Tab. 4). Only the periods matching this precondition are
used for estimating the Gaussian mixture model, while
$\text{Pr}_{\text{obs}}(\text{foehn})=0$ is set for all remaining
observations.

The underlying concept involves that two unobservable Gaussian components (or
clusters) exist, one describing ‘foehn’ conditions and the other ‘no foehn’
conditions. To distinguish between these two components, a main covariate is
required. In this study, the potential temperature difference ($\Delta\theta$)
between the valley station ($t_{\text{valley}}$) and the crest station
($t_{\text{crest}}$) is used, employing a simple dry adiabatic temperature
reduction $\Delta\theta=t_{\text{valley}}-t_{\text{crest}}-0.01\Delta h$ based
on the difference in altitude ($\Delta h$; cf.Tab. 1 & 4). During foehn
events, the air descends on the leeward side of the mountains, resulting in a
well-mixed atmosphere with neutral stratification, hence a potential
temperature difference close to zero.

However, the potential temperature difference alone might not be sufficient to
adequately separate the two states. Therefore, an additional concomitant model
is employed to weigh the two components conditional on additional covariates.
In this study, binary logistic regression is used for the concomitant model,
employing relative humidity and mean wind speed as additional covariates.
Models with this specification have been shown to work well for stations in
the Alpine region (e.g., Plavcan _et al._ 2014; Plavcan and Mayr 2015). Figure
3 provides an illustration of this model, depicting the use of $\Delta\theta$
to separate the two components (‘foehn’/‘no foehn’) and the effect of the
concomitant model on the joint density. More details can be found in Appendix
B.

Figure 3: Illustration of the two-component Gaussian mixture model used for
classification. The potential temperature difference ($\Delta\theta$) is used
as the main variable to separate component 1 (foehn) and component 2 (no
foehn); the concomitant model includes relative humidity (rh) and mean wind
speed (ff) at the valley station. Solid lines: density of the mixture model
once with a probability $p_{2}=0.9$ to observe component 2, once with
$p_{2}=0.3$.

Once estimated, the mixture model provides the posterior probability
$\text{Pr}_{\text{obs}}(\text{foehn})$ for each 10 min interval where the
required in-situ measurements are available (Eqn. 1). As these foehn
‘observations’ need to be combined with ERA5 in the next step, the results are
upscaled to an hourly temporal resolution ($\text{Pr}_{\text{1h}}$).
Therefore, the following basic assumption is applied: Each hour is considered
a ‘foehn’ event if at least half of the 10 min posterior probabilities within
that hour (e.g., 00:10–01:00 UTC for 01:00 UTC) are $\geq 0.5$, otherwise the
hour is considered a ‘no foehn’ event. If fewer than 4 out of 6 individual 10
min probabilities are available, the hour is excluded. This is similar to
Gutermann _et al._ (2012) and Mony _et al._ (2021) who employ a ‘four out of
six rule’.

$\text{Pr}_{\text{1h}}=\begin{cases}\text{missing}&\text{if}~{}~{}~{}\sum\text{Pr}_{\text{obs}}\in[0-1]<4\\\
\text{foehn}&\text{if}~{}~{}~{}\frac{1}{N}\sum\big{(}\text{Pr}_{\text{obs}}\geq
0.5\big{)}\geq 0.5\\\ \text{no foehn}&\text{else}\end{cases}$ (2)

###  3.2 Supervised learning: Modeling foehn probability

The outcome of the previous section is a binary time series, with each hour
labeled as either a ‘foehn’ or ‘no foehn’ observation (Eqn. 2). This serves as
the response variable in the supervised machine learning model, using the ERA5
data (Sec. 2.2) as input variables. The goal is to capture the relationship
between the two data sets with model of the form

$\text{Pr}_{\text{1h}}(\text{foehn})=f_{\text{learn}}(\text{ERA5}),$ (3)

where the probability $\text{Pr}_{\text{1h}}$ of the binary response is
modeled as a function $f_{\text{learn}}()$ of the available information
extracted from ERA5. $f_{\text{learn}}()$ can be any learner suitable for a
binary response such as logistic regression, decision trees, random forests,
or neural networks, to mention a few. For this study, three different
learners/models are employed:

lasso

Logistic regression with lasso (L1) regularization (Friedman _et al._ , 2010;
Tay _et al._ , 2023).

stabsel

Logistic regression with lasso-based stability selection (Meinshausen and
Bühlmann, 2010).

xgboost

Extreme gradient boosting (Chen and Guestrin, 2016).

In order to investigate the possible benefits of incorporating large-scale
information from the stations’ neighborhood, two variations of each learner
are considered: one utilizing the ‘full’ set of 497 variables and one only
using the ‘direct’ set of 155 variables (Sec. 2.2).

To account for location and time of day, separate models are estimated for
each of the six stations (Tab. 4) for each hour of the day ($0000$ UTC, $0100$
UTC, …, $2300$ UTC), resulting in a total of 864 models. Depending on the
station, the training data for these models include 10–18 years of data (see
Tab. 1). In addition, a six-fold cross-validation (CV) is performed using a
fixed period of 12 years (2011–2022) where, in each fold, two consecutive
years are left out as test data. Ellbögen and Innsbruck are missing one fold
(with test data 2013–2014) where no measurements from the crest station are
available and thus the classification is not possible.

###  3.3 Long-term foehn reconstruction

Once the models from the previous section are estimated, they can be applied
to the entire ERA5 period available. Although this is a prediction from a
statistical perspective, it is termed ‘reconstruction’ in this article, as
these predictions are applied backwards in time. The result is an hourly time
series of foehn probabilities $\widehat{\text{Pr}}_{\text{1h}}$ from January
1, 1940 to December 31, 2022 (83 years).

This reconstruction can fill possible gaps in historic records during the
lifetime of an AWS and extend ‘foehn observations’ to times prior to or beyond
the installation period of the AWS. The results may also serve as valuable
input for investigating the effects of foehn, e.g., in ecology where the
warming and drying effects of frequent foehn events can significantly impact
flora and fauna or favor the occurrence of forest fires.

Moreover, the reconstruction offers the possibility to analyze foehn
occurrence from a climatological perspective: Did the occurrence of foehn
increase/decrease along with the changing climate over the decades? Are there
changes in the seasonal or diurnal patterns? These questions are investigated
in more detail in the next section.

###  3.4 Season-trend decomposition

The comprehensive reconstructed data set allows to study foehn occurrence in a
climatological context. For this analysis, the hourly probability (Eqn. 3) is
aggregated by (i) taking the highest probability
$\widehat{\text{Pr}}_{\text{1h}}$ per day (0000 UTC–0000 UTC), before (ii)
calculating monthly averages. The resulting time series contains “monthly
means of the daily maxima”, which are then modeled using a season-trend
decomposition.

Due to the nature of the data there is a large year-by-year but also within-
year variability depending on the prevailing weather situation. To decompose
the signal, a season-trend decomposition is employed separating the signal
into long-term changes and a remainder component containing the residual
variability.

In this study, the regression-based decomposition of Dokumentov and Hyndman
(2022) is used, which also provides confidence intervals for the estimated
season and trend components. The model for the monthly mean foehn
probabilities $y_{t}$ assumes an additive decomposition into a smoothly
changing long-term trend $T_{t}$, a smoothly changing seasonal component
$S_{t}^{(m)}$, and a remainder $R_{t}$:

$y_{t}=T_{t}+\sum_{m=1}^{12}S_{t}^{(m)}+R_{t}~{}~{}\text{with}~{}~{}t\in
1,\dots,J,$ (4)

where $m\in\\{1,12\\}$ is the seasonal frequency (i.e., monthly) and $J$ is
the sample size (83 years $\times$ 12 months = 996). The model is estimated
via the _R_ package `stR` (Dokumentov and Hyndman, 2023).

##  4 Results

First, this section investigates which insights can be gained from the
reconstruction about the foehn occurrence at the six different target
stations. Different temporal scales are considered for this, namely: Inter-
annual changes in the foehn probabilities in Section 4.1, long-term trends and
seasonal patterns in Section 4.2, as well as changes in the diurnal patterns
across decades in Section 4.3. All of these results are based on the
reconstruction using the ‘lasso’ learner with the ‘full’ covariate set for the
full time period (without cross-validation).

Second, the performance of the supervised learning model is assessed under
different model specifications, namely: Using the ‘full’ set of all 497
variables vs. the 155 ‘direct’ variables only in Section 4.4 and comparing the
performance of the three supervised learners (lasso, stability selection,
extreme gradient boosting) in Section 4.5. All of these results are based on
out-of-sample Brier scores obtained in a six-fold cross-validation.

###  4.1 Temporal changes: Average annual foehn probabilities

The primary outcome of this study is the complete reconstruction of hourly
foehn probabilities over $83$ years (see Sections 3.1–3.3), yielding time
series with $N\sim 7.27\cdot 10^{5}$ observations. To carve out inter-annual
variations in the foehn probabilities, we aggregate the reconstructed data by
taking the daily maximum of $\widehat{\text{Pr}}_{\text{1h}}$ (Eqn. 3) before
calculating annual means. This can be interpreted as the average probability
of observing a foehn event on any given day within that year.

Figure 4 contains the result for all six stations, with Ellbögen exhibiting
the highest mean annual probabilities (on average $0.334$), while Altdorf and
Innsbruck show the lowest (on average $0.130$ and $0.133$ respectively).
Additionally, the annual mean of daily maxima from the classification is shown
for years with at least 80% of measurements available at the AWSs. The results
show an overall good agreement between the two signals from the reconstruction
and from the classification, with some larger gaps due to data availability as
well as some noticeable differences for specific stations in particular years.

Figure 4:  Annual mean of highest daily foehn probability (solid lines) for
all six stations from 1940 to 2022 based on the ‘lasso’ model with ‘full’
covariate set. Additionally, annual means of daily maxima from the foehn
classification using AWS data (dashed lines) and the availability of in-situ
measurements (straight bands at the bottom) are shown. Appendix D shows a
comparison to all other models and variants.

The reconstruction reveals a pronounced inter-annual variability, with certain
years exhibiting a much higher annual mean than the long-term average, whilst
others distinctly fall below. This variability is not random, as one can see
similar patterns among the four south-foehn stations. For instance, all four
stations show unusually high mean probabilities for 1951 and 1972. Similarly,
the two north-foehn stations exhibit a similar temporal behavior over time.

Moreover, the figure suggests a possible increase in the annual mean foehn
probability for south-foehn stations between 1940 and 1980. Hence this
question is investigated in more detail in the next section.

###  4.2 Temporal changes: Climatological trends and seasonal patterns

In this section, the analysis from the previous section is taken a step
further. Rather than focusing on the inter-annual variation the goal is to
bring out the long-term climatological trends and changes in the seasonal
patterns. Hence the reconstructed hourly time series are again aggregated but
to monthly (rather than annual) means of the daily maxima of
$\widehat{\text{Pr}}_{\text{1h}}$. Based on the season-trend decomposition
outlined in Section 3 Figure 5(a) illustrates the resulting smooth trends
($T_{t}$) along with the corresponding $95\%$ confidence intervals. Figure
5(b) depicts the corresponding smoothly varying seasonal signals
($S_{t}^{(m)}$) averaged over decades for visual purposes (where decade 1940
corresponds to the years 1940–1949 etc.).

The estimated trends (Fig. 5a) show a significant increase for the two
stations in Western Austria (Innsbruck, Ellbögen) between 1940 and 1980, and a
plateau thereafter. All other stations show a linear increase over the study
period. For four of the six stations (all except Montana and Comprovasco) the
smooth trend differs significantly from a constant.

Figure 5(b) shows the analysis of the seasonal changes which reveals the
different characteristics between north-foehn stations and south-foehn
stations. The two stations located south of the main Alpine ridge, Comprovasco
and Lugano, show one distinct maximum in spring and a minimum during autumn.
This pattern is stable over the entire study period and no changes in the
seasonal pattern ($S_{t}^{(m)}$) are found. The picture looks different when
focusing on the four south-foehn stations which all show two maxima in spring
and autumn with lower probability of foehn occurrence during summer and
winter. Although not significant, the season-trend decomposition indicates an
increase in the probability of foehn in spring (April, Mai) as well as in
autumn (October, November) with a slight decrease in late summer (August,
September).

Figure 5:  Results of the season-trend decomposition. (a) Estimated long-term
trends with estimated 95% confidence interval (blue) on top of annual mean
probabilities (black), the horizontal dashed line shows the mean of the trend
for better orientation. (b) Estimated seasonal pattern averaged over the whole
period 1940–2022 (dashed) and separately by decade (solid, label indicates
first year of each decade).

###  4.3 Temporal changes: Intra-daily patterns

For certain applications, information about intra-daily patterns and their
changes over time can be of great interest. With the hourly temporal
resolution of the reconstruction, such insights are now possible across
several decades.

Figure 6 shows Hovmøller diagrams for Ellbögen, depicting the decadal mean
probability per time of day and month. Despite pronounced variability between
the decades, the plot supports the previous findings showing an overall
increase over time with strongest increase in spring (April, May) and in
autumn (October, November). In addition, this visualization gives insights
into the diurnal pattern. Generally, foehn occurrence in Ellbögen is more
likely during the day (around 1000 UTC–2200 UTC) in spring and autumn, the
period where indication for a certain increase was found (Sec. 4.2). The
minimum average foehn probability is in the early morning and tied to the
length of the night. The time of the minimum varies somewhat interdecadally.

Figure 6:  Hovmøller diagrams showing the average foehn probability for
Ellbögen over the course of the day (y-axis) for all months (x-axis) for the
decades 1940–2019.

###  4.4 Model performance: Benefit of full covariate set

As described in Section 3.2, two variants of the supervised learning methods
are estimated: one using only the 155 ‘direct’ variables and one using the 497
‘full’ variables including large-scale atmospheric conditions (Sec. 2.2). For
both variants, a six-fold CV is performed (Sec. 3.2).

To investigate the benefit of the ‘full’ covariate set, Figure 7 shows the
Brier scores (BSs) for the ‘lasso’ model for all six stations. For each
station and each variant, BSs are shown for the test data set (out-of-sample)
as well as for the training data set (in-sample). This shows that the models
based on the ‘full’ variable set clearly outperform the models based on the
‘direct’ variables only. Although the overall performance of the less complex
models based on the ‘direct’ variable set is still decent, including the
additional large-scale spatio-temporal information substantially improves the
overall model performance.

Figure 7:  Comparison of ‘lasso’ model performance using the ‘direct’ variable
set (blueish) versus the larger ‘full’ variable set (orange with gray filling)
for all stations. Shown is the average Brier score from the six-fold CV for
the test (out-of-sample) and training (ins-sample) period. Lower is better.

In addition to the predictive skill, Figure 7 also shows the stability of the
model with largest variance in the BSs visible on the test data sets in
Ellbögen. On the training data sets the scores barely vary due to the large
sample size (10 years, hourly data).

###  4.5 Model performance: Comparison of supervised learners

While the previous section focuses on the benefits of using more input data,
this section compares the three different supervised learners described in
Section 3.2. For simplicity, only the results for models based on the ‘full’
variable set are shown as they have been shown to outperform those only using
the ‘direct’ variable set (Sec. 4.4).

Figure 8 illustrates that the BSs from the six-fold CV on the test set (out-
of-sample) is comparable for all three supervised learners with only minor
differences. The average BS is slightly lower for ‘lasso’ ($0.0261$), followed
by ‘stabsel’ ($0.0276$), and ‘xgboost’ ($0.0283$).

Figure 8:  Comparison of all supervised learning models using the ‘full’
variable set. Brier scores for training (in-sample; green) and test (out-of-
sample; violet) data set based on six-fold cross-validation. Lower is better.

On the training data (in-sample) the picture is similar for ‘lasso’ and
‘stabsel’ but ‘xgboost’ has much lower BSs. This indicates that ‘xgboost’
might be subject to some overfitting, despite careful tuning of the
hyperparameters (App. C.3).

##  5 Discussion and outlook

Using a novel combination of unsupervised and supervised learning we are able
to accurately reconstruct long-term foehn time series (starting from 1940) at
hourly resolution. More specifically, foehn classification is first
accomplished by (unsupervised) Gaussian finite mixture models based on AWS
measurements and then linked to ERA5 data using binary supervised learners
such as lasso, stability selection, or extreme gradient boosting. The
resulting foehn reconstruction enables novel analyses, exemplified here by
investigating long-term changes in trends, seasonal patterns, and diurnal
cycles of foehn occurrence.

The season-trend decomposition based on the period 1940–2022 reveals that all
six stations considered have either experienced a linear increase in foehn
occurrence (probability) over the entire study period or an increase between
1940 and the early 1980s that leveled off afterwards. For four of the six
stations it can be shown that this increase is significant. Although no
significant change in the seasonality is found, the results for all south-
foehn stations (Altdorf, Montana, Ellbögen, and Innsbruck) indicate a slight
increase in the occurrence of foehn in spring and autumn, with a slight
decrease in late summer.

The high quality of the foehn reconstruction is partially due to using a large
set of predictor covariates that not only contains information at the target
location but also includes additional large-scale atmospheric information from
the stations’ surrounding. The benefits of this full set of covariates are
similar for all three supervised learners considered: logistic regression with
lasso regularization (‘lasso’), logistic-regression-based stability selection
(‘stabsel’), and extreme gradient boosting (‘xgboost’). Lasso performs best in
our application, closely followed by the other two learners. Some further
improvements might be gained for xgboost if overfitting on the training data
can be further reduced, e.g., by a different hyperparameter tuning strategy.

For a comparison to existing publications, we complement the Brier scores
shown in the main paper by other popular scores for binary outcomes, such as
the false negative rate (FNR, also known as miss rate), false positive rate
(FPR, also known as false alarm rate), and percent correct (PC, also known as
accuracy). Based on the best model (lasso with full covariate set), we obtain
the following performances for Altdorf: $15.7\%$/$0.4\%$/$98.8\%$
(FNR/FPR/PC). These align well with the existing literature and stand out for
an exceptionally low FPR. Sprenger _et al._ (2017) report
$11.8\%$/$33.8\%$/$96.5\%$ and Mony _et al._ (2021) report
$21.4\%$/$21.4\%$/$97.7\%$. Similarly for Lugano we obtain
$15.3\%$/$0.7\%$/$98.2\%$, while Mony _et al._ (2021) report
$22.1\%$/$22.1\%$/$97.1\%$. More details are included in the appendix (App.
E).

The same holds for the foehn classification when compared to existing
literature and the Swiss foehn index (SFI) operationally used at MeteoSwiss in
terms of “average foehn hours per year”. The results from the Gaussian mixture
model (Sec. 3.1) for Altdorf show an average of 482.4 h/year which aligns well
with the SFI (458.4 h/year), as well as the results reported by Dürr (2008)
(478 h/year), Jansing _et al._ (2022) (465.8 h/year), and MeteoSwiss (2024)
(477 h/year). Montana exhibits 1007.4 h/year (SFI 904.7 h/year111Average based
on the years 2009–2016 only.), while Lugano and Comprovasco show 644.7 h/year
(SFI 563.0 h/year, MeteoSwiss (2024) 551 h/year) and 1077.4 h/year (SFI 953.4
h/year), respectively.

It should be emphasized that the excellent performance of the combined
approach depends crucially on the high quality of the automatic foehn
classification in the Gaussian finite mixture model. Thus, if the separation
of the AWS measurements into ‘foehn’ or ‘no foehn’ components works well (as
for the six stations presented), then different binary supervised learners are
able to link this reliably to the ERA5 data. However, there may also be
stations where a distinction into just two components might not be sufficient.
For example, we have found this to be the case for Aigle, Switzerland, where a
three-component mixture model appears to be necessary to separate light down-
valley winds, strong humid (catabatic) outflows, and actual foehn situations
(details not shown).

## Computational details

The results in this paper were obtained using R 4.2+. The majority of data
preparation and handling is done using the R packages stars 0.6.4, sf 1.0.15
and zoo 1.8.12. foehnix 0.1.6 is used for foehn classification, the supervised
learning is based on glmnet 4.1.7 and xgboost 1.7.5.1. The season-trend
decomposition is based on the R package stR 0.6.

## Acknowledgments

The study is partly based on the preliminary work of Morgenstern (2020). The
computational results presented here have been achieved (in part) using the
LEO HPC infrastructure of the University of Innsbruck.

## References

  * Armi and Mayr (2007) Armi L, Mayr GJ (2007).  “Continuously Stratified Flows Across an Alpine Crest with a Pass: Shallow and Deep Föhn.”  _Quarterly Journal of the Royal Meteorological Society_ , 133(623), 459–477.  10.1002/qj.22.
  * Armi and Mayr (2011) Armi L, Mayr GJ (2011).  “The Descending Stratified Flow and Internal Hydraulic Jump in the Lee of the Sierras.”  _Journal of Applied Meteorology and Climatology_ , 50(10), 1995–2011.  10.1175/JAMC-D-10-05005.1.
  * Cape _et al._ (2015) Cape MR, Vernet M, Skvarca P, Marinsek S, Scambos T, Domack E (2015).  “Foehn Winds Link Climate-Driven Warming to Ice Shelf Evolution in Antarctica.”  _Journal of Geophysical Research: Atmospheres_ , 120(21), 11037–11057.  10.1002/2015JD023465.
  * Chen and Guestrin (2016) Chen T, Guestrin C (2016).  “XGBoost: A Scalable Tree Boosting System.”  In _Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , KDD ’16, pp. 785–794. Association for Computing Machinery, New York, NY, USA.  ISBN 9781450342322.  10.1145/2939672.2939785.
  * Chen _et al._ (2023) Chen T, He T, Benesty M, Khotilovich V, Tang Y, Cho H, Chen K, Mitchell R, Cano I, Zhou T, Li M, Xie J, Lin M, Geng Y, Li Y, Yuan J (2023).  _xgboost: Extreme Gradient Boosting_.  R package version 1.7.5.1, URL https://CRAN.R-project.org/package=xgboost.
  * Courvoisier and Gutermann (1971) Courvoisier H, Gutermann T (1971).  “Zur praktischen Anwendung des Föhntests von Widmer.”  21, 7pp.  Accessed 2024-03-21, URL http://www.agfoehn.org/doc/Courvoisier_1971.pdf.
  * Datta _et al._ (2019) Datta RT, Tedesco M, Fettweis X, Agosta C, Lhermitte S, Lenaerts JTM, Wever N (2019).  “The Effect of Foehn-Induced Surface Melt on Firn Evolution over the Northeast Antarctic Peninsula.”  _Geophysical Research Letters_ , 46(7), 3822–3831.  10.1029/2018GL080845.
  * Dokumentov and Hyndman (2022) Dokumentov A, Hyndman RJ (2022).  “STR: Seasonal-Trend Decomposition Using Regression.”  _INFORMS Journal on Data Science_ , 1(1), 50–62.  10.1287/ijds.2021.0004.
  * Dokumentov and Hyndman (2023) Dokumentov A, Hyndman RJ (2023).  _stR: STR Decomposition_.  R package version 0.6, URL https://CRAN.R-project.org/package=stR.
  * Dürr (2008) Dürr B (2008).  “Automatisiertes Verfahren zur Bestimmung von Föhn in Alpentälern.”  _Arbeitsbericht_ , MeteoSchweiz.  Accessed 2024-03-21, URL https://www.meteoschweiz.admin.ch/dam/jcr:1fbe0d4e-fc7c-4121-b89e-1429a931272f/ab223.pdf.
  * Elvidge _et al._ (2020) Elvidge AD, Kuipers Munneke P, King JC, Renfrew IA, Gilbert E (2020).  “Atmospheric Drivers of Melt on Larsen C Ice Shelf: Surface Energy Budget Regimes and the Impact of Foehn.”  _Journal of Geophysical Research: Atmospheres_ , 125(17), e2020JD032463.  10.1029/2020JD032463.
  * Francis _et al._ (2023) Francis D, Fonseca R, Mattingly KS, Lhermitte S, Walker C (2023).  “Foehn Winds at Pine Island Glacier and Their Role in Ice Changes.”  _The Cryosphere_ , 17(7), 3041–3062.  10.5194/tc-17-3041-2023.
  * Friedman _et al._ (2010) Friedman JH, Hastie T, Tibshirani R (2010).  “Regularization Paths for Generalized Linear Models via Coordinate Descent.”  _Journal of Statistical Software_ , 33(1), 1–22.  10.18637/jss.v033.i01.
  * Grün and Leisch (2008) Grün B, Leisch F (2008).  “FlexMix Version 2: Finite Mixtures with Concomitant Variables and Varying and Constant Parameters.”  _Journal of Statistical Software_ , 28(4), 1–35.  10.18637/jss.v028.i04.
  * Gutermann _et al._ (2012) Gutermann T, Dürr B, Richner H, Bader S (2012).  “Föhnklimatologie Altdorf: Die lange Reihe (1864–2008) und ihre Weiterführung, Vergleich mit anderen Stationen.”  _Fachbericht MeteoSchweiz_ , 241, 1–53.  Accessed 2024-03-21, URL https://www.meteoswiss.admin.ch/services-and-publications/publications/reports-and-bulletins/2012/foehnklimatologie-altdorf--die-lange-reihe--1864-2008--und-ihre-.html.
  * Hersbach _et al._ (2023a) Hersbach H, Bell B, Berrisford P, Biavati G, Horányi A, Muñoz Sabater J, Nicolas J, Peubey C, Radu R, Rozum I, Schepers D, Simmons A, Soci C, Dee D, Thépaut JN (2023a).  “ERA5 Hourly Data on Pressure Levels from 1940 to Present.”  10.24381/cds.bd0915c6.  Accessed 2024-03-21.
  * Hersbach _et al._ (2023b) Hersbach H, Bell B, Berrisford P, Biavati G, Horányi A, Muñoz Sabater J, Nicolas J, Peubey C, Radu R, Rozum I, Schepers D, Simmons A, Soci C, Dee D, Thépaut JN (2023b).  “ERA5 Hourly Data on Single Levels from 1940 to Present.”  10.24381/cds.adbb2d47.  Accessed 2024-03-21.
  * Jansing _et al._ (2022) Jansing L, Papritz L, Dürr B, Gerstgrasser D, Sprenger M (2022).  “Classification of Alpine South Foehn based on 5 Years of Kilometre-Scale Analysis Data.”  _Weather and Climate Dynamics_ , 3(3), 1113–1138.  10.5194/wcd-3-1113-2022.
  * Laffin _et al._ (2021) Laffin MK, Zender CS, Singh S, Van Wessem JM, Smeets CJPP, Reijmer CH (2021).  “Climatology and Evolution of the Antarctic Aeninsula Föhn Wind-Induced Melt Regime From 1979–2018.”  _Journal of Geophysical Research: Atmospheres_ , 126(4), 1–19.  10.1029/2020JD033682.
  * Mayr and Armi (2008) Mayr GJ, Armi L (2008).  “Föhn as a Response to Changing Upstream and Downstream Air Masses.”  _Quarterly Journal of the Royal Meteorological Society_ , 134(635), 1357–1369.  10.1002/qj.295.
  * Meinshausen and Bühlmann (2010) Meinshausen N, Bühlmann P (2010).  “Stability Selection.”  _Journal of the Royal Statistical Society B_ , 72(4), 417–473.  10.1111/j.1467-9868.2010.00740.x.
  * MeteoSwiss (2024) MeteoSwiss (2024).  “Wetter und Klima A bis Z: Föhnhäufigkeit.”  Accessed 2024-03-21, URL https://www.meteoschweiz.admin.ch/wetter/wetter-und-klima-von-a-bis-z/foehnhaeufigkeit.html.
  * Mony _et al._ (2021) Mony C, Jansing L, Sprenger M (2021).  “Evaluating Foehn Occurrence in a Changing Climate Based on Reanalysis and Climate Model Data Using Machine Learning.”  _Weather and Forecasting_ , 36(6), 2039–2055.  10.1175/WAF-D-21-0036.1.
  * Morgenstern (2020) Morgenstern DS (2020).  _Multidecadal Foehn Time Series Reconstruction Using Machine Learning and ERA5 Reanalysis Data_.  Master’s thesis, Universität Innsbruck.  URL https://bibsearch.uibk.ac.at/AC15677542.
  * Plavcan and Mayr (2015) Plavcan D, Mayr GJ (2015).  “Towards an Alpine Foehn Climatology.”  33rd International Conference on Alpine Meteorology, Innsbruck. Accessed 2024-03-21, URL https://www.uibk.ac.at/congress/icam2015/abstracts_oral_presentations.htm#O14.4.
  * Plavcan _et al._ (2014) Plavcan D, Mayr GJ, Zeileis A (2014).  “Automatic and Probabilistic Foehn Diagnosis with a Statistical Mixture Model.”  _Journal of Applied Meteorology and Climatology_ , 53(3), 652–659.  10.1175/JAMC-D-13-0267.1.
  * Reinhard _et al._ (2005) Reinhard M, Rebetez M, Schlaepfer R (2005).  “Recent Climate Change: Rethinking Drought in the Context of Forest Fire Research in Ticino, South of Switzerland.”  _Theoretical and Applied Climatology_ , 82(1), 17–25.  10.1007/s00704-005-0123-6.
  * Richner and Hächler (2013) Richner H, Hächler P (2013).  _Understanding and Forecasting Alpine Foehn_ , chapter 4, pp. 219–260.  Springer-Verlag, Dordrecht.  ISBN 978-94-007-4098-3.  10.1007/978-94-007-4098-3_4.
  * Schoennagel _et al._ (2004) Schoennagel T, Veblen TT, Romme WH (2004).  “The Interaction of Fire, Fuels, and Climate across Rocky Mountain Forests.”  _BioScience_ , 54(7), 661–676.  10.1641/0006-3568(2004)054[0661:TIOFFA]2.0.CO;2.
  * Speirs _et al._ (2013) Speirs JC, McGowan HA, Steinhoff DF, Bromwich DH (2013).  “Regional Climate Variability Driven by Foehn Winds in the McMurdo Dry Valleys, Antarctica.”  _International Journal of Climatology_ , 33(4), 945–958.  10.1002/joc.3481.
  * Sprenger _et al._ (2017) Sprenger M, Schemm S, Oechslin R, Jenkner J (2017).  “Nowcasting Foehn Wind Events Using the AdaBoost Machine Learning Algorithm.”  _Weather and Forecasting_ , 32(3), 1079–1099.  10.1175/WAF-D-16-0208.1.
  * Stauffer (2023) Stauffer R (2023).  “foehnix: A Toolbox for Automated Foehn Classification based on Mixture Models.”  R package version 0.1.6, URL https://retostauffer.github.io/Rfoehnix/.
  * Tay _et al._ (2023) Tay JK, Narasimhan B, Hastie T (2023).  “Elastic Net Regularization Paths for All Generalized Linear Models.”  _Journal of Statistical Software_ , 106(1), 1–31.  10.18637/jss.v106.i01.
  * Turton _et al._ (2018) Turton JV, Kirchgaessner A, Ross AN, King JC (2018).  “The Spatial Distribution and Temporal Variability of Föhn Winds over the Larsen C Ice Shelf, Antarctica.”  _Quarterly Journal of the Royal Meteorological Society_ , 144(713), 1169–1178.  10.1002/qj.3284.
  * Wessen and Laffin (2022) Wessen JMV, Laffin MK (2022).  “Regional Atmospheric Climate Model 2 (RACMO2), version 2.3p2.”  10.5281/zenodo.3677641.  Accessed 2023-08-23.
  * Widmer (1966) Widmer R (1966).  “Statistische Untersuchungen über den Föhn im Reusstal und Versuch einer objektiven Föhnprognose für die Station Altdorf.”  _Vierteljahresschrift der Naturforschenden Gesellschaft Zürich_ , 111, 331–375.  Accessed 2024-03-21, URL https://www.ngzh.ch/archiv/1966_111/111_3-4/111_20.pdf.
  * Zumbrunnen _et al._ (2009) Zumbrunnen T, Bugmann H, Conedera M, Bürgi M (2009).  “Linking Forest Fire Regimes and Climate – A Historical Analysis in a Dry Inner Alpine Valley.”  _Ecosystems_ , 12(1), 73–86.  10.1007/s10021-008-9207-3.

##  Appendix A ERA5 fields

As described in Section 2.2 information from 90 different ERA5 fields is
bilinearely interpolated to (i) the target location of the different stations
used (cf. 2.1) as well as to a series of neighboring points arranged in a
‘star’ formation to include large-scale information. The fields used are
listed in Table 2 and 3.

CDS ERA5 fields name
---
`100m_u_component_of_wind` | `100m_v_component_of_wind`
`10m_u_component_of_wind` | `10m_v_component_of_wind`
`10m_wind_gust_since_previous_post_processing` | `instantaneous_10m_wind_gust`
`2m_dewpoint_temperature` | `2m_temperature`
`surface_pressure` | `mean_sea_level_pressure`
`eastward_gravity_wave_surface_stress` | `northward_gravity_wave_surface_stress`
`convective_precipitation` | `large_scale_precipitation`
`total_precipitation` | `boundary_layer_dissipation`
`boundary_layer_height` | `charnock`
`friction_velocity` | `gravity_wave_dissipation`
`high_cloud_cover` | `low_cloud_cover`
`medium_cloud_cover` | `total_cloud_cover`
`total_column_cloud_liquid_water` | `surface_net_solar_radiation`
`surface_net_thermal_radiation` | `surface_sensible_heat_flux`
`surface_solar_radiation_downwards` | `surface_thermal_radiation_downwards`
Table 2: List of ERA5 single level field retrieved. CDS ERA5 fields name
---
`divergence` | `geopotential`
`potential_vorticity` | `specific_humidity`
`temperature` | `u_component_of_wind`
`v_component_of_wind` | `vertical_velocity`
`specific_cloud_liquid_water_content` | `vorticity`
Table 3: List of ERA5 pressure level fields retrieved; each on the standard
pressure levels 500 hPa, 700 hPa, 750 hPa, 800 hPa, 850 hPa and 900 hPa.

To keep the approach generic, these neighbouring points are ‘arbitrary’ points
not linked to stations or cities, but geometrically placed around the target
station relative to its location to the main mountain range. Figure 2 (main
article) shows these ‘stars’ which are constructed as follows: First, the
orthogonal from the center location (station location; `C`) to the main
mountain range is computed. Along this orthogonal axis, two additional points
upwind of the foehn wind direction (`U`, `UU`) and two points in the downwind
direction (`D`, `DD`) are defined with a radius of 1∘ and 2∘ from the center
point. Second, four points with a radius of 2∘ from the center point deviating
by 45∘ from the orthogonal are added, defining points on the right hand side
(`*R`) and left hand side (`*L`) of the main foehn direction in the upwind
direction (`UR`, `UL`) and downwind direction (`DR`, `DL`). The downwind
direction for stations located north of the main alpine range is towards north
(south-foehn stations), for those located south of the main alpine range the
downwind direction is toward south (north-foehn stations).

Based on these points, the ‘direct’ variable set (155 variables) as well as
the ‘full’ variable set (497 variables) are calculated as described in Section
2.2 (main article). The extended (‘full’) variable set contains a series of
spatial and spatio-temporal information such as:

  * •

surface pressure differences between `C` to `U`, `U` to `C`, `UL` to `DR` and
others,

  * •

potential temperature differences between `C` to `U`, `U` to `C`, `UL` to `DR`
and others,

  * •

differences in vertical temperature differences between `C` to `U`, `U` to
`C`, `UL` to `DR` and others,

  * •

sum of precipitation on the upwind side (sum of `U`, `UU`, `UR`, `UL`) as well
as on the downwind side (sum of `D`, `DD`, `DR`, `DL`) as well and the
difference between these two sums,

  * •

mean cloud cover on the upwind side (sum of `U`, `UU`, `UR`, `UL`) as well as
on the downwind side (sum of `D`, `DD`, `DR`, `DL`) as well and the difference
between these two means.

Moreover, temporal variables are calculated including e.g., changes in the
geopotential height, temperature, or changes in humidity on different levels
over the past 3 hours or over the next 3 hours and variations thereof.

##  Appendix B Unsupervised learning: Two-component mixture model

For all eight AWSs used (Tab. 1; main article) observations are available for
the most recent 10–18 years on a 10 min temporal resolution. Since foehn winds
show a characteristic wind direction, only observations meeting a specific
precondition are used for classification as described in Section 3.1. Table 4
shows the wind sectors used for each target site, as well as the amount of
data used for classification.

|  Stations, conditional wind direction, data used
---|---
Altdorf |  Valley station: Altdorf (60∘–240∘)
Crest station: Gütsch, Andermatt (105∘–285∘)
20.0 % within sector, 57.1 % outside sector, 22.9 % removed
Innsbruck |  Valley station: Universität Innsbruck (130∘–230∘)
Crest station: Sattelberg (90∘–270∘)
5.8 % within sector, 49.3 % outside sector, 44.9 % removed
Ellbögen |  Valley station: Ellbögen (45∘–225∘)
Crest station: Sattelberg (90∘–270∘)
22.6 % within sector, 49.0 % outside sector, 28.4 % removed
Montana |  Valley station: Montana (30∘–100∘)
Crest station: Gütsch, Andermatt (120∘–180∘)
23.2 % within sector, 52.9 % outside sector, 23.9 % removed
Comprovasco |  Valley station: Comprovasco (270∘–30∘)
Crest station: Gütsch, Andermatt (290∘–89∘)
33.5 % within sector, 51.2 % outside sector, 15.4 % removed
Lugano |  Valley station: Lugano (270∘–80∘)
Crest station: Gütsch, Andermatt (290∘–89∘)
26.9 % within sector, 62.1 % outside sector, 10.9 % removed
Table 4: Overview of the stations and defined wind sectors used as
precondition for classification. For each target location (left column) the
corresponding valley station and crest station and their required wind sectors
are listed. In addition, the percentage of observations used for
classification (within sector), not used in classification (outside sector)
and removed due to missing data is shown (based on all available data). More
station details cf. Table 1.

The (conditional) two-component Gaussian mixture model with concomitant
variables Grün and Leisch (2008); Plavcan _et al._ (2014); Stauffer (2023) is
defined as follows:

$h(\mathbf{y},\mathbf{X},\mathbf{\theta},\mathbf{\alpha})=\underbrace{\big{(}1-\pi(\mathbf{X},\mathbf{\alpha})\big{)}\cdot\mathcal{N}(\mathbf{y},\mathbf{\theta}_{1})}_{\text{first
component}}+\underbrace{\pi(\mathbf{X},\mathbf{\alpha})\cdot\mathcal{N}(\mathbf{y},\mathbf{\theta}_{2})}_{\text{second
component}},$

where the joint (mixed) density $h()$ is the sum of two Gaussian components
$\mathcal{N}()$ times the probability $\pi()$ that an observation falls into
the second component. The potential temperature $\mathbf{y}$ is used to
separate the two components with $\theta_{\bullet}$ representing the
distribution parameters of $\mathcal{N}()$. The concomitant model modeling
$\pi$ is a logistic regression model employing an intercept as well as
relative humidity and wind speed ($\mathbf{X}$) and can be written as:

$\log\big{(}\frac{\pi}{1-\pi}\big{)}=\mathbf{X}^{\top}\mathbf{\alpha};~{}~{}~{}\pi=\frac{\exp(\mathbf{X}^{\top}\alpha)}{1+\exp(\mathbf{X}^{\top}\alpha)}.$

Once the required parameters $\theta$,$\alpha$ are estimated, the final
a-posteriori probability $\hat{p}\in[0,1]$ can be calculated using

$\hat{p}(\mathbf{y},\mathbf{X},\mathbf{\theta},\mathbf{\alpha})=\frac{\pi(\mathbf{X},\mathbf{\alpha})\cdot\mathcal{N}(\mathbf{y},\mathbf{\theta}_{2})}{\big{(}1-\pi(\mathbf{X},\mathbf{\alpha})\big{)}\cdot\mathcal{N}(\mathbf{y},\mathbf{\theta}_{1})+\pi(\mathbf{X},\mathbf{\alpha})\cdot\mathcal{N}(\mathbf{y},\mathbf{\theta}_{2})}$

The result is the probability
$\hat{p}(\mathbf{y},\mathbf{X},\theta,\alpha)=Pr_{\text{obs}}(\text{foehn})$
(Eqn. 1) with the same temporal resolution as the in-situ observations. After
the temporal upscaling (Sec. 3.1) this result used to label the data as
‘foehn’ or ‘no foehn’ event ($\text{Pr}_{\text{1h}}$; Eqn. 2).

##  Appendix C Supervised learning: Reconstruction

The hourly labeled data (Sec. 3.1, App. B) are are combined with the hourly
data derived from ERA5 (Sec. 2.2, App. A), serving the training data for the
supervised learning.

###  C.1 lasso: Logistic regression with lasso (L1) regularization

Due to the large number of possible covariates from ERA5 and the fact that
many of these variables are highly correlated, penalization is required. This
study employs lasso (least absolute shrinkage and selection operator) with L1
regularization, optimizing the following penalized log-likelihood:

$(\hat{\beta_{0}},\hat{\beta})=\underset{(\beta_{0},\beta)\in\mathcal{R}^{k}}{\text{argmin}}\Big{[}\frac{1}{2N}\sum_{i=1}^{N}(\text{Pr}_{\text{obs},i}-\beta_{0}-\text{ERA5}^{\top}\beta)^{2}+\lambda||\beta||_{1}\Big{]},$
(5)

with $\text{ERA5}\in\mathcal{R}^{N\times k}$, where $N$ is the number of
hourly foehn probabilities (Sec. 3.1) and $k$ the number of available
covariates (Sec. 2.2).

To find the optimal tuning parameter $\hat{\lambda}$, a 30-fold cross-
validation is performed which yields the optimal (regularized) set of
regression coefficients $\hat{\beta_{0}},\hat{\beta}$ using the _R_ package
`glmnet` (with s = "lambda.min").

###  C.2 stabsel: Logistic regression with lasso-based stability selection

As the lasso penalty shrinks more and more parameters to $0$ with increasing
$\lambda$, this model can also be used for stability selection (Meinshausen
and Bühlmann, 2010). The penalized regression model is estimated $M$ times for
a series of values $\lambda\in[0,\infty]$. The first $K$ parameters which
enter the model ($\neq 0$) will be memorized.

In this study, this is repeated $M=200$ times, estimating the model on a
random subset with only 50% of the available data. As the number of ‘foehn’
events can be rather low (0.6%–29% depending on station/time of day) a
stratified random subset (bagging) is drawn each time containing the same
proportion of ‘foehn’ and ‘no foehn’ events as the full data set. In each of
the $M=200$ iterations, the name of the first $K=40$ covariates entering the
model is stored. Those covariates which have been selected more than $60\%$ of
the time (at least 121/200 times) are then used to estimate an additional
unregularized logistic regression model.

  1. 1.

Whilst iteration is $\leq 200$:

    1. (a)

Draw stratified subsample of size $N/2$

    2. (b)

Estimate regularized coefficients $(\hat{\beta_{0}},\hat{\beta})$ for
different $\lambda$s

    3. (c)

Extract and keep the name of the first $K=40$ covariates entering the model

  2. 2.

Select the covariates which a selection frequency $>0.6$

  3. 3.

Estimate final (unregularized) logistic regression model with the covariates
from Step 2

Stability selection allows to automatically identify the most important
covariates among the ones available and to strongly reduce the complexity of
the final model (Step 3) without losing much predictive performance.

###  C.3 xgboost: Extreme gradient boosting

Extreme gradient boosting (XGBoost; Chen and Guestrin 2016) is a supervised
machine learning algorithm using a series of decision trees to best predict
the labeled data (response) using a series of explanatory variables
(covariates). Similar to random forests, XGBoost uses estimates and combines a
series of decision trees (a weak learner). The main difference is that XGBoost
uses a gradient decent method (thus the name) by optimizing the model based on
the outcome of the previous tree by finding additional split points to improve
the objective function.

Overly simplified, XGBoost starts with estimating one decision tree and
evaluates its residuals (i.e., the yet unexplained part/misclassified data).
These residuals are then used as input for the next iteration to estimate a
new tree to improve the misclassification from the previous step.

A series of tuning parameters exist such as the number of iterations to be
performed, the maximum depth of an individual tree, or the learning rate to
update the parameters in each iteration. Choosing these tuning parameters
wrong can lead to overly conservative models (not flexible enough) or models
which overfit the data, both of which can result in a poor (out-of-sample)
predictive performance. One way to avoid these scenarios is to perform proper
cross-validation. We performed a (limited) grid search using the _R_ package
`xgboost` (Chen _et al._ , 2023) with the following parameter space:

  * •

`eta`, learning rate: $\\{0.1,0.125,0.15,0.2\\}$

  * •

`max_depth`, maximum depth of each weak learner (tree): $\\{5,10,20\\}$

  * •

`min_child_weight`, minimum sum of instance weight needed in a child:
$\\{2,4,6\\}$

  * •

`gamma`, minimum loss reduction required: $\\{2,5,10\\}$

This yields an overall number of $108$ unique combinations. For each
combination, a 10-fold cross-validation using random stratified subsamples and
a maximum of $20$ boosting iterations is performed. It is then evaluated which
parameter set and boosting iteration resulting in the lowest average out-of-
sample error to estimate one final model using all data. A summary of the
selected tuning parameters is shown in Table 5. All models are estimated using
the ‘gbtree’ booster and the binomial logistic objective function.

| mean | median | range
---|---|---|---
subsample | 0.500 | 0.5 | 0.5-0.5
eta | 0.200 | 0.2 | 0.2-0.2
max_depth | 11.602 | 10 | 5-20
min_child_weight | 3.423 | 4 | 2-6
gamma | 2.868 | 2 | 2-10
nround | 19.997 | 20 | 19-20
Table 5: Mean, median and range of the selected XGBoost parameters via
stratified 10-fold cross-validation over all stations and times of the day
using the ‘full’ variable set.

##  Appendix D Long-term reconstruction

Figure 4 (main article, Sec 4.1) only shows results of the reconstruction
based on the ‘lasso’ model with the ‘full’ variable set. Figure 9 shows the
same information for all six stations including the results of all three
models (lasso, stabsel, xgboost) and both variants using either the ‘direct’
variable set or ‘full’ variable set.

Figure 9: Annual mean of daily maxima as in Figure 4 for all stations, all
three models in both variants.

The results of these six versions (three models, two variants each) show
differences but agree in large parts, showing similar features across the
years.

##  Appendix E Model performance/model comparison

Analogous to Figures 7 and 8 (main article; Sec,4.4 & 4.5) Figure 10 shows
average Brier scores from six-fold CV for all three models and model variants
(using the ‘direct’ or ‘full’ variable set). As discussed in the main article,
the ‘xgboost’ model might be subject to some overfitting (training) as it
always shows lowest scores across all stations and both model variants.

Figure 10: Comparison of average Brier score of all models (lasso, stabsel,
xgboost) and variants (‘direct’/‘full’) for all six stations. Average scores
of the six-fold CV are displayed for the test (red) and training (green)
period.

In addition, event-based metrics are used in the discussion (Sec 5) to compare
the results of this study to existing literature. Table 6 contains the results
for all six stations and both model variants (‘direct’/‘full’ variable set)
for all three supervised methods (‘lasso’/‘stabsel’/‘xgboost’).

Station | Variable set | FNR | FPR | PC
---|---|---|---|---
Altdorf | direct |  $20.0$/$18.9$/$13.6$ |  $0.6$/$0.7$/$0.3$ |  $98.4$/$98.3$/$99.0$
| full |  $15.7$/$17.1$/$10.7$ |  $0.4$/$0.6$/$0.2$ |  $98.8$/$98.5$/$99.2$
Innsbruck | direct |  $42.9$/$40.6$/$31.7$ |  $0.6$/$0.8$/$0.3$ |  $97.6$/$97.6$/$98.4$
| full |  $37.2$/$36.1$/$24.0$ |  $0.5$/$0.8$/$0.2$ |  $98.0$/$97.8$/$98.8$
Ellbögen | direct |  $19.4$/$19.9$/$12.9$ |  $2.9$/$3.0$/$1.6$ |  $94.2$/$94.0$/$96.4$
| full |  $13.9$/$15.9$/$7.9$ |  $1.9$/$2.4$/$1.0$ |  $96.0$/$95.3$/$97.8$
Montana | direct |  $24.1$/$24.4$/$14.5$ |  $1.7$/$1.9$/$0.8$ |  $95.7$/$95.5$/$97.6$
| full |  $19.4$/$21.7$/$12.1$ |  $1.3$/$1.8$/$0.7$ |  $96.6$/$95.9$/$98.0$
Comprovasco | direct |  $17.4$/$17.6$/$10.3$ |  $1.7$/$1.8$/$0.9$ |  $96.4$/$96.3$/$97.9$
| full |  $13.7$/$15.3$/$8.4$ |  $1.3$/$1.6$/$0.7$ |  $97.2$/$96.7$/$98.4$
Lugano | direct |  $17.2$/$17.3$/$11.8$ |  $0.9$/$1.0$/$0.5$ |  $97.9$/$97.8$/$98.6$
| full |  $15.3$/$16.4$/$9.5$ |  $0.7$/$0.9$/$0.4$ |  $98.2$/$97.9$/$99.0$
Table 6: Miss rate (FNR), false alarm rate (FPR) and percent correct (PC) in
percent based on the entire period (in-sample) for model comparison. Scores
are shown for the ‘lasso’/‘stabsel’/‘xgboost’ models in this order.