On the nonstationary identification of climate-influenced loads for the semi-probabilistic approach using measured and projected data

Abstract A safe and economic structural design based on the semi-probabilistic concept requires statistically representative safety elements, such as characteristic values, design values, and partial safety factors. Regarding climate loads, the safety levels of current design codes strongly reflect experiences based on former measurements and investigations assuming stationary conditions, i.e. involving constant frequencies and intensities. However, due to climate change, occurrence of corresponding extreme weather events is expected to alter in the future influencing the reliability and safety of structures and their components. Based on established approaches, a systematically refined data-driven methodology for the determination of design parameters considering nonstationarity as well as standardized targets of structural reliability or safety, respectively, is therefore proposed. The presented procedure picks up fundamentals of European standardization and extends them with respect to nonstationarity by applying a shifting time window method. Taking projected snow loads into account, the application of the method is exemplarily demonstrated and various influencing parameters are discussed.


Introduction
Changes in natural processes as well as anthropogenic drivers, such as greenhouse concentration and reflectivity or absorption of the sun's energy, can lead to changes in climate conditions (National Academy of Sciences (2020), Fahey et al. (2017), and Hayhoe et al. (2018)). The existence of the corresponding global climate change is supported by numerous research from various scientific disciplines and is generally accepted, see, Pachauri and Meyer (2014). Institutions, such as the Climate Computing Centre Potsdam in Germany, are continuously contributing within an international scientific community (Intergovernmental Panel on Climate Change-IPCC) to model future climate scenarios and to enhance projections of global and regional climate developments. In general, significant changes in intensity and frequency of extreme weather events are predicted (Pachauri and Meyer (2014), Hübener et al. (2017a), and Perkins et al. (2012)). and its components is studied at the University of Weimar. Amongst others, the aim is to provide methods to trace and adjust loads within the structural analysis of individual structures/members based on available climate data, such as for wind, snow or temperature. Projections of the aforementioned climate research are used as data basis, which allow an outlook on future developments of structural loads and enable an assessment of the associated structural safety (see, e.g., Vanem (2015), Croce et al. (2021), Croce et al. (2019), Ren et al. (2019), Yilmaz (2017), and Aziz and Yucel (2021)). The corresponding data is usually provided as time series with high temporal (daily values) and spatial (12 km) resolution for periods up to the year 2100, see "Regional Climate Projections Ensemble for Germany" (2022), Jakob et al. (2014), Warrach-Sagi et al. (2018). A general overview is given in Climate Projections (C3S; 2022). However, it should be noted that precise projections, i.e. regarding the temporal and spatial occurrence of extreme events, are very complex and subjected to comparably high uncertainties from today's perspective, see, Salas et al. (2018), Blöschl and Montanari (2010), Kundzewicz and Stakhiv (2010), Hübener et al. (2017b), and Scafetta (2021).
Next to the projections, records of climatic parameters may also be taken into account for statistical analyses considering associated climate loads. Following extreme value theory, probability density and cumulative distribution functions can be developed from the corresponding time series, allowing to identify probabilities of exceedance and non-exceedance of decisive load events within the service life of a structure (Coles (2001)). The current design standards presuppose stationary conditions, in which all extreme value events to be expected in the service life of a structure are statistically described by a density function, see, e.g., CEN (European Committee for Standardization; 1990), Ragno et al. (2019), and Fischer (1999). The concept of stationarity does not consider, whether the probability of occurrence of a certain extreme event and the corresponding probability distribution are connected to temporal changes (Coles (2001), Jacob (2013)). Due to the long service life of structures of 50 or 100 years, the validity of stationary assumptions may therefore be scrutinized, see, e.g., Coles (2001), Fischer et al. (2012), Soukissian and Tsalis (2015), and Milly et al. (2008). This is supported by numerous studies, illustrating fundamental changes in frequencies and intensities of extreme events connected to trends and shifts in corresponding data (Melillo et al. (2014), Coumou and Rahmstorf (2012), Mazdiyasni et al. (2017), Mallakpour and Villarini (2017), Wahl et al. (2015), Vahedifard et al. (2016), andJongman et al. (2014)). Accordingly, Coles (2001), Gilleland andKatz (2011), andCooley (2012), as well as Salas and Obeysekera (2014) have proposed concepts of nonstationarity to capture the temporal changes of climate variables. Based on these approaches, many scientific studies have been carried out to determine and assess the impact of climate change on structural loads, e.g., for temperature, wind, snow, or rain (Hamdi et al. (2018), Jeong et al. (2020), Qin et al. (2022), Croce et al. (2021), Le Roux et al. (2020), Aziz and Yucel (2021), Ren et al. (2019), Tabari (2020), Tian et al. (2021), Tramblay et al. (2013), and Yilmaz (2017)). The entire data set is represented with a density function providing a certain 98% quantile in the upper half of the figure. In this stationary view, the value will not be exceeded in the observation period with a probability of 98%, defining a certain failure probability for a structure.
The subdivision of the data set into two periods illustrates a time-dependent change in the density functions and a variation, in the present example of Figure 1 an increase of the 98% quantile. The associated time-dependent occurrence probability of an extreme value event influences the failure probability of the structure, which increases from the first to the second observation period in the example of Figure 1. This effect can be attributed to the temporal trend in the data leading to approaches of so-called nonstationary extreme value analysis. Since climate developments are generally accompanied by such trends (cf., Figure 2), the nonstationary approach may be of particular importance for the derivation of reliable design loads and for the approximation of the structural reliability in the considered service life (Coles (2001), Gilleland and Katz (2011), Cooley (2012), and Salas and Obeysekera (2014)). Respective comparisons of stationary and nonstationary approaches are made by Coles and Pericchi (2003), Serenaldi and Kilsby (2015), Montanari and Koutsoyiannis (2014), and Obeysekera and Salas (2014), as well as Madsen et al. (2013), considering sources of uncertainty in nonstationary models which may lead to physically inconsistent results in case of misspecification.
The treatment of nonstationary time series in extreme value analysis is considered in numerous studies and publications using different methods, see for instance, Coles (2001), Hennemuth et al. (2013), Read and Vogel (2015), Salas et al. (2018), andSlater et al. (2021). A typical approach for corresponding nonstationary models is to include time dependency of distribution parameters or associated physical processes as a covariate, as for instance, shown by Coles (2001), Vanem (2015), and Tian et al. (2021), accounting for linear trends and shifts between different periods. In the context of design tasks, it is not expedient to analyse time varying quantiles, leading to approaches of the expected waiting time (EWT; Olsen et al. (1998), Cooley (2012, and Salas and Obeysekera (2014)) and the expected number of exceedance (Parey et al. (2010), Cooley (2012). A different approach is presented by Rootzén and Katz (2013) as well as Mudersbach and Bender (2017) taking time-dependent non-exceedance probabilities and the intended service life of a structure into account. Mudersbach and Bender (2017) define a shifting time window where a certain section of a time series or length of time, respectively, is considered in the sense of a section, in which stationarity is assumed (Führböter (1976), Mudersbach (2009)). By continuously shifting the window in the direction of the time axis and applying statistical evaluations, nonstationarity in environmental events can be modeled. This non-parametric procedure allows determining the temporal development of probability distribution parameters directly from corresponding measured data, such that there is no need to specify a corresponding prior model, which is common in the sense of a parametric consideration of nonstationarity, see, Mudersbach (2009). In addition, a main advantage of this approach  is the straightforward inclusion of the intended service life of a structure and the possibility to follow the temporal development of parameters.
In this study, the procedure of the shifting time window is enhanced with respect to directly determining climate loads for the structural design, i.e. design values of loads. A probabilistic design procedure is applied to estimate the structural reliability, which is a probabilistic measure of structural safety and defined as the complement of the probability of failure, given that one or more of the structures' material properties, geometric dimensions or load properties are random or incompletely known (Ditlevsen and Madsen (2007), Nikolaidis et al. (2004), andSpäthe (1992)). In this paper, methods based on the First Order Reliability Method (FORM) considered in European standardization, namely EN 1990 (CEN (European Committee for Standardization), 1990), are extended by the nonstationary approach incorporating the reliability index and the sensitivity factors. In conjunction with climate projections, a procedure for predicting future probabilities of extreme weather events emerges. Results of the proposed method will continuously improve with the enhancement of the abovementioned climate models and projection data, respectively. In addition, monitoring data from measuring systems and weather stations may be incorporated in the analyses as well.
This paper is structured as follows. First, required relationships of classical stationary extreme value representation as well as safety concepts of European standardization are presented. Subsequently, the nonstationary representation based on shifting time window is enhanced as mentioned above and the methodology is outlined. Finally, the application of the procedure is demonstrated with an example of snow load concluding with a discussion and summarizing of results.

Classical stationary representation of extreme values
For the safe design of structures, frequencies and intensities of extreme load events (maximum and minimum values) during the service life of a structure are of particular importance. In case of climate actions and their natural variations, loads will always be time-dependent. Considering longer periods, the occurrence probability of an extreme load will obviously increase. Corresponding parameters therefore depend on the observed time duration, the so-called reference period t (e.g., service life of a structure).
Statistical modelling of extreme values is treated in many fundamental publications; see for instance, Coles (2001) or Späthe (1992). Assuming the random variable X to be a climate load, the considered reference period t is divided into N ¼ t=T observation periods T. By singling out the extreme value (maximum or minimum, respectively) from each observation period, a probability density function f X ðxÞ and a cumulative distribution function (CDF) F X ðxÞ of these extreme values is obtained. However, a statistical independence of neighbouring extreme values is to be assured.
In the following, an observation period of T ¼ 1 year is assumed, leading to the sample of observations x 1 ; . . . ;x N . These are associated to the random variables X i and to distribution functions F X i x ð Þ with i ¼ 1; . . . ; N. Considering maximum extreme values, a new random variable "upper value" X U ¼ max X 1 ; . . . ;X N ð Þ with the distribution function F X U x ð Þ is specified, expressing the probability of the random variable X U not exceeding a defined value x in N intervals of time T (i.e. the reference period t), see, Späthe (1992) or Fischer (1999): If the random variable of the "lower value" X L ¼ min X 1 ; . . . ;X N ð Þ is of interest, as for instance, when analysing minimum temperatures, comparable relations may in general be developed based on the distribution function of minimum extreme values, see, Späthe (1992). However, by transforming minimum extremes with a multiplication of (−1), correspondences of maximum extreme values can be maintained according to Rootzén and Katz (2013). Assuming stationarity, all X i are independent and identically distributed, such that the corresponding F X i x ð Þ are equal for each i ¼ 1; . . . ; N and thus can be treated as one CDF F X i x ð Þ¼F X x ð Þ. With a threshold x ¼ r and the observation period of T ¼ 1 year, Eq. (1) may be expressed for a reference period of t years: With regard to structural analyses, the reference period depends on the service life of a structure, which is often defined to be t ¼ 50 years in accordance to (CEN (European Committee for Standardization), 1990). Depending on the individual design situation and the connected safety requirements, EN 1990 defines an annual non-exceedance probability of F X ðrÞ ¼ 0:98 (characteristic value) for climate loads due to temperature, wind or snow. Considering Eq.
(2), this leads to and the corresponding quantile of The non-exceedance probability P X U � r ½ � represents the "safety level for the entire lifetime t". The 98% quantile of the annual extreme value distribution equals the 36.4% quantile of the 50-year extreme value distribution, as can be seen in Figure 3. Within the service life of 50 years, r is exceeded at least once at mean, as it is 1 À 0:364 ð Þ¼ 0:636>0:364 and thus the exceedance probability greater than the non-exceedance probability.

Basics of FORM
In the course of the design of structural components, it is to be ensured that the effects of actions (E) caused by loads do not exceed the component resistance (R) with a defined reliability, i.e. the probability of failure is limited to a desired magnitude (Ditlevsen and Madsen (2007), Nikolaidis et al. (2004), and Späthe (1992)). Here, R and E are assumed to be random and independent quantities (random variables) with normally distributed densities f E e ð Þ (expected value μ E , standard deviation σ E ) as well as f R r ð Þ (expected value μ R , standard deviation σ R ). The associated two-dimensional joint probability density function  Fig. 4a reflects realized pairs of resistance and action values. Failure occurs when R<E, leading to the following linear limit state function: Eq. (5) delimits the density of the realizations f R;E r; e ð Þ into the ranges "failure" and "no failure". The smaller the distance between the limit state line and the centre of the two-dimensional density function, the larger the failure range and thus the probability of a failure P f . According to Hasofer and Lind (1974) as well as Rackwitz and Fiessler (1978), the distance corresponds to the so-called reliability index β, see Fig. 4b. The associated design point on the limit state function (r � ;e � ) or (r � ;ẽ � ), respectively, i.e. the point where the limit state is exceeded with greatest probability, can be described by sensitivity factors α E and α R as shown in Fig. 4c. With respect to the normative semi-probabilistic concept, fractiles r p and e p of the distribution densities are defined as characteristic quantities of resistance and action. On the basis of Eq. (5), a nominal safety factor is introduced by r p ¼γ p � e p and decomposed into partial safety factors with γ p ¼γ R � γ E . Considering the quantiles r p and e p as well as the values of the design point r � and e � , they are described according to Späthe (1992) Introducing the design values of Eq. (6) and (7) into the condition E d � R d (cf., Eq. (5)) yields in the following basic equation for design: Based on Fig. 4, the design values or γ E and γ R , respectively, are calibrated or determined such that they are associated to a specified β Target . In doing so, partial safety factors take into account statistical uncertainties of the scattering action and resistance quantities, which are mainly composed of the random properties of load, material and geometry as well as uncertainties of modelling (Ditlevsen and Madsen (2007), Nikolaidis et al. (2004), and Späthe (1992)).
As can be seen, the safety index β is directly related to the probability of failure. The failure occurs for Z ¼ R À E<0. Since R and E are scattering, normally distributed quantities, this also applies to Z according to variance propagation with μ Z ¼μ R À μ E and σ 2 Z ¼ σ 2 R þσ 2 E according to Hasofer and Lind (1974). The associated distribution function reflects the failure probability P f with F Z ðz ¼ 0Þ ¼ P½R À E<0�, which corresponds to the fractile of the density function f Z z ð Þ at z ¼ 0, see, Figure 5. Based on the densitiy function shown in the figure, the corresponding reliability index can be However, using the standard normal distribution Φ(.), the probability associated with the random variable Z being smaller than the fractile value at z ¼ 0 may be determined as shown in Figure 5 on the right, which corresponds to the failure probability: The presented relations correspond to fundamentals of the First Order Reliability Method (FORM), which is to be applied for the calibration of design values according to EN 1990. Methodological background can for instance, be found in Ditlevsen (1981), Beacher and Christian (2003), Nikolaidis et al. (2004), and Melchers and Beck (2018). Inaccuracies of the FORM are treated by Der Kiureghian et al. (1987) or Zhao and Ono (2001) and advantages are discussed by Rackwitz (2001). The associated R-E model was firstly validated under the assumption of normally distributed and mutually independent basic variables with a linear limit state function (Hasofer and Lind (1974)), cf. Figures 4 and 5. If non-normal distributions are to be considered, such as the extreme value distributions, they are transformed at the design point to normal ones, so that P f � Φ À β ð Þ is a good approximation, see, Rackwitz and Fiessler (1978), Späthe (1992), as well as Melchers and Beck (2018).

Calibration of the design values based on EN 1990
With reference to the definition of the design point based on sensitivity factors in Fig. 4c as well as Eq. (3) and (9), the minimum exceedance or non-exceedance probabilities for design values of effects of actions E d and assigned resistances R d result to (CEN (European Committee for Standardization; 1990)): To avoid costly reliability analyses, FORM sensitivity factors α E and α R as well as target values for the safety index β Target are defined in (CEN (European Committee for Standardization), 1990) for calibrating design values. Since the present work essentially concentrates on the load side, the focus in the following will be on the determination of design values E d .
The target value of the reliability index is specified in (CEN (European Committee for Standardization), 1990) for different limit states, reliability classes and reference periods. If reliability class RC 2 and a reference period of t ¼ 1 year are selected in the context of the ultimate limit state, β 1;Target ¼ 4:7 is defined. For the extrapolation to an arbitrary reference period of t years, the following formula given in (CEN (European Committee for Standardization), 1990) and based on Eq. (2) and (9) can be applied: Thus, for a reference period of t ¼ 50 years, the target reliability index is determined to be When considering e.g., "lower values", the formulation is adjusted according to Herbrand et al. (2021) considering negative β values: According to (CEN (European Committee for Standardization), 1990), the FORM sensitivity factor of an action α E;t is defined to be α E;50 ¼ À 0:7 for a reference period of 50 years, providing that the ratio of standard deviations of effects E and resistances R is within the limits 0:16 � σ E;50 =σ R;50 � 7:6. If the considered action does not represent the leading action but an additional accompanying one, α E;50 ¼ 0:4 � ðÀ 0:7Þ ¼ À 0:28 may be assumed (for the origin of the values see, e.g., Späthe (1992)). For other reference periods, α E;t can be approximated according to Herbrand et al. (2021)

Representation of extreme values
Effects of nonstationarity described in section 1 are considered according to Mudersbach and Bender (2017). Based on Eq.
(1), nonstationary quantiles may be determined by defining areference year (e.g., year of manufacture; index i ¼ 1) as "entry into a trended data set" and a target safety level using the probability p t;Target for adverse exceedance or non-exceedance within the service life t. According to Mudersbach and Bender (2017), the following equation describes the non-exceedance probability p t;Target considering "upper values": Nonstationarity is reflected by expressing the probability of the random variable X U not exceeding a defined value r in N intervals of time T, where in the individual intervals i (time windows) the probability distribution functions F X i ðrÞ are not assumed to be the same, respectively, cf., Eq.
(2). For each interval, an individual probability distribution function is determined, i.e. the time-dependent probability variation of the random variable not exceeding the defined quantile r is considered. In Eq. (15), the determination of the quantile r leading to p t;Target is basically an optimization task or iteration, respectively, in which the combination F X 1 r ð Þ; F X 2 r ð Þ; . . . ; F X N r ð Þ � � of a pre-set r-value fulfilling Eq. (15) is searched.
According to Eq. (3), p t;Target ¼ 0:364 is defined to maintain the safety level of the standard.
Considering "lower values", the previously (in section 2) mentioned transformation with a factor of (−1) can be applied in conjunction with Eq. (15).
To cover the temporal changes within the data series, the distribution functions F X 1 x ð Þ; . . . ;F X N x ð Þ are determined within the framework of frequentist analysis with a shifting time window model according to Figure 6 (see, Mudersbach (2009) or Mudersbach and Bender (2017)). Based on the period t, an interval is extracted from the data series of the random variable X and divided into N time windows of the length k. The first window is arranged such that the above-mentioned reference year corresponds to the last value of this window. Finally, for each time window i ¼ 1; . . . ; N the distribution function F X i x ð Þ can be determined.
The time window length k can in general be arbitrarily chosen, however, the following restrictions ought to be considered (Mudersbach (2009)): • Stationarity is presupposed within a time window. Too large window lengths make this assumption no longer justifiable.
• The window length k should not be chosen too small, since the statistical power with regard to the population to be represented is also reduced with the number of available sample elements.
Thus, the exact choice of k has influence on the calculation result. With regard to standardization, however, data series with minimum observation periods of 20 years are required and more than 30 years are recommended in order to reliably cover the natural range of variation of climatic action (CEN (European Committee for Standardization), 1991a). This also corresponds to definitions of the World Meteorological Organization (WMO) and the German Weather Service (DWD), which combine Applying the presented method and selecting a window length k, advantageously only a certain and for the intended service life of a structure decisive part of the available time series is considered. The procedure allows incorporating trends in different time fractions and directly capturing the variability of temporal distribution parameters based on the predicted data, i.e. avoiding to review the suitability of different trend models. A disadvantage compared to further nonstationary approaches mentioned in the first section is the necessity of the choice of an appropriate time window length as well as a lower flexibility of considering different trend models for individual distribution parameters. Furthermore, a larger number of regression analyses may result in additional uncertainties in the context of error propagation.

Inclusion of the reliability index and the FORM sensitivity factors
In principle, design values considering nonstationary conditions can be determined using the relationships of the previous section. For reasons of simplicity, novel formulations based on the reliability index and the FORM sensitivity factors to directly incorporate specifications of EN 1990 are presented here. Considering positive β E;t -values (e.g., in case of "lower values") and making use of Eq. (1) and (9) gives the following correspondence based for the entire service life: Accordingly, with reference to Eq. (13) for the analysis of negative β E;t -values, which is the case for "upper values", it is: The variation of the distribution functions F X i x ð Þ due to nonstationarity is covered by varying parameters α E;1;i � β 1;i , with i ¼ 1; . . . ; N. The individual values of α E;1;i � β 1;i for the time windows depicted in Figure 6 are determined by rearranging Eq. (10), taking into account the specific distribution functions F X i x ð Þ¼F E i e ð Þ: Comparable to the previous section, an iterative process is performed for calibrating the design value E d based on Eq. (16) to (18). The value α E;t � β t;Target is to be understood as the target value applying to nonstationary conditions. It is defined such that the target safety level of the stationary case is reached according to section 3.2 at the end of the service life, i.e. α E;t � β t;Target ¼ À 0:7 � 3:83 ¼ À 2:68 after a period of 50 years. Eq. (18) gives values of α E;1;i � β 1;i for all time windows and a pre-set value of e. Introducing these into Eq. (16) or (17), respectively, the nonstationary value of α E;t � β t for the pre-set value of e is determined. In case it corresponds to the target value of α E;t � β t;Target , the corresponding quantile e � ¼ e is obtained representing the calibrated design value of the effect E d .

Uncertainties
An essential aspect for assessing the results obtained from the presented procedure is the awareness of uncertainties associated with the assumptions. Details on treating uncertainties in the context of reliability analyses can for instance, be found in Späthe (1992), Melchers and Beck (2018), and Most (2012), while uncertainties of nonstationary extreme value statistics are dealt with in Coles (2001), Serenaldi and Kilsby (2015), Read and Vogel (2015), Salas et al. (2018), and Slater et al. (2021), as well as Rootzén and Katz (2013). Next to the approximations and assumptions already mentioned in the previous sections, the choice of probability distribution functions and their adjustment to the given data has major influence on results of the time window model, potentially producing relatively large uncertainties. In addition, design events with sufficiently large return periods are mostly uncertain, as they do not allow to physically justify statements on their occurrence probabilities. In combination with uncertainties in data series of climate projections, substantial influences may have impact on the overall results. Consequently, effects of climate change on the safety elements of the semi-probabilistic concept can scientifically be substantiated with the presented methodology, however, a claim that predicted and quantified values will in fact occur cannot be deduced. The results obtained should always be critically scrutinized and validated by further investigations, see, e.g., Coles (2001). or the Kolmogorov-Smirnov test (Massey (1951)), can be used equally, however, this was not pursued further in the context of this study. For reviewing the significance of the trend in the time series introducing nonstationarity, the Mann-Kendall test is selected (see, Van Gelder et al. (2008)). Finally, the design loads for different future periods can be specified and comparisons to the stationary representation according to section 2 are possible. The methodology is illustrated in the following section in detail considering a steel girder with a snow load. Standard errors (SE) of the considered statistics and design values are estimated in the course of parametric bootstrapping procedures (Panagoulia et al. (2014)).

Example
The relations presented in the previous sections have been implemented in a tailored software based on the environment of Visual Basic for Applications (VBA). They are exemplified in the following using a simple design example. A single-span beam of structural steel S 235 as part of a roof is stressed by a uniformly distributed snow load q s according to Figure 8. Stability failure of the beam with respect to lateral torsional buckling is prevented by means of constructional provisions. For a better illustration, further loads (for example, the deadweight of the structure) and uncertainties of the load model associated with the snow are neglected.
According to (CEN (European Committee for Standardization), 1991a), q s depends on the snow load s on the ground, the shape parameter μ, the environmental coefficient C e and the temperature coefficient C t by: The parameter b denotes the width of the roof influencing the load introduction into the beam. Assuming elastic behaviour, the effect of action expressed by the maximum bending moment at midspan yields to: As can be seen, the random variable M y is related to the ground snow load s in terms of a linear transformation, which is also valid for corresponding quantiles of the probability density function. The non-exceedance probability related to a certain quantile s p is equal to the one of the corresponding moment ðM y Þ p (representing the quantile of the effect of action). Density or distribution functions are therefore developed based on the ground snow load s in the following.
For that purpose, snow depths h s are adapted from the regional climate projections (Climate Projections (C3S; 2022)) referring to the city of Bad Reichenhall in Germany, see, Figure 8. Data of the worst emission scenario RCP 8.5 is considered, as it presumably provides the highest time-dependent impact on the corresponding loads. However, uncertainty in the data is not incorporated in the following analyses. The data series (Hübener et al. (2018)) taken from (ReKliEs-De-"Regional Climate Projections Ensemble for Germany" (2022)) provides daily values of snow depths for the period from 1.1.1950 to 31.12.2100. It was bias-corrected, processed identifying annual extremes, and provided by Uzair and Abrahamczyk (Forthcoming). Within this processing, quantiles of the simulated data are compared to available quantiles of recorded data in the considered period. Recorded data are provided by the "Climate Data Center" (CDC) of the German Weather Service (DWD) in a daily resolution with a high data quality (see, DWD Climate Data Center (CDC; 2022), Behrendt et al. (2011), and Kaspar et al. (2013)). The differences in the quantiles can be traced either additively or multiplicatively and subsequently they are deployed for the adaption of the projected data. In this specific case, the additive adjustment is applied to determine the bias-corrected time series as it produces a more realistic forecast according to Uzair and Abrahamczyk (Forthcoming). The procedure preserves absolute changes of the data as well (Cannon et al. (2015)). However, it should be noted that time series considering the extremes of periods from e.g., 1 st of July to 31 st of June may be more productive for snow loads, aiming to avoid that two consecutive annual extremes are generated by the same environmental event. Depending on statistical independencies of events, this can of course not be generalized. As the methodology of the nonstationary approach is to be the main focus of this example, the abovementioned available data incorporating the annual extremes is taken as a basis in this study.
Considering the density of the snow ρ s and taking into account g E ¼ 9:81m � s À 2 , snow loads s may be determined as follows: The German Weather Service describes the relationship between the snow density in [kg � m À 3 ] and the snow depth in [m] using an empirically derived function, see, Caspar and Krebs (1974): Actual snow densities can in general be approximated by Eq. (22) with a certain degree of accuracy. It would be preferable and more precise to rely on values of the water equivalent, for which predictions are rather difficult and currently not available in ReKliEs-De. Since the present paper strongly concentrates on the methodology as mentioned above, the inaccuracy regarding the considered load data is for sure acceptable. However, if data series from weather stations (instead of projections) are fed into the procedures, values of the water equivalent could of course be taken into account without difficulty.
Based on Eq. (21) and (22), the annual maxima (T ¼ 1 year) of the snow load s shown in Figure 9 are determined. The time series is analysed using the nonparametric statistical Mann-Kendall test, which is applied in many studies for determining monotonic trends (see, Van Gelder et al. (2008), Ren et al. (2019), and Tian et al. (2021)). In comparison to other procedures, no prior assumptions need to be made regarding the course of the trend, as the increase or decrease of ratios between values is analysed (Faulkner et al. (2019)). The Mann-Kendall test shows a highly significant trend at a 5 % significance level for the considered time series, which introduces nonstationarity. However, it is to be pointed out that such tests are subjected to uncertainties. The plausibility of the results should be physically justified by additional preliminary information related to the underlying stochastic process (Serenaldi et al. (2018)). In the present case, climate warming in Europe can in principle lead to less snow events resulting in lower snow loads in the future, such that a significant downwards trend of the corresponding extremes is credible (Hübener et al. (2017a)).
Starting from the reference year 2021, a service life of t ¼ 50 years is assumed. Referring to a selected time window length of k ¼ 40 years, the relevant measures are given in the time interval ½2021 À k þ 1; 2021 þ t À 1� ¼ ½1982; 2070� (see, Figure 9). Based on this data, the nonstationary extreme value analysis is performed according to section 4.1 by determining the distribution functions F S i s ð Þ for the time windows i ¼ 1; 2; . . . ; N ¼ 50 considering the different probability distribution types of Table 1, see, Späthe (1992). While the Gumbel and GEV distribution are often preferably used for time-variant climate actions, the Lognormal and the Weibull distributions offer the advantage of defining a lower limit value of s ¼ 0 for the considered snow loads. The parameters of the distribution functions (a, u, A, k, m u , σ u , f 1 , f 2 , τ) illustrated in Figure 10 are estimated based on the method of moments, see, e.g., Späthe (1992) and Diburg et al. (2008). The corresponding time histories of mean m and standard deviation σ derived from the individual time windows are plotted in Figure 9. Standard errors (SE) for the mean values, standard deviations, as well as for design values are estimated via parametric bootstrapping according to Panagoulia et al. (2014). In the context of shifting time windows bootstrap samples for the statistics of interest are determined for each window by producing 1000 observations, each including k values, and considering estimates as distribution parameters fitted to the original data. As a result 1000 bootstrap samples, each consisting of 50 distribution parameter sets are available to estimate the SEs.
Considering the non-exceedance probability of p 50;Target ¼ p t 1;Target ¼ 0:364 according to Eq. (3), the characteristic value r ¼s k is determined using Eq. (15). In the present example, the characteristic ground snow load of s k ¼ 1:65kN � m À 2 (SE = 0.05 kN � m À 2 ) is obtained when applying a Gumbel distribution. Figure 11 displays the development of the non-exceedance probability within the service life of the structure, where in Figure 11b the annual probabilities p 1;i ¼F S i s k ð Þ and in Figure 11a cumulative values p n ¼F S 1 s k ð Þ � F S 2 s k ð Þ � � � � � F Sn s k ð Þ up to the observed time of n years are shown. As can be seen in Figure 11b, the annual target probabilities of 0.98 defined by the code are partially undercut or exceeded, successively, provided that the specified value of 0.364 is obtained for the total service life of t ¼ 50 years. Since in the first years the annual probabilities are lower than the annual target value of 0.98, the cumulative probabilities in Figure 11a continuously fall below the corresponding target values in this example.

Table 1. Considered probability distribution functions for the description of snow load
Distribution type

Definition range
Gumbel Next, the design value of the ground snow load is determined based on the nonstationary approach in section 4.2. According to Section 3.2, the target reliability index is set to β 50;Target ¼ 3:83 for the service life of 50 years and the FORM sensitivity factor to α E;50 ¼ À 0:7.
Using the probability distribution functions of the individual time windows, Eq. (18) is applied to determine the corresponding non-exceedance probabilities F E i e ð Þ¼F S i s d ð Þ (i ¼ 1; 2; . . . ; N ¼ 50) and the associated annual parameters α E;1;i � β 1;i for a pre-set value of e ¼ s. Introducing them into Eq.
(17) yields a value of α E;50 � β 50 for the entire service life of the structure. Provided it meets α E;50 � β 50;Target ¼ À 0:7 � 3:83 ¼ À 2:68, the pre-set s-value reflects the design value of the ground snow load s d . It is calculated to s d ¼ 3:20kN � m À 2 (SE = 0.20 kN � m À 2 ) with the Gumbel distribution in the present example.
Comparable to the non-exceedance probabilities shown in Figure 11, the individual reliability indices are shown in Figure 12 in terms of the annual values α E;1;i � β 1;i and the cumulative value of α E;n � β n at an observation time of n years. The aforementioned target values are also included for reasons of comparability. In addition, Figure 13 depicts the development of corresponding design snow loads depending on the reference period as well as the annual loads related to the annual Figure 10. Illustration of the probability distribution functions of Table 1 for a mean value of m = 5 and a standard deviation of σ = 1.5. Figure 11. Non-exceedance probabilities considering s k ¼ 1:65kN � m À 2 (SE = 0.05 kN � m À 2 ) for the target safety level of EN 1990 with a service life of t = 50 years (p 50;Target ¼ 0:364) and comparison to corresponding target values, a) cumulative probabilities p n and b) annual probabilities p 1 .
target values of α E;1 � β 1;Target ¼ À 3:79. For comparison, characteristic values are included in the diagrams of Figure 13 as well (cf., Figure 11). The decreasing trend of the time series in Figure 9 is reflected in the course of the loads. In principle, the annual loads follow a combined trend of mean and standard deviation shown in Figure 9, however, strongly influenced by the standard deviation. Since the sliding time window approach is nonparametric, short-termed developments in the time series contrary to the global trend are taken into account. Therefore, combined values representing the load development in dependency of the time (reference period) also cover short termed changes of climatic developments during the service life.
If the stationary approach according to section 2 is considered, the characteristic value as 0.364fractile value of the density function for the period [1982,2070] is determined with s k ¼ 1:71kN � m À 2 (SE = 0.21 kN � m À 2 ) and is thus 3.6 % higher than the respective value considering nonstationarity. Regarding the design value, the fractile related to α E;50 � β 50;Target gives s d ¼ 3:35kN � m À 2 (SE = 0.50 kN � m À 2 ), a value 4.7% larger than with the nonstationary approach. The difference is manageable in the present example, although the considered data shows a significant trend.
Based on the nonstationary design value and characteristic value of the ground snow load, Eq. (7) gives the following partial safety factor: Compared to the specification γ E ¼ 1:5 of the code, this difference appears to be quite substantial. However, when determining the snow load according to (CEN (European Committee for Standardization), 1991b), a design value of s d ¼γ E � s k ¼ 1:5 � 2:2kN � m À 2 ¼ 3:3kN � m À 2 results, i.e.
a value close to s d ¼ 3:20kN � m À 2 . The difference in the partial safety factor is therefore put into perspective from a practical design point of view. This would indeed have to be considered when  approaches of introducing an additional safety factor to account for changing climate, i.e. next to the current partial safety factor, are to be followed.
Since the design value of the material strength exceeds the "guaranteed" nominal minimum yield strength of f y ¼ 23:5kN � cm À 2 , it is state of the art for members not susceptible to stability to assume f y;d ¼ 23:5kN � cm À 2 and at the same time consider nominal cross-section dimensions for the calculation of the plastic moment resistance. The uncertainties caused by geometric tolerances are "covered" by the lower material strength. According to Kindmann et al. (2017), the plastic moment resistance of the IPE 220 section is: Comparing with the design value of the acting bending moment in Eq. (24), the profile provides sufficient capacity and is verified based on Eq. (8).
To exemplify the influences of the parameter k as well as the choice of the probability distribution function type, the example is additionally evaluated for time window lengths of k ¼ 20; 25; 30; 35; 40; 45; and 50 years considering the distribution functions specified in Table 1. A summary of the results is given in Table 2. Following the AICc criterion (see, Laio et al. (2009)), better conformity to the data set is obtained for Weibull and Gumbel distributions (in this order) compared to the Lognormal and GEV distributions (see, Table 2). With regard to the criterion, distributions with lowest AICc-values, which are determined based on the Likelihood function, explain the dependent variable best (here the snow load), taking the complexity of the distribution function via the number of estimated parameters as a penalty term into account. More complex models are therefore not consistently ranked better (Tian et al., 2021).
As can be seen in Table 2, the choice of a distribution function for both, the stationary as well as the nonstationary approach, is less important regarding characteristic values compared to design values. However, with quite similar AICcs of the Gumbel and Weibull distributions, large differences in design values illustrate the strong impact of the distribution type. With exception of the GEV, the choice of the window length has no significant effect on the design values in the present example. The comparably strong changes with the application of the GEV distribution are related to the curvature parameter τ, which is considered as a time-varying value in this analysis. It complies to the "shape parameter" according to Coles (2001), which may also be estimated to be constant over time for reasons of simplicity, what is not followed here. However, in the present example it is observed that larger window lengths are connected to increasing differences between the design values of the stationary and nonstationary approaches based on the Gumbel, Lognormal and Weibull distribution, whereas partial safety factors, and thus the safety distances between characteristic and design values, decrease. Particularly in case k is larger than 30 years, the results in Table 2 exhibit an expected lower snow load in the nonstationary analysis due to the decreasing trend in the time series (cf., Figure 9), though differences to the stationary approach are relatively small. It should be noted that all findings relate to the present example and are not generally applicable.
For reasons of comparison, the Gumbel distribution and the data series connected to window length of 40 years is analysed by the common nonstationary approach of including time as covariate. Here, the distribution parameters with mðtÞ and σðtÞ are chosen to be linear functions of the time. Two nonstationary models are considered, i.e. M 1 with mðtÞ ¼m 0 þm 1 � t and σðtÞ ¼ σ as well as M 2 including mðtÞ ¼m 0 þm 1 � t with σðtÞ ¼σ 0 þσ 1 � t. A maximum likelihood estimation (MLE) leads to the distribution parameters and design values presented in Table 3. The deviance statistic according to Coles (2001) gives D = 0.53 connected to a maximized log-likelihood of lðM 1 Þ ¼ À 12:58 and lðM 2 Þ ¼ À 12:31. In conjunction with a significance level of α ¼ 0:05, the X 2 1 distribution yields 3.84 > 0.53, such that a trend in σ cannot be confirmed and the model M 1 is to be preferred. As can be seen in Table 3, exhibiting smaller standard errors the shifting time window approach leads to slightly lower loads than the preferred model M 1 and to slightly higher loads compared to M 2 . Advantageously, a definition of functions for the temporal parameters mðtÞ and σðtÞ is not needed with the time window approach.

Summary and outlook
Aspects of nonstationarity may be of particular importance considering loads on structures due to global climate change. Corresponding projections and respective climate data are available as time series, which are picked up in a framework of extreme value statistics for predicting the development of structural loads within the service life of a structure. In case of a pronounced trend in the data series of an action, which is reviewed by the Mann-Kendall test in the present work, the importance of nonstationary approaches increases.
Based on approaches of Rootzén and Katz (2013) as well as Mudersbach and Bender (2017) including shifting time windows, the present work focuses on the identification of safety elements for nonstationary time series regarding the semi-probabilistic structural design of European standardization, i.e. EN 1990 (CEN (European Committee for Standardization), 1990). The procedure of the shifting time window is enhanced to directly determine design values of environmental loads. For that purpose, methods based on the First Order Reliability Method (FORM) are extended by the nonstationary approach incorporating the reliability index and the sensitivity factors. Climate projection data is represented by probability density functions, Table 3. Estimated nonstationary Gumbel-models and associated design value estimates in [kN/m]2 for a model with a linear function for the mean mðtÞ ¼m 0 þm 1 � t and a constant standard deviation σ in comparison with a model with a linear function for the mean mðtÞ ¼m 0 þm 1 � t and standard deviation σðtÞ ¼σ 0 þσ 1 � t, with standard errors ( where different typical distributions of extreme value statistics are considered and the quality of conformity using the Akaike Information Criterion with small sample size correction (AICc) is analysed. The design values and partial safety factors derived from the procedure may be applied for quantitative statements on the temporal development of the structural reliability.
To satisfy basic requirements of EN 1990, values are calibrated to ensure the required reliability level at the end of the structural service life considering presupposed parameters of the standard. In principle, presuming the availability of different climate projections or measuring data, respectively, safety elements for the semi-probabilistic concept and the reliability of a structure can be determined for relevant climatic actions.
Applied to the example of a single-span beam under snow action, it is illustrated how different probability distribution functions and assumptions of time window sizes effect the predictions of loads. The present study exemplifies that the proposed procedure allows to trace the design loads within the service life of a structure. The development of abovementioned parameters related to the structural reliability are depicted as well. For reasons of comparison, results obtained by the common nonstationary approach of considering time as covariate are also included revealing good agreement.
However, it has to be considered that climate projections are known to contain uncertainties, which, in conjunction with further assumptions of the presented procedure, as for instance, the size of the selected time window, may lead to a considerable overall uncertainty in the results. In the present work, the influence of different time window lengths on the design loads estimated by the nonstationary procedure is exemplarily illustrated.
In the future, it is expected that the overall uncertainty reduces: On one hand, the reliability of the database will continuously enhance taking new information into account, on the other hand, improvements with regard to uncertainty quantification and the choice of appropriate time window lengths will also lead to more certainty in the results. In current studies, First Order Reliability Algorithms for determining the sensitivity factors and reliability indices are incorporated in the presented nonstationary processes for that purpose as well. In addition, the presented methodology could serve for a probabilistic framework considering multiple projection scenarios with associated occurrence probabilities for estimating climate contributions.