Compound Flooding in Halmstad: Common Causes, Interannual Variability and the Effects of Climate Change

2.1 Models

The model data used for this study are derived from dynamical downscalings of global coupled climate models from the coupled model intercomparison project phase 5 (CMIP5) (Taylor, Stouffer and Meehl, 2012). The atmospheric variables from the global models are further downscaled using the regional atmosphere model named RCA4 (Berg et al., 2013; Samuelsson et al., 2011). The model has a grid covering Europe and large parts of the North Atlantic, with a horizontal resolution of 25 km and 40 vertical levels. The global model is used as a boundary condition to drive the regional model.

The ocean variables are similarly dynamically downscaled using the ocean-sea ice model NEMO-Nordic (Hordoir et al., 2019; Madec and the NEMO Team, 2016). The regional ocean model covers the North Sea and the Baltic Sea with a horizontal resolution of 3.7 km and has 56 levels in the vertical. The regional ocean model is forced with downscaled atmospheric fields from the RCA4 model as its upper boundary condition. On the open boundary in the North Sea it uses a combination of monthly temperature and salinity data from the global models and hourly sea level and currents data that has been derived from the RCA4 downscalings using a machine learning approach (Hieronymus, Hieronymus and Hieronymus, 2019). NEMO-Nordic is also forced with river run-off from the E-HYPE model

The process-based hydrological model HYPE (HYdrological Predictions for the Environment; Lindström et al., 2010) simulates components of the catchment water cycle and water quality at a daily time step. The model is semi-distributed, in which a river basin may be subdivided into multiple subcatchments, which are further subdivided into hydrological response units (HRUs) based on soil type and land use classes, and is set up for a pan-European domain based on irregular polygonal catchment delineations as presented in Donnelly et al. (2016) and Hundecha et al. (2016). The E-HYPE simulations that are used as inputs to NEMO-Nordic are taken from the C3S data set (Berg et al., 2021a), where specifically the time series of river discharge from the E-HYPEcatchM00 model are used. The model simulations are driven with meteorological parameters from the accompanying data set in C3S (Berg et al., 2021b).

The complete model system, with its ocean, atmosphere and river components, is composed of standalone models forced with the same downscaled atmospheric forcing. That is, the hydrological model and the ocean model are forced using the same CMIP5 projections downscaled with the same RCA4 model. It is thus not a fully coupled system, but impacts in form of sea level, precipitation and river streamflow extremes still co-vary in a consistent manner, as all these impacts result from the same modelled atmospheric states.

Table 1 shows the climate scenarios and global climate models used in the downscalings. All scenario runs are concentration based (i.e., they are forced with atmospheric concentrations of greenhouse gases and aerosols). Our historical simulations start in 1971 and are forced with greenhouse gases and aerosols that are consistent with the observed evolution of these quantities until the year 2005. The historical simulations that were ran with global models in CMIP5 actually start in 1850, but only the period from 1971 has been downscaled, which is why our historical simulations start in 1971. The future scenarios start in year 2006, and range from low concentrations under RCP2.6 to very high concentrations under RCP8.5, while RCP4.5 is a middle of the road scenario. In terms of global mean surface temperature (GMST) increase at the end of the century compared the preindustrial; RCP2.6 gives 1.6°C, RCP4.5 gives 2.4°C and RCP8.5 gives 4.3°C (IPCC, 2013). These values refer to ensemble medians for the whole CMIP5 ensemble. GMST is, however, an emergent property of climate models so values differ between models.

The RCP scenarios have not been assigned probabilities of occurrence by their makers, but some work using integrated assessment models (IAMs) has been aiming to assign such probabilities. Using information about fossil fuel resource availability, economic modelling and assuming current policies; Capellán-Pérez et al. (2016) and Huard et al. (2022) both found RCP4.5 to be the most likely scenario for the current century. The probability of having emissions as high as RCP8.5 decreased with time and reached zero in several of the IAMs probed by Huard et al. (2022) toward the end of the current century. Similarly, both Capellán-Pérez et al. (2016) and Huard et al. (2022) found the probability of having emission as low as those in RCP2.6 to be very low under current policies toward the end of the century. To sum up, future emission trajectories remain highly uncertain, but at least under current policies RCP4.5 appears to be the most likely pathway out of the three probed here.

2.2 Observations and validation

The model validation is done using daily means. This is partly because it is the available temporal resolution of the river flow data. Modelled and observed sea levels, meanwhile, are available with hourly frequency. However, since our interest is in multiple extremes, it is the possibility of these hazards to interact with each other that is our primary interest. Since the river Nissan is regulated, it is possible to, to some extent, prevent short-duration sea level extremes from coinciding with peak river discharge. Moreover, given that our focus is on long time scales, trends and interannual variability, daily means appears as a good resolution choice.

Figure 2 shows a quantile-quantile plot of daily mean sea level from our three modelled historical simulations plotted against sea level observations collected at Halmstad harbour between 2009 and 2022. The model shows a considerable underestimation of the sea level variability compared to the observations. Multiple reasons may help explain the bias. The atmospheric model very likely underestimates the surface winds, which is the primary driver of sea level variability in the area (Hieronymus, Hieronymus and Arneborg, 2017). Moreover, it is not just the wind at the time of a particular storm that is important, but also the wind the weeks and months before as it sets a baseline from which the storm surge acts. Such preconditioning has a considerable influence on sea level extremes inside the Danish Straits just to the south of our location in the semi-enclosed Baltic Sea (Andrée et al., 2023; Bellinghausen et al., 2023; Hieronymus and Hieronymus, 2021), and could potentially be important also in the Kattegat. However, it has been demonstrated that the modelling system performs much better in the Baltic Sea than on the Swedish west coast (Dieterich et al., 2019; Hieronymus, Hieronymus and Arneborg, 2017), even though the wind is underestimated over the Baltic Sea as well, suggesting that lack of preconditioning is unlikely to be a primary cause of our bias. The sea level dynamics are much more constrained by basin geometry in the semi-enclosed Baltic Sea than on the more open Swedish west coast, making the latter case harder to model. Complex basin geometry, poorly resolved by our model, is thus also likely contributing to the underestimation of sea level variability. A two-way nesting scheme giving higher resolution locally could be a way forward (Andrée et al., 2023). However, for Halmstad, there is also another potential reason for the underestimation of the modelled sea level variability. That is, the presence of a local effect that occasionally raises the sea level in Halmstad way beyond what is seen at neighbouring stations that otherwise have very similar sea level variations.

Figure 2

Figure 3 shows some sea level and wave parameters in Halmstad and at neighbouring tide gauges in Ringhals to the north and Viken to the south, during the hours when the sea level in Halmstad exceeds 1 m above the mean. The 1 m threshold is of course rather arbitrary, but the choice of threshold is not important. It merely serves the purpose of demonstrating that on some, but not all, occasions when sea levels in the area are very high they are much higher in Halmstad than in neighbouring locations. The hourly sea level measurements in Halmstad for the period 2009–2022 contain in excess of 150 observations that exceed 1 m above the mean. For a large-scale storm surge one would expect similar sea levels to be reached at these neighbouring stations. The panels showing sea level difference show that this is generally the case. However, we see at least three peaks at around t = 25, 87 and 113 when sea levels in Halmstad are significantly higher than at the neighbouring gauges. It is also evident that these peaks coincide with peaks in significant wave height. Moreover, all these incidents occur when waves are coming from the southwest. In fact, the wave direction during the times when we see a large local effect are slightly more from the south than the typical wave direction during the times when sea levels are 1 m above the mean or more in Halmstad. This might be important, as it is clear from Figure 1b that waves can only enter the small bay inside the Laholm Bay where Nissan has its river mouth from the south and southwest. This small bay is also where Halmstad harbour and our tide gauge is located. Wave set-up is a change in the mean sea level at the shoreline that is due to the presence of breaking surface waves, which in contrast to short-lived swash bores (uprush of waves collapsing on a beach) can last for significant time (Longuet-Higgins and Stewart, 1963), and may thus influence the sea level measured by a tide gauge. Moreover, an investigation from the neighbouring Baltic Sea suggests that wave set-up could give rise to considerable local water level variations in some locations (Pindsoo and Soomere, 2015).

Figure 3

Other potential explanations could be wind set-up, which, in contrast to wave set-up, is part of the model’s dynamics but may be under-resolved. Another possibility is some interaction with the out-flowing river plume. We note, however, that the river discharge was not particularly high during any of these events, so such an interaction appears unlikely to be a sole or sufficient cause for the local effect. While these events present intriguing aspects, they do not constitute the main focus of our current investigation. This is primarily due to the fact that a comprehensive investigation of this phenomenon would necessitate observational data that is currently unavailable. Most prominently, one would need more tide gauges inside the Laholm Bay so the scale of the phenomena could be determined accurately.

An implication of the phenomenon’s small spatial scale and possible relation to surface wave dynamics, however, is that such events are unlikely to be well represented by the ocean model. Surface waves are not part of the model’s dynamics and a horizontal resolution of 3.7 km may well be too coarse to capture the effect. Given that wave set-up can last for significant time and have a magnitude reaching 1/4 of the wave height at the breaking line (Longuet-Higgins and Stewart, 1963; Pindsoo and Soomere, 2015), the phenomenon could explain the local effect and perhaps, at least, partly the bias seen in Figure 2. However, since that figure is based on daily means and the phenomenon is intermittent, it is likely not the sole or even primary cause of the model bias, which we believe underestimated surface winds and challenging basin geometry to be. Nevertheless, the model bias, although unfortunate, is unlikely to affect the conclusions of this study. Here our focus is on the co-variability of different hazards, which should be well captured owing to all models using the same atmospheric forcing. Moreover, the ocean model shows a high correlation to sea level observations also on the west coast when high-quality forcing is used, so the timing of extremes is generally well represented (Hordoir et al., 2019). Long-term trends in sea level are not based on the regional ocean model, but on global climate, ice sheet and glacier models (Fox-Kemper et al., 2021), so these trends are unaffected by this bias.

The E-HYPE model is calibrated on a set of calibration stations in different hydrological response classes around Europe. This approach allows parameter settings to be applied to ungauged catchments with generally good results, as evaluated with independent validation stations for each class. The model is then evaluated with separate validation stations. The performance for the validation stations around the Baltic Sea are presented in Figure 4 for the Kling-Gupta Efficiency score (KGE; Gupta et al., 2009), the Nash-Sutcliffe Efficency (NSE; Nash and Sutcliffe, 1970), and the percentage bias (Pbias) for river discharge. NSE is a measure of predictive skill of hydrological models, with values ranging from negative infinity to one, which is a perfect score. NSE = 0 is equivalent to using the mean flow as a predictor. KGE is a goodness-of-fit measure, where a value of KGE = –0.41 has a similar interpretation as NSE = 0 (Knoben et al., 2019). The model is here evaluated for the catchments with validation stations that contribute to the region of the ocean model. There was no gauge with sufficient data records for the modelled period, nor sufficiently close to the Nissan river mouth for performing a local validation. The E-HYPE model performs overall well on all three scores in this evaluation, and in an earlier evaluation for all of Europe presented in Donnelly et al. (2016) and Hundecha et al. (2020).

Figure 4

2.3 Methods

As discussed in the introduction, one may view compound flooding as the joint action of two or more random variables. Here we let X(t) represent time varying sea level and Y(t) time varying river discharge. Important metrics such as flooded area or property damage are then given by some functions f(X,Y). The most simple model, that might apply reasonably well to for example flooded area, would be to assume $f\text{=}\tilde{X}\text{+}\tilde{Y}$, where $\tilde{X}$ is the area flooded by the ocean and $\tilde{Y}$ the area flooded by the river. In this simplified view it is easy to link the variability in the two random variables to the variability in their sum through

where σ is the standard deviation and r the Pearson correlation coefficient. The expected value, E, of the flooded area is, of course, simply

where the ≤ sign is there instead of an equal sign because the same area could be affected by both riverine and coastal flooding and should not be counted twice. The above decomposition clearly highlights that the joint flooding risk depends on the variance and magnitude of the individual variables and the correlation between them. These are therefore the main parameters investigated in the study. It should, however, be noted that our decomposition is not a full model of the dependence structure between two random variables. Such dependencies are often modelled using copulas that enables an isolation of the dependence structure in multivariate distributions (Dupuis, 2007). Some applications of copulas modelling to compound flooding is found in Bevacqua et al. (2017) and Phillips et al. (2022).

From our different downscaled climate scenarios, see Table 1. We have extracted time series of daily mean sea level, daily mean precipitation and daily mean river discharge in Halmstad. These time series are linearly detrended to remove long-term trends. This is done to study the co-variance of extremes, while trends are treated separately. Extremes are studied using time series of daily maxima for each year and calendar month extracted from our different model projections. The reason for using annual maxima time series for different months is twofold. Firstly, time series of block maxima can only converge in distribution to the Generalized Extreme Value (GEV) distribution by the Fisher-Tippett-Gnedenko theorem, so the limiting distribution is known. This extreme value version of the central limit theorem is often used with annual blocks, to infer return level curves for different hydrometeorological and oceanographic parameters. Annual maxima time series are thus a natural starting point when investigating extremes. Secondly, the division into monthly maxima gives us a way to infer the seasonality of the different extremes. Moreover, the window is short enough that multiple events occurring in the same month can be expected to lead to worse consequences than if they had been separated longer in time. That is, they can be considered as compounding events, at least in an approximate sense. From the model we also extract a normalized station based (Lisbon-Rejkavik) monthly mean North Atlantic Oscillation (NAO) index. Variability in annual sea level maxima has been connected to the NAO in observational studies on the Swedish west coast (Hieronymus and Kalén, 2020) and warrants investigation also in the model. A positive NAO typically means mild and wet winters for Sweden, whereas a negative NAO brings cold and dry air masses from the Arctic. The river flood response to these changes depends on the river, and largely on the interaction between rainfall and snow melt.

3.1 Hazard co-variance and its trend

A fundamental prerequisite for having compounding hazards is that the hazards involved have similar seasonality. In Figure 5 we have computed the 100-year return levels for river streamflow, precipitation and sea surface height from a long annual maximum time series created by concatenating data from all our different climate projections, leading to the creation of a time series more than 1000 years long. The return levels in Figure 5 are normalized through division by the return level from the month with the highest values for the respective hazard. This gives return level values between zero and one, where one is the highest monthly return level. It is evident from the figure that river discharge and sea surface height extremes are winter intensified, while extreme precipitation occurs predominantly in summer. A natural assumption is thus that compound flooding in Halmstad is most likely to occur in winter as a consequence of co-occurring sea level and river streamflow extremes. Compound pluvial and coastal flooding will thus be our primary focus in this investigation.

Figure 5

Equation 1 showed that the correlation coefficient between, and the standard deviations of, the different hazards are important parameters for determining compound flooding risk. In our case with three different parameters we have three different correlations to compute. River discharge to sea surface height, river discharge to precipitation and precipitation to sea surface height. These correlation coefficients for the different months are shown in Figure 6. All three panels clearly show that the annual maxima of these hazards are quite strongly correlated in winter months and uncorrelated in summer months, although with a considerable spread. Correlation coefficients exceeding |0.2| are significant at the 0.05 level in the longer RCP scenarios, while correlations coefficients exceeding |0.33| are significant in the shorter historical simulations. In this paper we focus primarily on Pearson correlation coefficients as those are the ones the that feature in our Equation 1. However, many studies of compound weather events instead use rank correlations that evaluate the monotonic rather than the linear relationship. In Figure 7 we show the Pearson, Kendall and Spearman correlation coefficients between annual maxima of sea surface height and river discharge. The results are qualitatively similar for all correlation coefficients, indicating that the relationship between the two variables is approximately linear.

Figure 6

Figure 7

Low-pressure systems bringing both high precipitation and storm surges to the Swedish west coast appear as an obvious cause for the co-variability of our three hazards. A natural assumption is therefore that this co-variant interannual variability could be linked to the NAO, a connection also made for several other hazards in the Baltic Sea area by Rutgersson et al. (2022). Of particular interest for our study is that surface wave energy was recently found to correlate with the NAO in the area, making the case for compound flooding even stronger (Adell et al., 2023), as surface waves can also affect sea levels.

Correlation plots between the NAO index and our three hazards are shown in Figure 8. It is clear from the different panels that our three hazards are all correlated to the NAO in winter months. Note that the NAO index here derived is based on monthly means rather than monthly maxima (as in the hazard time series). This choice is made because we expect a connection between the atmospheric mean state and the prevalence of low-pressure systems, rather than a direct connection between extreme NAO values and extreme hazards at the site. The correlation persists, however, but is somewhat weaker if a block maxima of the NAO is used (not shown). Also noteworthy is that there seems to be no consistent difference between the scenarios or GCMs in terms of the strength of the correlations.

Figure 8

To investigate further, we calculated the running correlation coefficient between the different hazards over 30 years periods for all months. That is, correlation coefficients are calculated for sliding 30 year periods, like a moving average but for the correlation coefficient. Figure 9 shows the result for the correlation between river discharge and sea surface height. Correlation coefficients exceeding |0.36| are significant at the 0.05 level for these 30-year period data slices. The correlation coefficients between the other hazards show similar behaviour (not shown). In all cases it is evident that the natural variability dwarfs potential warming driven trends in the correlation coefficients. A very large ensemble would be needed to get a good enough signal to noise ratio, so that a trend could be reliably estimated in the presence of this considerable natural variability. Consequently, a useful first approximation is to take these correlations to be independent of climate scenarios.

Figure 9

The correlations computed from the annual maxima time series indicate coherent interannual variability between our different hazards. Next we consider the strength of the correlation on different time scales. Figure 10 shows the correlation coefficient between the low pass filtered daily mean river streamflow and sea surface height as a function of the length of the filter window. It is clear from the figure that the highest correlations are found on time scales between about 100 days and a year. These results, once again suggest an important role of interannual variability for compound flooding in the area.

Figure 10

Having found strong, and to a first approximation emission scenario independent, correlations between our different hazards and the NAO we proceed to investigate how the standard deviation of the different hazards and the NAO index might change in a warming world. Running standard deviations over 30-year periods are shown in Figures 11, 12, 13, 14. Given the very considerable natural variability, it is clearly not possible to discern a climate driven trend in any of these standard deviations with the amount of data we have at our disposal. However, absence of evidence does not mean evidence of absence. Thus, just as for the correlation coefficients we may conclude only that if there is a warming-driven trend, it is currently hidden within a sizeable natural variability.

Figure 11

Figure 12

Figure 13

Figure 14

As a first approximation, it thus seems fair to infer that trends in compound flooding will primarily be driven by the trends in the means of river discharge and sea level, rather than by trends in the extremes. Hieronymus and Hieronymus (2023) reached a similar conclusion for sea level extremes in the area.

3.2 Trends in individual hazards

Figure 15 shows mean sea level projections for Halmstad for a number of different shared socio-economic pathway (SSP) radiative forcing combinations based on the IPCC’s AR6 mean sea level projections but corrected with more accurate land uplift data (Fox-Kemper et al., 2021 Vestøl et al., 2019). The radiative forcing, the number after the dash given in Wm^–2, is the net change in the radiative balance at the top of the atmosphere in the year 2100 compared to a pre-industrial background. A higher number implies greater warming. The mean sea level, in contrast, to the standard deviation of the daily mean sea level shows a very distinct trend under all scenarios. Coastal flood risk simulations from a number of locations around the Swedish coast have shown that flood risk is driven by the risk of high extremes in the next few decades, and by mean sea level rise in a longer perspective (Hieronymus, 2021; Hieronymus and Kalén, 2022; Hieronymus, 2023). The same is true in Halmstad, where it is clear that the land uplift of about 2 mm per year is not enough to offset significant relative sea level rise even under the very low emission scenario SSP1-1.9.

Figure 15

On the other end of the spectrum we have the very high emission scenario SSP5-8.5 and its alternative projection SSP5-8.5 low confidence. In the latter, the sea level contributions from Antarctica and Greenland that are usually taken from large model intercomparison programs are instead taken from some of the highest estimates in the published scientific literature (Bamber et al., 2019; DeConto et al., 2021). In these two unlikely scenarios (Hausfather and Peters, 2020), the projected mean sea level rise toward the end of the century, is so large that it may exceed even the largest daily sea levels in the observations. In all emission scenarios, it is extremely likely that sea levels far outside of the historical range will be seen in the current century. Moreover, while the SSP5-8.5 projections represent improbable high emissions, having mean sea level rise in the meters is expected also for much lower warming scenarios albeit significantly later (Fox-Kemper et al., 2021). Therefore, when considering coastal flooding, whether as a single or multiple hazard, the most concerning trend is evidently mean sea level rise.

River discharge is affected by the competing processes of a general increase in precipitation and the increasingly intense evapotranspiration driven by temperature. The mean river discharge from Nissan shows clear evidence of increases in winter (January and February), and decreases in late spring (April and May) toward the end of the century. The median ensemble changes are rather weak for RCP2.6 (10 to 20% in January/February, and –20% in April/May), but rather strong in RCP8.5 (55 to 60% in January/February, and –35 to 45% in April/May). The annual maximum discharge shows no significant increase for RCP2.6, but some increase of 7% in RCP4.5 and a strong increase of 22% for RCP8.5, at the end of the century.