Simulating Sundowner Winds in Coastal Santa Barbara : Model Validation and Sensitivity

Gert-Jan Duine 1,* , Charles Jones 1,2 , Leila M.V. Carvalho 1,2 and Robert G. Fovell 3 1 Earth Research Institute, University of California, Santa Barbara, CA 93106, USA; cjones@eri.ucsb.edu (C.J.); leila@eri.ucsb.edu (L.M.V.C.) 2 Department of Geography, University of California, Santa Barbara, CA 93106, USA 3 Department of Atmospheric and Environmental Sciences, University at Albany, SUNY, Albany, NY 12222, USA; rfovell@albany.edu * Correspondence: duine@eri.ucsb.edu


Introduction
Downslope windstorms are frequently observed on the south-facing slopes of the Santa Ynez Mountains (SYM) in coastal Santa Barbara (SBC), CA [1,2].The SYM range is about 100 km long and exhibits a predominant east-west orientation with peaks reaching between 800-1200 m.Despite its relatively narrow width (5-10 km), this mountain range plays a critical role in the local meteorology of coastal Santa Barbara, notably in the development of downslope winds.Locally, they are known as Sundowner winds (or "Sundowners"), due to the typical onset timing around sunset [2][3][4][5].All major wildfires in coastal Santa Barbara have intensified under the influence of Sundowner winds.Despite its importance for firefighting purposes, evacuation strategies, aviation and navigation hazards, the existing network of ground stations is insufficient to examine mechanisms associated with the mesoscale evolution of these winds.Therefore, to advance understanding of Sundowner winds, it is imperative to assess the skill of mesoscale models in reproducing the main observed features, such as onset and end of events, timing of maximum intensity and spatial variability of the winds.Downslope windstorms occur when the atmospheric flow interacts with sufficiently large mountain ranges, in particular when strong winds blow with a significant vector component perpendicular to the barrier in a stable environment.Examples of downslope windstorms include the Bora in Croatia, Foehn in the Alps, Chinook in the Rockies, and Santa Ana winds in the Los Angeles basin [6][7][8][9][10].Mechanisms to explain the development of downslope windstorms are generally associated with the reflection of mountain waves, represented by strong positive shear or a two-layer stability structure [11,12], an elevated inversion close to the mountaintop [13], a critical layer associated with either a mean state [14] or a self-induced critical layer caused by mountain wave breaking [15].Although much has been investigated about windstorms in Alpine mountains, for example, [7,12], less is known about windstorms affected by mountains with moderate elevations that are influenced by a shallow, cool and stable marine boundary layer such as observed in coastal Santa Barbara.
Sundowner events can occur under a variety of synoptic forcings and year-round, although with greater frequency during Spring [1,4,16].More importantly, in recent years the wildfires that significantly affected SBC under the influence of Sundowners have occurred in all seasons: summer (Sherpa fire, June 2016 and the Whittier fire, July 2017), late spring (Jesusita fire, May 2009), and winter (Tea fire, November 2008 and the Thomas fire, December 2017).These events further emphasize the importance of understanding mesoscale features of Sundowner winds year-round.Recent studies found that the intensification of downslope winds in the lee of the SYM are attributable to mean-state critical layers (from wind reversal) and self-induced critical layers associated with gravity wave-breaking [4], with the strongest events occurring under influence of strong synoptic pressure gradients [5,17].Although these studies have used mesoscale model simulations and indicated biases in surface windspeed in some locations [4], they have not provided adequate sensitivity analysis to properly identify the best choice of schemes to simulate Sundowners.
Given the deficient network of surface stations, understanding mesoscale features associated with Sundowners depend on accurate simulations of these events.The objective of this paper is therefore to evaluate the skill of high-resolution (1 km horizontal grid spacing) simulations of significant Sundowner events that are seasonally spread, and to enhance our understanding of the mechanisms explaining the development of these downslope windstorms.Here, simulations are performed with the Weather Research & Forecasting (WRF) model [18].For model validation, we focus on simulated wind speed, relative humidity, and temperature at the surface.Model sensitivity is further regarded considering vertical profiles of wind and temperature.Simulations of these variables, including mountain flows, turbulent zones and subsequent wind circulations are affected by model physics options [19] and are mostly sensitive to roughness length z 0 [10,[20][21][22] and planetary boundary layer (PBL) schemes [23,24], while some studies indicate that the PBL scheme choice is of secondary importance [10,22,25].Therefore, the focus of this investigation is on model sensitivity to these parametrizations, including a multi-physics ensemble approach, that appeared a necessity in regions with a low availability of surface stations.
The paper is organized as follows.Section 2 presents the observational dataset, verification strategy and model set-up.Section 3 presents the case studies.Section 4 examines model sensitivity and validation from a statistical perspective, whereas Section 5 discusses a more in-depth model sensitivity to choice of parametrization schemes typically hidden in regions when validating with low availability of observations.Section 6 provides the main conclusions.

Observations
Observations include surface station data retrieved from the Mesowest website [26] maintained by either government agencies or public utilities (Table S1).The NWS stations are located in airports over flat terrain.They measure winds at 10 m above ground level (agl), temperature and relative humidity at 2 m agl; reports at 10 min intervals are used for statistical analyses.The Remote Automated Weather Stations (RAWS, [27]) are the most abundant and are primarily installed in regions prone to fire hazards and high winds (e.g., valleys, ridges, canyons, and sloping terrain).RAWS provide 10 min averages for winds (at 6.1 m agl), temperature and relative humidity (at 2 m agl), and an hourly maximum gust.The SB County Air Pollution Control District (SBCAPCD) station, located on the lee of SYM, observes winds at 6.1 m agl, and temperature at 2 m agl, and reports at minute intervals.In total, 30 surface stations were available in the study area, eight of which are located in the lee of the SYM (Figure 1b).
Additionally, this study uses the National Oceanic & Atmospheric Administration (NOAA) Earth System Research Laboratory (ESRL) 449 MHz wind profiler installed in August 2016 at the Santa Barbara Municipal Airport (KSBA, number 11 in Figure 1b) that provides hourly averages of winds on a vertical resolution of 100 m, ranging from 200 m to about 5 km.

Model Setup and Sensitivity Experiments
This study uses WRF model version 3.9.1 configured with four two-way nested grids with 27 km, 9 km, 3 km and 1 km grid spacing (Figure 1a) and 55 vertical levels, with seven sigma-levels in the lowest 1 km agl.The innermost domain (see Figure 1b) covers the full extent of the SYM, the Santa Ynez Valley and San Rafael Mountains, which may modulate upstream conditions for the development of Sundowner winds [5].
ERA-interim reanalysis [28] at ∼0.7 • latitude/longitude resolution was used for model initialization and boundary conditions; grid-nudging was applied in the outermost domain.The USGS database was used for land cover.The first 18 h was considered spin-up time.Model output was saved on an hourly basis, which is the same time interval as most of the surface station observations.Subgrid physics parametrization schemes included the Thompson scheme for microphysics [29], RRTMG for long-and short-wave radiation [30], and Modified Tiedtke for cumulus in the domains run on 27 and 9 km [31,32].
Of particular importance for our sensitivity experiments are PBL schemes and land surface models (LSMs).PBL schemes compute the vertical transport of heat, momentum, and moisture between the land surface and the lowest troposphere due to turbulent mixing, while LSMs play a key role in partitioning the available energy at the surface and determining the values of z 0 .
Reference [22] showed the relative importance of LSMs over PBL schemes in their simulations of Santa Ana winds, and related this to the different values for z 0 in the LSMs.The influence of z 0 on the development of mountain waves and downslope windstorms has been shown before in idealized simulations [13,20,21].The values for z 0 depend on the roughness of the vegetation, which are provided by land-use and vegetation parameter lookup tables.Depending on the land-use in consideration, the selected LSM can provide different values for z 0 .For example, for open and closed shrublands, which comprise almost 40% of the land-use in the innermost domain (this fraction is even higher on the south side of the SYM), the differences for z 0 range from 0.15 in Pleim-Xiu (PX) LSM [38] to 0.01-0.06m in Noah [39] LSM (see [22], their Figure 8).As the study area is dominated by open and closed shrublands, we expect a similar sensitivity for the wind speed bias as in [10,22].To investigate this issue, higher values of z 0 from the PX LSM were implemented into the Noah LSM (hereafter referred to as Noah-z 0 ) following the strategy of [22].The updated Noah LSM was tested for all selected combinations with PBL schemes considered here.
The selected PBL schemes were run with two additional LSMs: the five-layer thermal diffusion (TD) scheme [40], and Rapid Update Cycle (RUC) [41].We found inconsistencies in the simulations of relative humidity and temperature using the TD scheme, which have been related to soil moisture issues [42].Also, the RUC scheme performed worse for relative humidity than the selected LSM schemes.Since both temperature and relative humidity at the surface are crucial during the fire season, and surface winds were not consistently better, we decided to exclude these configurations using the RUC and TD schemes from the analysis.
Other configuration settings included the hybrid vertical coordinate system, and the inclusion of subgrid-scale orography parameterizations [43,44], and the increase of vertical levels (75 instead of 55).These settings, however, have shown relatively minor impacts on our results and are not considered in our analysis.

Verification Strategy & Limitations
Model output is compared with observations on an hourly basis, starting 18 hours after model initialization.Of specific importance are sustained (i.e.10-min averaged) winds, temperature and relative humidity.The model statistical validation was assessed based on the event-averaged bias (MBE), mean absolute error (MAE), root mean square error (RMSE), and correlation.
Several adjustments have been applied to optimize the comparison between model and observations.Since RAWS measure winds at 6.1 m agl, the standard output of WRF was corrected from 10 m to 6.1 m following [10].This procedure was applied to all configurations, improving the potential bias in the wind speed.A more important consideration concerning the model set-up was the fact that the QNSE and the MYJ schemes specify smaller roughness length z 0 [10] that could introduce bias in simulated 10 m winds.After modifying their respective codes (this was done by commenting the lines U10E = UMFLX/EKMS10 + UZ0 and V10E = VMFLX/EKMS10 + VZ0 in both MYJ and QNSE surface layer scheme codes), the results for the MYJ and QNSE PBL schemes considerably improved (not shown).
Validation with station data were carried out by bilinearly interpolating simulated variables in the four grid cells surrounding the station location.Other tested approaches to compare grid cells with point observations included the choice based on the closest grid cell in vertical or horizontal space (as suggested by [44]), did not improve results or influenced the final conclusions of the manuscript.

Case Studies
To investigate the importance of synoptic forcing and seasonality, this study examines multiple Sundowner events associated with major wildfires that have affected the region: Jesusita Fire (May 2009), Tea fire (November 2008), Whittier fire (July 2017).We also investigated the Sundowner event on 11 March 2017 (which had no association with wildfire) for which data from the NOAA wind profiler were available.The case studies were selected following NWS LOX criteria for Sundowner events: sustained wind speeds of 13.4 m s −1 with a northerly wind direction, or wind gusts above 15.6 m s −1 based on the surface station observations in the foothills of the SYM [2].Notice that multiple surface stations exceeded the NWS LOX-criteria in most Sundowner events.
The Jesusita fire was among the most destructive wildfires affecting populated areas in SBC in the last decade.The fire ignited on 5 May 2009, around 13:00 Pacific Standard Time (PST).Two Sundowner events occurred during the next two days.The strongest one occurred during the evening of 7 May, and drove the wildfire downslope.In the Montecito hills (station 18, Figure 1b), wind speeds of up to 20.1 m s −1 with gusts reaching 30 m s −1 were observed (Figure 2a).This event lasted from 17:00 PST on 7 May to about 06:00 PST on May 8.The period of simulations was between 4 May 2009 10:00 PST and 9 May 2009 16:00 PST.Data from 17 stations were used for comparison with model output.
The Tea fire was ignited on 13 November 2008 around 17:50 PST with strong Sundowner winds on the eastern side of the SYM.Northerly winds meeting the NWS LOX Sundowner criteria were observed at MTIC1 (station 18, Figure 1b) from 15:00 PST until 23:00 PST reaching around 24.6 m s −1 with gusts over 30 m s −1 (Figure 2b).Model simulations started on 12 November 2008 10:00 PST and ran for 48 h.A total of 22 stations were available for model comparison.
The Whittier fire occurred on the western side of the SYM and was ignited in the Santa Ynez Valley on 7 July 2017.Once it reached the mountaintop of the SYM, Sundowners pushed it downslope during the night of 14 to 15 July 2017.RHWC1 (station 21, Figure 1b) observed winds of 13.4 m s −1 with gusts up to 18.3 m s −1 , meeting the NWS-LOX criteria between 16:00 PST and 03:00 PST.No other available nearby stations in the lee of the SYM met the criteria, making this a local event (Figure 2c).Simulations started on 13 July 2017 10:00 PST and ran for 48 h. 25 stations and the KSBA NOAA/ESRL wind profiler were available for model comparison.
For the purpose of seasonality, we complement with a Sundowner event that occurred on the evening and night of 11 to 12 March 2017.Northerly winds were strongest in the western part of the SYM, where the NWS criteria were met at the RHWC1 station (station 21, Figure 1b) with wind speeds up to 14.8 m s −1 and gusts up to 19.7 m s −1 (Figure 2d).Other stations recorded 4.5 m s −1 winds with gusts up to 13.5 m s −1 (MOIC1-station 16, Figure 1b), 12.8 m s −1 winds with gusts up to 16 m s −1 at MPWC1 (station 17, Figure 1b), and 10 m s −1 winds with gusts up to 18 m s −1 (SBVC1-station 25 in Figure 1b), highlighting the spatial variability of Sundowners.Simulations started on 10 March 2017 10:00 PST and ran for 72 h.Data from 27 stations and the KSBA NOAA/ESRL wind profiler were available for model comparison.

Winds, Temperature, and Humidity
We start our analysis with a statistical comparison of wind speed, temperature, and relative humidity (Figure 3).Correlations differ per event and range from around 0.4 for the Tea and 11-12 March cases, to about 0.6 for the Jesusita and Whittier cases (Figure 3a).RMSEs also differ per case from about 3.5-4 m s −1 (Tea, Jesusita) to lower than 3 m s −1 (11-12 March), and about 1.5 m s −1 (Whittier), see Figure 3d.
WRF simulations are robust for all events, but differences between the selected configurations are worth highlighting, especially when focusing on the dependency on z 0 .Recall that the PX and Noah-z 0 configurations have the same higher values for z 0 than the configurations using the standard values in Noah LSM.Higher correlations (Figure 3a) and lower RMSEs (Figure 3d) are found for all respective runs using the modified Noah-z 0 or PX LSMs in comparison with the configurations using the standard Noah LSM.Combining different LSMs into an average PBL-physics ensemble leads to slight improvements.For Jesusita, which encompasses the longest period of comparison among all case studies, the improvement is well-reflected in the values of RMSE, which go from 4.0 m s −1 (standard Noah) to 3.5 m s −1 (both modified Noah-z 0 and PX).The improvement of using higher z 0 also emerges when focusing on single PBL schemes (Figure 3a,d).Larger scatter among configurations is found for temperature and relative humidity than for wind speed (Figure 3b,c,e,f).However, using higher values for z 0 in LSM also improves comparisons for temperature and relative humidity.Highest correlations and lowest RMSEs are found for most configurations using the PX LSM with small differences among the configurations using different PBL schemes.The Noah LSM (both standard and z 0 ) appears relatively more sensitive to a particular PBL scheme than the PX LSM.The QNSE and MYNN schemes show better comparisons with observations (higher correlations and lower RMSEs) especially for the Tea, Whittier and 11-12 March cases.The MYJ scheme underperforms in the Tea and Whittier case studies relative to other PBL schemes.We further note that the Noah-z 0 + MYNN combination performs within the error range of the PX runs, but that the Noah-z 0 + YSU configuration shows a relatively poor performance when compared to the Noah-z 0 + MYNN and PX runs.
Figure 4 shows the statistical analysis for wind speeds for the four Sundowner cases, further generalizing [22] results.As discussed in that study, the different configurations show a large dependency on the use of z 0 : simulations using the standard Noah LSM with lower z 0 (blue dots) show higher wind speed biases than simulations with the higher z 0 (dark green and red dots).The skill of WRF in simulating winds during the Whittier and 11-12 March cases increases for the PX and Noah-z 0 LSMs (Figure 4c,d).The Tea fire simulations show an underestimation for the PX and Noah-z 0 LSMs (red and dark green dots), and a relatively good estimation for the standard Noah LSM.This is related to the general underestimation of winds at most stations (not shown).

Spatial Distribution of Winds
Figure 5 addresses wind biases in a spatial context.For simplicity, we focus on one PBL scheme (MYNN) that proved a relatively good simulation for all cases.Network-averaged statistics are comparable or improved for each event for the simulations using the higher z 0 (compare right with left panels in Figure 5).
In most of the investigated Sundowner events, WRF overestimates winds in the lee of the SYM, which is partly corrected for by using higher z 0 .Although overestimating winds leeward of mountains is not uncommon in downslope windstorm simulations [19], differences are case-dependent.The next section further investigates model sensitivity taking into account the high spatial variability of winds in the lee of the SYM by using one particular case study.

Comparison with Surface Wind Observations
We examine WRF skill and sensitivity in further detail using the 11-12 March 2017 event.This Sundowner event had a relatively high availability of surface stations in the lee of the SYM and included data from the wind profiler.Comparing WRF with surface observations revealed a high variability along the lee of the SYM for all case studies (Figure 5).For example, for the 11-12 March case, MBEs ranged from −4 m s −1 on the western side to more than 3.5 m s −1 on the eastern side (Figure 5g,h).
The RHWC1 station (station 21, Figure 1b) is located on the western side of the SYM at 447 m above mean sea level (msl), which is halfway up the slope of the SYM.Despite a very good skill in capturing wind direction (not shown), WRF greatly underestimated wind speed, leading to a negative MBE of 4.2 m s −1 (Figure 6a).The RHWC1 station is notorious for measuring very high wind speeds with a northerly component throughout the year (NWS LOX, personal communication).The onset and cessation timing of the stronger northerly winds is fairly correct however, leading to a relatively high correlation (0.81).Further eastward, the MOIC1 station is installed at 87 m msl in the foothills of the SYM (station 16, Figure 1b).At this location, wind speeds with a northerly direction reached only 5 m s −1 at peak intensity, while gusts reached 13 m s −1 (Figure 6b).The timing of the model is two hours early in onset, and two hours late in cessation of the strong NE downslope wind, decreasing correlation (0.70) comparatively to RHWC1.All model configurations simulated a northerly wind direction (not shown).However, WRF overestimates wind speed (black solid line), leading to a positive MBE of 3.5 m s −1 (MOIC1).We further note an interesting coherence between simulated winds and observed gusts (solid magenta line).For this event the gust factor for this station is ∼3, which is much higher than the value of 1.7 usually found for denser observation networks, indicating possible anemometer exposure issues for this site [22,45].
The airport (KSBA) is located about 3 km south of the SYM foothills (station 11 in Figure 1b), and some studies have indicated that Sundowners reach the airport only occasionally [1,16].Observed winds were much weaker, but occasionally stronger winds with a northerly direction appear between 20:00 PST (11 March) and 01:00 PST (12 March); these were absent in the simulations (Figure 6c).The weak winds in both simulations and observations lead to relatively low values for MAE, RMSE and MBE, but also decreased correlations.We further note that both configurations using Noah LSM and YSU PBL show much higher values for wind speed throughout the night than the other configurations.Owing to the low availability of the observations, it is hard to discern whether this particular combination performs better or worse than others.We address this particular issue in the next subsections.

Variability of Surface Winds and PBL Height
To further understand model sensitivity and the potential influence of conditions, we focus on north-south transects along the airport longitude (station 11, Figure 1b).We show 10 m meridional wind components and the PBL top, the latter being detected by the bulk Richardson number computed with respect to the surface layer [33] with a critical number of 0.25.Using a consistent strategy for the PBL top, rather than relying on heights reported by the PBL schemes (they all use different methods) enhances the comparability among PBL schemes, and solidifies conclusions drawn regarding upstream profiles.
At 15:00 PST, before the Sundowner onset, winds were calm (between 34.4 and 34.5 • N) in the lee of the SYM (Figure 7a).In the Santa Ynez valley, north of the SYM (between 34.5 and 34.7 • N), the PBL heights ranged from 1000 m to 2500 m msl (Figure 7b).Another indicator of the capping inversion for convective PBLs is the addition of 1.5 K to the minimum mixed-layer temperature [46].Using this approach, we found that most configurations simulated a capping inversion at around 2 km msl (not shown).Both PBL top detection approaches show a general lower PBL top for local TKE schemes (dashed lines) than nonlocal schemes, a difference commonly observed [47].South of the SYM, PBL heights were only 500 m msl or lower, which indicates the influence of a cold and shallow marine layer, emphasizing the very distinct boundary layers north and south of the SYM.
At 20:00 PST, strong northerly winds on the southern-facing slopes of the SYM were simulated in all configurations (Figure 7c).The strongest winds were located approximately halfway of the SYM slopes with magnitudes around 10 m s −1 , and quickly weakened in the foothills.Configurations using PX LSM (red lines) simulated slightly higher wind speeds than the simulations using the Noah LSM (green and blue lines).The PBL tops for the different configurations show very similar values at this time (Figure 7d).
A few hours later, at 02:00 PST, the differences among the configurations over the slope were still small, with similarly strong winds as simulated at 20:00 PST (Figure 7e).However, further south in the foothills (at x = 6 km), differences among configurations emerged.The foothills region appears as a transition zone, where most PX configurations (red lines) simulated near-zero winds (with the exception of YSU), whereas configurations using Noah LSM (green and blue lines) simulated winds of about 5 m s −1 .Further south (at x = 9 km), all configurations simulated near-zero winds, except the standard Noah + YSU and Noah-z0 + YSU configurations, which continued to simulate higher wind speeds on the foothills, the plains and over the ocean.The pattern of stronger winds over a larger area for the simulations using Noah and YSU started around 22:00 PST and lasted up to around 05:00 PST, as was indicated in Figure 6c.
Another relevant issue is the inversion near mountaintops that can play a major role in the development of downslope windstorms.Although we find some sensitivity to PBL schemes in the lee of the SYM for winds, these appear unrelated to the profiles upstream of the SYM, based on the capping inversion.However, all PBL schemes show a typical diurnal PBL evolution, which can be an important factor in the development of Sundowners.A subsequent study is needed to further investigate the potential influence of upstream profiles on Sundowner development, preferably by using a combination of radiosoundings up-and downstream of the SYM, and model experiments.,d,f) vertical cross sections of the PBL height for all configurations (see Figure 6 for the legend) at the indicated times (PST) for the 11-12 March case.The four vertical dashed lines in (e) are used in Figure 8, where the left dashed line indicates the airport location.

Along-Slope Vertical Wind Profiles and Jet Behavior
During the 11-12 March Sundowner event, some of the configurations stood out as outliers from the physics ensemble, notably those using either Noah LSM version in combination with the YSU PBL scheme.These differences are further examined by comparing vertical wind profiles from simulations and observations based on the NOAA wind profiler located ∼9 km south of the SYM crest (station 11 in Figure 1b; left dashed line in Figure 7e).
Winds at the airport at 20:00 PST, the approximate onset time of the Sundowner (Figure 6c), were northerly with peak in intensity (9 m s −1 ) at around 450 m agl (Figure 8a).As the wind profiler does not provide observations below 200 m, the low-level jet cannot be completely characterized based on wind shear at low elevations.Although none of the configurations simulated the northerly wind at the surface at this time (Figure 6c), a clear low-level jet was simulated by all configurations.However, the low-level jet was better characterized in configurations using the PX LSM (except PX + MYJ and PX + MYNN), and these configurations best resemble the vertical profile observed at this time.6 for legend) for the (a,e,i) wind profiler, (b,f,j) foothill, (c,g,k) slope, and (d,h,l) mountaintop locations, see also Figure 7e for the exact locations of each profile.Times (in PST) are indicated in all panels.
Upstream, in the foothills (x = 6 km), the configurations are more coherent in their simulation of the jet structure (Figure 8b).At this location, all configurations simulated a stronger low-level jet with its maximum closer to the surface.Further upstream on the slope (x = 3 km, Figure 8c), the jet's height above ground level is the lowest while wind speeds (15 m s −1 ) are the strongest from all along-slope profiles at this time.At the mountaintop (x = 0 km, Figure 8d), the configurations exhibited the most coherent profiles among all locations.We further note that all configurations at this location show a meridional wind component reaching zero at a similar altitude of around 2 km msl.
Two hours later at 22:00 PST at the mountaintop, the low-level jet is at similar intensity than at 20:00 PST (Figure 8h).Over the slopes, configurations were more coherent than at 20:00 PST (Figure 8g).They diverge in the foothills, where the weakest winds were simulated with PX + ACM2 (red solid lines marked with circles), and the strongest winds were simulated by configurations using the Noah LSM and YSU PBL schemes (solid green and blue lines in Figure 8f).The configurations with weaker and higher jets have a near-zero cross-mountain wind around 1250 m msl, whereas the strongest and lowest jets have zero cross-mountain winds at about 500-750 m msl.Further downstream at the airport, all configurations show better agreement, and are relatively close to observations for both wind speeds and low-level jet altitude (Figure 8e).Simulations with Noah LSM and YSU PBL show a distinct pattern from the other configurations.
Configurations using Noah LSMs and YSU PBL simulated higher wind speeds at the airport location until 02:00 PST (Figure 8i).The deviation of these configurations is evident in the foothills (Figure 8j), but not as much in the slopes (Figure 8k) or at the mountaintop (Figure 8l).This explains why configurations using Noah with YSU showed higher winds at the surface throughout the night compared to the other configurations (Figures 6c and 7e).

Erosion of the Marine Layer
To further understand the differences among the configurations during the 11-12 March case, we discuss the horizontal cross sections of θ v and meridional wind for a selection of configurations (Figure 9) at a time when largest differences between configurations were present (02:00 PST).The figure also highlights regions with a bulk Richardson number (computed based on adjacent model levels) lower than 1 (green contours) and lower than 0.25 (magenta contours).These indicate the presence of self-induced wave-breaking layers [48], and can be a crucial component for the development of downslope windstorms in the lee of the SYM [4].
All presented configurations simulated a wave-breaking layer at 02:00 PST as indicated by the profile of winds and θ v (Figure 9).However, some configurations have developed a clear hydraulic jump in the foothills located around 34.45 • N (Figure 9c-f) whereas in the Noah + YSU combinations this jump has not clearly developed (Figure 9a,b).The θ v profiles further indicate that configurations that have developed the hydraulic jump have not eroded the marine layer over the ocean (Figure 9c-f).In contrast, the Noah + YSU runs which did not develop a hydraulic jump have eroded the cold and stable marine layer (Figure 9a,b).
However, all configurations started with a stable marine boundary layer at the beginning of the downslope windstorm event at around 20:00 PST (this is shown by the PBL height in Figure 7d).Over time this was eroded by the Noah + YSU configurations (indicated by larger PBL heights over the ocean in Figure 7f and in Figure 9a,b).Interestingly, the marine layer has not been eroded in the PX + YSU configuration (Figure 9c).We recall that PX and Noah-z 0 share the same values for z 0 .However, the use of a different LSM leads to a different treatment of surface sensible heat flux H, which is intimately connected to friction velocity u * , especially during the night [49].
To further investigate this particular influence on the erosion of the stable marine boundary layer, it is instructive to explore cross sections of u * and H at times when the downslope windstorm developed (Figure 10).At 20:00 PST, when the configurations were still similar in their simulated surface winds (see Figures 6c, 7 and 8), differences in horizontal extents of u * and H among configurations also appeared small (Figure 10a,b).Two hours later, when configurations diverged (see Figures 6-8), the horizontal extents for u * and H significantly differed, with larger extents for the Noah + YSU configurations (Figure 10c,d), and became even larger at 02:00 PST (Figure 10e,f).Note the relative importance of negative values for H (more mixing) over those being close to zero (less mixing).
Although not investigated in this study, it is worth noting that a similar variability could be expected within a single PBL scheme, by using different constants or mixing lengths specifications.Nevertheless, this mechanism explains the erosion of the marine boundary layer later in the night as seen in Figure 9. Reference [2] hypothesized on an anecdotal basis that the erosion of a marine layer by a Sundowner may lead to a stronger Sundowner event the next day, warranting the importance of model sensitivity in the lee of coastal mountain ranges.

Surface Wind Fields along the SYM
Although the configurations were compared in a north-south transect across KSBA, similar differences can be expected in other transects along the lee of the SYM.The period between 20:00 PST and 02:00 PST is when the largest differences among configurations occurred.The horizontal extent of the strong winds (about 10-12 m s −1 ) was mainly constricted to the slopes and foothills of the SYM, while wind speeds decreased over the ocean, as indicated by the arbitrarily chosen reference configuration Noah-z 0 + MYNN (Figure 11a).Comparing the other configurations against the reference configuration reveals a very similar behavior over the slopes of the SYM (Figure 11b-f  As can be expected from the previous section, further south from the slopes and foothills, differences among the configurations emerge.We find the largest deviations for the configurations that use the standard Noah (Figure 11b,d) vs. those that use the Noah-z 0 (Figure 11a,c).The larger extent of strong winds for models using lower values for z 0 has been shown in [22], and can be expected by the lesser damping due to the lower values of z 0 [21].We find highest sensitivity for the YSU scheme, which shows the largest extent of strong winds (compare Figure 11c,d).Compared to the (arbitrarily chosen) reference scheme Noah-z 0 + MYNN, the smallest deviations are found for the configurations using PX + ACM2 (Figure 11f).In fact, most of the configurations using PX and the Noah-z 0 (with the exception of YSU) are very close to the ensemble mean (not shown).
The 10 m wind plots of Figure 11 spatially confirm what has been shown by the horizontal cross sections (Figure 7) and the vertical profiles (Figures 8 and 9): the configurations using YSU have more difficulty in slowing down the winds over the foothills, coastal plains, and the ocean.We showed in Section 5.4 that this was influenced by different mixing, and so a sensitivity to PBL schemes and LSMs, leading to the erosion of the stable marine boundary layer.Based on a multi-physics ensemble perspective, the 11-12 March case is more likely to be correctly simulated with configurations using PX or Noah-z 0 , excluding the YSU members.These configurations further show better agreement with the wind profiler on the coastal plains at times when wave-breaking occurs (Figure 8).The combination of the horizontal cross sections and vertical profiles of winds from the previous section suggests that the Noah + YSU does not properly simulate the 11-12 March case.Indeed, configurations using the YSU showed poorer skill for this particular case study (see Figure 3a,d).However, we must acknowledge that observations of turbulent properties in the vertical profile in the lee of the SYM and over the ocean are necessary for future model assessment.This could be typically done in an extensive field campaign.

Conclusions
We have investigated the sensitivity of the WRF model to simulate Sundowners.This investigation focused on the WRF model sensitivity to PBL schemes and LSMs during four selected Sundowner events that occurred in different seasons under different synoptic situations.By comparing simulations of winds with surface stations, a wind profiler, and using a multi-physics ensemble approach, we found that the model wind speed bias for the four events is explained by z 0 , and is therefore directly related to the choice of the LSM: higher z 0 leads, in general, to better simulations of Sundowners.Altering the values of z 0 in the Noah LSM (which by default has lower values for z 0 than PX LSM) leads to comparable results to the PX LSM scheme.Even though the availability of surface stations in the lee of SYM was low (depending on the case study, only three to five stations were available), and one may question the representativeness of some of the surface stations, these results agree with [22].Additionally, the simulations of relative humidity and temperature improve using LSMs with higher z 0 .
In general, all PBL schemes show similar timing for the onset and cessation of the Sundowner cases investigated here.The latter is most likely a consequence of the similar behavior and timing of PBL schemes for transitions from unstable to stable regimes.However, regarding the potential influence of upstream profiles on Sundowner development, a subsequent study is needed, preferably by using a combination of radiosoundings up-and downstream of the SYM, and model experiments.In a spatial sense, little sensitivity among PBL schemes is found over mountaintops and over the slopes south of the SYM.However, differences emerge from the foothills and further downstream over the coastal plains and the ocean.Configurations that simulated different locations of hydraulic jumps and consequent larger horizontal extents of u * and H in the lee of the SYM were also in favor of eroding the marine boundary layer.This sensitivity was not revealed when compared to the observations, owing to a low availability of surface stations.Nevertheless, this particular sensitivity is important for simulations of downslope windstorms in coastal mountain ranges in general.
The multi-physics ensemble approach showed that future simulations in this particular area will be more realistic when performed with LSMs that use relatively high values for z 0 , and with PBL schemes that adequately treat the interaction between downslope winds and the marine boundary layer.Proper simulations of the horizontal extent of strong surface winds in the lee of the SYM is critical for both firefighting purposes and evacuation strategies in an event of a wildfire during Sundowner winds.In a broader context, this study contributes to the improvement of understanding of downslope windstorms in coastal areas with similar topography and land cover that are strongly influenced by the diurnal cycle of the marine boundary layer.

Figure 1 .
Figure 1.(a) Domain configuration for the WRF simulations, and (b) innermost domain and the available stations for model evaluation.Station locations are indicated by their number corresponding to Table S1.The red line indicates the location for the vertical cross sections used in Figures 7, 9and 10.

Figure 2 .
Figure 2. Observed wind speeds (solid lines) and gusts (dashed lines) for the selected Sundowner events during (a) the Jesusita fire, (b) the Tea fire, (c) the Whittier fire, and (d) the 11/12 March 2017 case study.Station location in each of the panel titles are given in Table S1 and Figure 1b.Horizontal lines indicate the thresholds for Sundowner wind criteria as used by NWS LOX.Vertical dashed lines denote astronomical sunset and sunrise times.

Figure 3 .
Figure 3. Statistical comparison for all investigated Sundowner events (indicated in the legend in (a)) for wind speed (a,d), temperature (b,e), and relative humidity (c,f) and the model configurations indicated in each panel.

Figure 4 .
Figure 4. Scatterplot of network-and event-averaged sustained wind speed MAE vs. bias (MBE) from 13 physics ensemble members for the (a) Jesusita, (b) Tea, (c) Whittier, and (d) 11/12 March 2017 case studies.The legend in (a) applies on all panels.Data periods used for comparison are indicated in the text of Section 3.

Figure 5 .
Figure 5. Event-averaged MBE of wind speed (m s −1 ) for the (a,c,e,g) standard Noah and (b,d,f,h) modified Noah-z 0 configurations for the Sundowner events indicated in the panel titles.Each panel also shows the event-and network-averaged values for correlation, MAE, MBE, and RMSE.

Figure 6 .
Figure 6.Observed and simulated wind speed at the (a) RHWC1, (b) MOIC1, and (c) KSBA stations (see Figure 1b for their location) in the lee of the SYM during the 11-12 March case study.The statistics in each panel represent physics-averages of all ensemble members.The black (sustained wind) and magenta (gust) dashed lines show the threshold for Sundowner events.

Figure 7 .
Figure 7. (a,c,e) Horizontal cross sections of meridional wind component at 10 m agl and (b,d,f) vertical cross sections of the PBL height for all configurations (see Figure 6 for the legend) at the indicated times (PST) for the 11-12 March case.The four vertical dashed lines in (e) are used in Figure 8, where the left dashed line indicates the airport location.

Figure 8 .
Figure 8.Time evolution of vertical profiles of the meridional wind during the 11-12 March Sundowner case study from observations (black line, only at the wind profiler location) and model sensitivity experiments (see Figure6for legend) for the (a,e,i) wind profiler, (b,f,j) foothill, (c,g,k) slope, and (d,h,l) mountaintop locations, see also Figure7efor the exact locations of each profile.Times (in PST) are indicated in all panels.

Figure 10 .
Figure 10.Same as Figure 7 but for (a,c,e) friction velocity and (b,d,f) surface sensible heat flux.Note that the times and axis differ from Figure 7.
) indicating that different configurations show relatively low sensitivity in the 11-12 March 2017 case along the southern slopes of the SYM.

Figure 11 .
Figure 11.(a) Simulated 10 m winds using Noah-z 0 + MYNN configuration as a reference, and (b-f) differences between the indicated configuration and reference configuration after averaging between 11 March 2017 20:00 PST and 12 March 2017 02:00 PST.