Quantifying causal teleconnections to drought and fire risks in Indonesian Borneo

Fires occurring over the peatlands in Indonesian Borneo accompanied by droughts have posed devastating impacts on human health, livelihoods, economy and the natural environment, and their prevention requires comprehensive understanding of climate‐associated risk. Although it is widely known that the droughts are associated with El Niño events, the onset process of El Niño and thus the drought precursors and their possible changes under the future climate are not clearly understood. Here, we use a causal network approach to quantify the strength of teleconnections to droughts at a seasonal timescale shown in observations and climate models. We portray two drivers of June‐July‐August (JJA) droughts identified through literature review and causal analysis, namely Niño 3.4 sea surface temperature (SST) in JJA (El Niño Southern Oscillation [abbreviated as ENSO]) and SST anomaly over the eastern North Pacific to the east of the Hawaiian Islands (abbreviated as Pacific SST) in March‐April‐May (MAM) period. We argue that the droughts are strongly linked to ENSO variability, with drier years corresponding to El Niño conditions. The droughts can be predicted with a lead time of 3 months based on their associations with Pacific SST, with higher SST preceding drier conditions. We find that under the SSP585 scenario, the Coupled Model Intercomparison Phase 6 (CMIP6) multi‐model ensembles show significant increase in both the maximum number of consecutive dry days in the Indonesian Borneo region in JJA (p = 0.006) and its linear association with Pacific SST in MAM (p = 0.001) from year 2061 to 2100 compared with the historical baseline. Some models are showing unrealistic amounts of JJA rainfall and underestimate drought risks in Indonesian Borneo and their teleconnections, owing to the underestimation of ENSO amplitude and overestimation of local convections. Our study strengthens the possibility of early warning triggers of fires and stresses the need for taking enhanced climate risk into consideration when formulating long‐term policies to eliminate fires.

and thus the drought precursors and their possible changes under the future climate are not clearly understood.Here, we use a causal network approach to quantify the strength of teleconnections to droughts at a seasonal timescale shown in observations and climate models.We portray two drivers of June-July-August (JJA) droughts identified through literature review and causal analysis, namely Niño 3.4 sea surface temperature (SST) in JJA (El Niño Southern Oscillation [abbreviated as ENSO]) and SST anomaly over the eastern North Pacific to the east of the Hawaiian Islands (abbreviated as Pacific SST) in March-April-May (MAM) period.We argue that the droughts are strongly linked to ENSO variability, with drier years corresponding to El Niño conditions.The droughts can be predicted with a lead time of 3 months based on their associations with Pacific SST, with higher SST preceding drier conditions.We find that under the SSP585 scenario, the Coupled Model Intercomparison Phase 6 (CMIP6) multi-model ensembles show significant increase in both the maximum number of consecutive dry days in the Indonesian Borneo region in JJA (p = 0.006) and its linear association with Pacific SST in MAM (p = 0.001) from year 2061 to 2100 compared with the historical baseline.Some models are showing unrealistic amounts of JJA rainfall and underestimate drought risks in Indonesian Borneo and their teleconnections, owing to the underestimation of ENSO amplitude and overestimation of local convections.Our study strengthens the possibility of early warning triggers of fires and stresses the need for taking enhanced climate risk into consideration when formulating long-term policies to eliminate fires.

| INTRODUCTION
Indonesia on the island of Borneo in Southeast Asia is home to extensive peatlands of approximately 5.8 Mha (Wahyunto & Suryadiputra, 2008).These areas are susceptible to widespread fires which last for months, commonly from July to October following a lack of rainfall during the dry season in boreal summer (Harrison et al., 2020;Yin et al., 2016).Drought conditions associated with the variability of dry season rainfall sparked a number of large-scale fire events over the last few decades, notably in 1982, 1997, 2006and 2015(Chen et al., 2016;;Najib et al., 2022;Sloan et al., 2017).In 2015, peat fires have caused an estimated 100,000 premature deaths, major economic disruption with a cost of 16.1 billion US dollars to the Indonesian economy and, in 3 months, emitted more carbon than the entire European Union (EU) region (World Bank, 2016).Effective action to prevent the multi-hazard of drought and fire and their resulting impacts requires a more detailed understanding of climate-associated risk, which is underpinned by accurate representation of enabling conditions of fires and the associated climatic regimes.
The variability of dry season rainfall in Indonesian Borneo is linked to teleconnections, defined as anomalies in large-scale ocean-atmosphere circulations, which serve as sources of predictability (Allan et al., 1996;Diaz et al., 2001;Kretschmer et al., 2021;Wang, 2006).The relevant teleconnections of drought for the Maritime Continent region, including El Niño Southern Oscillation (ENSO) and Indian Ocean Dipole (IOD), are well-known and easily quantified by climate indices (Field et al., 2015;Nikonovas et al., 2022;Spessa et al., 2015).Droughts in the Maritime Continent are associated with El Niño as part of the anomalous subsidence arm of the Walker Circulation (Hendon, 2003) and the positive phase of IOD amid the cooling of eastern Indian Ocean (Ashok et al., 2001;Saji et al., 1999).In particular, they tend to co-occur with the onset phase of El Niño (Chen et al., 2017;Susilo et al., 2013), thus their predictability relies on ENSO forecasting (Wooster et al., 2010).The onset process of an ENSO cycle, characterized by rising sea surface temperatures (SSTs) over equatorial Eastern Pacific, is driven by multiple factors (Alexander et al., 2010;Jin, 1997;Solomon et al., 2008;Vecchi & Harrison, 2000) and the associated mechanism remains an active area of research (Chao et al., 2021;Yu & Fang, 2018).The most prominent mechanism is the elevated SST over eastern North Pacific east of Hawaii during boreal spring, which through wind-evaporation-SST feedback, weakens or even reverses the trade winds over equatorial Western Pacific, in turn triggering anomalous subsidence over the Maritime Continent, which contributes to droughts (Alexander et al., 2010;Vimont et al., 2001).There is a lack of research with regards to how these preceding teleconnections affect dry season rainfall, drought and fire risk specifically in Indonesian Borneo.A solid understanding of the physical processes combined with the data of specific observable variables that represent them is a necessary first step from which forecasters may be able to give a good prediction of droughts in Indonesian Borneo with a lead time of several months and enable early actions to prepare against fire onset.
It is important to understand how these teleconnections will change in the future and the subsequent effects on the likelihood of drought occurrence (Shepherd, 2021;Stott, 2015).Whilst uncertainties remain regarding the effect of anthropogenic climate change on the variability of ENSO and IOD (Beobide-Arsuaga et al., 2021; Intergovernmental Panel on Climate Change [IPCC], 2021), there is an emerging consensus that the occurrence of strong El Niño and positive IOD events may increase in the future (Cai et al., 2022;Li et al., 2021;McKenna et al., 2020;Power et al., 2013;Stevenson, 2012).It is even more likely that rainfall variability associated with ENSO will amplify in the future under a warming climate (Chung & Power, 2016;IPCC, 2021;Power et al., 2013).For instance, Coupled Model Intercomparison Phase 6 (CMIP6) models (Eyring et al., 2016) show that ENSOassociated rainfall variability will strengthen over the Central and Western parts of the Maritime Continent region (Chen et al., 2023), implying that dry episodes will become more strongly associated with El Niño.How a warming future climate may modify the linkage between drought and fire risk specifically in Indonesian Borneo and its large-scale climate drivers, particularly the preceding teleconnections, remains to be understood.
As the nature of teleconnections are basically causalities between remote climate features, we explore and quantify them using causal inference and causal network (Kretschmer et al., 2021).First introduced by Ebert-Uphoff and Deng (2012) to the climate science community, causal inference has been a powerful tool to detect the influence of teleconnections on climate and weather events, mostly at seasonal timescale, across the globe.It makes good use of climate time series data to address the subseasonal to seasonal prediction gap (Mariotti et al., 2018).Prominent methods of causal inference in climate science include linear structural equation models combined with hypotheses based on our understanding of the physical system (Kretschmer et al., 2016(Kretschmer et al., , 2021;;Pearl, 2013), Peter and Clark Momentary Conditional Independence (Runge, 2018;Runge, Nowack, et al., 2019) and Response-guided Causal Precursor Detection (Di Capua et al., 2019;Kretschmer et al., 2017).These methods are applied to detect drivers of climatic events with huge societal impacts, including severe winters in North America and Europe (Kretschmer et al., 2016(Kretschmer et al., , 2017)), Indian summer monsoon (Di Capua et al., 2019), crop yield in Morocco (Lehmann et al., 2020) and Atlantic hurricane activity (Pfleiderer et al., 2020).Despite the efficacy in detecting climate drivers, substantial challenges remain mainly due to the complexity of climate systems, autocorrelation among drivers and target variables, nonstationarity, non-linear dependency and noise (Runge, Nowack, et al., 2019).Nevertheless, given that the causal inference approach is transparent and informed by actual data, and could be conveniently interpreted in a geophysical sense, it could serve as a convincing means of evaluating physical climate models and narratives (Runge, Bathiany, et al., 2019).We explore the causal approach to quantify the effects of known observable predictive variables to droughts and fire risks in Indonesian Borneo in the past and future, represented by a number of indicators including number of dry days (denoted as NDD) (Novitasari et al., 2019;Polade et al., 2014) and maximum number of consecutive dry days (denoted as MCDD; Frich et al., 2002, Kurniadi et al., 2021) for droughts, and Canadian Forest Fire Weather Index (denoted as FWI) (Turner & Lawson, 1978), which correlates with actual fire occurrence strongly in various climatological regions across the globe (Dimitrakopoulos et al., 2011;Dymond et al., 2005;Ullah et al., 2013).
In summary, this paper aims to shed a light on the following problems:

| Observational and reanalysis datasets
To represent drought conditions and their teleconnections in the present climate, we used a large range of datasets as inputs.For rainfall and drought indicators, we used Climate Hazards InfraRed Precipitation with Station data (CHIRPS) daily precipitation (Funk et al., 2015) having tested other datasets (Table S2), with the reasons as follows: (1) CHIRPS covers a climatological time range of over 30 years and provides a daily temporal resolution which enables the computation of drought indicators; (2) we believe that CHIRPS can best represent rainfall in our region of study, with the in situ measurement gaps filled by satellite estimates such that the local rainfall maxima over the central parts of Borneo are properly accounted for (Figure S1); (3) CHIRPS is extensively used in rainfall studies over the Asian tropics (Guo et al., 2017;Narulita et al., 2021;Sulugodu & Deka, 2019).For estimating fire risk using Canadian Forest FWI, we used the European Centre for Medium range Weather Forecasts Reanalysis (ERA5) daily reanalysis data (Hersbach et al., 2020) including rainfall, air temperature at 2 m, relative humidity at 1000 hPa as well as zonal wind and meridional wind at 10 m, thanks to the availability of these variables at a high resolution of 0.25 online.We also used the ERA5 data for estimating zonal winds over equatorial Western Pacific, as it has a good track record of capturing ENSO and its related teleconnections as in previous studies despite being a reanalysis product (Capotondi & Ricciardulli, 2021;Nikonovas et al., 2022;Nurdiati et al., 2021).For sea surface temperature over eastern North Pacific, we chose monthly HadISST2.0as it is embedded with satellite measurements (Kennedy et al., 2013).To represent ENSO and IOD, Niño 3.4 SST and Dipole Mode Index were retrieved from the NOAA ESRL Physical Sciences Laboratory website, respectively (Rayner et al., 2003;Saji & Yamagata, 2003).Readers are referred to Table S1 for more details of the datasets including the regions of interest, available years and their sources.Based on the availability of data and the historical runs of climate model experiments (refer to Section 2.3 for details), data from 1981 to 2014 were extracted.

| Coupled models
For the model evaluation and projection of possible future teleconnections, we employed the ensemble of CMIP6 models (Eyring et al., 2016).From these models, we extracted precipitation, near-surface zonal wind, meridional wind, air temperature and relative humidity at daily timescale, and monthly sea surface temperature.We used historical, SSP245 and SSP585 experiments for all models listed in Table 1 (selected according to the

| Determining enabling conditions of fires
To represent the enabling weather conditions for fires in Borneo, we calculated a number of drought and fire risk indicators for our causal analysis.They include seasonal rainfall (in mm/day, denoted as BorneoP), NDD (Novitasari et al., 2019), MCDD (Frich et al., 2002;Kurniadi et al., 2021), and daily Canadian Forest FWI (Turner & Lawson, 1978), in each June-July-August (JJA) period.For indicators involving 'Dry days,' these are defined as daily precipitation <1 mm (Nastos & Zerefos, 2009;Pal & Al-Tabbaa, 2009;Polade et al., 2014).Observed drought indicators were calculated using CHIRPS daily precipitation data at a resolution of 0.25 × 0.25 as input, whereas historical and future projections of drought indicators were computed using CMIP6 daily precipitation data on their native grid (Table 1).
The FWI was developed by Canadian Fire Weather Service based on precipitation, surface air temperature, wind speed and relative humidity.Based on these inputs, the moisture content in organic matter is evaluated, represented by Drought Code, Duff Moisture Code and Fine Fuel Moisture Code.These give immediate indices of FWI, namely Buildup Index and Initial Spread Index, which reflects the potential of fire onset and possible speed of fire spread respectively.Details of the index can be found in its documentation (Turner & Lawson, 1978).As well as evidence of its efficacy from the literature support (Ullah et al., 2013), the capability of FWI to capture actual fire occurrence in the study region is backed by an exploratory analysis, where we used a database of fire occurrences between Year 2001 and 2021 from Fire Information for Resource Management System (FIRMS) (https://firms2.modaps.eosdis.nasa.gov).Figure S2 shows that the accumulated fire radiative power correlates well with FWI in each year with R 2 of 0.751 ( p < 10 −6 ).We employed ERA5 daily reanalysis data at a resolution of 0.25 as inputs to compute the observational estimate of FWI, and the daily CMIP6 outputs mentioned in Section 2.3 to calculate the FWI projections.Due to the slight variation of grid points among the variables in some of the CMIP6 models (ACCESS-CM2, ACCESS-ESM1-5, HadGEM3-GC31-LL and UKESM1-0-LL), variables were regridded with the area-weighted scheme to match the precipitation grids before the FWI was calculated.
We first computed the indicators in each grid cell and generated time series based on the area average within the designated area.For FWI, the values for each day were averaged to obtain a JJA mean for each year.

| Quantifying teleconnections
We attempt to extract causal information from the seasonal time series of various drought and fire risk indicators mentioned above and their known drivers, which are generated from the daily and monthly data collected (refer to Sections 2.2 and 2.3 for details).The time series are standardized so that strengths of associations can be directly compared against each other, and detrended to isolate variability from long-term trends.The standardized rainfall also means that BorneoP resembles the 3-month Standardized Precipitation Index every August (Figure S3), which fulfils the World Meteorological Organization's (WMO) recommendation for drought assessment (WMO, 2012).
We adopt the methodology of knowledge-guided causal inference described in Kretschmer et al. (2021) (abbreviated as Kretschmer method below).Kretschmer et al. (2021) suggests that under a pre-constructed causal hypothesis based on expert knowledge of physical mechanisms underlying the data (Alexander et al., 2010;Hendon, 2003;Saji et al., 1999) and the assumptions of linear dependence and the absence of confounders, linear regression can serve as a representation of causal effect (Pearl, 2009(Pearl, , 2013)).We recognize that correlation does not imply causation, and therefore we highlight the importance of physical knowledge about the processes involved to any conclusion of causality (Kretschmer et al., 2021).In our study, the causal hypotheses elicited in Figures 1 and 2 are validated by linear regression.It follows Kretschmer et al. (2021) who adopted a stepwise regression approach and drew conclusions based on a simple set of causal inference rules.The same method is used for revealing both present and future teleconnections, with observational and reanalysis data (see Table S1) and CMIP6 data (see Table 1) as inputs, respectively.The causal pathways are considered as follows: 1. ENSO and IOD in JJA (Figure 1) Based on our physical knowledge of the impact of ENSO on Borneo drought, a causal diagram (Figure 1) is constructed to showcase a 'collider structure' with two explanatory variables (ENSO and IOD) accounting for one target variable (Borneo drought).We quantify the causal effect from ENSO to Borneo drought (teleconnection 1a) by regressing Borneo drought (JJA, seasonal) on ENSO (JJA, seasonal), with the coefficient and p-value of correlation indicating the strength and significance of the relationship, respectively.The same approach is adopted for estimating the causal effect from IOD to Borneo drought.In mathematical terms, where x ENSO->Drought is the coefficient and c represents noise.
Similarly, to quantify teleconnection 1b, Considering both ENSO and IOD's influences on Borneo drought and the assumption that they propagate through separate causal chains, we also regressed Borneo drought (JJA, seasonal) on both ENSO and IOD (JJA, seasonal).This can be expressed as follows: We would expect that either coefficient x ENSO!Drought (MLR) or x IOD!Drought(MLR) is distorted from the initial relationship (compared with simple linear regression), with a larger distortion indicating a relatively weaker role of influence.This is known as collider bias (Kretschmer et al., 2021).
2. SST and zonal wind patterns over the Pacific in MAM (Figure 2) Based on our understanding of the physical processes involved in the onset of ENSO cycles as mentioned in Section 1 and Alexander et al. (2010), both the boreal spring (March-April-May, MAM) SST over eastern North Pacific east of Hawaii (denoted as SST) and the trade winds over the equatorial Western Pacific (denoted as uwind) are responsible for the variability of Borneo drought in JJA, with uwind being driven by SST, and driving Borneo drought eventually.This is illustrated using another causal diagram (Figure 2), which depicts the atmospheric circulation pattern over the months leading up to the Borneo drought.The causal diagram showcases a 'chain structure,' with two explanatory variables (SST and uwind) forming a chain and uwind acting as a mediator.We quantify the causal effects from SST and uwind to Borneo drought (teleconnections 2a, 2b and 2c) by regressing Borneo drought (JJA, seasonal) separately with SST (MAM, seasonal) and uwind (MAM, seasonal), and regressing uwind on SST, as follows: Given that uwind acts as a mediator in this causal chain, we would expect the product of coefficients of the causal links for stepwise correlations (x SST!uwind and x uwind!Drought ) approximately equals to the correlation coefficient between the two remote targets, SST and Borneo drought (x SST!Drought ).This is known as path-tracing rule for mediators (Kretschmer et al., 2021;Pearl, 2013), which informs the underlying physical process of teleconnections.Alternatively, we also carried out multiple linear regression, where a stack of SST and uwind can be regressed on Borneo drought, as follows: Here, given that uwind serves as a mediator in this causal chain and we are blocking the causal pathway from SST to Borneo drought via uwind, we expect that the coefficient of SST is reduced when uwind is included in the regression (Kretschmer et al., 2021).The difference in coefficients of SST (regressed on Borneo drought) between that determined over simple and multiple regression implies the degree of mediating effect of uwind.
The purpose of using two separate causal networks for quantification (instead of employing all potential drivers into the same model) is to demonstrate the drivers of Borneo droughts originating in different seasons.Quantifying the causal links shown in Figure 1 serves to compare the influences of the commonly known contemporaneous links, whereas those in Figure 2 inform the potential sources of predictability of Borneo droughts at seasonal timescale.
The outputs of the Kretschmer method are the correlation coefficients in Equations ( 1)-( 7).As the time series are standardized, the coefficients serve as β coefficients, which means that a change in the independent variable by a standard deviation leads to a change in the dependent variable by a standard deviation multiplied by the coefficient.Thus, an evaluation of the correlation coefficients for observed datasets gives insight into the strength of drivers of Borneo drought and fire risk.Comparing the correlation coefficients between the observations and the historical model simulation allows us to evaluate the representation of these teleconnections in models.Comparing the historical and future model simulations gives the overall change in the strength of the teleconnections and the range of different outcomes representing uncertainty.Paired t-tests were used to determine the significance of differences between future and historical runs.To eliminate the issue of multiple comparisons involved between the historical and future simulations of CMIP6 model members, the Bonferroni correction was applied to correct the p-values (Armstrong, 2014;Jafari & Ansari-Pour, 2019).We also evaluate whether there is intermodel consensus for future changes of drought and its teleconnections, defined by having more than two-thirds of models showing the same direction of change (Cai et al. 2021).The direction of changes in drought and its associated teleconnections are compared against each other for each model, which allowed us to determine the possible linkage between both.We acknowledge that the R 2 values of regression can also be used for assessing the strength of drivers but decide not to include all of them in this paper given their quadratic associations with their respective β coefficients (Joram, 2021).Having said that, we are reporting the R 2 values for Equation (7) as they can serve to represent the overall strength of drivers originating during MAM.

| Drivers of variability of rainfall, drought and fire weather in Indonesian Borneo region
First, we give an overview of drought conditions in Borneo.Based on CHIRPS precipitation data from 1981 to 2014, the daily mean rainfall in JJA is 6.54 mm/day, with a considerable interannual variability leading to a minimum of 3.57 mm/day in Year 1997 and a maximum of 10.15 mm/day in Year 2010.The average NDD in JJA in Borneo is 38.40, which is approximately 42% of the total number of days, with a range between 19.62 and 53.55 days.The average MCDD is 8.62, with its minimum and maximum being 3.21 and 14.19 days, respectively.As shown in Figure 3, the South of Borneo has on average more than twice the NDD and consecutive dry days than the north.We find some degree of discrepancy among these estimates by different datasets (Figures S1 and S4, Tables S3-S6), possibly due to data gaps (e.g., a sparse network of rain gauges affecting the robustness of rainfall estimation using APHRODITE; McAlpine et al., 2018) and modelling biases (particularly the 'drizzling bias' [Chen et al., 2021] leading to a lower NDD for ERA5).Based on ERA5 reanalysis data from 1981 to 2014, the highest mean FWI occurs over the south of Borneo with a value of around 5. FWI tends to be synchronous with rainfall and drought indicators, that is, maximum correlation between the monthly FWI and the monthly drought indicators occurs with a lag of zero months (Table S3).
In the following, we quantify the contribution of drivers as shown in Figure 1 to Borneo rainfall.Based on CHIRPS precipitation data and Niño 3.4 SST from 1981 to 2014 (referred to as 'observational estimate' below), the association of ENSO with BorneoP in JJA (teleconnection 1a) is strong and statistically significant.Using Equation (1), we fit x ENSO!Drought by regressing BorneoP against ENSO, and obtain a coefficient of −0.77, based on areal average rainfall.This statistically significant correlation is also evident based on other precipitation datasets (Table S8).By carrying out regression for each grid cell, we obtain the spatial distribution of correlation strength with ENSO, which shows slightly stronger correlations over the drier southern parts of Borneo compared to the wetter northern parts, as in Figure S5a.In contrast, the association between IOD and BorneoP (teleconnection 1b) is rather weak, insignificant and inconsistent.Using Equation (2), we fit x IOD!drought by regressing Bor-neoP against IOD and the coefficient is −0.16.By fitting both x ENSO!Drought(MLR) and x IOD!drought(MLR) into Equation (3), the coefficients of ENSO and IOD are −0.72 and +0.27, respectively, indicating the dominant role of ENSO relative to IOD in influencing BorneoP.Using ERA5 reanalysis precipitation data, we also uncover a possible multidecadal variation, with the statistically significant negative relationship between IOD and BorneoP suggested in literature (Saji et al., 1999) actually being present over mid-20th century  whilst gradually diminishing towards the end of 20th century (Table S9).Owing to the weak, insignificant and shifting relations between IOD and BorneoP and some reverse causality suggested in literature (Ashok et al., 2001), the pair is not taken forward for further analysis.
Next, we quantify the contribution of drivers as shown in Figure 2 to Borneo rainfall, based on CHIRPS precipitation, HadISST2.0and ERA5 wind data from 1981 to 2014 (referred to as 'observational estimate' below).We find that the associations between SST and equatorial uwind during MAM (teleconnection 2b), and the associations between uwind in MAM and Bor-neoP in JJA (teleconnection 2c) are strong and significant.Using Equation (6), we fit x SST!uwind by regressing uwind against SST, and the coefficient is 0.565.The positivity of x SST!uwind is because higher SST weakens or reverses the easterly trade winds, and uwind is positive when winds blow from west to east.Using Equation (5), we fit x uwind!BorneoP by regressing BorneoP against uwind and obtain a coefficient of −0.586.The negativity of x uwind!BorneoP is because weaker than usual or reversed easterly trade winds leads to low rainfall anomaly in Indonesian Borneo region.Based on Equation (4), the correlation coefficient between SST in MAM and Bor-neoP in JJA x SST!Drought is −0.260, which roughly equals to the product of x SST!uwind and x uwind!BorneoP , indicating the role of uwind as a mediator (Kretschmer et al., 2021).This is also evidenced by regressing BorneoP against both SST and uwind using Equation ( 7), which brings about the coefficient x SST!Drought(MLR) close to zero (0.103) and x uwind->Drought(MLR) strongly negative (−0.644).The R 2 value, which reflects the overall strength of this teleconnection pathway to BorneoP, is 0.350.By regressing BorneoP to SST for each grid cell, we find that the sign of correlation is negative across the study region (Figure S5b).The correlation coefficients of equations by regressing various drought and fire risk indicators are shown in Table 2. Owing to the interdependence of different target variables representing Borneo drought and fire risk, the values of coefficients are similar.
The same sets of correlation coefficients of equations using other datasets are shown in Tables S10-S14, leading to the same conclusions as the above.

| Evaluation of model simulations of rainfall and drought conditions in Indonesian Borneo region
We first evaluate the historical climatological BorneoP from the CMIP6 models in JJA in Figure 4a day and ensemble mean of 7.04 mm/day, compared with the observed daily mean of 6.54 mm/day) is due to the Institute for Numerical Mathematics (INM) and Institut Pierre Simon Laplace (IPSL) models showing severe wet bias, with means of over 10 mm/day, whereas the low bound is 3.53 mm/day.
We then consider another indicator-MCDD (Figure 4b)-which serves as a representation of drought conditions in Borneo.The observed value is 8.62 days, which is very close to the inter-model mean of 8.09 days.The model range is from 1.21 to 16.72 days.
Other drought indicators namely NDD (Figure S6) and FWI (Figure 4c) are not well represented by CMIP6 models.NDD is underestimated by the models by over 30% on average, and none of the models can capture the class of FWI appropriately.Therefore, we will mainly focus on MCDD as the selected drought indicator when discussing our results of Borneo drought teleconnections in the next section.

| Evaluation of model simulations of teleconnections
Based on our understanding of observed relationships between ENSO and BorneoP (teleconnection 1a), we next evaluate the associations seen in CMIP6 models against those shown in the observational estimate.Figure 5 shows the range of values of x ENSO!Drought for Equation (1) over different time periods.Compared with observations (x ENSO!Drought = −0.75),models on average underestimate the ENSO-BorneoP relationship (with a median of −0.496 and a range between −0.814 and +0.662).Most of the models underestimate or fail to capture this link.It is worth noting that some of the models which largely overestimate BorneoP (notably IPSL and INM models) also underestimate ENSO teleconnections with correlation coefficients much less negative than the observed value.
Next, we look at the teleconnection from the preceding season to BorneoP in CMIP6 models.Figures 6-8, S7 and S8 show the range of values of x SST!Drought , x uwind!Drought , x SST!uwind , x SST!Drought(MLR) and x uwind!Drought(MLR) , respectively, for Equations ( 4)-( 7), over different time periods.We find that models on average represent relationship between SST and BorneoP fairly well (with the median of x SST!Drought = −0.11)compared with the observational benchmark (x SST!Drought = −0.26),albeit slight underrepresentation by the majority of members, and substantial spread with a range of −0.629 and +0.279 (Figure 6).It is worth noting that members which correctly represent ENSO teleconnections with BorneoP tend to display a stronger relationship between SST and BorneoP, and vice versa.When the mediator uwind is considered, models underestimate the relationships compared with observations, as shown for the relationship between SST and uwind (median of x SST!uwind = 0.25 vs. 0.56 observed) (Figure 7) and that between uwind and BorneoP (median of x uwind!Drought = −0.40 vs. −0.59observed) (Figure 8).This may indicate some shortcomings of the representation of the wind fields over a large proportion of model members, which hinders proper estimation of the teleconnections.Because of this, the effect of uwind shown in models is generally weaker than that in observations, as reflected in multiple linear regression by a more negative coefficient of SST (median of x SST!Drought(MLR) = −0.05 vs. 0.1 observed) (Figure S7), and less negative coefficient of uwind (median of x uwind!Drought(MLR) = −0.4 vs. −0.64observed) (Figure S8).Models that correctly represent ENSO teleconnections tend to simulate a slightly stronger relationship between uwind and BorneoP, and vice versa.As per the R 2 values which reflect the overall strength of this teleconnection pathway to BorneoP, most of the models show a weaker relationship than observed (median of 0.21 vs. 0.35 observed) and the model spread is large with a range of 0.018 and 0.659.
We then look at the teleconnections for other drought indicators.We selected MCDD as the main indicator for analysis of drought teleconnections as the models generally represent it better compared with other indicators.Figure 9 shows the range of values of x ENSO!Drought for Equation (1) over different time periods by regressing MCDD over ENSO (with a median of 0.50 and a range between −0.56 and 0.811).Similar to what is seen for ENSO-BorneoP relationship, most models underrepresent the association between ENSO and Borneo drought.and ( 7), over different time periods, by regressing MCDD over drivers from the Pacific.Models as a whole pick up the signal between SSTs and MCDD fairly accurately, with a median of 0.174 (compared with 0.209 observed) and a spread between −0.223 and 0.441.However, the models tend to underestimate the relationships between uwind and MCDD (Figures S9-S11).The overall strength of the teleconnection pathway, as reflected by the R 2 , is well-represented by the models with a median of 0.187 (compared with 0.249 observed) and a spread between 0.009 and 0.738.
Model spread over JJA drought in Borneo and its associated ENSO teleconnections are pronounced in the historical runs.Some CMIP6 models are unable to realistically simulate the observed relationships whilst underestimating drought conditions to a large degree.Some sources of uncertainties such as ENSO amplitude (Beobide-Arsuaga et al., 2021) are well-documented in the literature, and they correspond to most of the models that have difficulties in simulating the ENSO-BorneoP relationship (INM-CM4-8, INM-CM5-0 and ACCESS-ESM1-5).However, some models with high ENSO amplitudes also show a lack of relationship between ENSO and Borneo drought, the most prominent one being IPSL-CM6A-LR, indicating some other sources of biases in the representation of teleconnections.
We conducted an exploratory analysis of the possible sources of model spread.We picked the strong El Niño years and created a composite of atmospheric circulation patterns in observations and the chosen CMIP6 models-IPSL-CM6A-LR as the model of interest and MIROC6 as the control being able to realistically simulate Borneo drought and its drivers.To identify El Niño events, we first calculated Niño 3.4 SST based on the SST outputs from the two chosen models followed by linear detrending (Kumar et al., 2016) to remove the warming trend of equatorial Eastern Pacific SSTs seen over recent years in both observation and CMIP6 historical runs.We then standardized the Niño 3.4 SST time series and select the years when Niño 3.4 SST is above 1 SD as a representation of strong El Niño events (King et al., 2019).The precipitation patterns during the simulated strong El Niño events were compared against those during the observed strong El Niño events.
We find that both models are able to simulate the large-scale anomalous dryness over the Maritime Continent during strong El Niño.However, IPSL-CM6A-LR simulates a local maximum rainfall anomaly in Borneo even during strong El Niño events (Figure 11,left).By checking the vertical velocity simulations (not shown), we find that this is a result of unrealistic simulation of strong local convection characterized by upward air motion, which overrides the large-scale teleconnection signal.This may be due to the model's reliance on novel convective parameterisation techniques (Couvreux et al., 2021).

| Model simulations of rainfall and drought conditions in Indonesian Borneo region and their teleconnections towards the future
There is no clear trend for mean rainfall towards the future, with none of the future experiments showing significant difference from the historical run at 95% confidence interval (CI) using t-test (paired between each member) (Figure 4a).There exists large and growing model spread towards the far future especially for SSP585 scenario.Whilst mean and median rainfall change by less than 10% (7.04 mm/day for the median and 7.31 mm/day for the ensemble mean), changes in individual model members cover up to plus or minus 30%, with a minimum of 3.89 mm/day and a maximum of 12.78 mm/day.It is worth noting that models with wet biases tend to get wetter, whilst those with dry biases show varied changes.For most models, the changes are more pronounced for high emission scenario.
Under the SSP585 scenario, Borneo drought (represented by MCDD) increases significantly (p = 0.006) in the far future by 4.48 days or 56% compared with historical baseline (Figure 4b), supported by inter-model consensus.There is inter-model consensus that FWI (Figure 4c) will increase towards the far future in both SSP245 and SSP585 scenarios, yet such increases are not statistically significant.Two models (HadGEM-GC31-LL and UKESM-1-0-LL) indicate large increases of FWI towards far for SSP585 scenario (Figure 4c).
The majority of models agree to a slight yet statistically insignificant strengthening of association between ENSO and Borneo drought, without inter-model consensus (Figure 9).As can be seen in Figure 12, models that project a strong increase in drought (indicated by MCDD) tend to produce a strengthening teleconnection from ENSO to MCDD, and vice versa.A statistically significant linear correlation ( p = 0.017) is shown by regressing the change of correlation coefficients between ENSO and MCDD (SSP585 far future minus baseline) against the percentage change of MCDD for each model.This provides strong evidence of the dependence of future drought changes on the ENSO signal, even though the models do not agree with each other over the future change of teleconnection strength.It is also noted that the direction of change of FWI is less aligned with its teleconnections (not shown), possibly implying a stronger influence of local drivers on the enabling conditions of fires.
Most model members show strengthening of SST-BorneoP and SST-uwind relationships in the future (Figures 6 and 7), and such changes are larger than changes to the ENSO-BorneoP relationship.In particular, for the SSP245 simulations, the strengthened relationship between SST and BorneoP over both near and far future compared with the baseline is statistically significant at 95% CI (p = 0.014) (Figure 6).We find that the changes of relationships between SST and uwind, and between uwind and BorneoP are not statistically significant (Figures 7 and 8).Based on Equation ( 7), most model members show slight weakening of the mediating effect of uwind towards the future, owing to strengthening of role of SST on BorneoP, yet such change is not statistically significant (Figure S11).The overall strength of the teleconnection pathway SST !uwind !BorneoP shows a sign of increase, as evidenced by an increase in multi-model median R 2 over both near and far future for both SSP245 and SSP585 scenarios (0.230-0.311, vs. 0.210 for historical simulations).There is multi-model consensus for such an increase in the SSP245 scenario over the near future.
Most model members show strengthening of SSTdrought relationships.Under both SSP245 and SSP585 scenario, the relationship between SST and MCDD strengthens significantly in the far future compared with historical baseline, with the mean of coefficient x SST!Drought increased by 0.156 and 0.161 for SSP245 and SSP585 scenarios respectively (p equals 0.043 and 0.001, respectively for the two scenarios), supported by inter-model consensus (Figure 10).There also exists inter-model consensus over the strengthening of this teleconnection for SSP585 in the near future, although the change is statistically insignificant.Based on Equation ( 7), a strengthening of the teleconnection from SST to MCDD through uwind is demonstrated by an increase in multi-model median R 2 (0.269-0.317 for the two respective scenarios and time periods, vs. 0.187 for historical F I G U R E 1 2 Δx ENSO!Drought versus ΔMCDD for SSP585 compared with historical baseline.simulations).This is supported by multi-model consensus for the SSP585 scenario.

| Discussion
This study contributes to ongoing efforts to identify precursors of fire risk in Indonesian Borneo (Field et al., 2015;Nikonovas et al., 2022;Spessa et al., 2015).Previous studies have put a particular focus on contributors to fire risk at provincial and local scale including land use and other non-climatic factors (Field et al., 2015;Harrison et al., 2020).We complement these studies by considering a large geographical area on the island of Borneo and demonstrate a strong and spatially coherent correlation between drought and fire risk indicators and ENSO across the region (Figure S6).We also find strong associations between these indicators and anomalies of SST and zonal winds over the Pacific in MAM.Over the short term, this enables data-informed regional-scale early warning triggers for the communities, as well as preventive measures such as imposing temporary bans of forest clearance and peatland drainage, when a drought is set to occur amid the emergence of its drivers with a lead time up to 3 months.We have conducted our study using state-of-the-art climate models under different emission scenarios, which facilitates consideration of enhanced climate risk when formulating policy reforms to eliminate human-induced fires in the long run.
We did not preprocess (apart from the necessary regridding for calculation of FWI) or downscale the CMIP6 model outputs as the aim of this study is to provide a general picture on the direction of change in drought over a relatively large area in Borneo and its associated large-scale climate drivers.As a pioneer of using CMIP6 data to analyse teleconnections with a causal approach, we also want to evaluate the capability and identify the issues of using raw data.We find that the capability of drought episodes to be represented in model outputs is sensitive to the metric used, with some (MCDD) being capable overall whilst others (such as FWI) are being considerably underestimated partly due to the 'drizzling bias' in climate model simulations (Chen et al., 2021).This could be addressed by preprocessing the rainfall inputs using suitable bias-correction techniques (Maraun, 2013).
Our study also contributes to the pool of literature on the large-scale climate drivers influencing the occurrence of hydrometeorological hazards in the Southeast Asia (Howard et al., 2022).We learn that compared with IOD, which is also a drought driver in the Maritime Continent (Ashok et al., 2001), ENSO exerts a much greater influence over Borneo, which paves the way for future teleconnection studies in the region to focus on drivers from the Pacific.In view of the occurrence of fires in boreal summer prior to the peak of El Niño, we propose an innovative attempt to traverse the spring predictability barrier of ENSO (Lai et al., 2018) where the predictable window spans.With the limited skill of fire risk prediction by state-of-the-art seasonal forecasting models being documented (Nikonovas et al., 2022, Spessa et al., 2015), we adopt a process-informed approach by focusing on particular teleconnections suggested by previous literature (Alexander et al., 2010;Vimont et al., 2001), and verify their presence using the available data and through a causal methodology.This is aided by the path-tracing rule suggested by Kretschmer et al. (2021), which serves to identify a robust signal by 'bridging' through the influence of zonal winds despite the weak direct association between SST and Borneo drought.Whilst we demonstrate that CMIP6 models are showing their skills overall in capturing this association, it would be interesting to explore whether the shortfall of seasonal forecasting models is caused by the underrepresentation of the causal pathways we identified.
We are aware that some known drivers of El Niño onset, including the oscillatory modification of the equatorial Pacific warm water volume associated with recharge and discharge mechanisms (Jin, 1997), trade wind-induced variations in evaporation associated with the Meridional Mode (Alexander et al., 2010), Rossby wave dynamics involving the initiation and westward propagation of subsurface ocean anomalies in Northwest Pacific (20 N-30 N, west of 170 E) (Solomon et al., 2008), and westerly wind bursts (Vecchi & Harrison, 2000), are not accounted for in this study.This is due to the known methodological challenges of the causal approach which makes it less suitable for datasets with fine temporal resolution where autocorrelation may arise (Runge, Bathiany, et al., 2019).The continued development of causal discovery algorithms, particularly when ensemble data inputs are enabled (Gerhardus & Runge, 2020, 2022;Runge, 2020), will hopefully serve to alleviate the issue.
Our study also sheds a light on the possible ENSO teleconnection changes under a warming future climate.By analysing the CMIP6 outputs, we find that the direction of change of drought occurrence depends on whether the ENSO influence strengthens, yet the overall strengthening of association suggested in literature (Chen et al., 2023) lacks multi-model consensus for the case of Indonesian Borneo region.Our unanticipated finding about the insignificant change of ENSO teleconnection yet a significant Pacific SST teleconnection raises questions of (1) changes of mechanism of ENSO onset under a warming climate (which is beyond the scope of this study but is relevant to future fire risks in the study region for long-term resilience building), and (2) quantification of ENSO under a warming climate.We note the limitation to represent ENSO solely using Niño 3.4 SST under a warming climate, as recent studies suggest that there might be better indicators such as rainfall over equatorial Central-Eastern Pacific (Cai et al. 2021;Ying et al., 2022), with high rainfall anomaly in the region under El Niño and vice versa.We conducted the analysis based on these matrices as a proxy of ENSO (including overall ENSO rainfall [8 S-8 N, 150 E-80 W], CP ENSO rainfall [8 S-3 N, 150 E-120 W] and EP ENSO rainfall [1 S-9 N, 140 W-80 W]) and came to the same conclusion, and therefore we focused this analysis on Niño 3.4 SST for simplicity.We appreciate the challenge to find a unified indicator and area to represent the diversity of ENSO, particularly in JJA due to the seasonal cycle of ENSO, which normally peaks in Austral summer.Amid the advancement in understanding the possible changes of ENSO and its teleconnections under a warming climate (Cai et al., 2021;Chung & Power, 2016;Power et al., 2013;Ying et al., 2022;Yun et al., 2021), a further step to develop an integrated multivariate climatic index to characterize ENSO state and strength would be desirable.This will open the door for ENSO teleconnection studies under global warming to uncover the possible changes of a wide range of atmospheric processes and high-impact weather events across the globe.

| Conclusions/main findings
Based on our study to extract causal information from various observational, reanalysis and CMIP6 model datasets, we argue that droughts occurring during JJA in Borneo are strongly linked to ENSO variability, with drought occurrence corresponding to El Niño conditions.Drought in Borneo could be predicted with a lead time of 3 months based on its association with SSTs over eastern North Pacific and zonal winds over equatorial Western Pacific in boreal spring.Through analysing the CMIP6 model outputs, we find that both Borneo drought in JJA and its association with SSTs in the eastern North Pacific in MAM may strengthen towards the end of the 21st century, although various sources of uncertainty remain as evidenced by the large spread between the models.This indicates that the worsening drought could be attributed to the strengthening of wind-evaporation-SST feedback in boreal spring.Forecasters may harness the enhanced relationships to improve seasonal drought prediction for the benefits of fire prevention.Specific findings supporting the main conclusion are as follows:

F
I G U R E 1 Causal hypothesis of synchronous Borneo drought drivers, with denoted teleconnections 1a and 1b.IOD, Indian Ocean Dipole.F I G U R E 2 Causal hypothesis of Borneo drought drivers with a lead time of 3 months, with denoted teleconnections 2a-2c.

F
I G U R E 3 Spatial distribution of various drought indicators over Indonesian Borneo.T A B L E 2 Causal effect of teleconnections on Borneo drought represented by different target variables.
, against the CHIRPS observational estimate.The slight wet bias of CMIP6 models on average (based on median of 7.00 mm/ F I G U R E 4 (a) Borneo Daily Mean Rainfall (BorneoP); (b) Borneo maximum number of consecutive dry days (MCDD); Borneo Canadian Forest Fire Weather Index (FWI) in JJA for Coupled Model Intercomparison Phase 6 experiments (historical: 1981-2014; near future: 2021-2060; far future: 2061-2100), compared with observational estimate (shown as blue line).

F
I G U R E 5 Coefficient of correlation between BorneoP and El Niño Southern Oscillation (ENSO) (x ENSO!Drought ) for Coupled Model Intercomparison Phase 6 experiments, compared with observational estimate (shown as blue line).

F
I G U R E 8 Coefficient of correlation between BorneoP in JJA and uwind in MAM (x uwind!Drought ) for Coupled Model Intercomparison Phase 6 experiments (historical: 1981-2014; near future: 2021-2060; far future: 2061-2100), compared with observational estimate (shown as blue line).F I G U R E 9 Coefficient of correlation between Borneo maximum number of consecutive dry days (MCDD) and El Niño Southern Oscillation (ENSO) (x ENSO!Drought ) for Coupled Model Intercomparison Phase 6 experiments, compared with observational estimate (shown as blue line).

F
I G U R E 1 1 Composite of rainfall anomaly of strong El Niño years from 1981 to 2014, from left to right: (IPSL-CM6A-LR simulation (left), MIROC6 simulation (middle) and CHIRPS observational estimate (right).

F
I G U R E 1 0 Coefficient of correlation between maximum number of consecutive dry days (MCDD) in JJA and SST in MAM (x SST!Drought ) for Coupled Model Intercomparison Phase 6 experiments (historical: 1981-2014; near future: 2021-2060; far future: 2061-2100), compared with observational estimate (shown as blue line).