Probabilistic forecasts of extreme heatwaves using convolutional neural networks in a regime of lack of data

Understanding extreme events and their probability is key for the study of climate change impacts, risk assessment, adaptation, and the protection of living beings. Extreme heatwaves are, and likely will be in the future, among the deadliest weather events. Forecasting their occurrence probability a few days, weeks, or months in advance is a primary challenge for risk assessment and attribution but also for fundamental studies about processes, dataset and model validation, and climate change studies. In this work we develop a methodology to build forecasting models which are based on convolutional neural networks, trained on extremely long 8000-year climate model outputs. This approach is parallel to weather model forecasting and has complementary scopes. Because the relation between extreme events is intrinsically probabilistic, we emphasize probabilistic forecast and validation. We demonstrate that neural networks have positive predictive skills, with respect to random climatological forecasts, for the occurrence of long-lasting 14-day heatwaves over France, up to 15 days ahead of time for fast dynamical drivers (500 hPa geopotential height fields), and also at much longer lead times for slow physical drivers (soil moisture). This forecast is made seamlessly in time and space, for fast hemispheric and slow local drivers. The method is easily implemented and versatile. We find that the neural network selects extreme heatwaves associated with a north hemisphere wave-number 3 pattern. We argue that this machine learning approach should be key in the future for quantitative process studies, model intercomparisons, and dataset studies. For instance, we find that the 2 meter temperature field does not contain any new useful statistical information for heatwave forecast, when added to the 500 hPa geopotential height and soil moisture fields. The main scientific message is that most of the times, training neural networks for predicting extreme heatwaves occurs in a regime of lack of data. We suggest that this is likely to be the case for most other applications to large-scale atmosphere and climate phenomena. Depending on the information to be learned, training might require dataset lengths as long as several thousands of years, or even more, for optimal forecasting skill. For instance, using 100-year-long training sets, a regime of drastic lack of data, leads to severely lower predictive skills and general inability to extract useful information available in the 500 hPa geopotential height field at a hemispheric scale in contrast to the dataset of several thousand years long. Even with several-thousand-year-long datasets, no convergence is observed in the predictive skills coming from hemispheric geopotential height fields. We discuss perspectives for dealing with the lack of data regime, for instance, rare event simulations and how transfer learning may play a role in this latter task.

Figure 1

(a) Snapshot of wind-speed velocity at the top of the troposphere, showing the jet stream over North America (from NASA). (b) Averaged horizontal kinetic energy at 500 hPa (midtroposphere) in the PlaSim model [28], showing the averaged northern hemisphere jet stream.

Figure 2

For all fields, we will use gridded data on the mid and high latitude northern hemisphere as represented by red meshlines. The figure also features in purple the area $\mathcal{D}$ (France). The North Atlantic Europe sector is represented in blue.

Figure 3

Snapshot of 2 m temperature anomaly and 500 hPa geopotential height anomalies. Temperature colormap is shown on the colorbar. The geopotential height anomalies are plotted via contour lines, colored green for positive and purple for negative with a separation of 20 m between the lines. The lowest level also corresponds to 20 m. This map displays a synoptic situation (daily averaged) on the first day of the strongest heatwave in the dataset.

Figure 4

Schematics of the CNN architecture. The numbers $a\times b\times c$ inside the boxes for Conv2D indicate the kernel sizes (a by b) and the number of filters $c$ in the convolutional layers. The box “Relu” actually consists of a sequence: batch normalization-activation (ReLu), and spatial dropout. Relu stands for rectified linear unit, it provides nonlinearity and is typically used to prevent vanishing gradients instead of sigmoid for the hidden layers.

Figure 5

Relevant climate fields for best prediction. All figures show the NLS [cf. Eq. (6)] vs the lead time $\tau$, for different combinations of the predictor fields. Fields are either integrated over the area of France (FI), masked over the France area (F) or over the northern hemisphere (NH). From $k\text{−}\mathrm{fold}$ cross-validation values, we plot the averaged NLS, plus or minus one standard deviation (shaded area and error bars): (a) single-field prediction: $T_{\text{FI}}$ (blue), $S_{\text{F}}$ (orange), $Z_{\text{NH}}$ (green), and $T_{\text{F}}$ (red); (b) fields combined pairwise for prediction: $(Z_{\text{NH}},S_{\text{F}})$ (blue), $(T_{\text{F}},S_{\text{F}})$ (orange), $(Z_{\text{NH}},T_{\text{F}})$ (green); (c) all three fields combined for prediction: $(T_{\text{F}},Z_{\text{NH}},S_{\text{F}})$ (orange).

Figure 6

Prediction skills when training the network with extra fields. The normalized logarithmic score NLS is on the $y$ axis. (a) Addition of stacked fields at extra time steps: $(T_{\text{F}},Z_{\text{NH}},S_{\text{F}})(t-\tau )$ (blue), $[(T_{\text{F}},Z_{\text{NH}},S_{\text{F}})(t-\tau ),(T_{\text{F}},Z_{\text{NH}},S_{\text{F}})(t-\tau -1)]$ (orange), $[(T_{\text{F}},Z_{\text{NH}},S_{\text{F}})(t-\tau ),(T_{\text{F}},Z_{\text{NH}},S_{\text{F}})(t-\tau -1),(T_{\text{F}},Z_{\text{NH}},S_{\text{F}})(t-\tau -2)]$ (green). The NLS features a very small improvement when adding the previous day information, but none if we add two previous days. (b) Addition of extra levels of geopotential height $[T_{\text{F}},Z_{\text{NH}}\left(500\text{hPa}\right),S_{\text{F}}]$ (blue), $[T_{\text{F}},Z_{\text{NH}}\left(500\text{hPa}\right),Z_{\text{NH}}\left(300\text{hPa}\right),S_{\text{F}}] [T_{\text{F}},Z_{\text{NH}}\left(850\text{mbar}\right),Z_{\text{NH}}\left(500\text{hPa}\right),Z_{\text{NH}}\left(300\text{hPa}\right),Z_{\text{NH}}\left(200\text{hPa}\right),S_{\text{F}}]$ (green). The NLS features a very small improvement, although not statistically significant, when adding more geopotential height fields.

Figure 7

Prediction skills versus dataset lengths and optimal geographical area. Panels (a) and (c) show PlaSim grid for the North Atlantic and Europe sector (NAE, a), and northern hemisphere mid and high latitude sector (NH, c), respectively. (b) Normalized logarithmic score versus lead time $\tau$, for neural networks trained with $(T_{\text{F}},Z_{\text{NH}},S_{\text{F}})$ predictors (blue) and $(T_{\text{F}},Z_{\text{NAE}},S_{\text{F}})$ predictors (red), and with datasets of length 7200 years (plain lines), 800 years (dashed lines) and 100 years (dashed-dotted lines). The results illustrate the lack of data regime, with very slow convergence of the prediction skill with the dataset length, and with a clear tradeoff between dataset length and size of optimal geographical area for best prediction.

Figure 8

Prediction skills versus dataset lengths and optimal geographical area for the 2 m temperature and soil moisture. Panels (a) and (c) show typical predictor fields, either masked over some restricted area [panel (a)], or over the whole mid and high latitude northern hemisphere [panel (c)]. Panel (b) Normalized logarithmic score versus lead time $\tau$, for neural networks trained with $(T_{\text{F}},Z_{\text{NH}},S_{\text{F}})$ predictors (blue) and $(T_{\text{NH}},Z_{\text{NH}},S_{\text{NH}})$ predictors (red), and with datasets of length 7200 years (plain lines), 800 years (dashed lines), and 100 years (dashed-dotted lines). The results illustrate the lack of data regime, with very slow convergence of the prediction skill with the dataset length, and with a clear tradeoff between dataset length and size of optimal geographical area for best prediction. The optimal area for 2 m temperature and soil moisture is the local one (France area).

Figure 9

Composites of the 500 hPa geopotential height maps $Z_{\text{NH}}$ (in meters), conditioned on committor values $p$ above the 99.9 percentile, at different lead time $\tau$. Panels (a), (b), (c), and (d) with values of $\tau =0,5,10,$ and 15, respectively. The regions of positive anomalies are indicated in red, while the negative anomalies are indicated in blue (we are using seismic colormap). The isolines are separated by a value of 20 m. The maximum geopotential height anomaly is indicated with a number colored in blue, while the minimum value is colored in red.

Figure 10

(a) Illustration of the $A\left(t\right)$(blue/red line) and $p\left[X\right(t),0]$ (in green) trajectories during summer of a specific year; red line corresponds to moments where the actual $A=A\left(t\right)>\alpha$ [see Eq. (1)], above 95 percentile threshold, while blue line denotes that $A=A\left(t\right)\le \alpha$. The filled green segments show the probability predicted by the neural network trained on $T_{\text{F}},Z_{\text{NH}},S_{\text{F}}$. This plot illustrates the general correlation between $A=A\left(t\right)\le \alpha$ and $T_{\text{F}},Z_{\text{NH}},S_{\text{F}}$. (b) Evolution of the committor function along a trajectory $p[X(t^{★}-\tau ),\tau ]$ as a function of lead time $\tau$, for a prediction of a potential heatwave at a prescribed physical time $t^{★}$. The orange points show the committor learned independently for each $\tau$, while the blue ones are obtained with transfer learning from one $\tau$ to the next. (c) Display of ${\sigma }_{\tau }\left(t\right)$, as defined in Eq. (10), for the same year as in panel (b), for the method with transfer learning (in blue) and without (in orange). We see that the latter is almost always above. This illustrates that variance between subsequent points is larger when transfer learning is not applied.

Figure 11

Normalized logarithmic score versus lead time $\tau$ for a neural network trained with the $(T_{\text{F}},Z_{\text{NH}},S_{\text{F}})$ predictors, with undersampling with $r=10$ (U10) (blue) and without undersampling (U1) (red). The two skill curves are within error bars of the experiment.

Figure 12

Benchmarks for the optimal normalized logarithmic score obtained by the proposed neural network for 5 percent heatwaves (blue) and 1 percent heatwaves (red).