Mesoscale modeling of a “Dunkelflaute” event

In the near future, wind and solar generation are projected to play an increasingly important role in Europe's energy sector. With such fast-growing renewable energy development, the presence of simultaneous calm wind and overcast conditions could cause significant shortfalls in production with potentially serious consequences for system operators. Such events are sometimes dubbed “Dunkelflaute” events and have occurred several times in recent history. The capabilities of contemporary mesoscale models to reliably simulate and/or forecast a Dunkelflaute event are not known in the literature. In this paper, a Dunkelflaute event near the coast of Belgium is simulated utilizing the Weather Research and Forecasting (WRF) model. Comprehensive validation using measured power production data and diverse sets of meteorological data (e.g., floating lidars, radiosondes, and weather stations) indicates the potential of WRF to reproduce and forecast the boundary layer evolution during the event. Extensive sensitivity experiments with respect to grid-size, wind farm parameterization, and forcing datasets provide further insights on the reliability of the WRF model in capturing the Dunkelflaute event.

In this paper, we focus on one such weather phenomenon called "Dunkelflaute" as it is rapidly becoming a major concern for the renewable energy community. 6 The word Dunkelflaute was coined by combining two German words "Dunkelheit" (darkness) and "Windflaute" (little wind) to describe heavy overcast skies and weak wind conditions. These meteorological events can last from a few hours to a few consecutive days. It is needless to say that under the influence of such a meteorological condition, little or no wind and solar energy can be produced.
On the 30th April 2018, an unexpected Dunkelflaute event occurred over the southern part of the North Sea and caused a large imbalance in renewable power generation and overall consumption. Given the acuteness of the situation, TenneT-the main transmission system operator for Germany and the Netherlands-had to issue an emergency alert in the Netherlands. 7,8 The crisis could not be avoided by simple load management or by making use of reserve power; instead, a substantial amount of electricity had to be imported from neighboring countries at high market price.
This Dunkelflaute event was not an isolated episode. As a matter of fact, over the past few years, several Dunkelflaute events occurred in Belgium, [9][10][11][12] Germany, 6,[13][14][15] and other neighboring countries. Some of them caused significant impacts on the power grids and electricity markets. There is no reason to believe that the occurrences of Dunkelflaute will subside in the future. Instead, with the ever increasing penetration of renewables in the power grid, the (negative) impacts of Dunkelflaute events will likely become more and more detrimental.
As a first step towards better forecasting of Dunkelflaute events, in this study, we investigate a recent Dunkelflaute event which occurred in Belgium. We analyze a diverse suite of observational datasets for detailed characterization. We evaluate the performance of a state-of-the-art mesoscale model, called the Weather Research and Forecasting (WRF) model, 16,17 in capturing this event. The organization of this paper is as follows. In Section 2, we briefly provide the meteorological background of the Dunkelflaute phenomenon. The selected case study is discussed in Section 3. The observational datasets and our atmospheric modeling framework are described in Sections 4 and 5, respectively. The simulated results along with in-depth analysis are documented in Section 6. At the end, we summarize our findings and elaborate on potential future research in this arena. A brief climatology of Dunkelflaute events in Belgium is documented in Appendix A.

| DUNKELFLAUTE: A METEOROLOGICAL PERSPECTIVE
The word Dunkelflaute does not exist in the vocabulary of meteorologists. Instead, they commonly refer to this phenomenon as the "anticyclonic gloom." 18,19 There are also localized names to describe this dull and drab weather phenomenon. For example, in the dialect of Lincolnshire, UK, it is known as moäky or moke. 20 One of the earliest references to anticyclonic gloom can be found in a sublime article by Captain C. K. M. Douglas. 21 Exactly 100 years ago, he conducted a comprehensive study on cloud characterization using aerial photography and stated: "It is not generally realised that when the sky is covered with a gloomy canopy of cloud, with the inevitable smoky haze over towns and for a considerable distance to leeward, one has only to ascend about a mile in order to enter a region with clear blue sky above, and a sea of white billowy cloud underneath, which stretches in all directions to a distant horizon which stands out sharply owing to the perfect visibility.

…
In winter the sky is very frequently overcast with a single sheet of low cloud varying in thickness from 500 to 2000 feet. This type of cloud is very characteristic of anticyclonic weather. There may be several days or a week of overcast sky, a state of affairs described by Mr. W. H. Dines as anticyclonic gloom." Over the years, our overall understanding of anticyclonic gloom has been steadily increasing. We now know that it is characterized by a high pressure system, extensive stratus and/or stratocumulus cloud cover, strong subsidence inversion, near-calm wind, low surface temperature, and possibly foggy nights. 22,23 We are also aware of the fact that anticyclonic gloom is predominantly a winter-time phenomenon. Although, it can happen during the summer months under appropriate synoptic meteorological conditions and sea state, for example, Galvin. 24 In spite of basic characterizations, more detailed knowledge (e.g., radiation and heat budgets, dynamical evolution, and dissipation) pertaining to anticyclonic gloom is severely lacking in the atmospheric science literature. For example, after the (inconclusive) studies by Priestley and Swinbank 25 and Robinson 26 around 1950, we have not come across any follow-up publication on heat budgets of anticyclonic gloom. In the era of advanced instrumentation (including remote-sensing) and cutting-edge numerical modeling (e.g., large-eddy simulation), such a critical knowledge-gap should not exist.
In this context, it is important to note that there have been several comprehensive modeling studies [27][28][29] probing both stratus-and stratocumulus-topped boundary layers. However, to the best of our knowledge, none of these studies focused on the weak wind (or calm) regimea necessary ingredient for the genesis of anticyclonic gloom. It is the sole purpose of the present study to investigate this specific regime using the WRF model.

| DESCRIPTION OF CASE STUDY
In January 2017, Belgium experienced a total of nine days of Dunkelflaute (anticyclonic gloom). [9][10][11][12] During these days, due to gentle breeze and overcast conditions, the energy production from wind and solar farms was far below their rated values. To further aggravate the situation, a nuclear power plant malfunctioned at the same time. Several measures, including more electricity generation from natural gas, additional power imports from neighboring countries, and other flexibility options, were undertaken to handle this acute problem. 30 Modeling by Elia of a similar extended period of low renewable energy generation (wind and PV) under a future 2040 scenario indicated that there could be an overall shortage of more than 1000 GWh of energy 10 which would be difficult to fill with existing storage technology.
The aggregated wind (onshore and offshore) and solar power production data for Belgium are publicly available from Elia (https://www.elia. be)-the transmission system operator in Belgium. The production data for the month of January 2017 are shown Figure 1. In lieu of any objective (quantitative) criterion in the literature, here, we define a period as Dunkelflaute when onshore wind, offshore wind, and solar power generation all fall below 10% of their respective nominal capacities. These periods are marked as lightly shaded gray regions in Figure 1. The weak wind conditions, occurring during nighttime hours (no solar power generation), are demarcated by a darker shade. In this paper, we focus on two of these Dunkelflaute periods: January 15-17 and January 22-25.
The surface analyses for January 15th and January 25th are shown on the top panel of Figure 2. It is clear that the Azores high was present in the eastern Atlantic on January 15th. At the same time, over the North Sea region, the gradients of surface pressure were relatively strong and were conducive to adequate wind power generation (see top-panel of Figure 1). In the following days, a stronger anticyclone (high pressure center) developed over continental Europe and moved slowly eastward. At its peak, it reached an impressive magnitude of 1042 hPa (not shown). On and around January 25th, the anticyclone was located over Germany. Given the vast expanse of the high pressure regions (refer to top-right panel of Figure 2), the pressure gradients were rather slack and led to extremely calm wind conditions. It is needless to say that neither the onshore nor the offshore wind farms in Belgium produced any energy on that day.
During both Dunkelflaute events of January 15th-17th and January 22nd-25th, the southern part of the North Sea was overcast with thick clouds (bottom panels of Figure 2), resulting in very little solar radiation reaching the surface. As shown in the bottom panel of Figure 1, the solar energy production was virtually negligible during these Dunkelflaute periods. Even though it is completely out of the scope of the present study, we would like to point out that the atmospheric circulation patterns over continental Europe were very unusual during January 2017 (refer to Appendix B). They caused extreme cold and violent storms at various places resulting in numerous deaths. 31 From that perspective, these Dunkelflaute events were far less sinister in nature.

| DESCRIPTION OF OBSERVATIONAL DATASETS
During January 2017, Belgium had only three operational offshore wind farms, C-Power, Northwind, and Belwind I, with a combined capacity of 712 MW. These wind farms and their associated wakes can be seen in the left panel of Figure 3. In this paper, we have used aggregated power production data from these wind farms in conjunction with solar farms to characterize the Dunkelflaute phenomenon. The temporal granularity of these datasets is 15 min. Note that these data have been post-processed by Elia to account for missing data or other data anomalies. 33,34 In addition to the power production data, we use meteorological data from several offshore and onshore locations. First of all, we make use of several lidar-based wind datasets which are publicly available via TNO (fomerly ECN, the Energy Research Center of the Netherlands).  The wind dataset from BWFZ has only a handful of missing samples. These measurements were validated against a cup-anemometer at Vlakte van de Raan wind station and were found to be of high quality. 35,36 In the HKZ region, two floating lidars were taking wind measurements at locations called HKZA and HKZB. Unfortunately, due to intermittent transmission failures, the wind dataset from HKZA has several gaps. The lidar at HKZB did not suffer from a similar data-loss problem. In fact, the wind data at HKZB were strongly correlated with conventional anemometer-based data from Lichteiland Goeree (LEG) and EuroPlatform (EPL) stations attesting to its high quality. 37,38 At LEG and EPL platforms, ZX wind lidars were also deployed and measured at several heights (from 90 to 315 m) every 10 min. 39,40 The locations of BWFZ, HKZ, LEG, EPL, and the offshore wind farms can be seen in the left panel of Figure 4.
Next, measurements from the Stabroek automated weather station (AWS) are utilized in this study. This AWS is operated by the Royal Meteorological Institute (RMI) of Belgium, and its location can be seen in the left panel of Figure 4. The recorded shortwave radiation data have a temporal resolution of 10 min. 41 Lastly, sounding datasets from three land-based stations (Herstmonceux, Norderney, and EDZE Essen, locations shown in the right panel of parameterizations; see Table 1 for details. The first five simulations start at 00 UTC on January 14, 2017, and last for a total of 14 days. The last simulation (called WRF−GFS+) is initialized at 00 UTC of January 21st. Model output are saved every 10 min.
All the simulations utilize 51 vertical levels with non-uniform grid spacing with the top of the model reaching approximately 16 km from the surface. The lowest grid level is approximately at 8 m from the surface, and there are 18 levels in the lowest 1 km of the model atmosphere. In every simulation (with one exception), grid nudging is applied above approximately 2 km in order to keep the simulations in sync with the largescale forcing data. In addition, the sea-surface temperature field is continuously updated throughout the simulations.
Numerous physical parameterization options are available in the WRF model to represent turbulence, land-atmosphere interactions, radiation, cloud microphysics, and other processes. Based on our past experience, we have used the following parameterizations: NOAH land surface model, 43 Rapid Radiative Transfer Model for Global Climate Models (RRTMG) for longwave radiation and shortwave radiation, 44 Table 2). The power and thrust curves of these turbines are described in Appendix C.
Three different large-scale forcing datasets are used for initial conditions (IC) and boundary conditions (BC). The ERA5 reanalysis dataset (horizontal grid size: 31 km; sampling rate: 1 hourly) is available from 1979 to present day from the European Centre for Medium-Range Weather Forecasts (ECMWF). 51 The ERA-Interim reanalysis dataset (horizontal grid size: 79 km; sampling rate: 6 hourly) is also available from ECMWF. 52 In addition to these reanalysis datasets, we also use an operational forecast dataset (horizontal grid size: 0.25 ; sampling rate: 3 hourly for 0 to 240-h period and 12 hourly for 240 to 384-h) from the Global Forecast System (GFS) to investigate if the Dunkelflaute periods could have been predicted in a real-time forecast scenario. The WRF−GFS and WRF−GFS+ simulations are initialized on January 14th and 21st, respectively, to investigate the impacts of prediction horizons on Dunkelflaute forecasting.
For the WRF runs involving ERA5, two nested domains are used (right panel of Figure 4). These domains are coupled in a one-way nesting mode. The outermost domain (d01) has a grid size of 9 km (domain size: 1890 km × 1674 km); whereas, the inner domain (d02) employs a grid size of 3 km (domain size: 819 km × 819 km). Given the coarser spatial resolution of the ERA-Interim and the GFS datasets, we deemed it necessary to use a three domain configuration. The largest domain (d01) for these runs uses a grid size of 27 km (domain size: 2673 km × 2943 km); it is not shown. The second domain (d02) and the third domain (d03) use grid sizes of 9 and 3 km, respectively. Their corresponding domain sizes exactly match the ones shown in the right panel of Figure 4.

| RESULTS
In this section, we compare the WRF model-based simulated results against various observational data. In addition, we investigate the sensitivities of the simulated results with respect to grid-size, IC/BC, and wind farm parameterization.

| Time series analysis
Measured wind speed time-series for the period of 14th-28th of January 2017 are shown in the top panel of Figure 5. It is evident that the wind speeds were mostly lower than 5 m s −1 during the Dunkelflaute periods. More interestingly, wind speed shears within the layer of 40-180 m were virtually absent for almost the entire two-week period; thus, it can be inferred that the lower part of the boundary layer was well-mixed. Starting January 27th, we do see some signatures of wind speed shear in the measurements. The simulated wind speed time-series from the WRF-ERA5 and WRF-ERA-I runs are shown in the bottom panel of Figure 5. The selected grid points are close to the BWFZ station. Output from both the runs show similar temporal evolution patterns and strongly agree with the measured data in terms of magnitudes of both wind speeds and wind shears. On January 21st, the simulated data exhibited a small level of shear which were not present in the measurements.
As illustrated by Figure 6, the observed and simulated wind direction time-series varied a lot during January 14th-28th. The WRF-ERA5 run has more-or-less captured the overall trend. Some deviations are however noticeable.
In the left panel of Figure 7, the WRF model-generated wind power production data are overlaid on top of the measured data by Elia. It is clear that the wind farm parameterization of the WRF model is able to accurately capture the magnitudes of the power production including the rapid ramp-down and ramp-up events. During the Dunkelflaute periods, the measured power production data had a few sporadic episodes of power generation; the simulated data were unable to capture such traits.
Given the spatiotemporally intermittent cloud patterns (see the bottom panel of Figure 2), rigorous validation of the WRF model-generated radiation data is a challenging task. As a poor man's choice, we have decided to compare the simulated (downwelling) shortwave radiation data against the observational data from the Stabroek station. From the right panel of Figure 7, it is clear that the WRF model significantly overestimated the magnitudes of shortwave radiation during the Dunkelflaute periods; in other words, the cloudiness was weaker in the model than in the reality. The measured and simulated air temperature time-series are compared in the left panel of Figure 8. Once again, we have selected the BWFZ location for comparison. The temperature data were measured at a height of 4 m. However, the simulated data are from 2-m height. In spite of the height difference, the measured and simulated data portray similar trends. As often happens in marine boundary layers, no sign of a diurnal cycle can be found in the temperature data during January 14th-28th. However, during this period, the air temperature dipped below freezing twice, due to the passages of cold fronts. Such temperature drops increased the energy demand, and in turn, the energy deficits during the Dunkelflaute event.
The simulated sea-surface temperature (SST) and surface sensible heat flux values from BWFZ are shown in the left and right panels of Figure 8, respectively. During the two-week simulation period, the SST decreased marginally. Because SST was almost always higher than the air temperature, there is expected to be positive sensible heat flux at the surface. Under such meteorological conditions, the surface layer is considered to be unstable (convective) and it promotes turbulent mixing. Wind shears decrease drastically due to mixing as depicted earlier in Figure 5.
During January 27th-28th, due to warm air advection, the air temperature became warmer than the sea-surface and created stably stratified conditions. Such atmospheric conditions are conducive to high wind shears as can be seen in Figure 5.

| Boundary layer structure
Measured profiles of potential temperature and wind speed spanning the entire boundary layer are not available over the southern part of the North Sea region. As a substitute, we utilize radiosondes launched from three neighboring locations over land (Herstmonceux, Norderney, and EDZE Essen; see right panel of Figure 4) to characterize the boundary layer structure. In Figures 9-11, we compare measured and simulated profiles for January 15th and 24th. For both days, we show profiles corresponding to 00 UTC (close to local midnight) and 12 UTC (close to local noon time). Simulated data from three runs (WRF-ERA5, WRF-ERA5#, and WRF-ERA-I) are considered. In addition, we also plot the vertical profiles which are directly extracted from the ERA5 reanalysis data.
Overall, all the simulations and ERA5 data more-or-less capture the magnitudes and shapes of the profiles. Among them, the ERA5 is a clear winner and the WRF-ERA-I performs the worst. The performances of the WRF-ERA5 and WRF-ERA5# runs are slightly poorer than ERA5. This reduction in performance is possible because the WRF runs only utilize ERA5 data during initialization and as boundary conditions. All the internal grid points of the WRF domains dynamically evolve without any data assimilation. The effects of grid nudging is felt only above 2 km.
Because the radiosondes were launched over land, the potential temperature profiles often portray the telltale signs of the circadian cycles. The observed wind speeds at higher altitudes were quite strong on January 15th and the simulated ones closely resembled them. However, the WRF results overestimated the turbine-layer wind speeds at Herstmonceux. On January 24th, the wind speeds subsided drastically. At all the locations, the observed and simulated wind speeds were approximately 5 m s −1 or even weaker in the lowest 1 km of the boundary layer.
In lieu of observed boundary layer profiles near the Belgian offshore wind farm region, we only document simulated time-height plots of wind speeds and potential temperatures in Figure 12. The estimated PBL heights are overlaid on these plots. From the temperature plot, it is clear that the PBL height was rather shallow (less than 500 m) for the large part of the 2-week simulation period. Furthermore, the temperature profiles are found to be uniform in height (i.e., well-mixed) within the boundary layer for most of the time. This result further corroborates our previous finding that the marine boundary layer near the BWFZ station was unstable during January 14th-27th of 2017.
The top-left panel of Figure 12 suggests that a sudden reduction of the wind speed happens on January 15th and persists until January 21st.
The gentle breezy condition returned again on January 22nd and lasted until January 27th. The depth of the weak flow field extends much higher than the PBL height. Thus, it is caused by a synoptic scale system (an anticyclone) and not modulated by boundary layer processes. The localized increase of wind speeds on January 21st could be due to an offshore low-level jet; it is a mere speculation and cannot be substantiated in this study.
In order to quantify the impact of the PBL schemes on the simulated results, we have computed the differences of simulated results from the WRF-ERA5# run (invokes the YSU PBL scheme) and the WRF-ERA5 run (utilizes the MYNN 2.5 level PBL scheme). The time-height plots of the differences of wind speeds and potential temperatures are shown in the bottom panel of Figure 12. In terms of potential temperature, the differences are not large. Most of the differences occur in the (subsidence) inversion zone overlying the boundary layer. The wind speed values also differ significantly within this zone. Within the PBL, the YSU scheme often produces slightly stronger wind speeds in comparison to the MYNN scheme.

| Sensitivity to grid-size and IC/BC
The observed and simulated hub-height (100 m) wind speeds from the BWFZ location are shown in the top panel of Figure 13. The simulated time series from WRF-ERA5 and WRF-ERA-I are quite similar. More interestingly, for both the runs, the simulated results are virtually insensitive to grid-sizes. The ERA5 reanalysis data match the observed data remarkably well, as found earlier in the context of vertical profiles.
In the context of (downwelling) shortwave radiation, the differences between different simulations are significant (see the middle and bottom panel of Figure 13). All the WRF-based results overestimate the magnitude of radiation during the Dunkelflaute periods. For other times, the simulated results are not too far off from from the observations.

| Sensitivity to wind farm parameterization
The numerical configuration and the physical parameterization options of the WRF-ERA5* run are identical to the default run WRF-ERA5.
Except, in the WRF-ERA5* run, the wind farm parameterization scheme is turned off. Whereas, in the WRF-ERA5 run, it is used in conjunction with 182 wind turbines. The differences in the simulated wind speeds and turbulent kinetic energy at the BWFZ station are shown in Figure 14.
The BWFZ station is located about 10 km from the Belgian offshore wind farms (refer to the left panel of Figure 4). It can experience wind farm wake effects when the prevailing wind direction is within 180-270 degrees. According to Figure 6, such wind directions happened a few times during Jan 14th-28th. During such periods (e.g., January 16th and 24th), the wake effects can be clearly seen in Figure 14. These differences are quite substantial in terms of both mean wind speeds and turbulent kinetic energy. For other periods, there are small differences as well.
This is due to the fact that in mesoscale simulations any perturbations (here imposed by the wind farm parameterization) in flow fields change future evolution (the so-called "butterfly effect").

| Reanalysis versus real-time
In addition to initial fields, the WRF model (or any other mesoscale model) requires boundary conditions spanning the entire simulation period.
For retrospective simulations, such boundary conditions can be extracted from the reanalysis datasets (e.g., ERA5 and ERA-Interim). Due to extensive data assimilation, such boundary conditions tend to be very accurate. However, in the context of real-time forecasting, such high-fidelity boundary conditions are not available. Rather, one has to use operational forecast data from a global model (e.g., GFS).
In this section, we compare the results from the WRF−GFS run against the default WRF−ERA5 run. We expect the accuracy of the WRF −GFS to be comparable with WRF−ERA5 run at the beginning of the simulation period; however, its performance will likely deteriorate with F I G U R E 1 4 Time-height plots of the differences between the WRF-ERA5 and WRF-ERA5* runs (d02 domains). The left and right panels represent differences of wind speeds and turbulent kinetic energy, respectively. The grid point closest to the BWFZ station is selected for all these plots. [Colour figure can be viewed at wileyonlinelibrary.com] increasing prediction horizon. The top-left panel of Figure 15 is in line with our expectation. During the period of January 14th-18th, in terms of wind power production, the performance of the WRF−GFS run is at par with WRF−ERA5. Afterwards, its results deviate significantly from the measured data from Elia. Even though the WRF−GFS run predicts the commencement of the second Dunkelflaute period rather accurately, it underestimates its duration. It predicts that the power output from the Belgian wind farms return to the nameplate capacity level by January 25th; in reality, such recovery happened after January 28th. In contrast, the WRF−GFS+ run, initialized on January 21st, significantly improves the quality of the wind power forecast for the second Dunkelflaute period (particularly during January 21st-26th) due to the reduction in prediction horizon (refer to bottom-left panel of Figure 15). In terms of (downwelling) shortwave radiation, the performance of the WRF−ERA5, WRF −GFS, and WRF−GFS+ runs are poor (right panel of Figure 15). More validation work, potentially involving satellite remote-sensing-based radiation data, is needed in this arena.

| CONCLUDING REMARKS
In this work, we simulate and characterize a Dunkelflaute (aka anticyclonic gloom) event via mesoscale modeling. An extensive suite of observational data assisted in the model validation. In addition to weak wind and cloudy conditions, we have found that the marine boundary layer was frequently well-mixed during this event. As a consequence, wind speed shears were negligible in measurements over the North Sea and also in corresponding simulated data. Over land, however, the observed and simulated profiles portray traits of stratification and wind speed shears.
For basic climatological characterizations of Dunkelflaute events, the ERA5 reanalysis dataset could be utilized owing to its global coverage, long-term availability, and high accuracy. However, this dataset does not include the wake effects of wind farms. In order to account for such effects, one should utilize the WRF model or other mesoscale models in conjunction with a suitable wind farm parameterization. Mesoscale simulations can also provide more advanced diagnostics (e.g., turbulent kinetic energy), and by virtue of their high spatial resolutions, they can resolve coastal effects.
In this study, we have noticed that in retrospective mode, some of the features (e.g., hub height wind speeds and power production) of the Dunkelflaute can be reliably simulated using coarse grid-sizes (e.g., 27 km). If this finding holds true for modeling of other Dunkelflaute events, then it will be possible to simulate these events with relatively low computational costs. In a real-time forecasting scenario, however, this specific event could not be predicted beyond four days. In our future work, we will find out if we can improve on the predictability of these events by coupling mesoscale modeling with deep learning approaches.
Before closing, we would like to mention that the renewable energy community is not the only stakeholder who is interested in a better understanding and forecasting capability of the Dunkelflaute phenomenon. It is also relevant for the air pollution community 53,54 and the astron-

APP E NDIX A: CLIMATOLOGY OF DUNKELFLAUTE EVENTS IN BELGIUM
We analyzed wind and solar power production data from Elia.be for the years 2013-2018. These data have a sampling rate of 15 min; thus, in a given year, we have approximately 35 000 samples. A particular sample is classified as a Dunkelflaute event if the production by onshore wind, offshore wind, as well as, solar farms fall below 10% of their respective nominal capacities during that particular 15 min period. We do not tag a sample as a Dunkelflaute event if solar power production is exactly zero and it happens outside the time-window of 09 UTC-15 UTC; this way, we effectively exclude nighttime conditions from the climatological analysis. It is clear from Figure A1 that every year all the three renewable power generation sources drop below 10% of their capacities for a substantial period of time (approximately ranging from 30% to 50%). Typically, offshore wind power production has a higher capacity factor in comparison to its onshore counterpart. On average, the Dunkelflaute events account for around 7%-8% of the time per year. These numbers do not vary much across the years.

APP E NDIX B : ANOMALY PATTERNS OVER EUROPE
The synoptic condition during the month of January of 2017 was quite unusual as can be seen in Figure B1; see also Dunstone et al. 56 These anomaly plots are created using the ERA5 data. First, a climatological mean is computed for the month of January from ERA5 data spanning the years 1981 to 2010. Then, this mean is subtracted from the monthly average of January 2017.
It is evident that during January 2017, over the North Sea, sea level pressure was a lot higher than the climatological mean over the North Sea region. Similar behavior was observed in the case of 500 mb geopotential height. Both zonal and meridional winds showed negative anomaly for the southern part of North Sea and surrounding areas. Extremely cold conditions prevailed over the southern part of continental Europe.

APP E NDIX C : POWER AND THRUST CURVES
The power and thrust curves of the wind turbines from the Belgian offshore wind farms are extracted from various sources (see Table C1) and are plotted in Figure C1. In some cases, in lieu of actual turbine data in the public domain, we had to utilize data from proxy turbines with similar characteristics.