Dynamical downscaling and data assimilation for a cold-air outbreak in the European Alps during the Year Without a Summer of 1816

The “Year Without a Summer” in 1816 was characterized by extraordinarily cold and wet periods in central Europe, and it was associated with severe crop failures, famine, and socio-economic disruptions. From a modern perspective, and beyond its tragic consequences, the summer of 1816 represents a rare opportunity to analyze the adverse weather (and its impacts) after a major volcanic eruption. However, given the distant past, obtaining the high-resolution data needed for such studies is a challenge. In our approach, we use dynamical downscaling, in combination with 3D variational data assimilation of early instrumental observations, for assessing a cold-air outbreak in early June 1816. We find that the cold spell is well represented in the coarse-resolution 20th Century Reanalysis product which is used for initializing the regional Weather Research and Forecasting Model. Our downscaling simulations (including a 19th century land use scheme) reproduce and explain meteorological processes well at regional to local scales, such as a foehn wind situation over the Alps with much lower temperatures on its northern side. Simulated weather variables, such as cloud cover or rainy days, are simulated in good agreement with (eye) observations and (independent) measurements, with small differences between the simulations with and without data assimilation. However, validations with partly independent station data show that simulations with assimilated pressure and temperature measurements are closer to the observations, e.g., regarding temperatures during the coldest night, for which snowfall as low as the Swiss Plateau was reported, followed by a rapid pressure increase thereafter. General improvements from data assimilation are also evident in simple quantitative analyses of temperature and pressure. In turn, data assimilation requires careful selection, preprocessing, and bias-adjustment of the underlying observations. Our findings underline the great value of digitizing efforts of early instrumental data and provide novel opportunities to learn from extreme weather and climate events as far back as 200 years or more.

1 Introduction

In central and western Europe, the extraordinary year of 1816 was referred to by historians as a “Year Without a Summer”, with particularly cold, wet, and cloudy conditions during the summer months; it is also known as “Eighteen Hundred and Froze to Death” for similar weather and climate in the northeastern USA (Auchmann et al., 2012; Briffa et al., 1998; Brönnimann and Krämer, 2016; Crowley et al., 2014; Stommel and Stommel, 1983; Wetter et al., 2011). The shifted precipitation patterns and summer cooling can partly be explained by the enormous and devastating eruption of Mount Tambora in Indonesia in April 1815 (Fischer et al., 2007; Harington, 1992; Oppenheimer, 2003; Raible et al., 2016; Robock, 2000, 2007; Schurer et al., 2019; Stothers, 1984; Wagner and Zorita, 2005) and, to a lesser extent, by random internal variability and low solar variability during the Dalton minimum (Anet et al., 2014). Besides diseases and a socio-economic depression after times of war, the adverse climatic and meteorological conditions led to delayed plant growth, crop failures, poor fruit harvests, rising food prices, and famine (Brázdil et al., 2016; Krämer, 2015; Luterbacher and Pfister, 2015; Trigo et al., 2009). Central Europe and Switzerland were among the most affected regions.

From a modern perspective and beyond its tragic consequences, the Year Without a Summer represents one of the few opportunities, or even the only opportunity, to not only understand the climatic and meteorological situations and developments after a major volcanic eruption but also to model the potential consequences of such adverse cold and rainy weather. However, 1816 was more than 200 years ago, and it is therefore not easy to obtain the necessary information.

The Year Without a Summer of 1816 in (central) Europe has been investigated by historians and historical climatologists using information from descriptive sources, weather observations, and early instrumental measurements (Auchmann et al., 2012; Brázdil et al., 2016; Brugnara et al., 2015; Harington, 1992; Pfister, 1999; Trigo et al., 2009). Most research uses monthly aggregated information to describe the summer of 1816, even if daily resolved observations or measurements are available. For instance, Brugnara et al. (2015; based on Dobrovolný et al., 2010) found that 1816 had the coldest summer in central Europe in instrumental records and the second-coldest since 1500 when including documentary evidence. Some studies have analyzed the weather in 1816 also at a daily resolution, mainly on the basis of instrumental observations. An analysis of twice-daily observations of temperature and cloudiness from Geneva, compared to derived weather types, revealed an increased number of cloudy days leading to a larger cooling effect in the afternoon than at sunrise and an association with more frequent precipitation during the summer of 1816 (Auchmann et al., 2012) than in a contemporary reference period. While many of these traditional reconstructions of historical extreme events rely on historical weather observations of a good quality and within a relatively dense spatiotemporal network (Maugeri et al., 1998; Wagner et al., 2016), modern applications in industry and science increasingly require numerical data from state-of-the-art products that typically provide gridded data.

One type of gridded data is atmospheric reanalyses which spans the full globe and up to multiple centuries, such as the current 20th Century Reanalysis version 3 (20CR; Slivinski et al., 2019). For the summer of 1816, the first dynamical atmospheric analysis with a scout version of 20CR was very successful in mapping pressure fields and the associated atmospheric dynamics over much of the eastern North Atlantic and Europe, despite a limited network of station barometer observations (Brohan et al., 2016). Most reanalyses provide subdaily gridded data of a three-dimensional atmosphere, but their spatial resolution is often too coarse to reproduce local processes in the atmosphere and at the surface. Particularly for modeling the impacts of weather conditions at the surface, data with high spatial resolution are required. For instance, weather reconstructions with horizontal grid sizes at a low kilometer scale were applied for modeling impacts of (extreme) events on agriculture (Blanco-Ward et al., 2019; Flückiger et al., 2017; Glotter et al., 2014), storms and related economic loss (Pinto et al., 2010; Stucki et al., 2015), or flood events and their impacts on infrastructure (Boé et al., 2007; Mahoney et al., 2022; Rössler and Brönnimann, 2018; Stucki et al., 2018), among others.

For Switzerland, gridded daily fields (e.g., of temperature and precipitation) with spatial resolutions as high as 1 km have been created with an analog resampling method (Flückiger et al., 2017; Imfeld et al., 2023; Pfister et al., 2020) and by other statistical approaches including principal component analyses (PCAs) of a modern dataset and interpolation of PCA scores from historical station data (Isotta et al., 2019; Stucki et al., 2020). Many of these statistical(–dynamical) reconstructions (or reanalyses) which numerically estimate the state of the atmosphere at a certain point in time involve data assimilation procedures (Carrassi et al., 2018). On a regional scale that has a finer temporal and spatial resolution, datasets have been derived from long-term reanalyses for several areas by statistical downscaling from global reanalysis products (Caillouet et al., 2016, 2019) and have been further refined by the assimilation of local and independent data such as historical surface observations of, e.g., temperature or pressure (Devers et al., 2020, 2021).

In this study, we aim to apply so-called dynamical downscaling, combined with the assimilation of pressure and temperature observations, to refine the information on the weather in summer 1816 from the coarsely gridded global 20CR reanalysis product. Dynamical downscaling procedures include nesting a limited-area domain from a weather forecast model into the global reanalysis product. This process can then be iterated to refine the global fields of atmospheric variables to local scales (Gómez-Navarro et al., 2018; Michaelis and Lackmann, 2013; Stucki et al., 2015, 2018, 2020). We use the Weather Research and Forecast Model (WRF; Skamarock et al., 2019) for downscaling 20CR to regional to local and (sub-)hourly scales, as well as for data assimilation. In contrast to statistical downscaling procedures, physical processes are inherent to the weather forecast model, leading to physically consistent simulations of the weather. Physical consistency of the high-resolution simulations is crucial, especially when they are used in a model chain (Maraun et al., 2010; Muerth et al., 2013). Furthermore, statistical downscaling is usually limited to a handful of atmospheric and surface variables, whereas dynamical downscaling simulations provide an encompassing set of variables and preserve their physical coherence (Fowler et al., 2007; Muerth et al., 2013). These properties make dynamical downscaling attractive despite the high computational costs associated with running a regional circulation model. For nowcasting and forecasting applications in modern periods, weather simulations from downscaling data of global circulation models or reanalyses have been further improved by assimilating additional local data from independent conventional observation sites, e.g., precipitation observations and Doppler radar or satellite data (Ban et al., 2017; Fatmasari et al., 2019; Gopalakrishnan and Chandrasekar, 2018; Thiruvengadam et al., 2020; Wang et al., 2013). Although such modern remote sensing techniques were not yet available in 1816, there is a considerable number of station observations that can be assimilated, particularly for central Europe. This opportunity motivates us to explore whether a combination of dynamical downscaling and assimilation of local temperature and air pressure observations can be at all feasible and successful for this region and for such a distant period like the Year Without a Summer, thus enabling us to analyze and learn from such historical extreme weather and climate events in detail.

Our simulations focus on an earlier cold spell over central Europe, which occurred approximately between 5 and 11 June 1816 and for which we have collected information from station measurements of temperature and pressure but also from weather diaries and records of observations with the naked eye regarding sunshine and cloudiness and the occurrence of precipitation, wind, and other variables. Being one of the most pronounced cold spells of this summer and with abundant data available, it serves as an excellent example for our analyses. We briefly discuss the representation of our case study period in the 20CR dataset, not least with respect to the feasibility of using 20CR data as atmospheric boundary conditions for the WRF simulations. Downscaling results are evaluated based on qualitative weather descriptions, as well as quantitative records of surface pressure and temperature. To provide an independent assessment, data from four stations are not assimilated but retained for validation purposes.

The article is organized as follows: the available station observations, the 20CR reanalysis product, the weather forecast model, and the data assimilation system are described in Sect. 2. In Sect. 3, we start with a brief description of the weather conditions in the summer of 1816 and, more concretely, during the cold spell in June 1816 from documentary evidence and observations, as well as its representation in the 20CR reanalysis product. In the second part, results from dynamical downscaling with the WRF model are shown. This includes (i) simulations without and (ii) with the data assimilation of early instrumental station observations. We illustrate the representation of the cold spell in early June 1816 in these simulations and compare them to independent station observations. A summary and the conclusions are given in Sect. 4.

2 Data and models

2.1 Observations

For the summer of 1816, a large number of surface station observations – at least by the standards of the early instrumental period – is available for Europe (see Tables 1 and 2, as well as Fig. S2 in the Supplement). For this study, a first source of pressure measurements is the International Surface Pressure Databank version 4.7 (ISPD; Compo et al., 2019; Cram et al., 2015; see Fig. S1 in the Supplement). A second source is daily to sub-daily measurements of temperature and pressure in Switzerland that come from the CHIMES project (Brugnara et al., 2020; Brugnara, 2022). The CHIMES dataset is based on the digitization of Swiss archival sources (Pfister et al., 2019). A total of 8 out of the 70 records (at 40 locations) in the CHIMES dataset cover the region of interest for the period from 5 to 12 June 1816 (see Table 1). As a special case, additional information from the records for Bern (observer Samuel Studer) was digitized by CHIMES and then quality-checked, adjusted, and made available by Hari (2021; see Fig. 1). A third source is observations of temperature and pressure from recently digitized stations in the Alpine area (see Brugnara et al., 2015, 2023; Fig. S2 and S3 in the Supplement and Table 1). Observations at locations within the innermost domain of the regional circulation model (see Sect. 2.3 and Fig. S2 in the Supplement) are summarized in Table 1. Additionally, there is a good number of European stations outside the innermost model domain (Table 2) for which pressure and also some temperature data are available (Brugnara et al., 2015). The pressure series compiled by Brugnara et al. (2015) have been supplied to ISPD (Compo et al., 2019; Cram et al., 2015).

Table 1 Weather stations with (assimilated) observational data of pressure and temperature between 5 and 12 June 1816. Independent records used for validation are marked with an asterisk. Irregular readings refer to potential missing information in the time series; parentheses indicate the time of the day for which missing information is most probable.

Table 2 Additional weather stations within outermost WRF domain with (assimilated) observational data of pressure and temperature between 5 and 12 June 1816.

Figure 1 Processed and original information from selected weather diaries in Switzerland. (a) Classes of cloud cover (grey symbols for cloud fraction from overcast (index 1) and mixed (0.5) to sunny (0) conditions) and amounts of precipitation ( filled blue circles; larger circles mean more precipitation) derived from twice-daily observation records in Geneva (Auchmann et al., 2012) for all months (day of month on x axis) of June between 1799 and 1821. The right margin summarizes instances of clear-sky conditions (where larger grey circles indicate brighter conditions) and of precipitation for each month. (b) As in panel (a) but with four classes of cloud cover (grey symbols) and three categories of precipitation (blue symbols), as derived from twice-daily observation records in St. Gall (Daniel Meyer; CHIMES) for all months in 1816. (c) As in panel (a) but with four classes of cloud cover (grey symbols) and observations of precipitation (filled blue circles) from twice-daily observation records in Bern (Studer; Hari, 2021) for June 1816.

For our data assimilation experiment, we used pressure and temperature data. Further available variables such as wind velocity and direction, precipitation type and occurrence, fresh snow, and cloud cover were used for validation purposes. As visible in Tables 1 and 2 and Figs. S2 and S3 in the Supplement, all major regions north and south of the Alps contribute data for assimilation, with the exception of eastern France. Note that we used pressure data for our regional assimilation that had already been assimilated in 20CR. This concerns the observations from ISPD; among them are three stations in the area of the European Alps (Table 1), namely Geneva, Turin, and Hohenpeißenberg, as well as all pressure observations shown in Table 2. We justify re-using these data by the fact that in the reanalysis, they serve to adjust the atmospheric state at low resolution, whereas, in our experiment, they represent a much smaller region with higher resolutions by a factor of about 2 to 20 for the outer- and innermost domains, respectively. They have less weight next to other stations, while they still provide very valuable information within our small network, especially on local effects that are resolved in the WRF model. Furthermore, individual measurements are often discarded from the assimilation if they show local effects not resolved in the reanalysis. For instance, the (bias-corrected) observation of pressure at Hohenpeißenberg for 6 June 1816 06:30 UTC was 885.70 hPa, which is the lowest but plausible value in our period of interest; it was judged to be too far from the first-guess pressure (−9.2 hPa; threshold in 20CR is 3.2 times the root of the sum of the squared variances of background and observation) by the algorithm and was thus not assimilated in 20CR (cf. Fig. S1 in the Supplement).

For an independent validation of pressure and temperature in our model simulations, we use the stations Delémont in the Swiss Jura and Augsburg in southern Germany. In addition, there are parallel time series from independent observers in Bern and Zurich, respectively (Bern Fueter vs. Bern Studer; Zurich Feer vs. Zurich Escher). In principle, these are also independent. Thus, we have two stations located very close to the assimilated station records and two stations that are located further away from the nearest assimilated station. On the one hand, a set of four stations is very small for relevant validation, but on the other hand, this allows us to use most of the rare and very valuable information for assimilation.

The quantitative station series were quality-checked manually prior to being used for regional data assimilation and validation. This included checking the metadata (e.g., coordinates, altitude, and observation time), converting data to modern units, and tracing gross errors in the data (e.g., temperature was too high by an order of magnitude). However, none of the used observation series are homogenized. Although we assume the raw records to be consistent over a period of only 10 d, potential errors must be considered and interpreted when showing the raw data. For assimilation, however, a simple correction approach was applied to correct for biases in the measurement series (see Sect. 2.4). In addition, the assimilation algorithm rejects values that are too far off from the first-guess simulated value. Hence, negative effects of potentially erroneous values on the assimilation can be considered to be rather small to negligible.

Furthermore, we use eye observations from selected weather diaries in Switzerland that were recorded in semi-standardized terminology (Auchmann et al., 2012; the St. Gall series from CHIMES; Hari, 2021). Adjustments included manual re-coding of the available information into classes (e.g., bright, partly cloudy, and overcast for cloud cover) or categories (e.g., spray, rain, and snow for precipitation types). Refer to Auchmann et al. (2012), Brugnara et al. (2015), or Brönnimann (2023) for information on how to best re-code the symbolic or word information, re-classify it, and attribute a plausible time to a record, among others. Given the inherent uncertainties, we consider the final data points as being of qualitative, complementary, and relative information which, when taken as a whole, may support or contradict our model outputs.

2.2 Reanalysis product

The NOAA–CIRES–DOE 20th Century Reanalysis version 3 (20CR; Slivinski et al., 2019) is used for synoptic analyses and as initial and boundary conditions for the downscaling experiments. 20CR is a four-dimensional global dataset that provides eight-times-daily fields of atmospheric variables on a ∼ 75 km horizontal grid from 1836 to 2015; a publicly available experimental extension with all 80 ensemble members goes back to 1806. The latest version of 20CR represents a substantial refinement of former versions which provided four-times-daily fields with ∼ 100 km horizontal resolution and temporal coverage from the mid-19th century onwards (Compo et al., 2011).

Here, we use the mean of the 80 ensemble members with exception of the analyses of the 20CR members presented in Sect. 3.1. Regarding initial and boundary conditions, previous studies for the same region found some deviations in variables such as the maximum wind speed from using the ensemble mean versus the individual members, but small smoothing effects in the pressure fields and overall limited benefit from applying an ensemble approach when comparing to station observations (Stucki et al., 2015, 2016). Regarding extreme events, smoothing effects may be more pronounced (e.g., Mahoney et al., 2022); the ensemble mean was found to provide accurate initial and boundary conditions, even if possibly less accurate than individual ensemble members (Michaelis and Lackmann, 2013). The variables taken as initial and boundary conditions encompass three-dimensional fields of temperature, humidity, geopotential height, pressure, and the horizontal components of wind speed, as well as two-dimensional fields of 2 m temperature and humidity, 10 m wind components, surface and sea level pressure, snow depth, skin temperature, sea surface temperature, a land/sea mask, and a sea ice flag. Furthermore, four layers of soil temperature and soil moisture from 20CR are used to initialize the regional model.

2.3 Regional circulation model

The non-hydrostatic Advanced Research Weather Research and Forecast Model version 4.1.2 (WRF-ARW; WRF hereafter; Skamarock et al., 2019) is used for dynamical downscaling from 20CR. The three nested limited-area domains have cell sizes of 27, 9, and 3 km, and the grid sizes are 127 × 109, 211 × 184, and 256 × 220 cells. The innermost domain is relatively large to avoid complex mountainous terrain at the boundaries, where possible. As for the downscaling procedure, the process of providing initial and lateral boundary conditions for each nested domain is as follows. For model initialization, the nested smaller domain receives the information from the coarser domain at the horizontal and vertical coordinates that both domains share. For all other coordinates (and for the outermost domain which typically does not share exact coordinates with the reanalysis), this information is spatially interpolated to the finer grid cells. The simulations are then incremented going forward in time until new information from the coarser domain is available, which is then fed in at the lateral boundaries of the nested domain.

There are 60 eta levels in the vertical, with a top level of 50 hPa. The model calibration builds upon previous WRF downscaling applications over the same region (Dierer et al., 2014; Gómez-Navarro et al., 2015; Stucki et al., 2015, 2016, 2018, 2020). The Thompson microphysics scheme (Thompson et al., 2008) is used for bulk microphysical parameterization and the Yonsei University (YSU) scheme (Hong et al., 2006) for the planetary boundary layer. The Kain–Fritsch scheme was used for cumulus parameterization in the larger domains (Kain, 2004) and turned off in the innermost domain. Spectral nudging (corresponding to a wavelength of about 1000 km) is applied to temperature, wind, and geopotential fields above the planetary boundary layer in the 27 km domain for consistency with large-scale forcing (von Storch et al., 2000). The use of spectral nudging has been analyzed and recommended in multiple studies, particularly for terrain with marked orography (Feser et al., 2011; Liu et al., 2012; Ma et al., 2016; Spero et al., 2014, 2018). In order to not restrain the simulation too much towards the large-scale forcing, the nudging coefficients were set to 0.0001 (Stauffer and Seaman, 1990). Tests with spin-up times of 24 h and 1 month revealed no conclusive benefits from the longer model spin-up. The WRF model is thus initialized on 4 June 1816 00:00 UTC, allowing for approximately 24 h of model spin-up before the cold spell starts. The simulation datasets are stored in hourly resolution.

Most of the above settings mostly correspond to common standards or even operational specifications, with exceptions. One of our more elaborate tests addressed the effects of using the standard modern-time land use scheme from the United States Geological Survey (USGS) vs. a reclassification of Anthromes v2 (Anthropogenic Biomes version 2; Ellis et al., 2010), which provides land use categories for the year 1800. One interesting feature appeared in that nightly temperature drops were more moderate in places where the local land use changed to more urban conditions in modern times (Fig. S4 in the Supplement); the nights are simulated up to 5 °C warmer in the built-up areas compared to cropland and woodland. The median shift during the night is 1.8 °C for Bern (3.6 °C for Geneva and 2.6 °C for Aarau), while the afternoon hours are around 0.25 °C warmer. Apart from this, only small changes in variables such as albedo, latent heat flux, or total column cloud fraction occurred in all three domains. Based on this, the historical land use scheme was used as the standard configuration.

2.4 Data assimilation system

In addition to dynamical downscaling without data assimilation described above, WRF simulations were combined with a 3D-Var (3D variational) data assimilation system (WRFDA; Skamarock et al., 2019). Simulations without (with) data assimilation are called NODA (DA) hereafter. The WRF setup for the DA simulations is the same as for NODA simulations, including spectral nudging in the outermost domain. With this combination of spectral nudging and subsequent regional data assimilation, we follow a number of previous studies that have already successfully adopted this technique for dynamical downscaling. For their regional 15 km Arctic System Reanalysis model, Bromwich et al. (2018) implemented nudging on temperature, geopotential height, and wind at wavelengths > 1000 km before using WRFDA for conventional observations, among others. Lin et al. (2021) combined spectral nudging with a 3D-Var assimilation of radar data for precipitation forecasts, and Yao et al. (2021) used atmospheric and snow data assimilation for springtime temperature simulations. For further details on the WRFDA system, see Barker et al. (2004, 2012), Huang et al. (2009), and Skamarock et al. (2019; their Chap. 11 and references therein).

Our basic idea for the assimilation was to mimic modern surface synoptic observations (SYNOP) of air pressure and temperature. The conversion into SYNOP implied assigning specific observation times for each series. While a number of records include time indications on (sub-)hourly scale, indications of the time of the day like “noon” or “sunset” were transformed using the R package suncalc (Thieurmel and Elmarhraoui, 2022). To mitigate the heterogeneity of observations and ensure compatibility with the 20CR input dataset, we implemented a simple bias correction, using a second harmonic fit (temperature) and a running mean difference (pressure), using data from the nearest grid cell and a time step from 20CR that is vertically interpolated to the station elevation. This bias correction was done for the whole year of 1816. While, for temperature, each time of the day (e.g., 00:00 and 03:00 UTC) was corrected separately, pressure was corrected using measurements from a moving window of 15 time steps. Potential timing errors were compensated by the bias correction procedure.

Observations of both temperature and pressure were assimilated in all three domains to ensure consistency between the nested domains as far as possible. In order to tune the WRF data assimilation system, a variety of analyses has been performed (not shown), including the assimilation of raw data, the assimilation of temperature only in the innermost domain, and a range of values for the observation errors and perturbation scaling. However, sensitivity studies in a technical sense were not possible, given the exploratory, “proof-of-concept”, nature of this study. In the applied data assimilation system, measurement errors were set to 1 K for temperature and 1 hPa for pressure (given bias-corrected observations) after some experience with measurement rejections by the WRFDA system. Observations are rejected if the innovation exceeded the observation error by a factor of 10 or for an altitude difference of more than 200 m. The background error covariance matrix was calculated using differences from 12 and 24 h forecasts for the same time steps (see Parrish and Derber, 1992). In order to obtain robust estimates of the model error that do not depend on a short simulated period or diurnal cycles, the differences between 12 and 24 h forecasts were calculated for 00:00 and 12:00 UTC over the period May–July 1816. The analysis was calculated every 3 h with an assimilation window of ±1.5 h around the analysis time step. Note that observations for assimilation were not available for all time steps; in these cases, empty station files had to be fed to WRFDA. As stated above, a limited set of tests regarding error estimates, data assimilation parameters, and bias correction has been performed (not shown). In all, we consider the described WRFDA configuration a subjectively optimal trade-off between enabling freedom of the model and restraining heavy dependence on observations, which would be problematic due to possible errors in the data.

3 Results

3.1 The meteorological situation in June 1816

In a first step, we describe the summer of 1816 from available studies and observations and, in particular, the cold-air outbreak over central Europe and the Alps between 5 and 11 June 1816. Long-term weather records by observers in Geneva, St. Gall, and Bern in Switzerland (Fig. 1) show that June 1816 was among the five rainiest months and clearly the least sunny month between 1799 and 1821 (Auchmann et al., 2012; Brugnara et al., 2022; Hari, 2021). The only day with reports of “plenty of sunshine” was 3 June 1816. Temperatures during the month of June were far below the seasonal average; measured afternoon temperatures in Geneva were below 10 °C during the first 10 d of the month (Auchmann et al., 2012; see their Fig. 2), leaving them among the lowest between May and October. Furthermore, it was extremely cold in the southeastern Alps: in Rovereto (Brugnara et al., 2023), 7 and 8 June 1816 were by far the 2 coldest summer days of the 1800–1839 period, i.e., colder than the next coldest day by 3 and 2 °C, respectively. According to historical information compiled by Pfister (1999), among others, frequent northwesterly flow in Switzerland led to long rainy episodes and persistent cloudiness, especially from 4 to 21 June. In fact, subdaily reconstructions of temperature and (interpolated) pressure by Brugnara et al. (2015; Fig. S5 in the Supplement) depict the Alpine region to the southeast of a surface high over the British Isles and a surface low over southern Scandinavia during this episode (see also Brohan et al., 2016) and show a temperature and pressure gradient along the Alps. Accordingly, documentary information report intermittent snowfall throughout June, even at low elevations (Pfister, 1999). There were 4 d of snow below 1500 m a.s.l. observed between 6 and 10 June 1816, and snow as low as 500 m a.s.l. (at Weggis, Lake Lucerne) was seen on 6 June. Digitized records from Bern (observer Studer; Hari, 2021), Aarau (Heinrich Zschokke) and St. Gall (Meyer) for 1816 support this information: May, June, and July each had > 20 d with precipitation; 5 to 11 June 1816 were all rainy days; and snow and rain was observed on 6 June 1816 in St. Gall. Further reports from the CHIMES raw material (Brugnara et al., 2020) include snowfall in Bern on 8 June and hail on 11 June (observer Fueter) and snowfall in Delémont on 6 June 1816 (unknown observer). Note that during the same period, the northeastern USA and Canada were also hit by snowstorms and a cold wave (Chenoweth, 2009).

Next, we aim to assess how well 20CR can reproduce and plausibilize the synoptic weather situation that we have outlined from the traditional reconstructions. In the first place, Fig. S1 in the Supplement shows that the spatial coverage of the North Atlantic region with instances of pressure information from ISPD is good, i.e., compared to other regions of the world and this point in time. Note that from this pressure information, the assimilation into 20CR generates a set of 80 ensemble members as deviations from the ensemble mean. Figure S6 in the Supplement delineates the 1005 hPa isobars over the North Atlantic region for all 80 members and the ensemble mean in 20CR for 6 June 1816 at 18:00 UTC. Overall, the sea level pressure (SLP) contour line of the ensemble mean is well located within the range of ensemble members, and the pressure minimum over the Alps is also well within the bulk of the members. Similarly, the ensemble mean of a 2 m temperature time series, extracted for a grid cell over Switzerland (7.7° E, 47° N), runs with the bulk of the members. Although slightly smoothed with regards to the ensemble variability, it clearly reflects the cold episode for this region.

Figure 2 Analysis of (a–d) temperature at 2 m above ground (shading; degrees Celsius), 10 m wind (grey vectors; reference vector approx. 10 m s⁻¹), mean sea level pressure (black contours; hPa), and geopotential height at 500 hPa (dashed light grey contours; decameters) over Europe and the Alps on 6 June 1816 21:00 UTC, as calculated from 20CR (a) and the WRF domains 01 (b), 02 (c), and 03 (d). Station observations of temperature are indicated as filled circles. Rectangles in panel (a) mark the WRF domains D01, D02, and D03.

In fact, 20CR indicates that the cold spell over central Europe started with a deep-cored low over Scandinavia on 3 June (not shown). While the low slowly retrograded in a southwesterly direction, cold-air masses were trapped over the subpolar eastern North Atlantic and trajectories show a south(-east)ward transport of cold air between this cyclonic system and a ridge to the west of the British Isles (not shown). On 6 June, a marked pressure gradient formed over the Alps, and the coldest air masses (around or below freezing point) reached the north side after sunset, whereas temperatures remained significantly higher south of the Alps due to foehn effects Fig. 2a). In the following days, the steering low made one more cyclonic pivot over Scandinavia to re-establish itself west of Denmark on 8–9 June (Fig. S7 in the Supplement). This led to a shift in the advection to a southwesterly direction over central Europe, then to a rather calm situation, and then back to northerly on 10 June. From here, the position of the Scandinavian low shortly renewed the cold-air advection towards the Alps before weakening and giving way to an extended Azores High on 11 June (not shown).

These analyses with the current 20CR version 3 can be seen as a continuation of the pioneering pre-20CR reanalyses by Brohan et al. (2016). They were already able to indicate higher pressure west of the British Isles and a low-pressure system over southern Scandinavia. However, their assimilation was based on a very limited network of station barometer observations. Thus, they indicated limited skill for the eastern North Atlantic. Furthermore, the synoptic situation resembles the one described by Brugnara et al. (2015) for July 1816, which points to repeated patterns of a zonal pressure dipole over northern Europe causing cold-air advection over central Europe and other European regions.

Hence, the current 20CR version 3 represents a big step forward, although our synoptic analyses still reveal coarse spatial patterns, and the ensemble mean may not necessarily concur better with the real atmospheric state than some individual ensemble members. It clearly captures the dynamics of a cold-air outbreak in two episodes between 5 and 11 June 1816 with a peak in the night from 6 to 7 June. This underpins the plausibility of the analysis and hence the quality of 20CR for this region and the early 19th century. In all, this also means that the 20CR ensemble mean is an adequate basis for the next step – the dynamical downscaling.

3.2 The cold period of 5–11 June 1816 in the WRF NODA and DA simulations

In this section, we explore the potential of WRF to produce more detailed weather maps that plausibly reflect the cold-air outbreak. For this, we zoom in from global to regional and local scales, following the three nested WRF model domains from the outermost to the innermost. With the refinement of the global information to regional scales, areas with below-zero temperatures at night become more distinct along the Alpine mountain ranges, central and northeastern France, western Germany, the Pyrenees and Cantabrian Mountains, and even southern England (not shown). For instance, analyses of the time step on 6 June 1816 at 21:00 UTC in domains 01 (not shown) and 02 (Fig. 2b) from the NODA simulation reflect and refine the predominant northerly flow of cold-air masses from the North Sea to the Alps and a split of the flow over northeastern France. Finally, domain 03 (Fig. 2c) incorporates the full Alpine bow and adjacent regions to the north and south. At this scale, local weather conditions are simulated at a horizontal grid size of 3 km. With a first swath of cold air reaching the north side of the central Alps on 5 June 1816 (not shown), areas of very low temperature are simulated on elevated Alpine terrain (e.g., below −20 °C in Valais) but also on the Jura Mountains and along the Alpine foothills during the following nights (e.g., from 6 to 7 June 1816). In contrast, many valleys and areas south of the Alps were even warmer than 10 °C. The same pattern is also evident from the temperature observations (Fig. 2a–c), although some biases with respect to 20CR and the NODA simulations are visible. This refined pattern of a north foehn situation becomes apparent from the SLP field, which delineates a north–south gradient across the Alps. Associated were northerly winds (> 10 m s⁻¹) across the Alpine rim (note even some counterflow) and a marked meridional temperature difference with comparatively warmer conditions south of the Alps. Accordingly, the distinct cloud layers at low and mid-levels on the northern side of the Alps were dissipated, and bright skies appeared on the southern side (Fig. 3a). Over the following days, the northerly flow changed to westerly and southwesterly from 8 to 9 June 1816 (not shown) to the north and west of the Alps, while the secondary low took shape over northern Italy (Fig. 3b and d). This induced substantially warmer air on both sides of the Alps, especially during daytime. On 11 June, the weather situation changed again (Fig. 3c). The Adriatic surface low weakened while the Azores High extended, and an associated low established over the North Sea (not shown). On the mesoscale, this brought a shift from westerly to calm, then increasingly northeasterly, winds and airflow to the north side of the Alps.

Figure 3 Analysis of (a–c) cloud cover (cloud fraction from 0 to 1) at low atmospheric levels (L), mid-levels (M), and high levels (H), as well as (total column) precipitable water (dashed red contour lines at 25 kg m⁻²), sea level pressure (smoothed dashed lines; hPa), and wind (grey vectors; every eighth in each direction is shown) (d) as in Fig. 2c. Simulations are with WRF NODA domain 03 for instances at times indicated above each panel.

Most of these spatial patterns and weather features come out similarly from the NODA and DA simulations (not shown). In contrast, the cold spell and a daily cycle of very low temperatures for the season are reproduced somewhat differently in the two simulations. Figure 4a–c shows temperature maps for the coldest point in time near sunrise on 7 June. The temperature maps show a considerable difference between the NODA and DA simulation. The DA simulation has lower temperatures by up to 4 °C for most of the domain and up to 10 °C in the inner-Alpine region.

Figure 4 (a–c) 2 m temperature maps for 7 June 1816 at 04:00 UTC. Shown are results for (a) the NODA run, (b) the DA run, and (c) the DA increment. Light grey arrows indicate the wind field. Station data for the corresponding assimilation window are indicated as colored circles. The straight line indicates the extent of the cross section. (d–e) Cross section (model orography is grey) of air temperature (shade; °C) from the Po plain (8.25° E, 45° N) to the Vosges mountains (6.85° E, 48.5° N) and up to 4000 m a.s.l. (meters above sea level) for 7 June 1816 at 04:00 UTC. Shown are cross sections for (d) NODA and (e) DA simulations. Freezing level and the 3 °C isotherm are emphasized as bold black lines.

The significantly lower temperatures in the DA simulation are also reflected in a cross section of the temperature field perpendicular to the Alpine bow, i.e., from the Po plain near Turin, Italy, to the Vosges mountains near Colmar, France, for the same instance in time (Fig. 4d and e). For this early morning, freezing levels drop from around 2000 m a.s.l. (meters above sea level) just near the Alpine rim to around 1200 m a.s.l. further north in NODA. Analogously, the 3 °C isotherm (approximately equivalent to the transient snowline on the ground) lies at around 1200 m a.s.l. at the northern alpine flank and drops to values to below 1000 m a.s.l. further to the north. Generally, both simulations are in line with the reports of recurrent snowfall below 1500 m a.s.l., arguably including snow and rain in or near St. Gall (around 800 m a.s.l.). However, only the DA simulation indicates freezing below 1000 m a.s.l. and a 3 °C isotherm reaching the lowest elevations north of the Alps. Thus, only the DA simulation closely fits with reports of snow as low as 500 m a.s.l. or even lower, for instance, at Bern (Swiss Plateau), Delémont (Jura), or Lake Lucerne (Alps). From this, we assume that the vertical temperature profile might be a little too warm in the WRF NODA simulation. We also infer that DA results in a more realistic representation of temperature and that, in general, our WRF simulations are able to reflect the effects of the large-scale cold-air outbreak on a local level. This leads us to a more detailed comparison of the simulations with observations for specific locations in order to gain a more general view over the differences between the NODA and DA simulations.

3.3 Verification of NODA/DA simulations with systematic observations and measurements

We assimilated as much of the very sparse station data as possible to obtain the most plausible simulation in terms of spatial and temporal dimensions. For this reason, we retained only four independent stations for a point-by-point verification of temperature and pressure. Two of them are quasi-parallel measurements (around 1 km apart) in Zurich and Bern and two stations are at least 50 km away from the next assimilated station. Due to these limitations, our verification is performed as follows. First, we consider variables that are independent of the assimilation. Precipitation, cloud cover, radiation, and wind can be compared with eye observations of, e.g., cloud cover, occurrences of rain or sunshine, and wind direction. Second, we include temperature and air pressure with a particular emphasis on the independent station records. Finally, we summarize the results with a quantification of the potential improvements from DA. First, we show the comparisons of observations with simulations for four locations in Switzerland from southwest to northeast, which are Geneva, Bern, Zurich, and St. Gall (Fig. 5; see Fig. S8 in the Supplement for the complete set of information). Comparisons with the WRF NODA/DA simulations are done using the grid point closest to the respective station location. Simulated pressure and temperature values were thus converted to station elevation for the quantitative assessment (cf. Sect. 2.4).

Figure 5 Meteograms of station observations and measurements (red colors) for (a) Bern with WRF NODA (lighter colors, grey, blue, and orange) and DA (darker colors) simulation output for the nearest grid point for the period between 4 and 12 June 1816 (x axis). The top plot in panel (a) shows observed cloud cover (red squares from bright with no fill to mixed with crosses and covered with fill) vs. simulated low-, mid- and upper-level cloud fraction (larger bars indicate more cloudiness). The second plot in panel (a) shows simulations of downward shortwave flux at the ground surface (orange lines; W m⁻²) and relative humidity at 700 hPa (dashed blue lines; ‰). The third plot in panel (a) shows thrice-daily observations of precipitation (red dots; qualitative) vs. simulated precipitation (vertical blue bars; mm). The bottom plot in panel (a) shows observations of wind direction (red vectors at unit length; north is up) vs. simulated wind direction (grey wind vectors; north is up) and velocity (black line and vector length; m s⁻¹ at 10 m above ground). Plots in the right column are as in panel (a) but for selected variables for (b) Zurich, (c) Geneva, and (d) St. Gall. Vertical red crosses in panel (b) indicate the relative humidity measurements on the ground for Zurich (dark for observer Escher; light for observer Feer), and red dots indicate precipitation in Zurich (darker for observer Escher; lighter for Feer). Vertical red bars in panel (c) show precipitation measurements in millimeters for Geneva. The tilted red cross in panel (d) means rain and snow for St. Gall. Grey plus signs mean confirmed observations, and minus signs mean no record of precipitation found.

Cloud cover observations from the Swiss stations agree about fairly bright skies for 4 June only. From there, intermittent overcast days are reported, but they are not coincident across the four stations. Note that in this context, there are inconsistencies between cloud cover and reports of rain in the observations for some stations. The station of Zurich (Fig. 5b) has two measurement series of relative humidity. The series by observer Feer appears to have values that are too low, possibly because of their being measured indoors. The measurements by Escher appear more realistic, and his readings at noon are among the highest of the summer (average of 81 % between 4–12 June, with maximum of 88 % on 6 June) and support the observations of dark skies. The simulations produce more consistent cloud cover; 4 June was mostly bright, followed by thick clouds from around 6 to 12 June but with clearly lighter cloud cover into 8 June. This is mostly consistent with the simulated values of tropospheric relative humidity and the reduced shortwave radiation; they indicate that the darkest days were on 6, 9, and 10 June (extending into 11 June). Note also that the finding of increased cloudiness over Geneva in the afternoons by Auchmann et al. (2012) is partly reflected in the simulated cloud cover (Fig. S8 in the Supplement).

The only station with precipitation measurements is Geneva (Fig. 5c). Qualitatively, precipitation observations agree very well with the simulated patterns of cloudiness, humidity, and precipitation and support the notion of the abovementioned cloudiest days. The other stations show intermittent rainfall from around 6 June onwards, with snow and rain on that day in St. Gall (Fig. 5d) but with no break on 8 June as in the simulations. This leaves us with the understanding that the simulations reproduce the observed weather evolution with brighter days at the beginning and dark and rainy days from 6 June (with a break on 8 June). However, the correspondence is not evident in all details, and some of the observation entries may be inaccurate. For instance, observations of cloud cover were noted twice or three times a day but at varying times, depending on the observer, and their assignment to a certain hour and the chosen categories by us are not concurrent.

The comparisons of observed and simulated wind further illustrate the qualities and limitations of the observations. Most of the information appears only reliable for periods with stronger winds and is thus hardly exploitable for comparisons. For Bern (and similarly for Schaffhausen; not shown), however, the wind information supports the simulated evolution of the mesoscale circulation with a predominant northerly flow until 9 June, followed by stronger southwesterly winds for 2 d, and then a drop in wind speeds and easterly (bise) wind directions towards the end of the period (Fig. 5a). Hence, we find that these particular records are a good example of astonishingly accurate observations, given that measuring highly variable wind parameters in an adequate quality is a difficult task, even with modern infrastructure and thoughtful site selection.

Differences between NODA and DA are mostly small to negligible, with the exception of slightly higher humidity and more precipitation in the DA simulation. In fact, the accumulated precipitation over the investigated period and across the Swiss Plateau is higher in DA (by up to around 20 mm; Fig. S9 in the Supplement), while it is lower over the parts of the Alps. This might be attributed to the assimilation of (in tendency lower) temperature and pressure values over the Swiss Plateau and its possible effects on the formation of convective precipitation in the simulation. Due to the lack of quantitative station data, a thorough assessment of simulated precipitation values is not possible.

Figure 6 Comparison of independent station records from (a) Delémont, (b) Augsburg, (c) Bern (observer Fueter), and (d) Zurich (observer Feer) with WRF simulations for the period between 5 and 11 June 1816 (x axis). Shown are 2 m temperature (left) and surface pressure (right) from the WRF NODA and DA experiments (lines) taken from the nearest grid point to the observation site, as well as raw and bias-corrected observations. Mean bias and root mean squared error (RMSE) between observations (bias-corrected) and the corresponding simulated values are indicated in each graph both for NODA and DA simulations. For reasons of completeness, the assimilated (and thus dependent) station data from Bern (c) and Zurich (d) are also indicated (grey is raw; black is bias-corrected). Note that in panel (d), the bias-corrected pressure observations (observer Escher; black dots) start only on 8 June, as this series starts on 1 June 1816, and corrections are made using a 15 d running average.

In contrast to wind, cloud cover, and precipitation, records of surface temperature and pressure were obtained from instrumental measurements at all locations. Our assessments for these variables feature the independent station series from Delémont and Augsburg, as well as the two independent parallel records from Bern (observer Fueter) and Zurich (observer Feer; Fig. 6). The difference in temperature readings for the parallel records is between 1 and 2 °C. Note that the bias correction applied to the data observed by Feer is larger than the one for the Escher data. Note also that in the original observations, there is a marked bias compared to simulations for some stations, especially where the elevation of the model terrain and observation site differs considerably (i.e., by more than 50 to 70 m). The reason for this bias is not entirely known and may lie in both simulation or observations. Evidence from the bias-corrected series (with 20CR; see Sect. 2.4), which do not exhibit a large bias, indicate that the latter is more probable. Warm biases of the order of 1–2 °C were in fact not uncommon in early instrumental temperature measurements during summertime. The surface pressure records from Studer and Fueter in Bern (Fig. 6c) agree very well with each other. The measurements by Fueter are higher by as much as 2.5 hPa, mostly because his values were not reduced to 0 °C. This causes an overestimation of about 2 hPa, which is well within the range of the general uncertainty in the SLP at that time of about 5 hPa which can largely be attributed to uncertainties in station elevation (see Brugnara et al., 2015; their Fig. 8). In fact, Studer even read a second barometer to verify the calibration; the difference between his two barometers in June 1816 is about 0.5 hPa. For the two series from Zurich (Fig. 6d), deviations of pressure readings appear larger than for Bern but with a deviation still within common measurement errors.

Our assessment of air temperature shows that the general agreement between observed and simulated near-surface (2 m) temperature is very good for the shown independent stations (Fig. 6; note that there is no parallel temperature series available for Bern). The daily cycle of the simulation mostly touches the measured values, and no general bias is apparent, although there are differences of up to 5 °C at times. For some instances, there are biases between raw and bias-corrected observations, with distinctive patterns depending on the measurement time (e.g., evening measurement at Delémont in Fig. 6a). This points to an error in the assumed observation times which was leveled out with the bias correction. Comparing NODA and DA simulations, DA clearly shows lower nighttime temperatures, especially for the coldest nights of 7 and 8 June. For 11 June, the 2 m temperature from the DA simulation is generally lower compared to NODA, and DA agrees better with the measured temperatures at all stations. Compared to the climatological means for the period 1981 to 2010 (provided by the Swiss Federal Office for Meteorology and Climatology MeteoSwiss), the daily maxima in the historical episode are only near or slightly above the daily minima for the modern period. These results are similar for the other stations (not shown). We conclude from the comparison to the historical measurements, as well as the modern climatology, that temperature seems well simulated in general, but some minima and maxima were probably more pronounced than seen in the simulation. Also, the pressure series from all stations agree very well with the simulated values. The pressure information supports the idea of brighter skies at the beginning and end of the period and clearly unsettled weather on 6 and 10 June. However, there are substantial differences between observations and the NODA simulation at the two pressure minima on 6 and 10 June. The readings differ from each other in time (up to 6 h), and they are all between 2.5 and 8 hPa lower than in the NODA simulation. In addition, the subsequent pressure increase is markedly stronger in the observations. The same deviations can also be observed for the other available stations. In the DA simulation, the surface pressure still do not reach the observed minima on 6 and 10 June, but they are more pronounced for these days (differences between 0 and 5 hPa), and there is a stronger pressure increase after the two observed minima. In short, both simulations reproduce the general evolution of surface pressure, although the variability in the pressure over time is arguably too small in the simulations and the temporal evolution too smooth. Data assimilation led to an improvement in this respect, and the values coincide better with observations than NODA.

Figure 7 Quantitative evaluation of simulated pressure and temperature at all observation sites. Independent (parallel) records are marked with an asterisk. Shown are bias (a) and the root mean squared error (b) for temperature, as well as bias (c) and RMSE (d) for surface pressure compared against homogenized station series. Annotations on the right-hand side of the respective graphs show median biases, as well as the RMSE over all records and over independent records, respectively.

Figure 6 shows also that mean bias and RMSE are substantially reduced by DA for the four independent stations and for both temperature and pressure, with very few exceptions. To get more spatial context, we calculated these measures for all available station records (shown in Fig. 7). For the majority, 14 out of 18 of the station records, the median bias of temperature is again substantially lower for the DA simulations compared to NODA. It drops by approx. 0.5 °C for all stations and by 1 °C for the independent stations. Similarly, RMSE values become substantially lower. On average, they drop by 0.2 °C for all stations and by 0.4 °C for the independent stations (Fig. 7a–b). Surface pressure biases (Fig. 7c) indicate an increase in the median bias by approx. 0.5 hPa for all stations, and an enhancement (reduction) by around 0.3 hPa for the independent stations. The RMSE (Fig. 7d) decreases by approx. 0.4 hPa with DA for all stations and by approx. 0.6 hPa for the independent stations. Overall, both non-independent and independent stations show the same patterns of improvements from DA. This reinforces our confidence that the spatial context of the local improvements is correct and that more independent stations would arguably show similar improvements. Furthermore, the differences in the quasi-parallel measurements in Bern and Zurich show approximate margins of error, which seem highly acceptable in a qualitative respect.

Note that the RMSE and the associated deviations from simulations are particularly high for temperature and pressure values from Marschlins and Turin, as well as temperature readings from Hohenpeißenberg and pressure readings from Aarau and Rovereto, arguably pointing to quality issues for the mentioned station series. In the case of Rovereto, the large pressure differences were found to arise from an uncertainty in observation times (not shown). The clearest outlier is the station of Marschlins. The data from this station were difficult to use in several respects. On the one hand, the data seem to be questionable (e.g., there are measurements of almost 30 °C; possibly on a sunny wall?); on the other hand, conversions into today's units of measurement could be erroneous. Third, the station is located on a valley slope where the elevations for the station itself, as well as in the WRF topography, are very uncertain (errors of ≥ 200 m are possible).

From these analyses, we can deduce that although state-of-the-art forecast verification cannot be done in our historical context, the comparisons reveal very good agreement between observation and simulations. For pressure, this might be expected, as pressure is assimilated in 20CR. Although deviations between observations and simulations are within the observation error in the early instrumental readings, simulated pressure variability is too smooth compared to the observations, especially in the NODA simulation. One reason for this might be the use of the ensemble mean from 20CR as atmospheric boundary conditions; using individual ensemble members might improve the simulated surface pressure. Generally, temperature observations also agree well with the simulations. Our quantitative assessments show that data assimilation clearly improves the results regarding both temperature and pressure. Strongly erroneous records are likely rejected by the assimilation algorithms, such that the simulation accuracy is not decreased.

4 Summary and conclusions

The Year Without a Summer in 1816 was characterized by exceptionally cold spells in central Europe. First, our analyses describe the meteorological situation of a cold spell over the European Alps in early June 1816 based on traditional reconstructions and in the 20th Century Reanalysis. Then, we provide weather simulations on an hourly temporal and 3 km (horizontal) spatial scale from two experiments, namely dynamical downscaling (i) without data assimilation (NODA) and (ii) with 3D-Var assimilation (DA) of bias-corrected pressure and temperature observations from stations in Switzerland and (central) Europe. Last, simulations are qualitatively and quantitatively compared to available early instrumental measurements and eye observations.

The cold-air outbreak over central Europe and the associated large-scale northerly airflow between a marked depression over Scandinavia and a ridge of high pressure west of the British Isles is well captured in the 20CR ensemble mean. The ensemble mean stands for a middle scenario within the bulk of the 80 realizations from the ensemble members and was found adequate to use as atmospheric boundary conditions for our downscaling simulations. Among others, this is in line with previous tests for the same region and the distant past (e.g., Stucki et al., 2018, 2020).

The quality of 20CR opened the field for experiments with dynamical downscaling, i.e., repeated nesting of the regional weather model WRF into the global 20CR data. Experiments with a land use scheme representing the early 19th century derived from the Anthromes project (Ellis et al., 2010) found differences in nighttime temperature for locations that have been urbanized since, while other variables only showed negligible differences compared to a modern scheme. We thus used the 19th century land use scheme for our analyses. Our downscaling simulations reproduce regional- to local-scale meteorological processes such as the foehn wind situation across the Alps with much lower temperatures on its northern side.

In general, the downscaled cloud cover, shortwave radiation and relative humidity agree well with eye observations of cloudiness or sunshine conditions. This is less true, of course, for timing and absolute values. In addition, relative humidity and precipitation were only measured at Zurich and Geneva, respectively, which means that comparisons of simulated versus observed precipitation are challenging. In all, we find that the general observations of a particularly cloudy and rainy episode between 5 and 11 June apparent in qualitative rain observations is reflected by the simulations, although they may disagree on intermittent brighter sky (instead on 7 June in the observations and on 8 June in the simulations). Also, we find that, arguably, the precipitation rates may be simulated too conservatively for the case of the only quantitative rain observation in Geneva. A robust evaluation is, however, difficult due to uncertainties in the observations and the strong spatial variability in the precipitation. The simulated evolution of advection is well reflected in some of the available wind observations. Indeed, some observers were set to deliver astonishingly accurate and meticulous meteorological information. Our validation analyses with independent (i.e., not temperature and pressure, which were assimilated in DA) weather variables showed that differences between NODA and DA are mostly negligible or small (e.g., cloud cover, shortwave radiation, humidity, wind, and precipitation). For temperature and pressure (assimilated variables), however, DA simulations are clearly closer to the observations. For instance, lower-pressure minima and a sharper rise after frontal activity are simulated with DA, whereas pressure variability is too small in the NODA simulations. Colder nighttime temperatures, lower freezing levels, and reported snowfall as low as around 500 m a.s.l. are only reflected in DA simulations. The general improvements with DA are also found in simple quantitative analyses of stations with independent and dependent temperature and pressure series. A careful selection and bias-correction of the assimilated station records is nevertheless crucial, as their quality largely affects the DA results.

In all, the analyses show that numerical weather simulations for this region and the early 19th century provide realistic atmospheric properties and dynamics, at a local kilometer-scale resolution. Despite a relatively sparse observational network and a rather short simulation period that does not allow a thorough validation of the general impacts of data assimilation on downscaling reconstructions, we conclude from our assessments that dynamical downscaling results are successfully improved by our data assimilation experiment. The continual historical observations and descriptions, available through digitizing efforts like the CHIMES project, are a prerequisite and extremely valuable for numerical studies of extreme weather and climate events of the past and for many more scientific purposes and practical applications. In this sense, the aspect of mutual exploitation becomes ever more important because better numerical methods allow the inclusion of more observations and more variables (such as information on wind, cloud cover, and rain/non-rain, as shown here), and this, again, will lead to better (i.e., regional, long-term, and high-resolution) reanalyses. With an envisaged extension to climatological timescales, our approach provides novel opportunities for the scientific community to learn from extreme weather and climate events as far back as 200 years. Given that global or regional gridded datasets may soon appear for periods beyond 1816 and more early measurements will become available, the prospects of soon entering the 18th century with such four-dimensional studies are very good.