Natural Hazards in a Changing World: Methods for Analyzing Trends and Non‐Linear Changes

Estimating the frequency and magnitude of natural hazards largely hinges on stationary models, which do not account for changes in the climatological, hydrological, and geophysical baseline conditions. Using five diverse case studies encompassing various natural hazard types, we present advanced statistical and machine learning methods to analyze and model transient states from long‐term inventory data. A novel storminess metric reveals increasing European winter windstorm severity from 1950 to 2010. Non‐stationary extreme value models quantify trends, seasonal shifts, and regional differences in extreme precipitation for Germany between 1941 and 2021. Utilizing quantile sampling and empirical mode decomposition on 148 years of daily weather and discharge data in the European Alps, we assess the impacts of changing snow cover, precipitation, and anthropogenic river network modifications on river runoff. Moreover, a probabilistic framework estimates return periods of glacier lake outburst floods in the Himalayas, demonstrating large differences in 100‐year flood levels. Utilizing a Bayesian change point algorithm, we track the onset of increased seismicity in the southern central United States and find correlation with wastewater injections into deep wells. In conclusion, data science reveals transient states for very different natural hazard types, characterized by diverse forms of change, ranging from gradual trends to sudden change points and from altered seasonality to overall intensity variations. In synergy with the physical understanding of Earth science, we gain important new insights into the dynamics of the studied hazards and their possible mechanisms.


Natural Hazards in a Changing World
Global data bases on natural hazards show an increasing trend in the number of reported events, the number of affected people, and the amount of economic losses in past decades (CRED & UNISDR, 2020).For instance, the direct economic losses as a consequence of natural hazard events have almost doubled in the past four decades, from 1.63 trillion USD in the period 1980-1999 to 2.96 trillion USD in 2000-2019.The bulk of the increasing damages is due to a growing exposure of people and economic assets in disaster-prone areas (EEA, 2017; invested in mitigation measures to reduce their vulnerability to impacts from natural hazards.For example, Thieken et al. (2016Thieken et al. ( , 2022) ) show that adapted spatial planning, comprehensive property-level mitigation, effective flood warnings, and a targeted maintenance of flood defense systems are effective measures to reduce losses from floods over time.Further, due to improved early warning systems and evacuation plans, the number of fatalities from natural hazard events has declined and societies have become less vulnerable in recent decades (CRED & UNISDR, 2018;EEA, 2017;UNISDR, 2011).Nevertheless, events with devastating damage and high numbers of victims may still occur, even in high-income countries.For example, the disastrous flood in western Germany in July 2021 (Dietze et al., 2022;Kron et al., 2022;Vorogushyn et al., 2022) caused damages of 33 billion EUR and killed 189 people (Kron et al., 2022).Despite growing data bases, it remains debated how much human perturbation of natural systems, and atmospheric warming in particular, has contributed to the increase in weather-related events and associated losses (Bronstert et al., 2020;F. E. Otto et al., 2020).The Intergovernmental Panel on Climate Change (IPCC) reports changes in the rates, intensities, and spatial extent of natural hazards with growing evidence (IPCC, 2012;IPCC, 2014IPCC, , 2019aIPCC, , 2019b)).
The above described dynamics in the natural risk components point out a shortcoming in the widely-used risk equation R = H • E • V (B.Merz et al., 2010;Plate, 2002;UNISDR, 2011), which describes risk R as a product of hazard H, exposure E, and vulnerability V (Peduzzi et al., 2001).A core assumption in this equation is that all terms are stationary over time.Yet, we know that at least the exposure (Gueneralp et al., 2015;Hallegatte et al., 2013;Jongman et al., 2012;Rohat, Flacke, Dosio, Dao, & van Maarseveen, 2019) and vulnerability to several hazards (Byers et al., 2018;Jongman et al., 2015;I. M. Otto et al., 2017;Rohat, Flacke, Dosio, Pedde, et al., 2019;Thieken et al., 2022) have changed in past decades, as did likely the hazard itself (Blöschl et al., 2019;IPCC, 2012;Kemter et al., 2020;Kundzewicz et al., 2014;Mueller & Pfister, 2011).Consequently, a more adequate form of the risk equation needs to account for transient states, where process variables are subject to change and require a time dependent evaluation of natural hazards H(t) and societal factors E(t) and V(t):

R(t) = H(t) ⋅ E(t) ⋅ V(t).
However, identifying and describing transient hazards is challenging.At this point, it is important to note the difference between the triggering natural hazard event and the damage it causes.The event that is directly related to the damage caused is often colloquially referred to as a natural disaster.The focus of this manuscript is on the changing dynamics of the occurrence of natural hazard events, that is, not on the study of the associated losses.We address the limitations encountered in natural hazard analysis and propose strategies to overcome them.In addition to presenting novel data evaluation approaches, we emphasize the importance of integrating domain knowledge to ensure model credibility.Given the inherent complexity of natural hazard events, advanced evaluation techniques and a thorough understanding of the hazard under study are essential to distinguish transient states from natural variability across temporal or spatial scales.In addition, an understanding of the hazard is essential to identify environmental drivers and their potential interactions.Overlapping or counteracting processes, coupled with limitations in data quality, can mask or distort change signals.These issues are discussed in more detail in the following section.

Natural Variability on a Temporal Scale
Natural hazards are characterized by temporal variability and random processes that control the occurrence of individual events.For example, the occurrence of earthquakes is assumed to follow a stochastic Poisson process (van Stiphout et al., 2012), and climatological data can have an inter-annual to decadal periodicity (Frei et al., 2000;Scherrer et al., 2016).Further, temperature or runoff data include characteristics of long memory processes (Franzke, 2012;Rust, 2007;Rypdal et al., 2013), which means that values measured today show a long-range dependence on values measured in the past (Hurst, 1951).This self-similarity leads to trends in time series of climatic variables, which are difficult to distinguish from climate change signals.The differentiation of non-stationary behavior from stochastic processes requires a careful analysis of reliable long-term data.

Natural Variability on a Spatial Scale
Spatial variability can be observed on different spatial scales not only in the hazard potential, but also in the trends of natural hazards.Global studies reveal geographical patterns in the development of climate extremes, such as in temperature, storms, precipitation, or droughts (IPCC, 2014;EEA, 2017;Stott et al., 2016;J. Vogel et al., 2021).Floods along major European rivers both have increased and decreased in magnitude because of atmospheric warming and subsequent altered hydrological conditions (Blöschl et al., 2019;Kemter et al., 2020).Geographical variations limit the spatial transferability of time-dependent hazard models and usually require the underlying data to be restricted to a certain area.

Superimposing and Counteracting Factors
Several interacting factors can intensify the occurrence of potentially adverse natural processes.For instance, landslide activity might be affected by various superimposing impacts of climate change (Öztürk, 2018), such as altered frequencies of landslide-triggering precipitation (Gariano & Guzzetti, 2016), glacier retreat, and permafrost degradation (Knight & Harrison, 2009), as well as shorter snow cover duration (Beniston et al., 2003).Not only the occurrence, but also the impacts of natural hazards can be either eliminated or enhanced by the interplay of influencing factors.As an example, for flash-floods Bronstert et al. (2018) describe the interacting and amplifying effects of extreme rainfall intensity, runoff formation, soil and river bank erosion, and subsequent urban-area inundation and extreme damage.By analyzing large-scale drought hazard, J. Vogel et al. (2021) detect an increase in the joint occurrence of droughts and heat waves in the Mediterranean, which exacerbates the negative consequences of the individual events.The joint or sequential occurrence of hazardous processes complicates the quantification of changing impacts (Zscheischler et al., 2018).

Limitations in Data Quality and Availability
In many cases, observation periods fall short of robustly tracking changes in the frequency or intensity of natural hazards.Monitoring processes may be interrupted in regions with disputed borders, for example, in High Mountain Asia, where data are confidential and physical access to those regions can be limited for nonresidents.Those regions often provide inconsistently filled and poorly annotated databases, and international recognition of natural hazards remains scant (Leinss et al., 2020).The Himalayas, for example, still have no permanently operating weather station above 5,000 m a.s.l., where climate-driven hazards such as flash floods, landslides, debris flows, or glacier detachments can develop (Salerno et al., 2015).In other cases, extreme events can disrupt continuous measurements.Peak discharges, for instance, can damage stream gauges or exceed their recording capacity (Korup & Tweed, 2007).Historic inventories are usually characterized by reporting gaps and biases, since only the largest events or those with disastrous consequences are documented (Wirtz et al., 2014).Reinsurance companies report a more consistent record after 1980 in their catalogs of disastrous natural events, though such data bases are likely biased towards events with insured losses (Wirtz et al., 2014).Public interests and political developments may also impact the documentation of natural hazards.Climate-related data are studied more intensively nowadays, for example, in the framework of political initiatives, such as the United Nations Sendai Framework for Disaster Risk Reduction 2015-2030 (UNISDR, 2015) and the EM-DAT International Disaster Database (CRED & UNISDR, 2020).Rising awareness and more intense reporting may increase the reported number and related losses of natural hazard events (CRED & UNISDR, 2018;Gall et al., 2009).Further, technical developments improve the conditions for data acquisition, storage, and exchange.For example, modern instrumentation and changing station density impacts on the continuous measurement of atmospheric and hydrological variables, such as air temperatures, precipitation or river discharge.It is noteworthy that in some regions an increasing density is reported, while in others the density decreases, see, for example, (Fekete et al., 2012).Yet, changes in instrumentation may result in inconsistent recordings with changes in precision, resolution, and measurement position.For instance, changes in data collection of sea surface temperatures led to erroneous warming patterns in the early twentieth century (Chan et al., 2019).Due to these obstacles in natural hazard recordings, their analysis requires not only good knowledge of statistics, but also of data acquisition, data processing, and the studied process to avoid pitfalls in the data evaluation.

Data Science for Natural Hazard Assessments
Despite these limitations, past years have seen a major advance in homogenizing global earth science and climate data sets.For example, global collections of daily air temperature date back to the 1850s (Brugnara et al., 2019) and enable more precise analyses of global warming.Reconstructing past events and homogenizing event catalogs have helped to establish continental-scale databases, for example, for river floods or tropical cyclones (Blöschl et al., 2017;Gudmundsson et al., 2021;Kossin et al., 2020).Such catalogs have allowed researchers to trace changing exceedance probabilities, or return levels, over many decades.In addition, remote sensing can help to independently estimate the timing, spatial extent, and some of the impacts of geomorphological and hydrological natural hazards, at least for the past few decades.As the ability to obtain, store, and analyze data has increased rapidly in the 21st century (Reichstein et al., 2019), more and more natural hazard assessments can be based on data-intensive algorithms.More satellite missions, seismic stations, and automatic river gauges, for example, have been producing data in amounts that exceed manual or visual interpretation (Gorelick et al., 2017;Marone, 2018).
Data-driven algorithms can help to recognize patterns and identify trends of natural processes that would otherwise have eluded human notice.For instance, the increasing number of seismic networks in combination with machine learning classifiers enables monitoring systems to go beyond solely detecting Earthquake tremors.Hammer et al. (2017) apply hidden Markov models, a probability-based sequential data analysis, to seismic data for an automatic detection of wet-snow avalanches in the Swiss Alps and provide information about avalanche activity which is independent from visual observations.Hibert et al. (2019) apply an automated classifier based on Random Forests to continuous seismic recordings and detect >5,000 landslides over a 22-year period in Alaska that had eluded scientific documentation.Cross-recurrence plots applied to runoff data helped to determine the (dis)similarity of rainfall-runoff events, and to quantify changes in the hydrological regime along the Elbe River, Germany, between 1901 and 2010 (Wendi et al., 2019).Marwan (2019) tied recurrence plot techniques to other environmental data, for example, to study transitions in past climates or the influence of external factors on ecological or climatic systems.The field of data science has brought forth machine learning classifiers such as random forests, support vector machines or neural networks that efficiently mine for structures in large data sets at high accuracy.
Another lively branch in data science uses Bayesian methods that solve, and predict from, multi-dimensional models under full probabilistic reasoning (Korup, 2021;Viglione et al., 2013).For example, K. Vogel et al. (2018) apply probabilistic graphical models, such as Bayesian networks, to detailed flood damage survey data to identify damage-driving or -reducing factors in the hazard component, the exposure, and vulnerability.A Bayesian treatment could enrich the current practice of natural hazard appraisals, because it informs practitioners clearly about uncertainties tied to expert judgment or data quality.

Potentials for Time-Dependent Hazard Models
While data science techniques continue to improve the quality of natural hazard research, the development and use of time-dependent hazard models remains relatively rare.In this study, we explore the potential of advanced statistical and machine learning techniques to identify and quantify changing hazards.To this end, we present five case studies that use rich databases and modern evaluation methods to analyze and model the evolving intensities and frequencies of different types of natural hazard events.Our overall goal is to establish the ability to characterize changes in natural hazards and project their shifting intensities under a warming atmosphere.Achieving this goal requires careful consideration of the inherent challenges, including seasonal and regional variability, interacting processes, and limited data, as outlined in Section 1.2.Each case study in this paper covers one or more of these issues and provides tools to address them.The case studies are drawn from recent publications authored by our team and were selected to show previously unknown aspects of changing hazard across a wide spectrum of hazard types.For ease of reference, Table 1 provides a brief summary of the hazard types, study regions, expected changes, methodologies used, challenges encountered, and relevant publications.
For the study at hand the examples have been tailored to adhere to the following guidelines: specify the hazard type under investigation and its anticipated changes, describe the available data and the applied evaluation technique, characterize any detected changes in the hazard, and discuss the findings within the broader context of  Earth's Future 10.1029/2023EF003553 the underlying physics.The method descriptions are restricted to the essential aspects, aiming to provide a fundamental understanding of the utilized concepts.The methodological details, for which we refer to the original studies, are not the primary focus of this paper.Instead, we approach the case studies from a new perspective to support two hypothesis: 1. Natural hazards can be subject to dynamic changes that require a time-dependent assessment of their risk potential.These changes take many forms, including gradual or abrupt shifts in frequency, intensity or seasonality.2. Advanced statistical and machine learning techniques, coupled with increased data availability, are improving time-dependent hazard assessments.However, the diverse nature of evolving hazards and their specific requirements underscore the need for close collaboration between earth sciences and data sciences to develop innovative analytical methods for reliable and comprehensive knowledge acquisition.
The presented case studies demonstrate a dynamic feedback between data science and physically oriented research branches, involving disciplines such as statistics, machine learning, seismology, hydro-meteorology and geomorphology.This interdisciplinary approach significantly enhances the selection of appropriate evaluation techniques and ensures a robust and high quality study design.

Motivation
The European Environment Agency reports an averaged annual economic loss of 4 billion EUR caused by meteorological events between 1980 and 2020 in EU member states (EEA, 2022), based on the CATDAT data set.
In 1999 the winter storm Lothar alone inflicted insured losses of approximately 6 billion EUR and caused 110 fatalities (MunichRE, 2015).European winter storms form over the North Atlantic and attain maximum wind speeds of 140-200 km/hr and even up to 250 km/hr in exposed coastal areas and higher mountains.A sound understanding of the factors that affect their genesis, trajectory, and evolution can aid in efforts to predict them.Thus far, M. Donat et al. (2011) have shown a statistically significant positive trend in the frequency of gale days and extreme wind speed percentiles, but this signal is only present in select parts of Northern Europe.Additionally, Befort et al. (2016) have found positive trends in European windstorm counts, but not in extreme cyclones.With regard to storm magnitudes, Ulbrich et al. (2009) show that parts of Europe can expect to see an increase in the frequency of extreme storms in warmer climates in the future.Windstorm magnitudes in particular were first introduced by Leckebusch et al. (2008) and have been used for a number of specific European case studies such as Pirret et al. (2017), but to date windstorm severities have not been studied systematically across a large set of storms.Here, we study windstorms which are tracked as persistent synoptic features (Leckebusch et al., 2008).Although windstorms by definition are "extreme" systems, their severity varies significantly between storms (Moran, 2021).Thus, to gain a comprehensive grasp of evolving windstorm effects, considering only their frequency isn't enough.It's crucial to also integrate their strengths, particularly concerning potential losses.Therefore, our goal is to analyze persistent trends in European winter windstorm severity from 1950 to 2010, using ERA-20C reanalysis data.

Data and Method
We use surface wind speeds from the ERA-20C reanalysis data set (Poli et al., 2016) in the period 1950-2010.This period is chosen as it is long enough to capture the decadal variability of windstorms, and the methods used to track them generalize well between data sets, for example, NOAA-20CR from 1950 onwards (Befort et al., 2016).
The procedure used to build and study windstorm tracks is described in detail in our previous publication by Moran (2021).All variables are analyzed in the extended winter season of October to March (ONDJFM), and we only consider windstorms that passed Europe in a region between 35°S-65°N and 10°W-30°E for at least one day.Windstorms are tracked using the method developed by Leckebusch et al. (2008), which identifies windstorms as moving clusters of surface wind speed anomalies over time.This algorithm also offers diagnostics of windstorm properties, such as their average surface wind speed magnitudes, the location of the windstorm center, or the area of the windstorm.Together, these properties allow us to compute the dimensionless Storm Severity Index (SSI) of a storm over its lifetime, defined as Earth's Future 10.1029/2023EF003553 VOGEL ET AL.
where k is a rectangular area of interest, t is the 6-hourly time step, and A k is the area (in km 2 ) of the data's grid box in the T159 spherical-harmonic grid.A k is normalized by a reference area A 0 , which corresponds to an area of 150,000 km 2 .We only aggregate SSI throughout the windstorm's lifetime T and, at each time step, the total area covered K.The SSI accounts for scaled surface wind speeds v k,t that exceed the local 98th percentile v k,98 , where the percentile follows the assumption that storm damages occur on 2% of all ONDJFM days of the time period (Palutikof & Skellern, 1991).From this, we can construct the annual dimensionless storminess metric by summing up the SSIs of each windstorm in each extended winter season.We choose storminess over windstorm counts because storminess incorporates the severity of the storms in a given year.As storm severity is expected to change with a warming climate (Ulbrich et al., 2009), we compare linear trends and R 2 values between the counts of the windstorms and the total integrated storminess.
Storminess is correlated with windstorm counts, with a value of 0.71 (significant to p < 0.05 using a t-test), yet, years with peaks in storminess do not necessarily correspond with peaks in windstorm counts, which we observe, for example, in 1989, 2004, 2007, and 2008.Thus, storminess provides a different picture of the potential impacts of windstorms in a given season.While a year might appear to have lower than average numbers of storms, it could have a high storminess, as a few individual significant windstorms can have a much greater impact.This demonstrates that storminess is an advantageous metric to use over the conventional windstorm counts.For a comprehensive results analysis, please see Moran (2021).

Summary and Conclusion
This study analyzed extreme windstorm events over Europe.Such events are of high relevance in particular for built-up areas and forested areas, and have created an enormous economic loss as well as a considerable number of fatalities across Europe.In the warming climate, it is of high importance whether the number and severity of such events might increase.In this regard, we analyzed the development of the frequency and strength of such events.We applied a linear trend analysis, which is the least complex form that can capture gradual monotonic change.The underlying data set encompasses available meteorological re-analysis data for Europe and the Northern Atlantic for the period 1950-2010.We constrained those data to the winter period (October-March), because this is the windstorm season over Europe.In general, we found that the winter storms increased over Europe.This is particularly visible when analyzing -rather than just the storm counts -the storminess index, which includes the severity of the storms in a given year and integrates the counts with their surface wind speed anomalies.The storminess shows a stronger trend compared to the windstorm count.Thus, the attribution of potential windstorm hazards is clearer when using the storminess metric instead of windstorm frequency.This can be used to give stakeholders in reinsurance or public policy more accurate information in terms of how impactful windstorms will be today and in the near future.This method can also be adapted to seasonal or decadal timescales, as well as to other regions.The method is also applicable to wind information gained from climate model simulation and therefore applicable for information regarding wind conditions in a future climate.

Motivation
Extreme precipitation can have adverse impacts on societies, such as harvest losses or increased accident risk in road traffic.Furthermore, heavy rainfall can trigger cascading hazards, such as riverine floods, flash floods or landslides.With ongoing atmospheric warming, extreme precipitation intensity events could become more frequent and more intense (Bürger et al., 2019;Bürger et al., 2021;IPCC, 2012;M. G. Donat et al., 2016), given that a warmer atmosphere has a higher storage capacity for water vapor (Vergara-Temprado et al., 2021).Despite decades of rainfall measurements, notions of historic trends in extreme rainfall remain vague, possibly due to a lack of methods to quantify such changes.The IPCC AR5 (IPCC, 2013), for example, concludes that "are precipitation (such as the highest annual daily precipitation total) events were likely to have increased over regions with sufficient data since the late nineteenth century" (p.213).
Besides inter-annual changes in intensities and frequencies, extreme precipitation might also change seasonally (Konapala et al., 2020).Alexander et al. (2006) found that indices of extreme daily precipitation suggested a global increase in extreme precipitation from September to May.Yet, there was no trend detectable during summer (June-August).These results are consistent with the findings from Zolina et al. (2008), who analyzed seasonal trends in the 95th and 99th percentiles of daily precipitation totals in Germany.They found positive linear tendencies for winter, spring, and autumn, though summer had mostly negative trends.However, the trends were in part neither significant nor coherent in space.Looking at sub-daily time-steps, that is, hours or even shorter, analyses in Germany, Austria, and Switzerland have revealed that high-intensity but short rain storms have become more frequent in the past 35 years according to rain station data (Bürger et al., 2014(Bürger et al., , 2019;;Mueller & Pfister, 2011).The trends in intensity are increasing the most between July and September.
The frequency of such high-intensity rainstorms may have doubled or even tripled over a period of a few decades.A change in the seasonal cycle of extreme rainfall might change flood hazard, such that periods of extreme precipitation overlap with high water levels in rivers from snow melt (Berghuijs et al., 2019;Vormoor et al., 2015).Furthermore, a shift in seasonality could also affect crop yields, as plants are sensitive to changing precipitation during early growing stages (Derbile & Kasei, 2012;Rosenzweig et al., 2002).Understanding changes in the seasonal cycle of extreme precipitation is thus vital for a better adaptation to the impact of climate change.Here, we show an approach to quantify temporal and seasonal changes in extreme precipitation.

Data and Methods
We adopt the block maxima approach from extreme value theory (Coles, 2001) to describe the magnitude and occurrence probability of extreme precipitation.This method extracts the maximum precipitation within a given block.For this analysis, we use a block size of 1 month.The probability distribution of the block maxima can be described by the generalized extreme value distribution (GEV) 1/ξ } with {z: 1 + ξ(z μ)/σ > 0} and the location parameter μ, the scale parameter σ and the shape parameter ξ.Instead of fitting a GEV with independent parameters to each month, we model the seasonality using harmonic functions in the form of linear models within the parameters of the GEV (Maraun et al., 2009;Mínguez et al., 2010;Rust et al., 2009).This approach offers more reliable monthly return levels and reduces model complexity.Simulations showed that including the seasonal cycles also offers more robust estimates of annual return levels (Fischer et al., 2018).To address inter-annual changes of seasonality, we extend the approach by a seasonal-temporal component (Peter et al., 2024).Conceptually, the time-dependency of the three GEV parameters in our model is represented by the following equation: where θ year, month denotes the location, scale, or shape parameter for a specific year and month.The model includes a constant offset θ 0 , followed by three approximating functions, f 1 , f 2 , and f 3 , each containing a different set of predictors.To capture the seasonal cycle, we utilize a series of harmonic functions for f 1 .For example, the first harmonic order can be represented as , where ω = 2π/365.25 is the angular frequency, c t the center day of month t = 1, .., 12, and a and b are coefficients to be estimated.The terms sin(ωc t ) and cos(ωc t ) serve as the predictors.To account for the changes in extreme precipitation over time, we employ Legendre polynomials up to order five on the years for f 2 .The interaction term f 3 captures the change of the seasonal cycle with the years and is a product of the predictors for intra-and inter-annual variations.This term helps to account for any joint effects between the monthly and yearly variations.We refer readers to Peter et al. (2024) for detailed expressions and descriptions of the seasonal-temporal approach and restrict the paper at hand on presenting the overarching concept.Analogous to the seasonal-temporal model, we present a seasonalspatial model in Fischer et al. (2019) to incorporating spatial variability of extreme precipitation.The seasonalspatial approach enables return level estimation at ungauged sites and can further enhance accuracy.
We apply the seasonal-temporal model to monthly maxima from 519 long time series of daily precipitation sums in Germany obtained from the German Weather Service (Deutscher Wetterdienst, DWD (2022)).The stations are roughly equally distributed in space across the country and each of it covers at least the period 1941-2021.To select appropriate seasonal-interannual models for each station, we employ a step-wise forward regression approach using the Bayesian Information Criterion (BIC).This method allows us to identify the necessary predictors that best explain the variations in the data.After determining the final model setup, we estimate the coefficients for the model.Once the model is established, we use it to calculate return levels (quantiles of the GEV distribution for a given non-exceedance probability p) for different return periods.The return levels represent the expected precipitation magnitude that is exceeded, on average, once within a specific return period T, where T = 1/(1 p).

Results
Implementing the aforementioned method uncovers that 65% (n = 338) of all evaluated stations encountered a historic change in daily extreme precipitation, as we comprehensively elaborate in Peter et al. (2024).Most of these stations show a change in the seasonality ( f 3 ), in the form of an altered amplitude or a shifted phase, rather than a long-term change in heavy rain over time ( f 2 ) with a consistent change for all months.Analyzed linear trends are in general weak, only 15%-35% (depending on return period and month) show changes of more than 5% from 1941 to 2021.A rise in return levels dominates especially in June.The changed seasonality is mainly characterized by an increased amplitude, later extreme precipitation events for shorter return periods (<10 years), and earlier extreme events for longer return periods (>10 years).Although changes are spatially very different, regions with similar and more pronounced properties stand out: one in southern Germany exemplified by the station Rain am Lech and one in the center of Germany represented by the station Wesertal-Lippoldsberg.Those observations and 100-year return levels are shown in Figure 2 (left and mid).The station Mühlhausen/Oberpfalz-Weihersdorf, as well located in the south of Germany, illustrates a phase shift to earlier times in the year (Figure 2, right).
We refer to Peter et al. (2024) for a more comprehensive and detailed analysis of the data evaluation.For the current paper, we focus on the stations presented in Figure 2, to illustrate different forms of changing extreme precipitation in Germany.At Rain am Lech the threat due to rare extreme daily precipitation has increased in all months of the year.The 100-year return level has risen about 58% from 1899 to 2021.The seasonal cycle at the station Wesertal-Lippoldsberg shows a changing amplitude with intensified events in summer and moderated in winter.This record benefits from the inter-annual modeling (see Peter et al. (2024) for model verification), although summerly return levels for the recent years are overestimated.Altered return levels in winter could be attributed to changes in large-scale atmospheric circulations (Willems, 2013) or changed temperature and block situations (Fauer & Rust, 2023), while a rise of return levels in summer may indicate more intense convective thunderstorms (Zeder & Fischer, 2020), accompanied by lightning, hail, and high rainfall intensity.We call for adapting hydraulic structures in all regions with increasing return levels, regardless of whether these changes are only seasonal or occur throughout the entire year.The Mühlhausen/Oberpfalz-Weihersdorf station shows a shift in the peak of the seasonal cycle by about 35 days from early August (1931) to the end of June (2021).Shifts in seasonality to more vulnerable times (e.g., coincidence with snow melt or growing stages of plants) can lead to potentially increased impacts on society.

Summary and Conclusion
This study analyzed extreme precipitation data regarding precipitation intensity.High-intensity rainstorms are of great relevance for a variety of hazards, such as local-scale flash floods, where they cause great damage particularly in urban areas or destroy crops and erode fertile soils in agricultural fields.We have analyzed Earth's Future 10.1029/2023EF003553 temporal -including seasonal -changes of extreme precipitation in Germany.We used a non-stationary generalized extreme value model to analyze the magnitude and occurrence probability of daily extreme precipitation intensity over Germany.The seasonal cycle was modeled with harmonic functions.Based on daily rainfall data obtained from the German Weather Service (DWD) from 519 rainfall stations, we processed a comprehensive data set for Germany of daily rainfall data for 80 years  or more.To summarize, we found in 65% of all stations a historic change in extreme precipitation.Most of these stations show a change in the seasonality and an intensification dominates especially in summer.However, trends are seasonally and spatially different, which is why a detailed spatial and seasonal analysis is recommended.For the time being, the detected trends are still weak, but the method gives consistent results and can be applied in the same manner for data, including the observations of the current years, and for higher resolved duration (e.g., hourly sums).The given approach needs an adaptation so that future return levels can be predicted in a meaningful way, since those are highly relevant for adaptation measures to climate change in urban hydrology and agriculture.Additionally, implementing other atmospheric covariates such as the global mean temperature or the North Atlantic Oscillation could offer a more physically motivated explanation of the changes in extreme precipitation.

Motivation
Rivers spill over their banks occasionally, posing a threat to residential areas and cities.Floods account for large numbers of casualties and economic losses amounting to billions of US$ each year (UNDRR, 2019).The 2002 Elbe flood (Kreibich et al., 2005;Ulbrich et al., 2003), the rain-on-snow flood in the Bernese Alps, Switzerland, in October 2011 (Rössler et al., 2014), and the devastating flood in western Germany and Belgium in July 2021 (Kron et al., 2022;Thieken et al., 2022;Vorogushyn et al., 2022) are only three recent examples of extremely damaging and costly floods in Central Europe.In general, extreme precipitation, excess soil moisture, and snowmelt largely control the timing, magnitude, and frequency of riverine floods (Berghuijs et al., 2019;Frei et al., 2000;R. Merz & Blöschl, 2003;Vormoor et al., 2016;Wetter et al., 2011).Changes in characteristics of river floods have been studied intensively in the past years.For example, Petrow and Merz (2009) found that flood hazard along large German rivers increased in winter, though not in summer, between 1951 and 2002.A larger sample of 4,262 hydrometric stations installed in Europe between 1960 and 2010 showed that atmospheric warming has been affecting the timing and magnitudes of river floods (Blöschl et al., 2017(Blöschl et al., , 2019)).Regional warming causes earlier snowmelt-induced floods across northeastern Europe in spring, and increasing autumn and winter precipitation has generated increasing flood magnitudes in northwestern Europe.Vormoor et al. (2015) project that flood seasonality in Norway could change such that rainfall replaces snowmelt as the dominant flood-generating process.In general, rising temperatures have been depleting seasonal snowpacks in recent decades (Laternser & Schneebeli, 2003;Marty, 2008;Scherrer et al., 2004).Projected warming is expected to further diminish alpine snow cover and reduce the seasonal redistribution of water from winter to summer via the build-up and melt of seasonal snowpacks (Hanzer et al., 2018;Marty et al., 2017).Furthermore, thinner snowpacks lower the potential of meltwater floods (Musselman et al., 2017).In addition, changes in the seasonality and peak values of extreme precipitation (see Section 3) could cause seasonal shifts in flood hazard.
Besides the possible impacts of climate change, the consequences or land-use change, modifications of the river network, or water resource management may significantly alter river runoff and hydro-extremes (Bronstert et al., 2007(Bronstert et al., , 2023;;Niehoff et al., 2002;Pfister et al., 2004).Particularly in high mountain regions, reservoirs for hydropower production affect river runoff, redistributing river runoff from summer to winter and thus decreasing runoff seasonality (Arheimer et al., 2017;Bosshard et al., 2013;Meile et al., 2011;Pérez Ciria et al., 2019;Verbunt et al., 2005).Throughout the twentieth century, and since the 1950s in particular, most rivers spawning from the Swiss and Austrian Alps have been interrupted by hydropower dams (Wildenhahn & Klaholz, 1996).

Data and Methods
The existing knowledge on the timing and magnitudes of floods in river basins largely hinges on time series of gauged discharge.River discharge data offer the longest available hydro-climatological record, for some locations dating back to the nineteenth century.Distinguishing between natural variability and long-term climate trends demands very long and high-quality recordings (Auer et al., 2005;Begert et al., 2005;Rottler et  estimating climate-driven changes in flood magnitudes.However, we call for analytical tools that can resolve non-linear responses in flood discharges to understand long-term changes in river runoff and underlying mechanisms. We investigate discharge time series from four gauging stations located in Central Europe.The stations are part of the hydrometric observation networks managed by the water authorities in Germany and Switzerland.Their recordings are regularly checked to ensure high quality and reliability.The selected gauges offer a wealth of daily-resolution discharge data dating back to at least 1869, which can be obtained from the Global Runoff Data Centre (GRDC, 2020).The depicted gauges stand out for the exceptional length of their records and represent different types of flow regimes: nival, pluvial, and complex.Specifically, the Wasserburg gauge, situated on the Inn River (a Danube tributary) in Upper Bavaria, Germany, and three other gauges located in the Rhine River basin -Basel (Switzerland, representing the nival regime), Würzburg on the Main river (a Rhine tributary, typical of a pluvial regime), and Cologne (on the Lower Rhine, exhibiting complex regime).For the purpose of this paper, we exemplarily focus on the Basel gauge, which represents the alpine-dominated part of the Rhine catchment.This region experiences seasonal snow cover, making it particularly relevant for assessing climateinduced changes in snowmelt, which could significantly impact river runoff (Addor et al., 2014;Hanzer et al., 2018;Stewart, 2009;Viviroli et al., 2011).Comprehensive investigations of the other river gauges, along with high-quality temperature and precipitation data from nine Swiss climate stations can be found in our preceding study (Rottler, Francke, et al., 2020).The inclusion of temperature and precipitation records in our analysis not only supports the evaluation of the presented methodology but also aids in attributing detected changes in runoff to relevant climatic factors.
We provide a comprehensive description of the applied methodology in Rottler, Francke, et al. (2020), while in this paper we focus on presenting the general concept of the methods employed.Specifically, we utilize a complete ensemble empirical mode decomposition (EMD) with additive noise to our data.EMD decomposes a time series into oscillatory modes, effectively separating short-term signals from underlying (non-linear) trends (Huang et al., 1998;Luukko et al., 2016;Wu et al., 2007).This approach offers distinct advantages over commonly used linear methods, which often lack physical justification.By employing EMD, we attain a more flexible characterization of trend signals, as it does not require a predetermined basis function.In Rottler, Francke, et al. (2020) we provide a quantifiable comparison between the non-linear signals detected through EMD and the commonly used linear trend analysis.To ensure a highly resolved examination of longterm changes in daily resolution hydro-climatic data, we couple EMD to quantile sampling and moving time windows.Our investigation spans the time frame from 1869 to 2016.Extracting quantiles from discharge time series on a daily basis provides valuable insights into the seasonal distribution of river runoff (runoff seasonality).Runoff seasonality specifies the difference in mean runoff between seasons.Hence, low runoff seasonality describes a relatively even distribution of runoff throughout the year.High runoff seasonality points to significant variations, such as low runoff in winter and high runoff in summer.In order to assess the onset and evolution of changes, we apply EMD on a daily basis and use the resulting residuals to represent the underlying signal of change in the data (Luukko et al., 2016;Rottler, Francke, et al., 2020).To make results of different days comparable, each residual is centered by subtracting its mean.Consequently, negative/positive values indicate lower/higher runoff values compared to the average during the period from 1869 to 2016.We extend the approach to annual quantiles of discharge to investigate changes in low or high runoff characteristics (changes in quantiles).

Results
Results for the High Rhine demonstrate that in snow-dominated river basins, runoff is low in winter, when most of the precipitation is solid and stored in temporary snowpacks.Most runoff is generated in warm seasons, when snow melts and precipitation is liquid, even at rather high elevations (Figure 3a).In general, snow-dominated river basins show a strong runoff seasonality.In recent decades, runoff has increased during winter and spring and decreased from July to October (Figure 3b).Consequently, the runoff seasonality is decreasing.Accordingly, quantile values of low probability levels (<0.6) are increasing and quantiles of high probability levels (0.6-0.85) are decreasing.However, values from the highest quantiles (>0.85) have increased again in recent decades (Figure 3c).
One factor potentially altering the seasonality of alpine rivers is the construction and management of reservoirs for hydropower.In Switzerland, first hydropower stations were built in the nineteenth century and almost all storage volume currently in use in the High Rhine has been made available until the mid-1970s.In accordance with Pérez Ciria et al. ( 2019); Brunner and Naveau (2023), our results indicate that the construction of reservoirs for hydropower productions need to be considered to explain decreasing runoff seasonality in the High Rhine over the time frame investigated (Figure 3).Building upon the weather data analyzed in (Rottler, Francke, et al., 2020), we further suspect that the increase in high quantiles of river discharge in the High Rhine (Figure 3c) predominantly results from higher precipitation rates, and less so from reservoir constructions or changes in snowmelt rates.Recent studies show that heavy precipitation events have become more intense and that precipitation has increased on average, particularly during winter (Frei et al., 2000;Lehmann et al., 2015;Mueller & Pfister, 2011;Rottler, Francke, et al., 2020;Schmidli & Frei, 2005).In addition, an upward shift of the snowline suggests that a larger area of the basin receives liquid instead of solid precipitation (Allamano et al., 2009a(Allamano et al., , 2009b;;Barnett et al., 2005;Beniston et al., 2018;Gobiet et al., 2014;Rottler, Bronstert, et al., 2020;Zubler et al., 2014).

Summary and Conclusion
This study analyzed long-term changes in the runoff of rainfall-and snowmelt-dominated river basins.As both the seasonal timing and severity of either rainfall-or snowmelt-induced floods depends on the air temperature in the river basin, changes to floods in such hydro-climatological environments might occur as a consequence of regional warming.We used publicly available observations of river runoff, temperature and precipitation covering almost 150 years to quantify changes and assess underlying forcing mechanisms.In this chapter, we present results for gauge Basel.The full analysis and further details are available in (Rottler, Francke, et al., 2020).Our approach's novelty lies in the combination of non-linear trend detection, moving average trend statistics, and quantile sampling, enabling us to separate short time signals from underlying trends in highly resolved (daily), very long, and high-quality hydro-climatic data.A sequence of analytical steps extracts high-resolution signals of long-term changes and provides a comprehensive picture of existing trends.The consideration of non-linearity in changes over time facilitates the attribution of changes.
Drawing on the findings of Pérez Ciria et al. ( 2019); Brunner and Naveau (2023) we propose that the construction and operation of reservoirs for hydropower production might be a contributing factor to the altered streamflow characteristics in the High Rhine region, leading to a control effect on river runoff and a reduction in runoff seasonality.This aspect is discussed in more detail in Rottler, Francke, et al. (2020).The cited study also includes detailed investigations of precipitation patterns and indicates an increase in (intense) rainfall, which might contribute to the observed increase in high river discharge quantiles.Furthermore, we suspect that an upward movement of the snowline, driven by rising temperatures, results in larger areas of the basin to receive rainfall instead of snowfall.As a consequence, the importance of snowmelt as a flood-generating process will further decrease, while rainfall-triggered flood events could become more severe and frequent.The investigation of downstream river gauges (i.e., Cologne) suggests that alpine signals can propagate downstream and affect runoff in river segments far beyond the alpine area.

Motivation
The Himalayas host the largest reservoir of ice outside the polar regions.Yet, glaciers in the Himalayas have thinned, retreated, and lost mass since at least the 1970s (Bolch et al., 2019), and rates of ice loss have accelerated since the 2000s (Maurer et al., 2019).Ongoing glacier retreat has created space for new glacier lakes, and existing lakes have been growing from sustained meltwater supply.Glacier lake abundance and area in the Himalayas increased by 11% and 15%, respectively, between 199011% and 15%, respectively, between and 201811% and 15%, respectively, between (Wang et al., 2020)).Lakes dammed by abandoned moraine dams account for one quarter of all (n = 8,300) Himalayan lakes (Maharjan et al., 2018), yet they had some of the highest growth rates in this period (Nie et al., 2017).Moraine dams are prone to failure because they are made of loose, poorly sorted sediment, possibly accommodating ice cores (Clague & Evans, 2000;Richardson & Reynolds, 2000).Steep slopes adjacent to moraine-dammed lakes could also trigger ice avalanches, rockfalls, and landslides, which could impact lakes and produce splash waves that overtop the lakes (Haeberli et al., 2017).When glacier lakes fail, several million cubic meters of meltwater and sediment can drain in devastating glacier lake outburst floods (GLOFs) (Osti & Egashira, 2009;Watanbe & Rothacher, 1996;Xu, 1988).Given ongoing atmospheric warming, glacier retreat, and lake growth, GLOFs are projected to become more frequent (Harrison et al., 2018).Yet, trends in the historic activity of GLOFs remain controversial (Carrivick & Tweed, 2016;Harrison et al., 2018;Nie et al., 2018;Richardson & Reynolds, 2000).Many outbursts may have gone unnoticed in Himalayan headwaters, limiting our knowledge on GLOFs to larger, destructive cases (Komori et al., 2012).Such inconsistencies in GLOF databases call for a more complete GLOF inventory to estimate changes in GLOF frequency, and hazard, eventually.

Data and Methods
We developed a processing chain in Veh et al. (2018) to detect previously unrecorded GLOFs from time series of Landsat satellite images.To systematically fill the GLOF inventory in the Himalayas, we first trained a Random Forest classifier (Breiman, 2001) that predicts land cover for >8,000 post-monsoon images (September to November) between 1988 and 2017.Random Forest is a popular ensemble learning classifier in remote sensing because this method combines many hundreds of independent decision trees by using bootstrapped samples of land cover classes (Belgiu & Drăguţ, 2016).Random feature selection within the model reduces overfitting, while improving accuracy.We defined five land cover classes, manually labeled 12,500 training pixels with one of these classes, and trained the Random Forest model.We aggregated the votes of the individual trees and obtained an overall accuracy of 90% of the model.This high accuracy is suitable to predict the land cover for each pixel of each Landsat image in our time series.The Landsat mission has a short revisit time of 16 days, and thus allows us to detect lakes that shrank between two successive land cover maps.The partial or complete loss of water bodies is one of two key diagnostics for outbursts, judging from the 17 previously reported cases that we collated in this period.Based on this observation, we developed an automated change point algorithm that explores the stack of land cover maps for pixels that shifted from water to land (Veh et al., 2018).A second key diagnostic for GLOFs are sediment fans below shrunken lakes.We visually assessed all sites of predicted water loss for downstream sediment aggradation to complete our search for previously unrecorded GLOFs.
We then merged the inventory of previously reported cases with our newly detected GLOFs (see Tables S1 and S2 in Veh, Korup, von Specht, et al. (2019)).We tested whether this updated inventory shows a trend in GLOF frequency over time using a Bayesian linear regression.The noise model for our target variable annual GLOF frequency follows a student-t distribution.This distribution allows for heavier tails than the Gaussian distribution that is conventionally used as the likelihood function in linear regression models.The robustness is achieved by accounting for a low number in the degrees of freedom ν in the student-t distribution (here ν = 4), which controls the shape, and thus the probability mass, in the tails of the distribution (Kruschke, 2015).Our estimates of trends in mean GLOF rates can thus account for potential effects of outliers in the response, that is, balancing years that had many or no GLOFs.We also binned the Himalayan arc into swaths of 50-km width to assess whether regional GLOF abundance correlates with glacier area (Arendt et al., 2017), glacier mass change (Brun et al., 2017), or changes in surface water (Pekel et al., 2016).
We chose the 100-year return level of outburst volume V 0 or peak discharge Q p as an objective metric of GLOF hazard.Return levels can be interpreted as the exceedance probabilities in an extreme value distribution.To determine the 100-year GLOF discharge, we chose the Generalized Pareto (GP) model fitted to time series of GLOF discharges.We simulated such GLOF discharge time series by assuming randomly occurring events.The parameter of the respective Poisson distribution corresponds to the observed regional GLOF frequency λ (i.e., number of GLOFs per year).Outburst volume and peak discharge are based on empirical estimates of V 0 and Q p for each lake in a given region (Cook & Quincey, 2015;Maharjan et al., 2018;Walder & O'Connor, 1997).The GP model then takes as input all values of simulated V 0 or Q p above a high threshold in a given time series of regional discharges.We showed in Veh et al. (2020) that our GP models are robust for the 20% highest simulated discharges and used this threshold to eventually estimate GLOF return periods.Accordingly, the 100-year GLOF is the value of V 0 or Q p that has a 1% chance to be reached or exceeded in the GP model, which is a widely used characteristic in river flood hydrology (Katz et al., 2002).

Results
Our automated search in Landsat images reveals 22 newly detected GLOFs, as detailed in our publication by Veh, Korup, von Specht, et al. (2019).We thus more than doubled the previously known GLOF count in the Himalayas since 1988 (Figure 4a).Yet, the total of 38 GLOFs shows no change in the annual frequency (Veh, Korup, von Specht, et al., 2019): the trend in our study period is 0.0006 +0.0415 / 0.0394 GLOFs yr 1 (mean and 95% highest density interval).Given increases in lake size and abundance, the activity of GLOFs per unit of glacial lake area has even decreased in the past 30 years.This finding could point at some resilience of glacier lakes to sudden regions, and in the Eastern (Western) Himalayas, where regional GLOF hazard is highest (lowest).Blue numbers are the mean 100-year return level of peak discharge Q p including the 95% highest density interval; brown numbers are the mean GLOF rate in this region.Brown ticks are historical GLOFs.The GLOF from lake Zhangzangbo in 1981 had a peak discharge of ∼16,000 m 3 s 1 (Xu, 1988) and a 100-year return period, accordingly.outbursts.For example, lakes could have decoupled from their retreating parent glaciers and thus could have become less prone to impacts from their calving tongues.
A spatial cluster of GLOFs has formed in the Central and Eastern Himalayas (Bhutan and Eastern Nepal), in contrast to the less frequently affected regions in the northern Hindu Kush-Karakoram ranges.We provide a detailed presentation of our swath analysis in Veh, Korup, von Specht, et al. (2019).It suggests that the regional GLOF clusters correlate with both the abundance of meltwater areas (Spearman's correlation coefficient ρ = 0.87; Figure 4a) and the change in lake areas, calculated as the difference between two epochs before and after the year 2000 (ρ = 0.77).Glacier cover (ρ = 0.18) and changes in mass balance (ρ = 0.26) have little association with the number of observed GLOFs.We infer that a higher abundance and variability of glacier lakes warrants a higher GLOF count, but not necessarily an increasing rate of lake outbursts.
We chose a stationary extreme-value model to estimate GLOF return periods, given that the GLOF frequency has remained constant over our study period.As explicated in Veh et al. (2020), our investigation yielded a contemporary 100-year GLOF volume V 0 of 33.5 +3.7 / 3.7 × 10 6 m 3 (mean and 95% highest density interval), along with a 100-year GLOF peak discharge Q p of 15, 600 +2,000 / 1,800 m 3 s 1 across the entire Himalayas (Figure 4b).The regional variations of GLOF frequency and lake sizes are key determinants in our hazard appraisal.The 100-year return level of Q p in the Eastern Himalayas is at least triple that of any other Himalayan region, because the Eastern Himalayas have the largest lakes and the most historic GLOFs.Accordingly, GLOF hazard is smallest in the western Himalayas, where GLOFs are rare and lake cover is smallest (Figure 4b).

Summary and Conclusion
GLOFs are extreme floods triggered by a geophysical process, such as a landslide, earthquake or dam failure.These floods transport large amounts of water and debris, and may have severe consequences for the river reaches and population downstream.It has been argued that the occurrence of GLOFs may increase because of regional warming, due to a rising number of glacial lakes, reduced stability of slopes above such lakes and dams forming the lake, and more frequent ice avalanches discharging into such lakes (Haeberli et al., 2017).Yet, quantifying trends in GLOF occurrence has been challenging as these floods are rarely reported, while others may drain without reporting in uninhabited valleys.This study more than doubled the number of known GLOFs in the Himalayas using an automated processing chain of optical satellite images (more than 8,000 images between 1988 and 2017).Potential GLOFs were identified by applying a Random Forest classifier for the land cover prediction and a change point algorithm to detect shifts from water to land.Yet, a Bayesian robust linear regression applied to the total of 38 previously recorded and newly detected cases did not reveal a change in the annual GLOF frequency.GLOFs occurred in regional clusters, which correlate with the abundance of meltwater areas and the change in lake areas over the last decades .This spatial clustering leads to large regional differences in the calculated GLOF return periods.The three stages of GLOF analysis proposed here (from GLOF detection to analyzing their frequency and estimating regional GLOF return periods) present a suitable framework for modern GLOF hazard assessment.Given a rapidly growing population, infrastructure, and hydropower projects in the Himalayas, this approach complements local expert-based hazard and risk appraisals.

Motivation
Earthquakes and co-seismic events such as tsunamis or landslides rank among the most devastating natural hazards worldwide.Earthquakes have caused a cumulative economic loss of 661 billion USD and 747,000 fatalities between 1998 and 2017 (CRED & UNISDR, 2018).Seismicity, for example, at active faults, is assumed to be resilient against anthropogenic impacts.However, earthquakes may be triggered, delayed or advanced in response to changing boundary conditions: stress transfer from adjacent regions, for example, due to plate tectonics, movement of magma, or tides (C.H. Scholz et al., 2019), can bring a fault closer to, or away from, failure.Man-made stresses from mining, geothermal exploration, or crustal load from large hydropower dams may also initiate seismicity.The same applies to the presence of fluid flows in Earth's crust.Stress from fluids include natural causes such as increased pore pressure after rainfall (Hainzl et al., 2006), as well as man-made intrusions such as wastewater injections into deep wells.The majority of this wastewater originates from conventional hydrocarbon deposits and is separated during oil production.Moreover new drilling techniques for oil and gas exploitation have shed light on injection-caused earthquakes.Hydraulic fracturing ("fracking") routinely produces micro-earthquakes below M = 2, though an increase in earthquakes with higher, potentially damaging magnitudes could not be proven (Ellsworth, 2013).However, wastewater injection into deep wells has produced a notable increase in harmful earthquakes.For example, the US state of Oklahoma has attained public and political interest due to a rapid increase of seismicity since the onset of high-volume wastewater injections in 2009 (Gupta & Baker, 2017;Langenbruch et al., 2018;Montoya-Noguera & Wang, 2017), including five events of M ≥ 5 in Oklahoma (Ellsworth, 2013;Langenbruch & Zoback, 2016).To prevent further, possibly larger earthquakes, injection volumes had to be reduced by law in 2016 (Baker, 2017).
Changes in the stress regime may not only affect the number of earthquakes, but also the ratio of large-magnitude to small-magnitude earthquakes (Amitrano, 2003;Goebel et al., 2013;C. Scholz, 1968).A decreasing ratio of large to small magnitudes can indicate an increasing stress level in Earth's crust (Urbancic et al., 1992;Wyss, 1973), for example, as a consequence of human activity.A reliable detection of changed seismicity can be challenging due to small sample sizes, binned data, and natural fluctuation.Here, we present a Bayesian approach to identify change points in earthquake numbers and magnitude distributions based on seismic catalogs.The approach is applied to seismic catalogs from Oklahoma and southernmost Kansas (US) to substantiate a rapid increase in seismicity in connection with the growing number of industrial projects and associated fluid injections.

Data and Methods
We apply a change point model to examine the footprint of fluid injections in seismic catalogs from Oklahoma and southernmost Kansas (US).The studied earthquake catalog comprises events with magnitude M ≥ 2 between 1 January 2005 and 31 December 2019, obtained from the Oklahoma Geological Survey (OGS, 2021).The distribution of earthquakes is assumed to follow a Poisson process.The count N(M) of earthquakes of magnitude M or greater is typically described by the Gutenberg-Richter relationship: log here, M c is a threshold above which the earthquake catalog is considered complete.The parameter a represents the total seismicity rate, and b characterizes the proportion of smaller to larger earthquakes within a specific region.
We focus on the temporal change of the two key parameters a and b and compare two model setups to investigate earthquake occurrences in a given catalog: M 0 assumes that a and b remain constant over time.In contrast, M 1 expects one change point with credible differences in a or b before and after the change point.We test the hypothesis of no versus one change point using the Bayes factor, B 01 .The Bayes factor relates the probability of observing the recorded data under the assumption of no change point, p(d|M 0 ) , to the probability of observing the data under the assumption of one change point, p(d|M 1 ) , Due to the fact that the Bayes factor quantifies the evidence of the supported model and strongly depends on the choice of the prior distributions, the Bayes factor threshold is different for the parameters a and b.Synthetic tests for a suggest a threshold for the Bayes factor at B = 0.05, whereas B = 0.5 is a reasonable choice for the b-value.This means, that for example, for the parameter a, B 01 < 0.05 can be interpreted as a decisive evidence against a model without a change point.The concept can be extended to models with several change points.Here, the Bayes factor B lk = p(d|M l ) p(d|M k ) is used, where M l and M k denote a model with l and m change points, respectively.We provide a detailed description of the extended approach in Fiedler, Hainzl, et al. (2018).
In the second step, we use the likelihood function to determine the most likely position of the change point and check its significance (α = 0.05).This approach avoids smearing effects of the commonly applied moving window approach (Kulhanek, 2005;Nuannin et al., 2005), where the b-value is estimated from sliding time windows with predefined window length and step size.Subsequently, the detected change points are visually compared with the onset of fluid injection wells (OCC, 2021), and a correlation analysis is performed.
The refined location of changes in the magnitude-frequency distribution can be used not only to analyze humaninduced earthquakes, but also to detect natural transitions in the stress regime.For example, in Fiedler, Hainzl, et al. (2018) we study possible precursor signals of large magnitude events.Decreasing b-values have been observed before major earthquakes (Schorlemmer et al., 2005;Spada et al., 2013;C. H. Scholz, 2015), but the value of changing b-values to predict large earthquakes is debated (Imoto, 1991;Kamer & Hiemer, 2013;Nakaya, 2006;Nanjo et al., 2012;Smith, 1981).

Results
Figure 5 shows a considerable increase in high-volume injection wells in Oklahoma after 2010.This onset is accompanied by change points toward an increasing seismicity rate, starting with a few change point locations in 2010 and leading to a high concentration of change points toward increased seismicity between 2013 and 2015.We analyze the correlation between the change points and the onset of injection wells in Fiedler, Zöller, et al. (2018) and find it to be especially pronounced for high-volume injection wells.This correlation is consistent with the results of Langenbruch and Zoback (2016), as well as the findings of transients between 2008 and the end of 2015 by Gupta and Baker (2017) and Montoya-Noguera and Wang (2017).About 25% of all time series suggest two or more change points.The later change points are usually associated with a decreased seismic activity, possibly related to reduced injection volumes (Langenbruch et al., 2018).From around 2016, the majority of change points reflects a decreasing seismicity rate.This indicates a positive effect of the prohibition of large injection volumes in 2016.

Summary and Conclusion
This study analyzed earthquake occurrences in Oklahoma (US) with regard to suspected changes in seismicity as a consequence of the growing number of industrial projects.Earthquakes belong to the most damaging natural hazards worldwide.Their occurrence can be favored by man-made stresses, for example, from mining, geothermal exploration or hydraulic fracking.While most human-induced earthquakes are limited to smallmagnitude events, the US state of Oklahoma reported a rapid increase in seismicity, including events of magnitude five or higher, since the onset of high-volume wastewater injections in 2009.We applied a change point detection algorithm to pinpoint transitions in the magnitude-frequency distribution of earthquakes recorded in comprehensive catalogs from Oklahoma and southernmost Kansas from 1 January 2005-31 December 2019.A Bayes factor that compares the performance of models with no, one or several change points is used to derive the number of transitions.Subsequently, the likelihood function is used to locate the most likely change point(s).A correlation analysis of the identified change points with the onset of injection wells confirmed an increase of seismicity as a consequence of high injection volumes.We further found a cluster of change points with decreasing seismicity around 2016, which indicate a positive effect of the injection reduction mandate from 2016.In contrast to trend analyses or the commonly applied moving window approach, the change point algorithm has no smearing effects and is suitable for detecting rapid modifications of the studied system.The precise determination of changed seismicity simplifies the detection of the seismic response to naturally or anthropogenically changed stress regimes and improves the root-cause investigation.

Need and Potential for Time-Dependent Hazard Models
Table 2 summarizes the five case studies presented, highlighting the specific challenges of the hazards studied and how these were addressed through tailored This process has led to new ways of evaluating data and the development of innovative model designs, providing new insights into changing hazard potentials.In summary, the individual studies support the hypotheses set out in Section 1.3: 1. Identified changes in hazard potential underscore the need for time-dependent hazard models to capture dynamics (e.g., in frequency, magnitude, space-time variability) of hazard events, and, thus, contribute to the dynamics of the associated risks.Significant differences in the manifestations of change are evident, whether hazards evolve gradually or abruptly, show local increases or decreases, or undergo seasonal shifts.2. Accordingly, data evaluation methods must be adapted to the specific nature of the hazard under investigation.
The choice of method used should be based on expected changes and problems encountered in data analysis, rather than on personal preference or popularity of the method.The studies presented are examples of how an efficient study design was developed through close collaboration between data and earth scientists.This facilitated the development of time-dependent hazard models and provided new insights into the hazard type studied.
In the following, we summarize how the presented hazard analyses pave the way for better quantification of transient hazards by addressing the aspects outlined in Section 1.2, namely natural variability on temporal and spatial scales, overlapping and counteracting factors, and limitations in data quality and availability.These approaches serve as showcases and inspiration for how data science techniques can be used to overcome existing challenges in natural hazard assessment.

Natural Variability on a Temporal Scale
Understanding the intrinsic temporal variability in natural hazards is key to identifying transient states.For instance, meteorological dynamics can be characterized by seasonal variations.Annual values, such as annual mean temperature or annual mean precipitation, ignore seasonal effects and may hide, for example, an increase or decrease in seasonal amplitudes or timing.Average values (e.g., annual rainfall values) may also mask transient features in shorter intervals, such as for example, peak rainfall intensity in a sub-hourly time scale.Sections 3 and 4 use oscillating functions or modes to capture seasonal behavior.By identifying and distinguishing seasonality from long-term trends, we receive more distinct information for both components and, in addition, are able to detect shifts in seasonality.Another often assumed form of temporal variability, namely the random occurrence of independent events, is typically modeled as a stochastic Poisson process.This is illustrated in Section 6 for modeling the temporal occurrence of earthquakes.Poisson processes can also be used to estimate return levels of outburst floods, given that associated triggers such as earthquakes, landslides, or ice avalanches can be considered as randomly occurring and independent events (Section 5).
Both the frequency and the intensities of events are prone to change.Hazard assessments are typically interested in the extreme or largest events in a given period, as those are associated with the largest impacts.Extreme events, such as severe wind storms (Section 2), heavy rainfall (Section 3), very high or low river runoff values (Section 4), or large glacier lake outburst floods (Section 5) may be defined as the highest (or lowest) quantiles of reported intensities.This allows us to capture the tails of the probability distribution, which are usually of high relevance with regard to the damage potential.The occurrence of extreme events can be tied to extreme value distributions in order to model the probabilistic behavior of rare events (Section 3 and 5).

Natural Variability on a Spatial Scale
Spatial variability and the peculiarities of different hazard types limit the transferability of hazard models and the joint evaluation of data from different regions.Section 3 emphasizes that return levels of extreme precipitation in Germany vary on local and regional scales.While some regions show a general increase in extreme precipitation intensities, others show a forward or backward shift in seasonality, an increase in seasonal amplitudes, or even decreasing intensities.Further, Section 5 shows a regionally contrasting pattern of glacier lake abundance, which controls the likelihood of observing a lake outburst flood.Both studies suggest models that can account for regional and local differences in magnitude and frequency of the underlying process.Extending Section 3, we refer to Fischer et al. (2019) and Ulrich et al. (2020), who show how the location can be included in modeling extreme precipitation.This approach is especially beneficial for rare events and small data sets, since data from neighboring stations is used to train the model.Section 4 showed that the overall/large scale hydrological runoff dynamics is resulting from a combination of different hydrographical regions and altitudes of the Rhine River basin.For understanding the whole region's dynamics it is essential to understand the mechanisms at the subregional scale.

Superimposing and Counteracting Factors
Superimposing or opposing effects complicate the interpretation of observed results.For instance, the hazard of flooding could increase due to seasonal shifts in precipitation, if extreme precipitation becomes more likely to coincide with snow melting.Rising temperatures can have an opposite effect (Section 4).In the European Alps, rising temperatures reduce snowpacks and consequently decrease the potential for meltwater floods.Moreover, human impacts, such as changing land use or reservoir construction, can hide less pronounced signals, such as changes in snowmelt-induced river runoff.This example illustrates that mechanistic process knowledge is essential to interpret the results and infer possible consequences.A purely data-driven analysis that ignores the interplay of contributing factors may fail to draw relevant conclusions.For example, an unchanged GLOF frequency (Section 5) may point at an unchanged hazard in an unchanged environmental system.Yet, this is in contrast to the observed growth of glacier lakes in number and size.At the same time, GLOF-triggering mechanisms might be affected by changed environmental conditions, for example, the number of lakes affected by ice avalanches may decrease due to glacier retreat.

Limitations in Data Quality and Availability
Depending on the hazard and region studied, aspects of scarce or inconsistent data can be of high relevance.
Especially remote regions, such as high mountain areas in general and the Himalayas in particular, suffer from poor instrumentation and limited data.For instance, considering the occurrence of GLOFs, data recordings are not only incomplete, but also biased towards large and destructive events.Remote sensing can help to systematically detect, and fill inventories of gravitational and hydrological hazard events without drawing on field work or eyewitness reports.In Section 5, the time series analysis of satellite imagery helped to compile a more complete database of GLOFs and to reduce the observation bias.As a result, it increased the reliability of the statistical trend analysis.In contrast, the remaining studies make use of permanent and continuous monitoring, which form the basis for robust statistics on seismic data (Section 6), meteorological data (Sections 2, 3 and 4), and discharge data (Section 4).The amount of data exceeds manual inspection and requires automated data processing with a careful selection of explanatory variables and model designs.Section 2 illustrates the benefit of including domain knowledge in the feature selection to retrieve more information when analyzing complex data.Compared to commonly used windstorm counts, the newly developed storminess metric shows a more distinct signal of the increasing hazard arising from European winter storms.

Data Science for Time-Dependent Hazard Assessments
In summary, to address aspects of natural variability, interfering impacts, and data quality, the analysis of natural hazards requires the integration of knowledge about the data acquisition and processing, a physical and mechanistic understanding of the studied process, and methodological expertise in applicable algorithms and potential pitfalls of data evaluation.As a consequence, natural hazard analyses benefit from a dynamic feedback between different branches of research, including statistics, machine learning, and physical models, which describe the seismological, hydro-meteorological, ecological or geomorphological processes.To exploit the available data  Note.This table summarizes the challenges in the presented hazard studies, and the approaches taken to address and overcome these challenges.It also outlines the innovative aspect of the study design and the main findings with respect to the changing hazard.
most effectively, their analysis needs a careful selection of suitable evaluation techniques with consideration of the research question and the underlying physics.Especially machine learning approaches may not be completely traceable and at some points resemble a black box.Moreover, random coincidences or artifacts in the data can lead to unexpected results.To diminish misinterpretation, the study should be designed to prove or disprove a certain hypothesis, such that "high wastewater injection rates increase the seismic activity" (Section 6).
Consequently, a useful study design demands physical understanding of the underlying process and knowledge about applicable algorithms.
In the presented hazard analyses, each study was designed according to the physical understanding of how frequencies, magnitudes, and return periods of the studied hazard might have changed.We applied linear models to express the belief that hazard frequencies or intensities might have changed more gradually, and less abruptly or in another more complex form.For instance, in Section 2 the linear model confirms an increasing trend in the number and intensity of winter storms in Europe.Surprisingly, in Section 5 the linear model shows no trend in GLOF frequency in the Himalayas, despite a growth of glacier lakes in number and size.This finding is difficult to explain, and its interpretation requires the consideration of trigger mechanisms and the underlying physical processes.The models in Sections 3 and 4 are not limited to linear trends, but also account for seasonal behavior and changes in seasonality.They identify local trends and seasonal shifts for extreme precipitation in Germany (Section 3) and river runoff in the European Alps (Section 4).In contrast to linear models, a change point model is preferred if the studied hazard is expected to be subject to more sudden changes.Unlike a trend analysis or a moving window approach, the change point model can capture sharp responses to abrupt system changes.In Section 6, a change point model helps to detect the onset of increased or decreased seismicity as a consequence of changed wastewater injection.In each studied research question, the custom-fit model selection is key to provide new findings from the available data.

Conclusion
We have presented five case studies to illustrate the multifaceted nature of changing natural hazards in recent decades.It is important to realize that such changes can become apparent in multiple data characteristics, for example, in frequency, magnitude, temporal or spatial pattern (either as an organized or random type of variability), or in a gradual/linear or sudden/abrupt manner.Our selection includes different hazard types and different manifestations of change: The increasing trend in European winter storms could be highlighted by showing not only an increase in the number of storms, but also in their severity (Section 2).Different forms of change could be observed for extreme precipitation in Germany.Increasing and decreasing trends in daily precipitation extremes could be identified, as well as shifts and changing amplitudes of the seasonal cycle.In addition, the prediction of return levels could be improved (Section 3).Decreasing snowpack, changes in precipitation and river engineering are leading to a change in the seasonality of river discharge in Alpine catchments, with a reduced risk of meltwater-induced floods but an increased frequency of rainfall-induced floods (Section 4).Despite the increasing number and area of glacial lakes in the Himalayas, no change in annual GLOF frequency has been observed over the last 30 years.However, spatial variations in GLOF hazard are significant and highly correlated with the abundance of meltwater areas and the change in lake areas over the past decades (Section 5).
The identification of change points in seismicity data from Oklahoma (USA) and their correlation with the onset of wastewater injection wells improves the attribution of increased seismicity to the onset of high-volume injection wells (Section 6).The synthesis of the case studies not only demonstrates the need for time-dependent hazard models to capture changing hazard potentials, but also illustrates the different ways in which change can manifest itself, such as gradual versus abrupt changes, altered amplitudes, or seasonal shifts.
Scientific advances in data-driven methods and the increasing availability of climate-related and environmental data have facilitated in-depth analysis of large data sets beyond human interpretation.This is paving the way for new opportunities to identify and quantify changes in natural hazards.However, each hazard assessment has different needs and shortcomings, requiring tailored methodologies.The case studies presented provide strategies to address specific shortcomings, such as variability on temporal and spatial scales, interacting processes, or limited data quality.They serve as examples for similar scenarios and suggest assessment techniques tailored to the specific nature and expected change of the hazard under consideration: Linear models serve to identify gradual trends, oscillatory functions adequately describe seasonal variations, extreme value distributions are effective in modeling changes in the most hazardous events, and change point detection techniques are used to identify sudden shifts and improve attribution to driving forces of change.
As a final conclusion, we state that 1.The presented combination of studies demonstrates the need for an adequate and customized analysis and model selection.The choice of an appropriate study design is critical for reliable and efficient hazard assessment under change.When applied correctly and supported by a mechanistic/physical understanding of the governing processes, this will enable to disentangle the dynamics about the hazard type under investigation.2. Time-dependent hazard analyzing methods and models enable to capture changing hazard potentials.They further illustrate the different ways in which change can manifest itself, such as gradual versus abrupt changes, altered amplitudes, or seasonal shifts.3. Based on our experiences, we think that an over-simplification and/or a broad averaging over large scales in time and space will result in severe loss of information and hide the inherent mechanisms/reasons for changes.
Consequently, this may yield to misleading system's analysis or even prognosis about the future development of such changes.4. The symbiotic alliance between data science and earth science is a key determinant of successful timedependent natural hazard assessments.

Figure 1 .
Figure1.Storminess (black)  and windstorm counts (blue) in the period 1950-2010, extracted from ERA-20C data.The dots correspond to the value of the extended ONDJFM winter in year/year + 1. Black labels are notable hazardous windstorms in the respective ONDJFM winter.Both trends (dashed lines) are significant to p < 0.05 using a t-test.

Figure 2 .
Figure 2. Monthly maxima of daily precipitation sums shown as Box-Whisker-Plots at three different stations in Germany and the 100-year return level of the first (dotted) and the last year (solid) of the observation period as well as for 1941 (dashed) (top row).The bottom row shows the temporal course of the monthly 100-year return levels.

Figure 3 .
Figure 3. Discharge of the High Rhine measured at the gauge in Basel, Switzerland, for the period 1869-2016: (a) seasonality of river runoff investigated by estimating discharge quantiles [m 3 s 1 ] on a daily basis, (b) onset and evolution of changes assessed by applying EMD on a daily basis and displaying the residuals [m 3 s 1 ] of the decomposition, and (c) changes in quantiles (determined on an annual level) using EMD and displaying the centered residuals [m 3 s 1 ] of the composition (adapted after Rottler, Francke, et al. (2020)).

Figure 4 .
Figure 4. Distribution and return periods from moraine-dammed glacier lake outburst floods (GLOFs) in the Himalayas.(a) Bubbles show the size and abundance of moraine-dammed lakes in 1°× 1°bins.Triangles are previously known GLOFs, and red dots are GLOFs detected from Landsat images between 1988 and 2017.(b) Estimated return periods of GLOFs in allregions, and in the Eastern (Western) Himalayas, where regional GLOF hazard is highest (lowest).Blue numbers are the mean 100-year return level of peak discharge Q p including the 95% highest density interval; brown numbers are the mean GLOF rate in this region.Brown ticks are historical GLOFs.The GLOF from lake Zhangzangbo in 1981 had a peak discharge of ∼16,000 m 3 s 1(Xu, 1988) and a 100-year return period, accordingly.Figure adapted from Veh et al. (2020).
Figure 4. Distribution and return periods from moraine-dammed glacier lake outburst floods (GLOFs) in the Himalayas.(a) Bubbles show the size and abundance of moraine-dammed lakes in 1°× 1°bins.Triangles are previously known GLOFs, and red dots are GLOFs detected from Landsat images between 1988 and 2017.(b) Estimated return periods of GLOFs in allregions, and in the Eastern (Western) Himalayas, where regional GLOF hazard is highest (lowest).Blue numbers are the mean 100-year return level of peak discharge Q p including the 95% highest density interval; brown numbers are the mean GLOF rate in this region.Brown ticks are historical GLOFs.The GLOF from lake Zhangzangbo in 1981 had a peak discharge of ∼16,000 m 3 s 1(Xu, 1988) and a 100-year return period, accordingly.Figure adapted from Veh et al. (2020).

Figure 5 .
Figure 5. Locations and occurrence times of transients for models with one, two or three change points, compared with approval dates of injection wells between 2005 and 2019 in Oklahoma.Wells that had an approved injection volume >10,000 barrels per day are black.(a) Map of the first occurrence of estimated transitions, (b) latitude-time plot, and (c) time-longitude plot with all estimated change points and injection wells.For all transients, it is illustrated whether an increasing or a decreasing rate is given.

Table 1
Description of the Case Studies Considered in the Individual Sections of This Paper

Table 2
Summary of the Case Studies Presented in Sections 2 to 6