Climatology of diablo winds in Northern California and their relationships with large-scale climate variabilities

The most destructive wildfires recently in Northern California have been linked to the occurrences of Diablo Winds (DWs). This study investigates the climatology of DWs during September–December 1979–2018, and their relationships with various climate modes using observations and two high-resolution reanalysis datasets. Our finding shows that DWs do not have a long-term trend in terms of the annual total number, total duration, and associated maximum wind speeds of DWs over the past 4 decades. However, their associated minimum relative humidity (RH) has decreased significantly, especially in October, which suggests that the dryness during DWs has become more severe with time, possibly leading to an increased chance of fires, and their destructive potential. We also find that the annual total number and duration of DWs have exhibited an quasi-periodic variation, with intervals ranging from 2 to 4 years. The periodic variability of DWs might be attributed to the teleconnections between DWs and climate oscillations, specifically the El Niño–Southern Oscillation (ENSO) and the quasi-biennial oscillation (QBO), through their modulation of pressure systems near California and the location of the Pacific jet stream. It is suggested that when La Niña and the QBO westerly phases co-occur in the spring, DWs in the following fall and winter tend to occur more frequently, and are associated with more intense high winds and dryness as compared to other springtime QBO and ENSO phases. This relationship may be used to predict the seasonal outlook of DWs.


Introduction
California with its Mediterranean climate is a fire-prone state. As the population grows and wildland and urban interface areas continue to increase (Williams et al. 2019), wildfires are becoming an ongoing threat to many of these communities as they adversely affect regional air quality and lead to loss of life and property. Wildfires are defined here as unpredicted and uncontrolled fires, which are different from prescribed fires that are planned and/or managed to reduce excess fuel buildup in forested areas. Among the top 20 most destructive wildfires in California history, six of them-the Camp Fire (2018), the Tubbs Fire (2017), the Tunnel Fire (1991), the Nuns Fire (2017), the Atlas Fire (2017), and the Redwood Valley Fire (2017)-have been linked to the occurrences of Diablo Winds (DWs). DWs play two important roles in wildfire activities: (1) causing ignition by pulling down power lines or poles, and (2) rapidly spreading fires and increasing the fire size (Mass and Ovens 2019;Monteverdi 1973;Pagni 1993;Smith et al. 2018). As a precautionary measure, several DW events in 2019 prompted the Pacific Gas and Electric Company to shut-off power in several areas of the San Francisco Bay area (SFBA), affecting the daily lives of many residents.
The DWs, which are a type of downslope wind, are characterized as warm, dry, and gusty northerly to northeasterly winds. They occur most frequently during the fall to early winter from the Northern Coast Ranges into the SFBA, and down the western slopes of the Sierra Nevada Mountain Range into the Great Central Valley . Occurrences of DWs are often associated with near-surface ridges that settle either along the coast of Oregon and Washington and/or over the Great Basin (Abatzoglou et al. 2013;Brewer and Clements 2020;Hughes and Hall 2010;Mass Electronic supplementary material The online version of this article (https ://doi.org/10.1007/s0038 2-020-05535 -5) contains supplementary material, which is available to authorized users. and Ovens 2019), which feature a stable low-level atmosphere. This synoptic pattern, along with a calm and stably stratified blocking low-level flow upstream of the mountain ridges, creates strong downslope winds near the surface on the lee side of the mountains in Northern California. The downslope winds descend adiabatically and isentropically, advecting warm and dry air from above the mountain tops downward toward the surface (Chen and Sun 2001;Werth et al. 2011). Warming and drying effects can be enhanced further by turbulent mixing under strong vertical wind shear over rough terrain (Elvidge and Renfrew 2016). In addition, when the air flows over the top of the Northern Coast Ranges or the Sierra Nevada mountains in a stably stratified environment, gravity waves can bring high momentum (energy) aloft down to the surface, also resulting in gusty downslope surface winds (Abatzoglou et al. 2013;Durran 1990). This combination of strong pressure gradients with dry and warm air forms the so-called DWs over Northern California.
As mentioned above, the development of DWs critically depends on large-scale weather patterns, especially the position and strength of low-level high-pressure systems. The presence of this particular synoptic pressure pattern over the local terrain produces strong offshore winds with low relative humidity (RH) values, leading to DW development. It is well-known that these synoptic-scale systems are modulated by atmospheric and oceanic modes that occur over the Pacific Ocean (Wang et al. , 2017. For example, the El Niño-Southern Oscillation (ENSO) has a significant impact on California weather and climate through its teleconnection with the Pacific jet stream and low-level pressure systems near California (Yang et al. 2018).
ENSO is an interannual phenomenon involving ocean-atmosphere interactions associated with changes in the Walker circulation and sea surface temperature (SST) over the tropical Pacific Ocean (Wang 2002). During El Niño phases, SST in the eastern Pacific Ocean becomes warmer than normal, accompanied by a weakening or reversal of the Walker circulation (i.e., weakened easterly trade winds, sometimes even turning easterly winds into westerlies). In contrast, during La Niña phases, SST in the eastern Pacific Ocean becomes colder than normal, accompanied by a strengthened Walker circulation (strengthened easterly trade winds).
The warmer-than-normal SST over the eastern Pacific Ocean during El Niño phases can result in an eastward shift of tropical convection associated with the upward branch of the Walker circulation, as well as a stronger and deeper Hadley circulation over the eastern Pacific, compared to during La Niña phases. This modulation of the Hadley circulation and relocation of tropical convection can cause at least two impacts over California. The first impact is a stronger Aleutian Low during El Niño phases compared to during La Niña phases. The second impact is a strengthened and southward-shifted Pacific jet stream during El Niño phases compared to during La Niña phases. Both impacts on California can lead to abnormally wet and cold winters during El Niño years, and abnormally dry and warm winters during La Niña years (Gershunov and Barnett 1998;Hoell et al. 2016;Liu et al. 2018;McCabe and Dettinger 1999).
The influence of ENSO on California weather and climate is one example of teleconnections. The underlying mechanism can be considered as a low-frequency Rossby wavetrain induced by ENSO (Su et al. 2005), propagating from the equator to extratropical regions over the North Pacific and North America, affecting weather and climate in these regions. The atmospheric variables (e.g., mean tropospheric temperature and geopotential height) tend to respond to the changes in SST associated with ENSO with a lag of 1-2 seasons on an interannual time scale (Su et al. 2005). In addition, this lagged relationship between ENSO-related indices and atmospheric variables (e.g. precipitation) was also observed over California (Liu et al. 2018;Mamalakis et al. 2018). Mason et al. (2017) examined wildland fire potential in the continental United States from 1979 to 2015, and found that in early summer, reduced (increased) wildland fire potential in the southwest tends to occur during El Niño (La Niña). In contrast, an opposite trend is observed for late summer. Inevitably, ENSO characteristics and its teleconnections with California's weather and climate have been shown to change under climate warming based on recent numerical modeling studies (Cai et al. 2014(Cai et al. , 2015aKao and Yu 2009;Wang et al. 2014Wang et al. , 2017Yoon et al. 2015a, b;Yu and Zou 2013;Yu et al. 2012;Zhou et al. 2014). Studies found that under climate warming, the wetting effect of El Niño on the California winter season weakens. In addition, increased water cycle extremes in California may occur in the future via increasingly extreme El Niño and La Niña events if climate warming continues (Yoon et al. 2015a, b;Yu and Zou 2013).
The quasi-biennial oscillation (QBO) could be another large-scale mode that is linked to the climate variabilities over California. While there has not been any study specifically focused on the impact of the QBO on California climate, the strong influence of the QBO on the Pacific jet stream location hints at a potential linkage between the QBO and California climate variabilities because of the strong connection between the Pacific jet stream and precipitation in California (Garfinkel and Hartmann 2011;Hansen et al. 2016;Wang et al. 2018).
The QBO is a stratospheric phenomenon and is characterized by alternating westerly and easterly wind regimes in the equatorial stratosphere with an average period of 28 months. Previous studies found that the QBO can influence weather and climate in the troposphere via at least two ways (Hansen et al. 2016). First, the QBO can directly modify temperature, convection, and vertical wind shear along the tropopause in the tropics, and then impose the consequent changes to the extratropics via the modification of atmospheric circulations, such as the Hadley circulation. For example, through the modification of convection and vertical wind shear along the tropical tropopause, stronger tropical deep convection, a stronger and westward-shifted Hadley circulation, and a stronger Walker circulation are observed during the QBO easterly phase (QBOE) compared to the QBO westerly (QBOW) in the North Pacific region (Garfinkel and Hartmann 2011;Huang et al. 2012;Liess and Geller 2012). These changes during QBOE can result in a strengthened and southward-shifted Pacific jet stream as well as a weakened Aleutian low, compared to those during QBOW (Wang et al. 2018). Second, the QBO can indirectly modify the properties of the stratospheric polar vortex by altering the background flow and the propagation of extratropical planetary waves in the polar stratosphere (Hansen et al. 2016). The stratospheric polar vortex is weaker and warmer during QBOE than during QBOW. The polar-vortexinduced changes can propagate down to the troposphere, affecting the surface weather and climate.
The QBO also shows a strong connection with ENSO, especially in the North Pacific Ocean (Sun et al. 2019;Wei et al. 2007). Hansen et al. (2016) found that the stratospheric equatorial QBO signals extend down to the troposphere over the North Pacific Ocean during the Northern Hemisphere winters only during La Niña but not during El Niño events. Moreover, the Aleutian low is deepened during QBOW (positive phases) as compared to QBOE (negative phases), and the Pacific jet stream is shifted northward during La Niña. However, there has still not been a consensus on how QBO has changed under climate warming. Some previous studies have suggested that QBO has a longer cycle and weaker amplitude as the climate warms (Kawatani and Hamilton 2013;Kawatani et al. 2012;Rind et al. 2014;Watanabe and Kawatani 2012), while other studies show a decreased or unchanged QBO periodicity and an increased QBO amplitude under a warmer climate (Giorgetta and Doege 2005;Schirber 2015).
Despite DWs' strong linkage to ignition and spread of destructive wildfires in Northern California, there have been very limited studies dedicated to the DWs, especially when compared to the Santa Ana Winds (SAWs) in Southern California (Abatzoglou et al. 2013;Guzman-Morales and Gershunov 2019;Guzman-Morales et al. 2016;Hughes and Hall 2010;Jones et al. 2010;Li et al. 2016;Monteverdi 1973;Rolinski et al. 2019). A long-term climatology study on DWs is still lacking and their relations with large-scale climate variabilities have not been elucidated. With these motivations, this study aims to (1) examine the DW climatological trends, (2) explore the DW's relations with the large-scale climate variabilities and the underlying mechanisms, and (3) investigate the predictability of DWs from a climatological perspective.

Datasets
In contrast to the previous investigations of DWs that were either case studies using numerical modeling or short-term climate studies using station data (i.e., remote automated weather stations, RAWS), this study focuses on the long-term climatology of DWs for the period of September-December 1979-2018, and the relationships between DWs and large-scale climate variabilities. In this study, in addition to RAWS data, we use two high-resolution (~ 30-32 km) reanalysis datasets: the National Centers for Environmental Prediction (NCEP) North American Regional Reanalysis (NARR) data, and the fifth-generation of European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis data (ERA5).
The RAWS data are hourly surface data acquired from the RAWS USA Climate Archive (https ://raws.dri.edu/). The winds are measured at 6.1 m above ground level (AGL) and the RH is measured between 1.2-2.5 m AGL. Wind and RH data are 10-min averages prior to the hourly record time (Brewer and Clements 2020).
NARR data are 3-hourly data with approximately 0.3° (32 km) grid spacing, and are obtained from the Earth System Research Laboratory (ESRL) Physical Sciences Division (PSD; https ://www.esrl.noaa.gov/psd/data/gridd ed/data.narr.html). They are created using the high-resolution Eta model, the Regional Data Assimilation System, and global observations. NARR data provide estimations of large-scale atmospheric conditions such as pressure and synoptic-scale winds, with the best performance in winter (Mesinger et al. 2006). Several studies have used NARR data to describe climatological characteristics of SAW events (Hughes and Hall 2010;Jones et al. 2010).
ERA5 data are hourly data with approximately 0.25° (~ 30 km) grid spacing, and are obtained from the Research Data Archive at the National Center for Atmospheric Research, Computational and Information Systems Laboratory (https ://doi.org/10.5065/BH6N-5N20). ERA5 incorporates a vast amount of historical observations into global estimates using advanced modeling and data assimilation systems. Several studies validated ERA5 with remote sensing and in-situ observations and found that ERA5 leads to a consistent improvement over its previous generation, the EC-interim reanalysis (Delhasse et al. 2020;Olauson 2018;Urraca et al. 2018).

Study Area
DWs often occur in and affect two areas of Northern California, namely the SFBA and the western slopes of the Sierra Nevada mountains (Brewer and Clements 2020;Smith et al. 2018). The Tunnel Fire, the Tubbs Fire, the Nuns Fire, the Atlas Fire, and the Redwood Valley Fire each occurred in SFBA, while the Camp Fire occurred in the western Sierra Nevada Mountains. Due to the vast topographic differences between SFBA and the western Sierra Nevada Mountains, this study only focuses on the DWs that affect SFBA. A future study will focus on the DWs over the western Sierra Nevada Mountains.

Definition of DWs
We identify DW events in a similar fashion as did Smith et al. (2018), but with adjustments to incorporate the reanalysis data that will be used for analysis. The DW events are identified independently for each dataset. We define a DW event when data from any of the chosen reanalysis grids or observation stations, which will be introduced in the next section, meet the following five criteria: (1) the surface wind direction is northerly to northeasterly (315°-90°), (2) surface RHs are lower than 30%, (3) surface wind speeds are stronger than 8 m s −1 for RAWS, 7 m s −1 for NARR, and 4 m s −1 for ERA5, (4) the first three conditions are met persistently for six or more consecutive hours, and (5) each event is at least 24 h apart. If two events are within 24 h of each other, they are considered as one event.

Reanalysis data validation
To ensure the capability of NARR and ERA5 data to resolve the features of DWs, we first compare NARR and ERA5 with RAWS data. Because there are strong inter-station variabilities in identifying DW-like events among the RAWSs (Smith et al. 2018), we select four stations (Table 1), which are located on the lee side of the Northern Coast Range and are close to the northern border of the SFBA, in order to focus on DW events that have an immediate impact on SFBA (Fig. 1). We choose three DW events for our reanalysis data validation: the Tunnel Fire in 1991, the Tubbs Fire in 2017, and the Camp Fire in 2018. Gridded values of wind speed, wind direction, and RH from NARR and ERA5 are compared against the measurements from the four RAWS stations (relative locations shown in Fig. 2). We note that there is a difference in the representativeness of observations and reanalysis data, in which the former are point values measured at 6.1 m AGL, while the latter are grid-cell averaged values at 10 m AGL for winds and 2 m AGL for RH. Figure 3 shows the time series of wind speed, wind direction, and RH of the DW event during the Tubbs Fire in October 2017. The time series plots cover at least one day before and after the fire event to show the performance of ERA5 and NARR through the entire period of the DW event. Three of the RAWS stations (caCKNO, caCHAW, and caCSRS) (left column of Fig. 3) consistently show a short period of strong northerly winds (blue lines), with rapid onset on October 8 and rapid decline on October 9. In contrast, the winds at the caCATL station, the most southeasterly located of the four RAWS, show a diurnal-like pattern with a decrease in wind magnitude during the day and an increase in wind magnitude overnight. At all four stations, the RH value rapidly decreases at the beginning of the DW event and remains low during the event, departing from its normal pattern of diurnal fluctuation.
For the 3-hourly NARR data (middle column), the time series plots of A and B grids show similar trends in wind speeds and wind directions to those from RAWS at the caCKNO, caCHAW, and caCSRS, while those trends at C and D grids are similar to those from the caCATL station. The slight differences in the location of the DW and higher wind speeds between the observations and the NARR data are reasonable, considering that the 32 km grid spacing of the NARR data cannot fully resolve the complex terrain of Northern California and locally strong winds influenced by such terrain. Even though the decrease in RH during the event is less dramatic in NARR as compared to RAWS observations, a decreasing trend is still present in association with the event, and the time RH reaches its minimum value, around 1 pm local time, which is similar in both datasets. For the 1-hourly ERA5 data, while the four grids consistently give weaker wind speeds than those from RAWS and NARR, they show similar variabilities in wind speed, wind direction, and RH as those from RAWS at caCKNO, caCHAW, and caCSRS. This suggests that although ERA5 data quantitatively underestimate the wind speeds, again due to different data resolutions, they are able to qualitatively capture the characteristics of this DW event. Similar results can be seen for the DW events during the Tunnel Fire in 1991 and the Camp Fire in 2018 [Figs. S1 and S2 in the supplementary information (SI)]. We note that in the comparison with RAWS, ERA5 data provide a better representation of wind patterns than NARR do, while NARR data provide a better representation of RH over the longer term (Fig. S3). In summary, the three RAWS stations, caCKNO, caCHAW, and caCSRS, and all four grids from ERA5 represent well the wind and RH features during the three DW events. However, only two grids from NARR reasonably capture the features of the DWs. Thus, for consistency, we utilize data over two grids (A and B) from NARR and ERA5 and two stations, caCKNO and caCHAW, from RAWS to identify DW events that have an immediate impact on SFBA. The caCKNO and caCHAW stations are selected instead of the caCSRS station based on the work of Smith et al. (2018), which suggests that observations at the caCKNO and caCHAW stations would capture Diablo-like wind conditions better than those at the caCSRS station.
In addition to the validations based on hourly time series of case studies, we also compare the annual variabilities of the total number and duration of DW events, as well as their associated maximum wind speed and minimum RH among the three datasets (Fig. 4). NARR and ERA5 present similar trends of DW climatology in terms of annual total number, total duration, associated maximum wind speeds, and associated minimum RH values with very high correlation coefficients (Table 2). After 1993 when observations became available, all three datasets present very similar trends of DW climatology. This again validates the use of NARR and ERA5 to study the long-term climatology of DW, and adds robustness to our subsequent results.
We would like to clarify that even though NARR and ERA5 data are capable of resolving the features of DWs and provide similar long-term climatological trends as those from RAWS, there are also notable differences in the number and timing of DWs when comparing the three different datasets. To quantify those differences, the contingency table and skill scores for the DW events identified using RAWS, NARR, and ERA5 are included in Tables S1 and S2. Because ERA5 data are better at representing wind patterns and NARR data are better at representing RH values based on the comparison with RAWS data (Fig.  S3), we also include the DW events identified using the combination of ERA5 and NARR data (RH from NARR and wind speed from ERA5) in Table S1 and Table S2 for comparison. The two tables show that (1) both NARR and ERA5 identify fewer DW events than those from RAWS for the period of 1993-2018, (2) for the period of 1979-2018, DWs from NARR have significant overlapping occurrences with DWs from ERA5, (3) DWs identified from ERA5 and RAWS have more overlapping occurrences than from NARR and RAWS, and (4) there are fewer DW events identified using the combination of NARR and ERA5 (winds from ERA5 and RH from NARR) compared to using either NARR or ERA5 individually. This might be because the RH in NARR and wind speed in ERA5 did not simultaneously meet the criteria of a DW, resulting in fewer DW events.

Relationship between Climate Index and DW Index
In order to relate DWs to climate variabilities, we construct a DW index (DWI) to normalize DWs' annual total duration as follows: where j is a given year between 1979 and 2018, X j is the total duration hours of DWs from September to December of year j, std is the standard deviation of X j , and − X is the average of the total hours of DWs for the past 40 years. Because the annual total duration and total number are highly correlated (r > 0.8, see Sect. 3), we only construct one index, the total duration instead of two for both total duration and total number of DWs.
To explore possible relationships between different climate variabilities and DWs, a large set of climate indices are investigated, including thirteen ENSO-related indices, four Atlantic-Ocean-related indices, two Pacific-Ocean-related indices, six teleconnection-related indices, and three atmosphere-related indices (Table 3) detailed description of each index, including data access and data characteristics, is listed in Table S3 (Liu et al. 2018). Because different climate modes may have different response times, the relations of the DWI to the 28 climate indices are investigated both contemporaneously and with various lead times. Up to three seasons lead-times with respect to DWI are considered for the climate indices, namely summer: June-August (− 1), spring: March-May (− 2), and winter: previous year December to present year February (− 3), where the values in parenthesis refer to the number of seasons in each lead time. The relations are examined using the Pearson's linear correlation coefficients, and are considered to be statistically significant when exceeding a 95% confidence interval.

Composites of different phases of climate indices
To explore possible atmospheric mechanisms that are responsible for the connections between DWs and the correlated climate variabilities, composite analyses of relevant atmospheric variables are conducted in response to different phases of each selected climate indices. The atmospheric variables are monthly anomalies, which are calculated by subtracting the monthly average based on the 1981-2010 period from each monthly data (Liu et al. 2018). The atmospheric variables are composited for both positive and negative values of the climate index to represent the circulation patterns for positive and negative phases of the climate indices. To ensure statistical significance of the differences between the positive and negative phases of each climate index, a two-tailed Student's t test at the 95% confidence interval is conducted. Note that when this study was conducted, the surface analysis from ERA5 was available from 1979 to 2018. However, ERA5s upper air data were only available for the most recent 10 years. Therefore, we use ECinterim instead of ERA5 for the upper-air composite analysis. EC-Interim is a global atmospheric reanalysis dataset available from 1 January 1979 through 31 August 2019. It has been superseded by the ERA5 reanalysis dataset. The spatial resolution of the EC-Interim dataset is approximately 80 km (T255 spectral) on 60 levels on the vertical axis, from surface pressure up to 0.1 hPa. Because the composite analysis focuses on large-scale environmental conditions, the resolution of 80 km of EC-interim is sufficient for this purpose.

Long-term trend of DWs
To understand how DW features have changed with time, we perform a long-term trend analysis using three correlation  Fig. 4a and b, and Table S4). The associated maximum surface wind speeds also show an insignificantly long-term trend (Fig. 4c). However, the associated minimum RH of September-December shows a significant downward trend verified by different statistical methods (Pearson's linear, Kendall's Tau, and Spearman's Rho coefficients) with a confidence level exceeding 95%. (Table 4a and Table S5) from both two reanalysis datasets. Further investigation of the monthly minimum RH trend over the 40 years reveals that according to ERA5 results, both October and November have a statistically significant downward trend (Tables 4b and S5). However, for NARR, only October has, but November does not. Throughout this paper, we only present a conclusion when both datasets show consistent results. Therefore, we conclude that October has a significant downward trend ( Fig. 5 and Table 4b).
Please note although September also show noticeable downward trends, it does not pass the significance test ( Fig. 5 and Table S5). This suggests that the dryness during DW events had become more severe with time, especially in October, possibly leading to an increased likelihood of fires and destructiveness of those fires that occur. It is worth mentioning that maximum surface wind speeds associated with DWs from ERA5 are consistently weaker than those from NARR as shown in Fig. 4c. The differences in maximum wind speeds associated with DWs between ERA5 and NARR might be caused by (1) different boundary layer and land surface parametrizations schemes employed in ERA5 and NARR models, which cause different representations of land-surface roughness and elevation between corresponding ERA5 and NARR grid points, and (2) different density and quality of assimilated data, as well as different data assimilation methods used to nudge wind fields towards observed values between ERA5 and NARR models (Ramon et al. 2019). While no significant long-term trend is found for the annual total number and duration of DWs, the interannual variability/oscillation is evident (Fig. 4). Because the annual total number and duration of DWs show wave-like patterns, their associated periodicities can be determined by using the spacing between peaks of the time series. For the total number of DWs per year, the period ranges from 2 to 4 years with an average of 2.8 years using NARR data, and 2-3 years with an average of 2.4 years using ERA5 data. For the total duration of DWs per year, the period ranges from 2 to 4 years with an average of 2.7 years using NARR data, and 2-3 years with an average of 2.2 years using ERA5 data. This quasi-periodic variation, with period lengths ranging from 2 to 4 years, implies possible linkages between DW development and variations in the atmospheric and oceanic modes. Note that the annual total number and duration both vary in a very similar manner, with correlation coefficients exceeding 0.8.

Correlation analysis
To investigate possible relationships between DWs and atmospheric and oceanic modes, we calculate Pearson's linear correlation coefficients between DWI and 28 climate indices with and without time leads. Since the DWI from RAWS is not available before 1993, our correlation analysis is conducted for the period of 1993-2018. Our results show that among the 28 climate indices, two climate indices, the Southern Oscillation Index (SOI) and the QBO index with two season's lead-times, referred to as SOI (− 2) and QBO (− 2), respectively, have strong positive correlations with DWI from all three datasets, NARR, ERA5, and RAWS (Table 5). In fact, results show that several ENSO-related indices with 1-3 seasons lead time have a significant correlation with DWI from one or more out of the three datasets (Table S6). Among them, SOI (− 2) has the highest correlation from all three datasets. Therefore, we chose SOI (− 2) to represent the impact of ENSO on DWI for our study. In comparison, SOI (− 2) has a slightly higher correlation coefficient than QBO (− 2) from NARR, while for RAWS and ERA5, QBO (− 2) has a higher coefficient than SOI (− 2). Note that only the most recent 26 years of data are used for the correlation analysis, which might not be long enough to identify the impact by the Pacific Decadal Oscillation (PDO) and the Atlantic multidecadal oscillation (AMO) whose oscillation periods are generally longer than 26 years. To visualize the relationship of DW features with SOI (− 2) and QBO (− 2) during the most recent 26 years, a scatterplot of DWI, the two climate indices, and the associated maximum wind speed/minimum RH are shown in Fig. 6. Figure 6 and Table S7 demonstrate three important results: (1) for DWI: QBO (− 2) and SOI (− 2) are both positively correlated with DWI (Table S7). Figure 6a shows that larger La Niña signals (larger positive values on the y-axis) and larger QBOW signals (larger positive values on the y-axis) are positively correlated to larger DWIs (larger positive values on x-axis). On the other hand, larger El Niño signals (larger magnitudes of negative values on the y-axis) and QBOE signals (larger magnitudes of negative values on the y-axis) are positively correlated with smaller DWIs (i.e., larger magnitudes of negative values on x-axis), (2) for maximum wind speeds of DWs: QBO (− 2) and SOI (− 2) also are positively correlated to associated maximum DW wind speeds. Larger La Niña signals (larger positive values on the y-axis) and larger QBOW signals (larger positive values on the y-axis) are positively correlated to larger maximum DW wind speeds (red-ish shading color); larger El Niño signals (larger magnitudes of negative values on the y-axis) and QBOE signals (larger magnitudes of negative values on the y-axis) are positively correlated to smaller maximum DW wind speeds (blue-ish shading color), and (3) for minimum RH values of DWs: QBO (− 2) and SOI (− 2) are negatively correlated with the minimum RH of DWs. Larger La Niña signals (larger positive values on the y-axis of Fig. 6b, c) and larger QBOW signals (larger positive values on the y-axis of Fig. 6b, c) are correlated with smaller minimum RH of DWs (red-ish shading color); larger El Niño signals (larger magnitudes of negative values on the y-axis of Fig. 6b, c) and QBOE signals (larger magnitudes of negative values on the y-axis of Fig. 6b, c) are correlated with larger minimum RH values (blue-ish shading color).
Based on the results of (1), (2), and (3), it is shown that during positive phases of QBO (− 2) (QBOW), the DW events tend to have positive DWI associated with higher maximum wind speeds and lower minimum RH values (Fig. 6). In contrast, during negative phases of QBO (− 2) (QBOE), the DW events tend to have negative DWI associated with lower maximum wind speeds and higher minimum RH values. This suggests that higher DW occurrences with more intense features (high winds and dryness) tend to occur during QBOW conditions than during QBOE conditions. Similar patterns are also observed for SOI (− 2): higher DW occurrences with more intense features (high winds and dryness) tend to occur during positive SOI (− 2) (La Niña) conditions than during negative SOI (− 2) (El Niño) conditions. Previous studies have found that the interaction between QBO and ENSO themselves can also influence the climate variabilities over the North Pacific Ocean, where the QBO effects can be enhanced during La Niña phases and diminished during El Niño phases (Hansen et al. 2016;Sun et al. 2019;Wei et al. 2007). To investigate whether the ENSO-QBO interaction is dominated by a linear interaction or a nonlinear interaction with the activities of DWs (Hansen et al. 2016), we conduct a multiple linear regression model analysis using the two climate indices serving as predictor variables and DWI serving as a response variable. It is assumed that the response variable, DWI, is a function of two linear terms and one nonlinear interaction term of the two climate indices (Hansen et al. 2016): Fig. 6 Scatter plots of DWI versus QBO (− 2) (the 1st and 3rd rows) and DWI versus SOI (− 2) (the 2nd and 4th rows). The colors represent associated a maximum wind speed (m s −1 ) and b minimum RH (%). The red-ish colors represent higher wind speeds and lower RH, while blue-ish colors represent lower wind speeds and higher RH where y represents the response variable (DWI), β o is the intercept term, β 1 and β 2 are the best-fit regression coefficients for ENSO and QBO, and β 3 is the best-fit regression coefficients for the nonlinear interaction term.
The three datasets consistently show that the coefficient (β 3 ) of the interaction term is close to zero (Table 6). Having a coefficient (β3) of the nonlinear interaction close to zero suggests that the potential nonlinear interaction of ENSO-QBO on the DW activities is mathematically insignificant, and therefore, we do not consider this effect in our analysis. Thus, only the linear terms of SOI and QBO are kept in the final regression equations. Please note that having β3 close to zero does not completely eliminate a possibility of the impact of nonlinear interaction between the QBO and ENSO on DW activities. The ENSO and QBO climate modes influence the DW activities mainly through linear interactions ( 1 and 2 > 0) , while their nonlinear interaction plays a minor role in DWs ( 3 ∼ 0).
Even though the nonlinear interaction between QBO-ENSO shows no significant contribution in the multiple linear regression, both modes have influences on the synoptic environment. Therefore, it is of interest to take a closer look at the combined relationship of QBO and ENSO with DWs at their various phase stages. Scatter plots of DW features, such as the annual total duration (DWI) and its associated maximum wind speeds and dryness (minimum RH), across different values of SOI (y-axis) and QBO (x-axis) are shown in Fig. 7. The figure shows that when SOI (− 2) and QBO (− 2) are both positive (La Niña/QBOW), the DW events in that quadrant generally have higher values for the associated DW features such as duration, maximum wind speeds and dryness. When SOI is positive but QBO is negative (La Niña/QBO easterly), there is a mix of higher and lower values for the associated DW features. On the other hand, when SOI and QBO are both negative (El Niño/QBO easterly), the DW events in that quadrant generally have lower values for their associated DW features. For the remaining quadrant where SOI is negative and QBO is positive (El Niño/ QBOW), there is again a mix of higher and lower values for associated DW features. (2) In summary, it appears that the QBO-ENSO interaction has no substantial impact on DW features. However, during La Niña and QBOW conditions, DWs do tend to occur more frequently and have more intense features (e.g., high winds and dryness) than the other ENSO and QBO conditions.
It is worth mentioning that the combined linear effect of ENSO and SOI explains about 50%, 52%, and 49% of DWI variation (quantified using the adjusted coefficient of determination, R 2 adj ) for RAWS, NARR, and ERA5, respectively. The high R 2 adj values from all three datasets support our hypothesis that the last 26-year periodicities in the annual total number and duration of DW events are likely caused by ENSO and QBO from a climatological point of view. Mesoscale and smaller-scale effects likely contribute to the remaining unexplained variations, which is beyond the scope of this study.

Predictability of DWs
The high correlation of ENSO and QBO teleconnections with DWs in the most recent 26 years provides an opportunity to predict DWs. Since the annual total duration of DWs (DWI) are highly correlated with the SOI and QBO two seasons ahead, SOI (− 2) and QBO (− 2) can potentially be used for the seasonal outlook of DWs. To examine the predictability of DWI using SOI (− 2) and QBO (− 2), we use stepwise multiple linear regression analyses to find the optimal regression equation using the data of 1993-2012 as the training period and the data of 2013-2018 as the validation period from each dataset (Table 7). A stepwise regression method is used here because it can reduce the set of predictor variables to those that are strictly necessary and can account for nearly as much of the variance as is accounted for by the entire variable set. The Pearson's linear correlation coefficients (r), root-mean-square error (RMSE), and bias are calculated to quantify the performance of the derived linear regression equations. Please note that we have tested several time ranges for training and validation periods. The values of β 1 and 2 can be different when a different training period is used. However, the overall performance of the linear regression equations is similar.
The results show that the linear prognostic equations from all three datasets consistently show reasonable estimates of DWI for both the training period and the validation period ( Fig. 8 and Table 8). This demonstrates the potential predictability of DWs during the wildfire season in Northern California using the previous spring's QBO or SOI as a predictor. Please note that the equations developed in this study serve as a potential outlook for DWs during wildfire seasons, rather than for an operational purpose.

Individual effects
A connection between DWs and the QBO-related as well as ENSO-related large-scale atmospheric circulation patterns likely exists through two possible mechanisms. First, they can affect the pressure systems near California, which , and c minimum RH (%). The red-ish colors represent higher wind speeds and lower RH, while blue-ish colors represent lower wind speeds and higher RH controls the low-level wind patterns over SFBA. Second, they can modulate the latitudinal location of the 350-hPa Pacific jet stream, which steers the path of winter cyclones and atmospheric rivers (ARs) Liu et al. 2018;Raphael and Mills 1996). The former will be referred to as the pressure anomaly mechanism, while the latter will be referred to as the jet stream displacement mechanism. Regarding the pressure anomaly mechanism, we correlate the March-May average QBO with the September-December average atmospheric variables (Fig. 9). We find that the QBO (− 2) demonstrates a strong positive correlation to 950-hPa geopotential height across the coast of Oregon and Washington as well as the Great Basin (Fig. 9a), which are the critical regions for high-pressure systems favorable to DW development. In addition, the QBO (− 2) shows a strong negative correlation with 2-m RH along California coastal regions (Fig. 9b). To further understand the correlation between the QBO (− 2) and the relevant atmospheric variables, composite analyses of these relevant atmospheric variables are conducted based on their associated positive and negative QBO phases. Figure 10a shows that a positive phase of the QBO (− 2) is characterized by large positive 950-hPa geopotential height anomalies (stabilization) and anticyclonic wind anomalies across the Gulf of Alaska and the Great Basin. This induces northerly to northeasterly 10-m wind anomalies over Northern California (Fig. 10d), which can transport low-RH air from inland and aloft, resulting in large negative 2-m RH anomalies over SFBA (Fig. 10d), favorable for DW development. In contrast, during the negative phases of the QBO (− 2), strong negative 950-hPa height anomalies with a cyclonic circulation anomaly are anchored over the northeast Pacific Ocean (Fig. 10b). The cyclonic circulation anomaly produces onshore flows over SFBA, transporting high-RH air from the ocean, and leading to an elongated band of positive 2-m RH anomalies along the California coast (Fig. 10e). This reduces the likelihood of DW development.
With regard to the jet stream displacement mechanism, our analysis shows that the QBO (− 2) has a negative correlation with 350-hPa wind speed anomalies over California and a positive correlation over Western Canada (Fig. 9c). This means that large negative 350-hPa wind anomalies (i.e.,  easterly anomalies) over the southwestern US and large positive wind anomalies (i.e., westerly anomalies) to the north exist during the positive QBO phases (Fig. 10g), suggesting a northward shift of the jet stream. This reduces the likelihood of winter storms and atmospheric rivers reaching and affecting California by deflecting them northward (Chang et al. 2015;Liu et al. 2018;Wang et al. 2014;Yoon et al. 2015b). This can be seen in the composite of the precipitation anomaly in which large positive precipitation anomalies are located over western Canada (Fig. 10j), while negative precipitation anomalies are located over California. Therefore, the drier conditions (lower RH) in California during the positive QBO phases favor the development of DWs and increase the fire hazard. In contrast, during the negative QBO phases, the opposite patterns are observed. The jet stream shifts southward as suggested by large positive 350-hPa zonal wind anomalies over the southwestern US and large negative 350-hPa wind anomalies to the north (Fig. 10h). This synoptic pattern favorably steers winter storms across the western US, which brings large positive precipitation anomalies in particular, to Northern California (Fig. 10k). The aforementioned differences in 950-hPa geopotential height and winds, 2-m RH, 10-m winds, 350-hPa zonal winds, and precipitation between the positive and negative SOI phases are all statistically significant based on a two-tailed Student's t-test at the 95% confidence interval as shown in Fig. 10c, f, i, and l. Similar to QBO, ENSO, represented by March-May average SOI, appears to influence the DWs through both the pressure anomaly mechanism and the jet stream displacement mechanism as shown in Fig. 11 and Fig. S6. The SOI (− 2) is positively correlated with the 950-hPa geopotential height across the western US (Fig. 11a), and negatively correlated with the 2-m RH along California coastal regions as well as along the Sierra Nevada Mountain Ranges (Fig. 11b), which is similar to those from QBO (− 2). In addition, SOI (− 2) is negatively correlated with 350-hPa wind speed anomalies over California and positively correlated with 350-hPa wind speed anomalies over Canada (Fig. 11c). This dipole correlation is similar to that from QBO (− 2).  Fig. 10 Composite of September-December average a-c 950-hPa geopotential height anomalies (shading; m) and wind anomalies (wind arrows; m s −1 ), d-f 2-m RH anomalies (%) and 10-m wind anomalies (m s −1 ), g-i 350-hPa zonal wind anomalies (shading; m s −1 ), and j-l precipitation anomalies (shading; mm day −1 ) for the March-May average positive (1st column) and negative QBO phases (2nd column) during 1993-2018 from NARR. The difference between the positive and negative QBO phase composites are shown in the 3rd column. The meshed areas represent areas above the 95% confidence level, and only wind vector differences exceeding the 95% confidence level are plotted

Combined effects
The next step is to analyze the combined QBO-ENSO signal (i.e., the anomalies that occur when both phases coexist). Figure 12 shows the composite of 950-hPa geopotential height and winds, 2-m RH, 10-m winds, 350-hPa zonal winds, and precipitation for the four possible combinations of QBO and ENSO phases. It shows that during co-occurrence of La Niña and QBOW, the strong positive 950-hPa geopotential height anomalies with anticyclonic wind anomalies are located across the Gulf of Alaska and the Great Basin (Fig. 12a). In contrast, the other three combinations show strong negative 950-hPa height anomalies with a cyclonic circulation anomaly anchored over the northeast Pacific Ocean (Fig. 12b, c, and d). In addition, negative 2-m RH anomalies over SFBA only appear during the La Niña/ QBOW phase combination (Fig. 12e), but not during any of the other three combinations (Fig. 12f, g, and h).
Similarly, only the composite for La Niña/QBOW shows a northward shift of the jet stream with large negative 350-hPa wind anomalies over the southwestern US and large positive wind anomalies over the northwestern US (Fig. 12i). This decreases the likelihood of winter storms and ARs reaching and affecting California. In contrast, the other three combinations show favorable wet conditions, during which the jet stream shifts southward as suggested by large positive wind anomalies over the southwestern US and large negative wind anomalies to the north (Fig. 12j, k, and l). This synoptic pattern favorably steers winter storms across the western US, which in particular brings large positive precipitation anomalies to Northern California.
In summary, the synoptic environment is the most favorable for DW development during La Niña and QBOW conditions rather than during the other phases of ENSO and QBO. This is because when La Niña and QBOW co-occur in the spring, synoptic environmental conditions in the following fall and early winter tend to be favorable for DW development. The favorable synoptic environment conditions include: (1) strengthened 950-hPa geopotential heights located across the Gulf of Alaska and the Great Basin providing northerly to northeasterly surface wind anomalies transporting low-RH air from inland and aloft, and (2) a northward shift of the Pacific jet stream, which decreases the likelihood of winter storms and ARs reaching and affecting California (low RH). This explains why during La Niña and QBOW conditions, DWs occur more frequently and

Concluding remarks
The most destructive wildfires in Northern California history have been linked to the occurrence of DWs. In spite of DWs' strong connection to ignition and spread of destructive wildfires in Northern California, a long-term climatology study on DWs is still lacking and their relations with largescale climate variabilities have not been explored. With these motivations, we investigate the climatology of DWs for the period of September to December 1979-2018, and their relations with various large-scale climate variabilities (28 climate indices) based on observational RAWS and two high-resolution (~ 30-32 km) reanalysis datasets: NARR and ERA5.
To validate the usage of NARR and ERA5 for this DWs study, we compare NARR and ERA5 winds and RH against RAWS observations during three infamous DW-driven wildfire events: the Tunnel Fire in 1991, the Tubbs Fire in 2017, and the Camp Fire in 2018. The comparisons show that NARR reasonably captures the characteristics of DWs with a rapid increase and decrease of strong northerly winds and dryness. As for ERA5, although it generally underestimates the surface wind speeds, ERA5 data are able to qualitatively capture the characteristic of DWs. We also compare the annual variabilities of total number and duration of DW events, as well as their associated maximum wind speeds and minimum RH values among the three datasets, and they present very similar trends. Our analysis reveals that there has not been a discernible long-term trend in the DWs' annual total number and duration, as well as associated maximum wind speeds over the past 4 decades. However, their associated minimum relative humidity (RH) has decreased significantly with time, especially in October. This suggests that the dryness during DWs has become more severe with time, especially in October, and possibly leading to an increase in the likelihood of fires and exacerbating the destructiveness of those fires that do occur.
We also find that the annual total number and duration of DWs have a quasi-periodic variation, with period lengths ranging from 2 to 4 years. The periodic variability of DWs might be attributed to the teleconnections between DWs and climate oscillations, specifically the El Niño-Southern Oscillation and the quasi-biennial oscillation (QBO), through their modulation of pressure systems near California and of the location of the mid-latitude jet stream over the Northern Pacific. It is suggested that when La Niña and QBO westerly conditions occur in the spring, the synoptic environment conditions in the following fall and early winter tend to be favorable for DW development. The favorable synoptic environment conditions include (1) strengthened 950-hPa geopotential heights located across the Gulf of Alaska and the Great Basin, which provide northerly to northeasterly surface winds advecting low-RH air from inland and aloft, and 2) a northward shift of the jet stream, which decreases the likelihood of winter storms and ARs reaching and affecting California. As a result, when La Niña and QBO westerly conditions occur in the spring, DWs tend to occur more frequently and the associated high winds and dryness become more intense in the following fall and winter. This relationship can be used to predict the seasonal outlook of DWs, which may help California design and implement preventive measures to reduce catastrophic DWdriven wildfires. However, there have been a few named fires that were attributed to DWs when QBOW and La Niña cooccurred in the spring. For example, among the top 20 most destructive wildfires in California history, four of them (the Tubbs Fire, the Nuns Fire, the Atlas Fire, and the Redwood Valley Fire) were contributed by DWs in 2017, following the co-occurrence of La Niña and QBOW in spring of that year.