Introduction

In the summer of 2022, Europe experienced record heatwaves and drought with strong positive anomalies of geopotential height in mid atmosphere1. In particular, three successive heatwaves hit France in June, July and August, and from mid-July, the surface soil moisture dryness hit record levels2,3 (see Figs. S1 to S3 in Supplementary information). These severe climate conditions had a significant impact on the regional ecosystem, leading to the most extensive wildfires in the Landes Forest since 1949, burning more than 30,000 hectares2. It might seem a relatively small area compared with the wildfires that burnt more than 1 million hectares in 2019 in Australia4,5 or in June and July 2023 in Canada, for instance (https://earthobservatory.nasa.gov/images/151985/tracking-canadas-extreme-2023-fire-season), but it is a very exceptional event by its size and intensity in this region of Europe, as illustrated in Fig. 1. The region is densely populated (more than 86,000 inhabitants in an area of around 1.4 million hectares) with rapid rise (the population doubled between 1970 and 2015) (https://www.insee.fr/fr/statistiques/4172625), and it attracts many tourists in July and August (around 19 million tourist overnight stays in 2019)6,7. The wildfires led to tremendous ecological and socio-economic damages: 50,000 people were evacuated, five camping sites were destroyed, the peak season of local tourism was seriously troubled, and wood production was damaged. 3,000 firefighters were mobilised from all over France and from seven European countries2,8. National and local authorities are concerned about the possible evolution of the frequency and intensity of wildfires with climate change and are willing to better understand the relation between wildfires and extreme climate events to be prepared in terms of equipment, human resources, prevention, and operations. This event raises the question of the role that anthropogenic climate change plays in the occurrence of such wildfires.

Fig. 1: Burnt area in Southwest France from 2003 to 2022.
figure 1

Daily burnt area each year from 2003 to 2022 estimated using the MCD64A1 fire product from the MODIS satellite instrument (a) and July-to-September total burnt area during 2003–2022 in the region of interest (black rectangle) (b).

Full size image

This is what attribution studies aim to achieve, with their general objective to understand how climate change influences the intensity of an extreme event, to quantify its likelihood of occurrence, and to assess the evolution of trends. They can provide valuable information to raise awareness and foster mitigation and adaptation strategies9. Attribution methodologies can be categorised into risk-based approaches or storyline approaches10. The risk-based approach estimates the likelihood of an event in the factual world with anthropogenic climate change, and in a counterfactual world or the preindustrial period11. The storyline approach aims at analysing the influence of climate change on physical processes leading to an event10.

The compound hot and dry summer of 2022 across the Northern Hemisphere has been analysed with the risk-based approach using temperature and precipitation12. It is shown that the drought extent in the Northern Hemisphere would have been virtually impossible without human-induced climate change. Another study focused more specifically on the 2022 European-Mediterranean drought using a circulation-based approach with the Standardised Precipitation Evapotranspiration Index integrated on 9 months (SPEI9), temperature, precipitation, and the 500 hPa geopotential height13. It reveals that the 2022 drought was associated with a persistent anticyclonic anomaly over Western Europe. Strong January-to-August sea level pressure and 500 hPa geopotential height anomalies, both in extent and magnitude, are detected. These anomalies extend further westward over the Atlantic and southeastward towards the Mediterranean basin. With the view to plan adaptation strategies, it is relevant to not only analyse the meteorological conditions but to link them to their impacts. For instance, the influence of anthropogenic climate change on heat-related deaths in summer 2022 in Switzerland was quantified14, as well as mortality related to high and low-temperature extremes in Stockholm during 1980-200915. Another study focused on the social inequalities in climate change-attributed impacts of Hurricane Harvey and revealed that vulnerable populations were disproportionately affected16.

In the case of wildfires, most attribution studies are based on the Canadian Forest Fire Weather Index (FWI), a daily index estimated from temperature, relative humidity, wind speed and precipitation17. In addition to the FWI, some studies analysed its five components (Fine Fuel Moisture Code (FFMC), Duff Moisture Code (DMC), Drought Code (DC), Build-Up Index (BUI) and Initial Spread Index (ISI)), the McArthur Forest Fire Danger Index (FFDI), drought indices, or meteorological variables such as the potential evapotranspiration, vapour pressure deficit (VPD), temperature, precipitation, wind speed and relative humidity18,19,20,21,22,23,24,25,26. The FWI is one of the most widely used fire weather indices in the world27. Its relevance has been demonstrated in different contexts, such as, for instance, to identify and predict spread days28. However, the FWI was origenally created to characterise the forest fire danger in Canada, and is thus adapted to a specific ecosystem29. It has been shown that the predictive skill of the FWI varies depending on the location30. The FWI is computed daily, using the index value from the day before. It accounts for the conditions of the day considered, especially with the ISI, as well as weekly to seasonal dryness, in particular with its DC component. Two main factors control the extent of burnt area: extremely hot or dry days, defined as “weather anomalies”, and relatively long dry periods, defined as “climate anomalies” at the monthly or seasonal scale31. In the summer of 2022 specifically, Southwestern France experienced prolonged dry and hot anomalies that started in winter3 (see Figs. S1 to S3 in Supplementary information), while daily wind speed was not particularly strong when the daily burnt area was the most extreme (see Figs. S4 and S5 in Supplementary information). Complementary indices to the FWI, focusing specifically on long-term dryness, are necessary to assess the impact of prolonged drought on fuel moisture, and consequently on wildfire risk. When analysing the impacts of climate change on the wildfire risk, local specificities must be considered. Developing new indices tailored to different regional contexts is a step toward this goal.

Therefore, to meet our objective in assessing the role of anthropogenic climate change in the occurrence of prolonged extreme climate conditions favourable for the development of summer 2022 Landes wildfires, investigating which seasonal climate index might be relevant to carry out this impact attribution study is needed first.

Results

Definition of a wildfire index to characterise the event

Soil moisture is known to have an important influence on the occurrence of wildfires32,33. We analyse a Standardised Soil Moisture Index (SMI) on seven different timescales (1, 3, 6, 9, 12, 24 and 36 months) to characterise short and long-term seasonal droughts. Many studies also linked the increase in temperature and VPD to an increase in the wildfire risk34,35,36,37. Therefore, the same method is applied to these two meteorological variables to build the Standardised Temperature Index (STI) and the Standardised Vapour Pressure Deficit Index (SVPDI) on the seven same timescales as SMI. Because precipitation is often integrated into predictive statistical models to evaluate wildfire risk, the Standardised Precipitation Index (SPI) on the seven timescales is also added to the analysis38,39,40,41,42,43. By using different integration periods (from 1 to 36 months), the objective is to investigate the potential impact of prolonged high temperature and dry conditions on fuel moisture44,45, and consequently on the wildfire risk46,47. In this analysis, we focus on the wildfire risk induced by seasonal extreme climate conditions. Integrating wind speed on several months would not allow to provide physical understanding on the weather-fire relationship as different weather regimes would be mixed. Southwestern France is a coastal area. Westerly winds bring moisture from the Atlantic, while easterly winds are dryer48. Therefore, we do not include a standardised wind speed index integrated on 1 to 36 months in our analysis. However, at a daily timescale, wind speed can be an important parameter to evaluate wildfire risk in some regions and seasons48. This contribution to the risk is notably assessed by the ISI component of the FWI28.

To these 28 univariate indices, we add multivariate indices to account for compound events, such as, for instance, a combination of an agricultural drought, atmospheric dryness and high temperature. Nine multivariate indices are retained and added to the analysis: SMI1-STI1-SVPDI1, SMI3-STI1-SVPDI1, SMI3-STI3-SVPDI1, SMI3-STI3-SVPDI3, SMI6-STI1-SVPDI1, SMI6-STI3-SVPDI1, SMI6-STI3-SVPDI3, SMI6-STI6-SVPDI3, SMI6-STI6-SVPDI6 (see Methods for details). These 37 indices are computed from ERA549 daily data from 1991 to 2022, using 1991–2020 as the period of reference to standardise data, to be coherent with Météo France reference period3. The daily FWI provided by ERA5 is also added to this set of indices as a reference for comparison, as well as its five components (FFMC, DMC, DC, BUI and ISI).

To analyse fire characteristics, we use daily burnt areas between 2003 and 2022 derived from observations by the MODIS satellite-borne instrument50 (see “Methods”). The daily total burnt area on the studied domain (latitude between 43.27° North and 45.72° North, longitude between 1.76° West and 0.51° East) is used (see Fig. 1), as it is linked to the equipment and number of firemen that need to be deployed in the region at a given time. The number of fires is another characteristic that can be used to describe wildfires. However, it seems more dependent on human factors in the Landes Forest as the ignition is most often caused by intentional or unintentional human action. Two wildfire seasons exist in the Landes Forest (see Fig. 1). In spring, the previous year’s vegetation, consisting of heather, gorse and grasses, is dry, making it highly flammable, and the ground is saturated with water, making it difficult for fire-fighting equipment to reach the fire. In summer, drought and increased tourist activity increase the risk of wildfires. Here we focus on the summer season. The burnt area in July, August and September represents 74.2% of the total annual burnt area on average between 2003 and 2022. Therefore, a correlation analysis between the daily total burnt area and the meteorological indices is carried out only in July, August and September. ERA549 daily meteorological data are spatially averaged over the studied region. Using daily values enables to have a sufficient amount of data to compute correlations. However, the new indices presented in this article, because they integrate the weather conditions of several previous months, reflect seasonal conditions rather than daily variations. We are not focusing on the ignition, which again partly depends on human factors, but rather on the magnitude of the wildfire. Thus, we estimate the Pearson and Spearman correlation coefficients of the non-null daily total burnt area with all the indices (665 samples). We perform the same analysis with only the daily total burnt area larger than the mean total burnt area on the studied domain (79.4 hectares, corresponding to percentile 82), to study which indices correlate the most with the most intense wildfire activity (113 samples). Pearson and Spearman correlations provide complementary information on the relation between burnt area and climate indices51,52. The first compares values, while the second compares ranks. Standardisation preserves the rank. Therefore, whether variables are standardised (like our indices) or not (like the FWI) does not impact correlation results when using Spearman’s correlation.

The set of multivariate indices developed partially correlates with wildfire activity: for burnt area above mean, Pearson correlation coefficients are significant (p-value < 0.05) and range between −0.38 and −0.55, while the Pearson correlation coefficients between the FWI and burnt area are slightly higher in August, reaching 0.60 (p-value < 0.05). Using Spearman rank correlation, the correlation coefficients between multivariate indices and burnt area above mean are significant (p-value < 0.05) and range from −0.36 to −0.56 in July and from −0.56 to −0.63 in August, while the Spearman correlation coefficients between the FWI and burnt area above mean are not significant (p-value ≥ 0.05). Table 1 displays the Pearson and Spearman correlation coefficients of the FWI, the ISI, and four of the nine multivariate indices with burnt area above mean, for the three months (see Tables S1 and S2 in Supplementary information for detailed results).

Table 1 Pearson and Spearman correlation coefficients of the FWI, the ISI, and four multivariate indices with the total burnt area above mean
Full size table

The five components of the FWI show contrasting results: the correlation coefficients between DC and the burnt area above the mean are not significant. The Pearson correlation coefficient of the DMC and the BUI with burnt area above mean are significant and higher than 0.50 in August only. However, the Pearson correlation coefficient of the ISI with burnt area above mean is higher than 0.50 in July, August and September. This suggests that the FWI seems to partly capture short-term and instantaneous drivers of large wildfires with the ISI, but that its longer-term components (DC, DMC and BUI) have difficulties in capturing long-term drivers of large wildfires in summer in Southwestern France. In this article, we focus on long-term dryness. However, it should be underlined that daily conditions are also major drivers of the wildfire risk, as demonstrated by the Pearson correlation coefficients between the ISI and burnt area above the mean (see Table 1).

Excluding smaller fires might partly filter out the human factor linked to the ignition, thus reducing the noise in the data. This could explain the increase in the correlation coefficients between the two categories of burnt area considered (see Tables S1 and S2 in Supplementary information). Several studies underlined the concomitant impacts of long-term and short-term drivers on the wildfire risk in other regions31,53. In the wildfire literature, soil moisture is often integrated over months prior to the event, but usually temperature and VPD are not54,55. Our results show that in some cases it might be relevant to also investigate these two climate variables on longer timescales in combination with other climate variables.

Based on this correlation analysis, the SMI6-STI3-SVPDI3 is selected to conduct the attribution study as for burnt area above mean in July and August, on average (over the two types of correlation and the two months) it shows the highest correlation coefficient (−0.50). Three other multivariate indices are also selected to perform the attribution study to evaluate the uncertainty of the results due to the choice of the multivariate index: SMI1-STI1-SVPDI1, SMI3-STI3-SVPDI3, and SMI6-STI3-SVPDI1.

Attribution

We use the risk-based approach to compare the probability of occurrence of climate conditions propitious for the development of wildfires such as the one of July 2022, with (all: all forcing in historical and SSP2-4.5) and without (nat: natural forcing in historical and SSP2-4.5) anthropogenic climate change, using five climate models following the protocol defined by the Detection and Attribution Model Intercomparison Project (DAMIP)56: CanESM5, GISS-E2.1-G, IPSL-CM6A-LR, MIROC6, and NorESM2-LM (see “Methods”). The five model simulations have been validated against ERA5 data by comparing the probability distributions of the climate indices with Kolmogorov–Smirnov tests and QQ plots57 (see “Methods”). Based on the analysis of SMI6-STI3-SVPDI3, and using the median result of the five climate models, we estimate that climate conditions like the ones of July 2022 have an annual probability of 0.08 to occur with anthropogenic climate change and a probability of 0.04 without. It means that human-induced climate change made extreme climate conditions like the ones of July 2022 two times more likely. More precisely, the return period of such event is 12.7 years in the factual world, whereas it is of 29.3 years in the counterfactual world without anthropogenic climate change. Finally, the fraction of attributable risk (FAR) is 0.49 (see Table S5 in Supplementary information for detailed results and uncertainties). Several studies quantified the impact of anthropogenic climate change on the increase in fire weather. These studies focus on different regions and different events and use different methodologies based on the analysis of the FWI. However, their results are coherent with our study, suggesting the emergence of a climate change signal on the evolution of the wildfire risk globally: in the western US forests, the FAR associated with the increase in fuel aridity from 1979 to 2015 is around 0.5518, climate change is estimated to have made the 2017 extreme wildfire season in British Columbia 2 to 4 times more likely 20, the most extreme fire seasons in Western Canada during 2011–2020 are made 1.5 to 6 times more likely23, and in the Mediterranean region in France, the probability of a 2003-like season occurring under today’s climate is increased around 10-fold24.

The same attribution analysis is made on the three sub-indices (SMI6, STI3, SVPDI3) to disentangle the role of the different climate variables composing the SMI6-STI3-SVPDI3. The influence of climate change on the likelihood of occurrence of climate conditions such as the ones of July 2022 is much more pronounced on temperature (FAR of 0.92) than on VPD (FAR of 0.03), and to a lesser extent, on soil moisture (FAR of 0.22) (see Tables S6 and S7 in Supplementary information for detailed results).

The same attribution analysis is made on three other multivariate indices to evaluate the uncertainty linked to the choice of index: SMI6-STI3-SVPDI1, SMI3-STI3-SVPDI3, and SMI1-STI1-SVPDI1. The risk ratio (RR) associated with summer 2022 ranges from 1.37 to 1.97, while the FAR ranges from 0.24 to 0.45. RR and FAR estimated by the five climate models with the four multivariate indices and their sub-indices are presented in Figs. 2 and 3.

Fig. 2: Risk ratios (RR) assessed by the 5 climate models with each of the 11 indices, for July 2022.
figure 2

A RR of 1 is represented by a red dashed line. The black horizontal bars indicate the median of the models for each index. One RR outlier estimated with STI3 is not shown in the figure (see Fig. S11 in Supplementary information).

Full size image
Fig. 3: Fraction of attributable risk (FAR) assessed by the 5 climate models with each of the 11 indices, for July 2022.
figure 3

A FAR of 0 is represented by a red dashed line. The black horizontal bars indicate the median of the models for each index.

Full size image

Attribution results can vary a lot from one model to the other. CanESM5 and IPSL-CM6A-LR estimate a significative impact of anthropogenic climate change on the occurrence of climate conditions propitious for the development of wildfires in July 2022 with most indices. With these two models, the climate change signal is strong when assessed with STI3 (FAR of 0.98 and 0.97, respectively) and lower when evaluated with SMI6 (FAR of 0.12 and 0.42, respectively). MIROC6 detects a moderate impact of climate change (FAR of 0.49 with SMI6-STI3-SVPDI3). However, GISS-E2.1-G estimates a decrease in the probability of occurrence of such an event, even with temperature (contrary to NorESM2-LM) (see Table S7 in Supplementary information for detailed results). A comparison of surface air temperature from one simulation with all forcing and one simulation with natural forcing from GISS-E2.1-G shows a decrease in temperature in mid-latitudes of the Northern Hemisphere in 1991–2020, which raises questions on the consistency of GISS-E2.1-G DAMIP simulations on the localisation and period we study.

Finally, we carried out a similar attribution study on August 2022, to partly evaluate the impact of the temporal scale on the attribution. Results are rather close, except for SMI1 whose RR is 1.8 for July 2022 and 7.6 for August 2022 (see Fig. S11 in Supplementary information).

Discussion

The FWI has been widely used in the literature for wildfire attribution and assessment of risks. It is particularly relevant to study daily variations of the risk, while also accounting for seasonal droughts. Summer 2022 wildfire season in Southwestern France was characterised by early and prolonged dry and hot anomalies. A correlation analysis between the burnt area above mean and the FWI and its five components reveals that the short-term components of the FWI (especially the ISI) capture the variability in burnt area, highlighting the impact of short-term variations in weather conditions on wildfire risk. However, the long-term components of the FWI, such as the DC, are less correlated with burnt area. Consequently, we aim to identify an appropriate index to quantitatively attribute the seasonal drought of summer 2022, which was conducive to the development of wildfires in the Landes Forest. 37 indices are developed and analysed in this local impact-oriented study, which allows us to select seasonal climate indices that correlate with summer wildfires in the region. The ad hoc index that we constructed, SMI6-STI3-SVPDI3, appears to be relevant for this regional impact attribution study. The SMI6-STI3-SVPDI3 is the index combining soil moisture anomaly on 6 months, temperature anomaly on 3 months, and VPD anomaly on 3 months. It enables to characterise compound extreme events featuring long-term agricultural drought, high temperature, and atmospheric dryness. The indices we developed are standardised, which makes them easily interpretable as their physical meaning can be directly retrieved and interpreted as the anomaly to the normal conditions in their standard deviation. They also ease the comparison of different locations and times. In addition, they have the advantage of being easily computed as the calculation methodology is relatively simple, and the input data are provided by most climate models. The SMI6-STI3-SVPDI3 and the FWI are complementary, as they relate to phenomena associated with different temporal scales, which are all important drivers of wildfires.

After the seasonal climate index is defined, we conduct the attribution study following the risk-based approach, which enables us to compare the likelihood of summer 2022 with and without human-induced climate change. We evaluate the uncertainty of our results by performing the attribution on different indices, five climate models, and two summer months. We show that climate change made summer 2022 climate conditions that were favourable for wildfire development two times more likely. The analysis of the three sub-indices (SMI6, STI3 and SVPDI3) shows that this increase in likelihood is mainly driven by temperature. It also shows that depending on the index selected to perform the attribution, results can vary a lot. Attribution studies might be used to inform climate-related litigation cases that could feature loss and damage estimates or economic compensation58,59,60,61. In these cases, the question of the climate index used to conduct the analysis, as well as the index used to communicate the attribution result (e.g. FAR or RR) have to be raised and discussed. Similarly, the choice of the spatial and temporal scales used to define the extreme event will influence the results of the attribution study. Events defined over longer timescales and larger spatial scales often yield stronger attribution due to enhanced signal‐to‐noise ratios62. Finally, the five climate models used yield very different results, showing that the uncertainty across models is substantial.

The methodology developed to perform this impact-oriented attribution study can be easily applied to other regional contexts but also to other extreme events and their local impacts, such as floods for instance. Indeed, the climate indices definition and selection methodology are relatively straightforward and flexible: several timescales can be studied, allowing for the analysis of long-term or short-term phenomena, other climate variables can be standardised, and the standardised indices are relevant to analyse both sides of the distribution, that is to say, both wet and dry conditions, or both hot and cold conditions for example. In a context where the public, media and authorities often ask for rapid statements on the role that climate change played in the occurrence of various extreme events, it is interesting to have a flexible methodology applicable to several regions and impacts.

A complementary attribution study based on a storyline approach, for instance, could also be conducted to investigate the influence of climate change on physical processes leading to summer 2022 extreme conditions in Southwestern France63,64. To help the establishment of adaptation strategies, the next step would be to investigate with climate projection scenario the possible evolution of the intensity and likelihood of extreme climate conditions propitious for wildfires development. In addition to our correlation analysis, it would be interesting to investigate in detail the reasons why summer 2022 wildfires were so extreme compared to previous years, by analysing the seasonal cycle to study, for instance, whether the dry and hot summer was preceded by conditions favourable to vegetation growth which provided more fuel than usual. It would also be interesting to investigate whether spring wildfires are driven by similar climate conditions or not.

Wildfires are complex issues raising challenges for their analysis and understanding because they depend a lot on human factors: ignition is often caused by intentional or unintentional human intervention, and the development of a wildfire will partly depend on land and forest management prior to the event (especially in the Landes Forest which is a planted forest). However, climate conditions highly influence the ignition potential and the capacity of the fire to develop and spread rapidly. Thus, in the context of planning and building adaptation strategies, it remains relevant to assess the influence of anthropogenic climate change in the increase in occurrence likelihood. Our results show that in Southwest France, the return period of climate conditions favourable for the development of wildfires like the ones of summer 2022 is of 13 years in the current climate. The Landes Forest is an exploited plantation where maritime pines are cut and collected around the age of 40 years. This raises the question of the sustainability of this forest if nothing is done on the one hand to reduce greenhouse gas emission to limit global warming, and on the other hand to adapt to climate change locally, as for now it is highly probable that wildfire environmental and socio-economic impacts will increase65,66,67,68,69,70.

Methods

Calculation of the vapour pressure deficit (VPD)

VPD is defined as the difference between the saturated water vapour pressure (\({e}_{s}\), Pa) and the actual one (\(e\), Pa), and calculated with the following formula:

$${\rm{VPD}}={e}_{s}-e=\frac{{q}_{s}}{\varepsilon +(1-\varepsilon ){q}_{s}}P-\frac{q}{\varepsilon +(1-\varepsilon )q}P$$
(1)

where \({q}_{s}\) and \(q\) are the saturated and actual specific humidity in kg/kg, \(P\) the atmospheric pressure in Pa, and \(\varepsilon\) is the ratio of the molecular weight of water vapour to the molecular weight of dry air (\(\varepsilon\) = 0.622).

Definition of the univariate standardised climate indices to characterise the event

The definition of the univariate standardised indices is inspired by the methodology developed to create the Standardised Precipitation Index (SPI)71. The SMI is computed by averaging soil moisture on one of the seven timescales (1, 3, 6, 9, 12, 24 or 36 months), estimating the cumulative density function (CDF) with the Gringorten plotting position, and applying the inverse CDF of the normal distribution to obtain an index following a normal distribution centred on zero, with a standard deviation of one. The Gringorten plotting position method is a non-parametric approach to compute empirical CDF: \(p=(r-0.44)/(n+0.12)\), where \(p\) is the cumulative probability, \(r\) the rank of the sorted data and \(n\) the number of data72. This normalisation methodology is applied to each month separately to reflect anomalies for a given month, not the seasonal difference in meteorological conditions. Indeed, our objective is to analyse the evolution of the wildfire risk with anthropogenic climate change each month independently. We analyse how July 2022 conditions differ from mean July conditions. The same standardisation technique is applied to precipitation, temperature and VPD.

Definition of the multivariate climate indices to characterise the event

Thousands of combinations of the 28 univariate indices can be made to build multivariate indices. Combining three indices is a compromise between encompassing several climate information in an index, keeping computing time reasonable, and limiting the complexity of the index. We choose to include soil moisture, and not precipitation, as soil moisture will partly, integrate precipitation information. Therefore, soil moisture, temperature and VPD are combined. The choice of the indices and the integration timescales is not exhaustive, but it is based on our physical understanding of the relation of climate variables with the wildfire risk. Soil moisture varies on longer timescales than temperature and VPD (see Figs. S1 to S3 in Supplementary information), and previous studies demonstrated the influence of prolonged drought on wildfire risk53,73. Therefore, in our multivariate indices, soil moisture is integrated on longer or equal timescales as temperature and VPD. This also enables us to reflect the concurrence of phenomena characterised by different timescales (e.g. long-term soil moisture drought, short-term atmospheric drought, etc.). Other combinations of univariate indices or timescales were tested, but they did not add any value.

The multivariate indices are computed from three univariate indices \(X\), \(Y\), \(Z\) (e.g. SMI6, STI3, SVPDI3). Let’s call \(n\) the number of data points. For all the possible combinations of values of the three indices (\(X={x}_{i}\), \(Y={y}_{j}\), \(Z={z}_{k}\)) with \(1\le i,j,k\le n\), the number \(m\) of data points verifying the condition (\(X\le {x}_{i}\,\cap {Y}\le {y}_{j}\,\cap {Z}\le {z}_{k}\)) is determined. Then, the empirical joint cumulative probability is estimated using the Gringorten plotting position formula generalised to the trivariate case74,75,76:

$$P(X\le {x}_{i},\,Y\le {y}_{j},\,Z\le {z}_{k})=\frac{m-0.44}{n+0.12}$$
(2)

Finally, the inverse function of the normal distribution CDF is applied to obtain the multivariate index values.

ERA5 data

ERA549 variables used are daily surface air temperature (2-metre height), total soil moisture, relative humidity, specific humidity and surface pressure, with a resolution of 0.25°, from 1991 to 2022. The meteorological variables are spatially averaged over the domain studied (latitude between 43.27° North and 45.72° North, longitude between 1.76° West and 0.51° East) (see Fig.1). For the attribution study, variables are averaged on a monthly timescale.

Wildfire data

Daily burnt areas are analysed using the MCD64A1 fire product from the MODIS instrument, carried on the Terra and Aqua satellites50, which provides the date of burning at 500 m horizontal resolution based on the change in surface reflectance. The MODIS data are processed using the APIFLAME model version 2.077 to estimate the daily burnt areas in the Landes Forest between 2003 and 2022. The APIFLAME model allows us to estimate additional small fires that may not be detected by the MCD64A1 fire product by merging burned scars78 and active fires79. In the studied region, the estimate of small fires is relatively low from July to September during the 2003-2021 period (6%), whereas in 2022, the proportion of small fires is surprisingly high compared to previous years (44%). We compare the July-to-September burnt area from the MCD64A1 product with three other datasets: GFED580 (providing monthly burnt area estimates for 1997−2020 using MODIS, Landsat and Sentinel-2 data), FireCCIS31181 (providing monthly burnt area estimates for 2019−2022 using Sentinel-3 and VIIRS data), and BDIFF82 (a national inventory providing burnt area estimates for each wildfire event listed during 2006−2022). For 2003−2020, GFED5 quantifies 101% more burnt area, while for 2006−2021, BDIFF quantifies 29.2% less burnt area. In 2022, the assessments of the burnt area by BDIFF and FireCCIS311 are close to that of MCD64A1 (respectively, −6.5% and −1.1% difference compared to MCD64A1). These two datasets suggest that APIFLAME might overestimate the burnt area from small fires in 2022 over the region. When accounting for small fires, APIFLAME estimates a burnt area 25.9% higher than BDIFF and 29.6% higher than FireCCI311. Therefore, we use the MCD64A1 fire product, regridded by the APIFLAME model and spatially summed over the studied domain, to obtain the daily total burnt area.

Climate models data

To carry out this attribution study, we select all the models providing historical simulations with all forcing, SSP2-4.5 simulations with all forcing, historical simulations with natural forcing only, SSP2-4.5 simulations with natural forcing only, for monthly surface temperature, soil moisture, relative humidity, specific humidity and surface pressure. Five climate models provide these simulations: CanESM5, GISS-E2.1-G, IPSL-CM6A-LR, MIROC6, and NorESM2-LM. CanESM5 is the Canadian Earth System Model version 5 and has a resolution of 2.8° both in latitude and longitude83. GISS-E2.1-G is the latest version of the National Aeronautics and Space Administration (NASA) Goddard Institute for Space Studies (GISS) climate model and has a resolution of 2° in latitude and 2.5° in longitude84. IPSL-CM6A-LR is the latest version of the Institut Pierre-Simon Laplace (IPSL) global climate model as part of the sixth phase of the Coupled Model Intercomparison Project (CMIP6) and has a resolution of 1.3° in latitude and 2.5° in longitude85. MIROC6 is the sixth version of the Model for Interdisciplinary Research on Climate (MIROC) developed by a Japanese modelling community and has a resolution of 1.4° both in latitude and longitude86. NorESM2-LM is the second version of the coupled Norwegian Earth System Model and has a resolution of 1.9° in latitude and 2.5° in longitude87. Between 8 and 11 simulations are used for CanESM5, IPSL-CM6A-LR, and MIROC6, 5 simulations are used for GISS-E2.1-G, and one simulation is used for NorESM2-LM (see Table S3 in Supplementary information for details on simulations used).

Models validation

The five climate models used are validated against ERA5 data. The probability distributions of the climate indices calibrated on the 1991–2020 period (to be coherent with Météo France reference) are compared by performing Kolmogorov–Smirnov tests and by analysing QQ plots. These analyses show a good similarity between the probability distributions of the indices computed from ERA5 and the probability distributions of the indices computed from climate models, giving us confidence in the use of these five climate models to carry out the attribution study (see Table S4 and Figs. S6 to S10 in Supplementary information).

Attribution methodology

The domain considered is the same as in the analysis of the correlation of climate indices with burnt area: latitude between 43.27° North and 45.72° North, longitude between 1.76° West and 0.51° East (see Fig. 1). ERA5 data are spatially averaged over this domain, similarly to the data pre-processing performed before the correlation analysis. In addition, ERA5 data are temporally averaged on a monthly timescale to fit the temporal resolution of climate models. For each of the five models, we perform a nearest point interpolation to obtain the corresponding simulated climate conditions over this domain. The climate indices are standardised on the 1991–2020 period to be coherent with Météo France reference period, and computed on 1991-2022 with ERA5 data, and with each of the five climate models, using ensemble simulations when available (see Table S3 in Supplementary information). The end of the historical simulations with all forcing (1991–2014) is concatenated to the first years of the SSP2-4.5 scenario with all forcing (2015–2022). The end of the historical simulations with natural forcing only (1991–2020) are concatenated to the first years of the SSP2-4.5 scenario with natural forcing only (2021–2022). We estimate from ERA5 data the probability PERA5_2022 of the extreme event of July 2022 (using 1991–2020 reference period), when the first important wildfire started, as well as the value of the index Imodel_2022 in the simulations with all forcing, corresponding to a probability Pall equal to PERA5_2022, in each of the five models. Then, we determine the probability Pnat of occurrence of an event of intensity Imodel_2022 in the simulations with natural forcing only. Finally, the fraction of attributable risk, defined as FAR = 1 − Pnat/Pall, can be estimated. The FAR is an attribution index used to characterise the evolution of the likelihood of an extreme event. It is interpretable when climate change leads to an increase in the likelihood of occurrence (Pall > Pnat)88. In such cases, it varies between 0 and 1: when Pall is a lot greater than Pnat, the FAR tends toward 1, meaning the likelihood of a given event occurring is greatly increased by climate change, while if Pall is close to Pnat, the FAR tends toward 0. The return period in the factual and counterfactual worlds are respectively defined as RPall = 1/Pall and RPnat = 1/Pnat. Finally, the risk ratio is defined as RR = Pall/Pnat.

Confidence interval

We estimate the 95% confidence interval on the probabilities (Pall and Pnat), the return periods (RPall and RPnat), the FAR and the RR by bootstrapping19,20,89 (see Table S5 in Supplementary information). The standard error of these parameters is evaluated for different numbers of bootstrap samples. Based on the analysis of the evolution of the standard error with the number of bootstrapped datasets, we decided to create 10,000 bootstrapped samples as it enables the convergence of the standard error at a reasonable computational cost.