Entrainment Rates and Their Synoptic Dependence on Wind Speed Aloft in California's Central Valley

Daytime atmospheric boundary layer (ABL) dynamics—including potential temperature budgets, water vapour budgets, and entrainment rates—are presented from in situ flight data taken on six afternoons near Fresno in the San Joaquin Valley (SJV) of California during July/August 2016. The flights took place as a part of the California Baseline Ozone Transport Study aimed at investigating transport pathways of air entering the Central Valley from offshore and mixing down to the surface. Midday entrainment velocity estimates ranged from 0.8 to 5.4 cm s−1 and were derived from a combination of continuously determined ABL heights during each flight and model-derived subsidence rates, which averaged -2.0 cm s−1 in the flight region. A strong correlation was found between entrainment velocity (normalized by the convective velocity scale) and an inverse bulk ABL Richardson number, suggesting that wind shear at the ABL top plays a significant role in driving entrainment. Similarly, we found a strong correlation between the entrainment efficiency (the ratio of entrainment to surface heat fluxes with an average of 0.23 ± 0.15) and the wind speed at the ABL top. We explore the synoptic conditions that generate higher winds near the ABL top and propose that warm anomalies in the southern Sierra Nevada mountains promote increased entrainment. Additionally, a method is outlined to estimate turbulence kinetic energy, convective velocity scale (w*), and the surface sensible heat flux in the ABL from a slow, airborne wind measurement system using mixed-layer similarity theory.


Introduction
The atmospheric boundary layer (ABL) is the layer of air in direct communication with the surface of the Earth via convective thermals and friction-induced turbulence, characterized by eddies that turnover on time scales of less than about one hour (Stull 1988). ABL entrainment is the spontaneous mixing that occurs at the interface between the top of the turbulent boundary layer and the lowest portion of the laminar, or less turbulent, free troposphere, and this process is known to significantly influence the heating, composition, and growth rate of the ABL. Thus, the nature of the ABL is, to first order, often assumed to be driven by the surface buoyancy and momentum fluxes that lead to a mix of convective and shear-driven turbulence, respectively. However, features of the broader meteorological (synoptic) setting that influence the top of the ABL are also known to be important, such as subsidence (Ball 1960;Driedonks and Tennekes 1984) and vertical wind shear in the lower free troposphere (Fedorovich and Conzemius 2008). Vertical wind shear changes the structure of turbulence and the mean flow in the ABL relative to purely convective boundary layers (Sykes and Henn 1989;Moeng and Sullivan 1994;Khanna and Brasseur 1998).
As a result, the shear at the top of the ABL, across the entrainment zone, is particularly responsible for enhancing entrainment fluxes and thus ABL growth (Kim et al. 2003;Conzemius and Fedorovich 2006;Pino and De Arellano 2008). A large-eddy simulation (LES) study conducted by Conzemius and Fedorovich (2006) demonstrated that while the surface heat flux is the strongest predictor of ABL growth, vertical wind shear across the ABL inversion also enhances it, and typically has a far greater impact on the entrainment flux than surface shear does. This potential for shear enhancement of the entrainment flux and ABL growth is determined by the atmospheric stability across the entrainment zone, as greater stability will lead to less vertical mixing of momentum and thus allow stronger shear to develop. This interplay between shear and stability has traditionally been relegated to some parametrization of the bulk Richardson number across the inversion (Fairall 1984;Sorbjan 2004;Conzemius and Fedorovich 2006;Traumner et al. 2011).
Early mixed-layer models parametrized boundary-layer growth by assuming that the downward sensible heat flux at the ABL height was some fixed fraction of the upward surface sensible heat flux (Ball 1960;Tennekes 1973). This was based on the assumption that the sensible heat flux at the top of the ABL is driven solely by the surface heat flux. The (negative) ratio of the two fluxes, the entrainment efficiency (A R ), was concluded to be~0.2 by a literature survey conducted by Stull (1976) with most reported values falling within the range of 0.1 to 0.3 in the absence of significant mechanical turbulence. However, individual field observations have reported values as large as 0.64 (Davis et al. 1997) and even 1.05 (Angevine et al. 1998). Several field studies estimating entrainment efficiency over a wide variety of surfaces present results that suggest a robust dependence on wind shear (or simply wind speed) at the ABL top, although this is by no means conclusively established. Barr and Betts (1997) studied A R using radiosonde budgets over a boreal forest and found an increase in A R with both increasing mixed-layer wind speed and increasing entrainment layer wind jump. Davis et al. (1997) found a "modest positive correlation" between A R and the jump in the mean wind across the ABL top using airborne budgets over a similar boreal forest. Flamant et al. (1997) report values ranging from 0.15 to 0.30 depending on the shear across the ABL over the Mediterranean Sea using airborne lidar measurements. In a study over agricultural flat terrain in Illinois using radar wind profilers and radiosondes, Angevine (1999) concluded that there was a "weak dependence" of entrainment efficiency on wind speed or shear at the ABL top. In this work, we present further data that evince a role of wind The San Joaquin Valley (SJV), the southern branch of California's Central Valley, is long and wide, running from the San Joaquin River Delta near Stockton (at sea level) over 400 km south-east to Arvin, CA (≈ 135 m above sea level (a.s.l.)). Its expansive width ranges from 60 to 100 km bound by asymmetric rims with the crestline of the Coast Range to the south and south-west rising near 1.5 km a.s.l., the Tehachapi Mountains at its south-eastern terminus rising to approximately 2 km a.s.l., and the high Southern Sierras to the east and north-east averaging around 3 km a.s.l. (see Fig. 1). The summertime Pacific high drives a persistent north-westerly boundary-layer wind off the California coast that runs roughly parallel to the coast south of Cape Mendocino (Dorman et al. 2000). At the shoreline, however, friction resulting from this flow interacting with the coastal plains causes it to align more directly with the pressure gradient, giving it an onshore component. This onshore wind is mostly blocked because the shallow coastal marine boundary layer typically lies below the elevation of the coastal mountains. Nevertheless, the flow is able to inundate the interior of the state at breaks in the orography such as the Petaluma Gap and especially the San Francisco Bay area, feeding through the Carquinez Strait into the Central Valley (Schultz et al. 1961;Frenzel 1962;Moore et al. 1987;Zaremba and Carroll 1999). This flow is augmented by the extreme land-ocean thermal contrasts found throughout the region. Upon entering the valley, the flow is diverted north into the Sacramento Valley and south-east into the SJV as it butts up against the Sierra Nevada mountains. The up-valley wind in the SJV, aligned with the background, north-westerly flow offshore, is enhanced by a pronounced horizontal temperature gradient that rises and falls diurnally from ≈ 3 K (late afternoon) to ≈ 6 K (near midnight) between Sacramento and Bakersfield (Zhong et al. 2004). This sustains a pressure gradient that drives up-valley winds that persist throughout the entire diurnal cycle except near the terminal and sidewall mountains where drainage flows develop overnight (Zhong et al. 2004;Bao et al. 2008;Beaver and Palazoglu 2009).
An observational study of ABL depths throughout the Central Valley by Bianco et al. (2011) highlighted the unique and unexpected occurrence of a minimum during the summer months when surface radiative forcing is strongest. The authors conclude that synopticscale subsidence cannot explain the reduction of the valley ABL heights, suggesting instead that cold air advection and/or surface irrigation, which is difficult to constrain in mesoscale models (Alexander et al. 2022), may be the more likely causes. Some evidence for the surface irrigation hypothesis, which indirectly lowers the ABL heights by reducing the surface sensible heat flux, was provided in a later modelling study by Li et al. (2016). Nonetheless, we suggest that the conclusions of Bianco et al. (2011) may reflect an underappreciated consequence of the thermally forced mountain-valley system in the SJV, which can result in a strong daytime subsidence aloft. A more recent airborne study estimated midday entrainment rates of the valley ABL in the winter and summer (Trousdell et al. 2016), but did not seriously consider the flow above the ABL. Faloona et al. (2020) discuss the general circulation across California and proposed a three-layered model of the atmosphere above the SJV based on observations from over two dozen flights from two different aircraft. The three-layered model includes the ABL and free troposphere, but adds a "buffer layer" between approximately 600 m and 2 km above the valley floor that acts as a storage basin for pollutants above the ABL but below the free troposphere. In this study, we analyse a subset of that dataset from a sequence of six summertime flights using ABL budgets, an abundance of deeper profiles, and the Weather Research and Forecasting model (WRF) to investigate how the layering of the atmosphere above the SJV and the synoptic conditions influence the ABL entrainment.

Aircraft Instrumentation
Aircraft data were collected by a Mooney Bravo and Mooney Ovation, which are fixed-wing single engine airplanes operated by Scientific Aviation Incorporated (Boulder, Colorado, USA). The wings are modified to continuously sample air via inlets that flow to the on-board instrumentation. Temperature was measured by a probe (HMP60, Vaisala, Finland). Methane, carbon dioxide, and water vapour are measured with a cavity ring down spectrometer (2301f, Picarro, Santa Clara, California, USA), which is operated in its 1 Hz precision mode. Winds are measured using a Duel-Hemisphere Global Positioning System (VS101, Hemisphere, Scottsdale, AZ, USA) combined with direct airspeed measurements obtained via an onboard flight computer (PFD1000, Aspen Avionics, Albuquerque, NM, USA), as described in Conley et al. (2014). In short, this low-cost horizontal wind measurement system employs a Global Positioning System (GPS) with twin antennas that allow for a high heading accuracy (± 0.11°), which improves the GPS-derived ground velocity vector. Combined with the true airspeed measurement, a simple vector subtraction is then performed to obtain horizontal wind speeds at a rate of 1 Hz.
The aircraft flight paths were designed to probe the regional ABL by a series of long flight legs that run roughly parallel to the SJV axis between Fresno and Visalia (Fig. 1). In general, level legs were flown at constant altitudes within the ABL spanning from 60 m tõ 1.5 km above ground level (a.g.l.), well above the ABL top, which is determined from in situ profiles of temperature and relative humidity. Periodically, the aircraft conducted deep vertical profiles from near the surface (≈ 10 m) to at least 3 km a.g.l. in order to capture the vertical structure of the measured scalars. The ABL heights were determined from these profiles by visual inspection of temperature and water vapour. Namely, from a given profile, the altitude where water vapour reached its background free-tropospheric quantity and the altitude where potential temperature began to increase by >~10% of its mixed-layer value were averaged to obtain the ABL height at the specific time and location of the profile that penetrated the ABL top. In total, approximately 10 penetrations of the ABL top took place on each 4-h flight (see Fig. 5). The flights specifically targeted the time of the day (~1200 to 1600 LT, local time = UTC-8 h) when the ABL is actively growing, but after the initial rapid growth phase through the near-neutrally stratified residual layer in the morning after sunrise.

Model Configuration
The WRF model version 3.8.1 (see Skamarock et al. 2008) was used to provide the mean vertical velocity component, which are essential to the study but not measured by aircraft during the flights (see Alexander et al. 2022). The model was configured to use two-way nested domains with 12-and 4-km horizontal resolutions. Much of the coarser domain covers the western United States, while the finer resolution domain is centred over California. In the vertical dimension, the model used 50 terrain-following levels with approximately thirty levels below 5 km and an increased resolution (~50 m vertical spacing) in the lowest 500 m, corresponding to the approximate daytime SJV boundary layer depth. Surface characteristics were prescribed by the United States Geological Survey (USGS) land use dataset, and the North American Regional Reanalysis (NARR) provided initial conditions and boundary conditions, including sea surface temperatures, at a 3-hourly timestep. In addition, observational nudging for wind speeds and surface-level temperature and water vapour was applied to the coarse domain with four-dimensional data assimilation (FDDA) using the 6-hourly datasets from the National Centers for Environmental Prediction (NCEP) Administrative Data Processing (ADP) Global Surface Observational Weather Data (ds461.0) and Upper Air Observational Weather Data (ds351.0). This simulation implemented the Yonsei University planetary boundary layer (YSU PBL) and the Rapid Update Cycle (RUC) land surface model (LSM) parametrizations, which were coupled using the MM5 surface-layer scheme (Benjamin et al. 2004;Hong et al. 2006;Jimenez et al. 2012) without top-down mixing. The YSU PBL parametrization is a nonlocal turbulence closure that uses an eddy-diffusion model that allows for counter gradient fluxes, and employs explicit treatment of the entrainment zone (Hong et al. 2006). The RUC LSM models bare soil and vegetated ground for both the heat and soil moisture diffusion over 9 soil layers (Smirnova et al. 1997). The RUC model imposes a "proxy" minimal irrigation for cropland designations during the growing season by forcing all soil to be at least 20% higher than the wilting point Benjamin et al. 2016). Other physics parametrizations include the rapid radiative transfer model for general circulation models (RRTM-G) shortwave and longwave radiation parametrizations, the Morrison double moment microphysics parametrization (Morrison et al. 2009), and the Kain Fritsch Cumulus scheme in the coarser model domain (Kain 2004).

Turbulence Measurement Methodology and Similarity Technique
Utilizing the horizontal 1 Hz wind measurement from the Mooney aircraft, we tested a novel, low-cost method to obtain aircraft turbulence measurements. To check the capability of this airborne wind measurement system to capture the most relevant turbulent eddy scales, a fast Fourier transform was applied to the data to create power spectra of the horizontal winds. The expected power spectra for the inertial subrange of turbulence have the form: where S(k) is the power contained at a given eddy size (inversely proportional to the eddy wavenumber, k), ε is the viscous dissipation rate of turbulence kinetic energy TKE, and α k is the Kolmogorov constant, which is~0.52 for the streamwise component of wind (Stull 1988;Sreenivasan 1995;Ni and Xia 2013). The wavenumber is related to the frequency (f ) of the wind measurement and the true airspeed of the aircraft (U TAS ) by the conversion factor k = 2πf /U TAS . This analysis was applied to 77,300-s samples of level flight data from the ABL to generate independent power spectra that were then averaged together within equal-log bins. No windowing was ultimately used to generate the spectra in this analysis, after testing a windowing function and observing that it did not discernably affect the results. The Nyquist frequency of the wind measurement system is 0.5 Hz, and the average true air speed was 75 m s −1 . Thus, the smallest possible turbulent eddies that could be measured have length scales of tens of metres. The limitation of the aircrafts ability to measure smaller eddies implies that our measurements of wind variance, σ u 2 and σ v 2 , will be underestimated. To compensate, we estimated how much of the total wind variance was not being captured by the measurement system by multiplying the power spectra (S) by the frequency contained in each interval. When the x-axis is transformed to a log scale, this procedure results in a graph where the area under the curve is proportional to the wind variance (Stull 1988). The inertial subrange of turbulence is expected to follow a -5/3 power law (Eq. 1), which becomes -2/3 (-5/3 + 1) when multiplying the y-axis by f . Applying a f −2/3 fit through the inertial subrange to extend the spectra out to the higher frequencies shows that about 18% of the total variance is unaccounted for (Fig. 2) by taking the ratio of the red shaded area to the total red area + blue area. Thus, we apply a correction factor of 1/(1-0.18) = 1.218 to our direct measurements of σ u 2 and σ v 2 . After applying this correction to the horizontal wind variances, we can estimate both the convective velocity scale (w * ) as well as TKE. The variances of the horizontal winds have been shown to relate to the convective velocity scale (w * ) in the convective mixed-layer as (Panofsky et al. 1977;Caughey and Palmer 1979): Using this relationship, and the observation from our data that σ u~σ v for a coordinate system that is aligned with true north, we can calculate TKE (per unit mass): by including another relationship from Lenschow et al. (1980) which has been rearranged to solve for σ 2 w : .
From Eq. 2, we estimate w * as: With Eq. 6 and the definition of w * ≡ [Q SH z i g/θ v ] 1/3 , the surface sensible heat flux (Q SH ) can be estimated as: where g/θ v is the buoyancy parameter. The kinematic flux units of K m s −1 are converted to W m −2 by multiplication of ρ c p , where ρ is the surface air density and c p is the specific heat capacity of air. The air density is determined from the Visalia Municipal Airport's Automated Weather Observing System (AWOS). It is worth noting that the above equations are framed in the convective-dominated mixedlayer perspective, while our ABL environment contains significant shear as exhibited by our vertical profiles. However, σ w is likely dominated by convection in our ABLs to enough of a degree where Eq. 4 remains relevant, as suggested by the fact that the highest ABL wind shear present in the Lenschow et al. (1980) study is about twice as great as the shear we estimate in our study by dividing the winds near the ABL top by the depth of the ABL (see Sect. 3.4). Analysing the average w * for all six flights as a function of σ u and σ v sampling times (Fig. 3), the value of w * asymptotically approaches ≈ 1.5 m s −1 at 300 s of sampling time. Thus, we conclude that 300-s samples adequately capture the dominant production-scale turbulent eddies of the ABL, providing a relatively accurate estimate of turbulence parameters from a low frequency (horizontal) wind measurement.

Scalar Budgeting
Our analysis here aims to quantify the source terms for enthalpy to the ABL, i.e. the entrainment and surface sensible heat fluxes, in the flight region. Outlined in the seminal work of Lenschow et al. (1981) are original applications of the scalar budgeting techniques (Telford and Warner 1964;Warner and Telford 1965;Lenschow 1970) originally designed to help aircraft validate the emerging technique of eddy covariance for measuring sensible heat fluxes. Lenschow et al. (1981) further describes the effectiveness of well-designed aircraft ABL studies in determining the net source or sink of ozone in a given region by careful measurement of the other dynamically controlled terms. The technique can be applied to any scalar in a turbulent medium (i.e. potential temperature, NO x , water vapour, SO 2 ) to enable the calculation of important residuals including source or sink terms for non-conserved species (Kawa and Pearson 1989;Conley et al. 2009;Faloona et al. 2009;Bandy et al. 2011;Caputi et al. 2019). Our scalar budget equation for potential temperature is the Reynolds-averaged conservation equation of a scalar in the turbulent mixed layer of the ABL, ignoring radiative heating/cooling, assuming horizontal homogeneity of turbulence and neglecting the molecular diffusion term: where θ is potential temperature and − w θ z i is the entrainment heat flux at the ABL top. The last term is a bulk approximation of the vertical flux divergence over the depth of the ABL. The advection term here is the mean streamwise horizontal wind (U) multiplied by the horizontal potential temperature gradient ( ∂θ ∂ x ) where x is the axis of mean wind, which typically aligns with the valley axis.
Without a high-frequency turbulence probe on the aircraft, the entrainment flux cannot be directly measured. By using the zeroth-order jump model derived by Lilly (1968), we parametrize the entrainment flux as: where θ is the difference in potential temperature between the free troposphere bottom and ABL top, and w e is the entrainment velocity, the measurement of which is discussed in Sect. 2.5. To obtain the time and space derivative terms in Eq. 8, we first determine the ABL heights from the vertical flight profiles, which are conducted at various locations and times throughout the flight. We then apply a linear regression to approximate a time-dependent ABL height, i.e. ∂z i /∂t, to determine when the aircraft was inside the ABL and parse the dataset accordingly. A multi-linear regression is then applied to the subset of data in the ABL to obtain all of the derivative terms found in Eq. 8 in a similar manner to the technique outlined in Conley et al. (2009) and Caputi et al. (2019). Our multi-linear regression was performed using QR (matrix) decomposition to obtain a linear least-squares fit. The mean winds for the horizontal advection term are computed by taking the vector mean winds of all data points inside each 20-m vertical bin, then averaging the mean winds from each bin to obtain a representative average boundary layer wind. The surface sensible heat flux (Q SH ) then becomes the one unknown variable in Eq. 8 which allows for closure of the budget.
The above process is repeated for a water vapour scalar budget analysis, using the measured water mixing ratio q which replaces the scalar θ in Eqs. 8 and 9. Likewise, the surface latent heat flux (Q LH instead of Q SH ) becomes the single unknown variable in the scalar budget equation, which is computed by closing the equation for each flight.

Entrainment Velocity and Boundary Layer Height
The mean vertical wind at the ABL top, which can be quite significant in mountain-valley systems (Rampanelli et al. 2004), is important to consider for quantifying boundary-layer growth. In the summertime afternoon in the SJV, subsidence of order 1-5 cm s −1 is estimated to occur based on nearby observations in the Sacramento Valley (Myrup et al. 1983), although there remains a lack of direct measurements of this quantity. To estimate the entrainment velocity, we expand the total derivative of the ABL height ( dz i dt ) into the Eulerian derivative, which adds an advection term (Albrecht et al. 2016;Trousdell et al. 2016). The resultant z i budget equation leads to a relationship between the entrainment velocity (w e ), the ABL growth rate ( ∂z i ∂t ) and the mean vertical velocity at the ABL height (W z i ): The advection and growth rate of ABL height (first two terms on right-hand side of Eq. 10) are directly observed by aircraft, using the wind speed from within 150 m of the ABL top (U zi ) as discussed in Sect. 3.4. The last term, W zi or vertical wind at the ABL top, has evaded careful measurement due to its relatively diminutive magnitude (Angevine 1997;Lenschow et al. 1999Lenschow et al. , 2007. Although recent attempts with airborne measurements to measure vertical velocity in the fair-weather atmosphere have been somewhat successful (e.g. De Wekker et al. 2012;Cooper et al. 2016), the expected magnitude of synoptic-scale subsidence in the SJV is several factors smaller than what has been demonstrated in those studies. With a lack of methodology to directly measure W zi , we instead rely on the WRF model described in Sect. 2.2. To extract this quantity from the WRF results, a polygon was generated that encompasses the general region of the flight area, which is shown in Fig. 4 along with the average 700 m a.s.l (i.e. near ABL top) subsidence values from 1200 to 1600 LT across all 6 flight days. As shown in Fig. 4, the model simulates areas of uplift all along the rim of the Central Valley from the coast range on the west side and the Sierra Nevada foothills on the east, representing the daytime up-slope winds along the foothill regions. In the valley, subsidence is generally observed, with the exception of a narrow band of weak uplift that runs through the valley centre. This is consistent with a water tank study by Reuten et al. (2007) that showed counter-directed horizontal rolls in a plain-plateau system, with low-level convergence in the centre of the plain region.
To determine W zi for a given day, the vertical velocities at the heights just below and just above the observed ABL height were taken at each gridbox within the polygon. Then, Fig. 4 Mean 700 m a.s.l. vertical velocity simulated by WRF for the afternoon hours (1200-1600 LT) of the six flights, showing subsidence throughout most of the SJV. Velocities over mountainous terrain that are likely orographically driven are not shown the average vertical velocity was computed for both the extracted plane of grid cells below z i and the one above. A linear interpolation was then used to determine the value of W at the exact ABL height. To investigate the sensitivity of these subsidence rates, we repeated the calculations over slightly different polygonal areas containing the majority of the flight tracks and shifting the time interval by an hour for eight distinct calculations. The standard deviations of these varying circumstances for all six flights averaged 0.5 cm s −1 , therefore we use this value as an estimate of uncertainty in WRF-derived subsidence (see Table 1) even though this does not address the accuracy in the model's vertical velocity.

Scalar Budget Results
The results of our scalar budget analysis for boundary-layer height and potential temperature are presented in Tables 1 and 2, respectively, where each computed term in the budget equation is shown for each flight. The flight durations were typically 3.5-4.0 h taking place between~1200 and~1600 LT. Of note is that the entrainment velocity ranges from 0.8 to 5.4 cm s −1 with an average of 3.0 cm s −1 , which is consistent with prior results obtained in this environment during the summertime (Karl et al. 2013;Trousdell et al. 2016). The average surface sensible heat flux is 170 W m −2 and the average ABL height is 540 m a.g.l., consistent with those reported in Bianco et al. (2011) andFaloona et al. (2020). The average advection term of the enthalpy budget is relatively small, -0.07 K h −1 , which is very similar to the values reported in Bianco et al. (2011) for the similar region and confirms findings from prior studies that use the scalar budget technique showing that the average advection is a relatively small contribution to the average total budget Conley et al. 2011;Trousdell et al. 2019;Caputi et al. 2019). However, in some of those aforementioned studies, the advection term on any given day can be quite large and contribute significantly to the scalar budget. In our present study, on no given day does the advection contribute more than 0.3 K h −1 of the average 1.15 K h −1 heating to the ABL, leaving the surface and entrainment heat fluxes to be the dominant sources of enthalpy to the ABL ( Table 2).
The results of the scalar budget analysis for water vapour are presented in Table 3, which again shows the contribution of each term in Eq. 8 by flight. The anomalous results for the 6 August flight show a surface latent heat flux of -279 W m −2 , which is obviously implausible. Looking at the budget, this may be associated with an observed strong drying trend of -1 g kg −1 h −1 during that flight. However, looking at archived meteorological data of specific humidity from sites in Visalia, Mendota, and Porterville (NOAA ESRL https://www.psl.noaa. gov/data/obs/datadisplay/) shows a transient, regionwide drop of 1.5-2.0 g kg −1 between 1500 and 1600 LT, which may have biased this trend measurement from the flight. Removal of the negative result on the 6 August flight results in an average latent heat flux of 200 W m −2 , about 15% larger than the sensible heat flux, resulting in a Bowen ratio of 0.87. This value is consistent with other reports in semi-arid regions with irrigation (Todd et al. 2000). The average advection is -0.07 g kg −1 h −1 , which is again negligible compared to the entrainment flux (-0.84 g kg −1 h −1 ) and time rate of change (-0.61 g kg −1 h −1 ).

Model Comparison
Here, we infer the entrainment velocities from the WRF model output over the same polygon used to capture subsidence (see Sect. 2.5). This is done by tracking the subsidence at the ABL top as well as the ABL derivative terms from Eq. 10, allowing for comparison with the semi-independent aircraft-observed entrainment rates. The surface fluxes and Richardson numbers in WRF are again averaged over the aforementioned polygon over the flight periods (1200-1600 LT) to generate the results in Table 4. From our scalar budgeting, we found an average surface sensible heat flux of 173 W m −2 (± 21). In the WRF model runs, the surface heat fluxes are on average about 65% greater, with a mean of 286 W m −2 . In the model proposed by Tennekes (1973), the ABL heights scale with the square root of the integrated surface heat fluxes and thus should be approximately 30% greater in the WRF model than what is observed. However, the entrainment velocities, even when accounting for the increased heat flux by scaling by w * , are much larger in the WRF model (mean w e /w * = 3.51 × 10 -2 ) compared to the observations (mean w e /w * = 2.29 × 10 -2 ). Consequently, the ABL heights are about 60% larger in the WRF simulations, with a mean z i = 891 m compared to the observed average z i = 541 m. A more in-depth analysis of simulated boundary-layer heights and surface fluxes, including permutations of different land surface models and PBL schemes, is presented in Alexander et al. (2022). Despite these differences, the day-to-day variability appears to be somewhat captured with the WRF model  for both the ABL heights (R 2 = 0.31, p = 0.13, one-tailed) and the surface sensible heat fluxes (R 2 = 0.53, p = 0.05, one-tailed), although only the latter to statistical significance. A visualization of the simulated versus observed boundary-layer heights for July and August can be seen in Fig. 5. The day-to-day latent heat fluxes, both those obtained from the scalar budgeting method and those from the WRF output, correlated with California Irrigation Management Information System (CIMIS) reference evapotranspiration data from the same region (R 2 = 0.43 and 0.89, respectively, p = 0.11 and 0.002, one-tailed), with 6 August removed from the measurement set because of its un-physical result. The mean latent heat flux derived from the WRF model for these 5 days was 17% higher than the scalar budget mean result (234 W m −2 and 200 W m −2 , respectively). Additionally, the latent heat fluxes in WRF did appear to capture the day-to-day variability of the aircraft-observed latent heat fluxes from the scalar budget technique, again with the non-physical negative result removed (R 2 = 0.73, p = 0.03, one-tailed).
The average reference evapotranspiration (ET O ) reported by the CIMIS stations for the 5 afternoons, when multiplied by the latent heat of vaporization to obtain a latent heat flux, was 499 W m −2 . By design, CIMIS data are obtained over grass at well-watered locations. Thus, the regional latent heat flux is expected to be significantly lower than what the CIMIS reports because the valley as a whole comprises a variety of land types with different water availability and crop coefficients (Trousdell et al. 2016). From our data, we observe that the regional latent heat flux from the scalar budget analysis is approximately 40% of the CIMIS values. This is slightly lower than what is reported in Trousdell et al. (2016), where regional latent heat fluxes obtained from applying the scalar budget equation to airborne data in the San Joaquin Valley were found to be 55% of the CIMIS values.

Similarity Method Comparison
The values of w * and Q SH for each 300-s period of flight within the ABL, as outlined in Sect. 2.3, were computed, and flight-averages of these values were calculated (see Table 5). We compared the heat fluxes and convective velocity scales obtained by our scalar budget technique to those obtained by the similarity technique used in Sect. 2.3. The convective velocity scale resulting from the scalar budget is inferred by taking the cube root of the surface heat flux's multiplication with z i and g/θ v . A mean convective velocity and surface sensible heat flux of 1.45 ± 0.10 m s −1 and 220 ± 40 Wm −2 , respectively, are obtained by the similarity method. Compared to the values obtained by the scalar budget technique (1.37 ± 0.15 m s −1 and 170 ± 50 W m −2 ), the convective velocities are within each other's envelope of uncertainty and the heat flux is 29% larger. Despite the larger heat fluxes and convective velocities obtained by the similarity method, they are still lower than that of the WRF model (1.95 m s −1 and 286 W m −2 ). When comparing the scalar budget and similarity techniques, the day-to-day variability is somewhat captured in the w * values (R 2 = 0.47, p = 0.07, one-tailed) but not the sensible heat flux (R 2 = 0.13, p = 0.24, one-tailed). Additionally, the normalized mean bias is only 0.09 for the convective velocity when comparing the two methods. This comparison is encouraging, as both methods are largely independent, and both have significant sources of uncertainty.

Error Analysis
The error for each derivative term in our multi-linear regressions is a root-mean-square (RMS) error derived in the fit. The kinematic entrainment fluxes are the product of the entrainment velocity and a scalar delta term (i.e. − w e θ for the kinematic entrainment sensible heat flux and − w e q for the water vapour entrainment flux). The delta term error was assigned to be 1 K for potential temperature and 1.3 g kg −1 for water vapour, which were based on variations observed between many vertical profiles that penetrate the ABL top. This error term for the scalar jump across the boundary-layer top is purposefully conservative to account for the ill-defined nature of the idealized discontinuity across the inversion base and associated limitations in the zero-order model (Conzemius and Fedorovich 2006) required for observational studies of this type. The entrainment velocity contains (1) derivatives of ABL height, whose errors were previously mentioned, (2) subsidence obtained from the WRF model, for which we estimate an error of 0.5 cm s −1 based on the observed sensitivity of the vertical velocity values when random adjustments to the polygon were made, and (3) the mean horizontal wind at ABL height, which assigned an error of 0.3 m s −1 based on the measurement capabilities of the instrument (Conley et al. 2014) and the standard error of our 5-min level legs. The same error for horizontal winds near the ABL height applies to the ABL horizontal winds used in calculating the advection terms. The surface heat fluxes (both sensible and latent, for the potential temperature and water vapour budgets, respectively) are residual terms from Eq. 8, whose errors are calculated by combining the individual errors for all of the other terms in a standard error propagation.

Wind Shear, Wind Speed, and Entrainment Efficiency
A direct relationship has been expressed in the literature between the non-dimensional form of the entrainment velocity (w e /w * ) and a negative power (between 1 and 1.5, see Traumner et al. 2011) of the Richardson number computed across the entrainment zone, which we will investigate here. From Wyngaard (2010), a generalized gradient Richardson number for the atmosphere is: where u and v are the components of horizontal wind in any coordinate system. In this equation, the numerator quantifies the degree of stratification which suppresses turbulence (when positive), while the denominator quantifies the wind shear which is proportional to the production of shear-induced turbulence in the turbulence kinetic energy budget equation (Stull 1988 Eq. 5.1a) under the K-theory approximation that substitutes bulk gradients for fluxes. Thus, lower Richardson numbers are associated with a higher likelihood of turbulent mixing.
In modelling studies of ABL entrainment, such as Conzemius and Fedorovich (2006), researchers are able to employ horizontal homogeneity in their models, making the quantification of entrainment zone shear (∂U/∂z| EZ ) and thermal structure (∂θ/∂z| EZ ) relatively straightforward. In the real ABL, however, the boundary-layer height, entrainment zone thickness, wind velocity, and thermal stability all vary in space. Even our frequent (≈ 10 per flight) measurements of ABL height for an aircraft study of this type can only provide information for budgeting over a large region, not a high-resolution picture of the spatial variability of boundary-layer heights across the SJV (see Fig. 5). Thus, high quality measurements of shear and stability inside a spatially varying and relatively thin entrainment zone were not possible. Furthermore, measuring the entrainment zone thickness would have required a turbulence probe to measure the buoyancy fluxes at multiple altitudes, which was not present during this field campaign. As a result, several assumptions and simplifications are required to estimate a Richardson number-like quantity across our entrainment zone.
First, we estimate an average entrainment zone thickness based on a study of lidar observations showing that entrainment velocities of our magnitude (see Table 1) are typically associated with an entrainment zone depth of approximately 200 to 300 m (Traumner et al. 2011, Fig. 6). Due to issues arising from using small vertical bins for aircraft data, we determine the entrainment zone winds at the ABL top, U zi , as the vector mean wind speed within 300 m (150 m in each vertical direction) of the flight-averaged ABL top, which is calculated by averaging every 1-s data point within this vertical region for both the u and v components of wind. As a matter of interest, our observations of U zi are highly correlated with the jump in wind speed across the entrainment zone in the WRF model output (R 2 = 0.60, p = 0.04, Fig. 6 Observed correlation between entrainment velocity divided by the convective velocity scale (w e /w * ) and the modified inverse Richardson number (Eq. 12) one-tailed). Second, due to our shallow ABLs as well as the finding from Conzemius and Fedorovich (2006) that surface shear can indirectly affect entrainment zone shear by inducing drag, we assume that horizontal winds approach zero at the surface and estimate the average shear across the full depth of the ABL as a proxy for entrainment zone shear, i.e. ∂U/∂z| EZ U zi /z i . Third, Fedorovich et al. (2004) define a stratification parameter, simplified here as G = (N FT /N EZ ) 2 , where N FT and N EZ are the Brunt-Väisälä frequencies ([g/θ ∂θ/∂z] 1/2 ) of the lower free troposphere and entrainment zones, respectively. If the approximation is made that the buoyancy parameters of the free troposphere (or in our case, buffer layer) and entrainment zone are not significantly different, the stratification parameter simplifies to G = (∂θ/dz) FT /(∂θ/dz) EZ , which the authors find is approximately constant at 1.2. Thus, for our purposes here we approximate the lapse rate within the entrainment zone as (∂θ/∂z) EZ ≈G −1 γ ϑ , where γ ϑ is the lapse rate just above the ABL in the lowest 500 m of the buffer layer, which does not vary much (4.2 ± 0.4 K km −1 ). Combining these three items, we arrive at a modified inverse Richardson number for the entrainment zone that conforms to our observational capabilities: The convective velocity scale used to non-dimensionalize the entrainment velocity (w e /w * ) was obtained from our estimated surface sensible heat fluxes from the scalar budget technique (see Sect. 2.4). The rationale for scaling the entrainment velocity by the convective velocity is that the turbulence kinetic energy driving entrainment is primarily provided by the integrated buoyancy production (Driedonks 1982), as such a scaling reveals the local controls on entrainment with the surface heat flux dependence removed. As shown in Fig. 6, we observe a strong correlation between the scaled entrainment velocity and inverse Richardson number for the 6 flights (R 2 = 0.73, p = 0.02, one-tailed). This suggests that, in line with our stated hypothesis and previous research outlined in the introduction, increased wind shear at the top of the ABL results in a greater entrainment velocity and entrainment efficiency.
Earlier mixed-layer models parametrized boundary-layer growth by assuming that the downward virtual temperature flux at the ABL height was some fixed fraction of the upward surface sensible heat flux (Tennekes 1973;Dreidonks 1982;Culf 1992). This was based on the assumption originally stated by Ball (1960) that the sensible heat flux at the top of the ABL is driven principally by the surface heat flux (to be more precise, the buoyancy flux, but in the case of continental boundary layers the buoyancy flux is comprised mostly of the sensible heat flux). The (negative) ratio of the two fluxes called the entrainment efficiency (A R ) was concluded to be 0.2 with a reported range of 0.1-0.3 in the absence of mechanical Fig. 7 Relationships between the entrainment ratios (A R ) and mean boundary-layer wind speeds (at the top for this study) from a variety of studies. The mean boundary-layer wind speeds in Betts and Barr (1996) and Barr and Betts (1997) are approximated as near the centres of the reported bins (< 5, 5-10, and > 10 m s −1 , and < 4, 4-6, > 6 m s −1 , respectively) turbulence (Stull 1988). A field study conducted by Betts and Ball (1994) studied A R and classified the different days from their study as those with low ABL wind speeds (average 4.8 m s −1 ) and those with high wind speeds (average 11.3 m s −1 ) and found that with higher wind speeds, A R was greater. They found the shear in both categories to be relatively small across the ABL height (although they reported slightly stronger shear for the high wind days), and they suggest that the difference in A R is due to greater mechanical turbulence on the high wind days driving entrainment.
Here, we compute our entrainment efficiencies for each flight from the surface sensible and entrainment heat fluxes calculated from our scalar budgets, and our resultant values of A R range from 0.04 to 0.50 with an average of 0.23. Figure 7 compares these values with U Zi , which again supports the hypothesis that shear enhancement, related to stronger mean wind speeds in the entrainment zone, results in a greater entrainment efficiency (R 2 = 0.68, p = 0.02, one-tailed).
Results from other observational studies are included in Fig. 7 for comparison. Above the grasslands of Kansas using a budgeting technique from radiosondes, Betts and Barr (1996) reported an average A R of 0.39 (± 0.19) and pointed out an increase with mean wind speed. Davis et al. (1997) used airborne eddy-flux measurements of buoyancy fluxes in the ABL over the boreal forests of central Canada to find an average A R = 0.08 with a modest positive correlation to the jump in wind across the ABL interface. Then applying the same technique to the same Canadian boreal forests of the Davis et al. (1997) study, Barr and Betts (1997) reported an A R of 0.21 and further noted a stronger correspondence between entrainment ratio and mean wind speed than to that of wind shear across the entrainment zone. Flamant et al. (1997) used airborne lidar and turbulent flux profiles over the Mediterranean coastline and found a range of 0.15-0.30 in A R , which mainly depended on the jump in wind across the ABL top. Angevine (1999) formed budgets using wind profilers and radiosondes over the US Great Plains to obtain A R = 0.35 (U < 5 m s −1 ) and A R = 0.57 (U > 5 m s −1 ), finding no distinction when comparing values with different jumps in wind across ABL top. In a modelling study, Kim et al. (2003) found that A R ranged from 0.13 to 0.30, with geostrophic winds ranging from 5 to 15 m s −1 . Lastly, Conzemius and Fedorovich (2006) found in their LES study that shear at the entrainment zone is much more important in enhancing entrainment than surface shear. A comparison between A R and the average wind speed near the top of the boundary layer from our study is shown in Fig. 7, along with a comparison to results from the aforementioned similar studies that looked at the relationship between wind speed and A R . While the ranges of A R among the different studies appear comparable, there appears to be a much more dramatic dependence on wind speed in the low wind environment of the SJV.
The interpretation of the modified Richardson number presented here, as well as the relationship between entrainment zone wind speed and wind shear in the ABL, can be debated.
In theory, enhanced entrainment should be the result of wind shear at the ABL top as simulated by Conzemius and Fedorovich (2006). At least in the zero-order model being used here, this theory has the advantage of being physically intuitive and having measurable quantities in our aircraft dataset. However, the fact that the previous studies, as well as the current one, found clearer relationships with mean wind than with observed wind shear could support the idea that the direct measurement of wind shear is generally challenging and ultimately too noisy to show a reliable correlation. Alternatively, it is possible that the consensus of LES models showing the ABL entrainment dependency on entrainment zone wind shear is less straightforward in the real atmosphere, where the coupling (or lack thereof) between surface and entrainment shear is more complex than in idealized simulations. We believe that the large variability in entrainment ratios observed here and in past studies should serve as a call to the community to study its dependence on wind (and shear) further in canonically convective boundary layers.
As a matter of interest, some researchers have speculated an inverse relationship between A R and the Bowen ratio (Angevine 1999;LeMone et al. 2019). While the theoretical significance of such a relationship is not clear given that both parameters are mathematically dependent on the surface sensible heat flux, we do find an inverse relationship between A R and the Bowen Ratio when removing the non-physical result of the 6 August flight (R 2 = 0.67, p = 0.05, one-tailed).

Synoptic Environment and Buffer Layer
To understand what might control the valley winds, and thus entrainment variability above the ABL from the synoptic perspective, we looked to the North American Regional Reanalysis (NARR) dataset (Kalnay et al. 1996) for the July flight dates and the August flight dates separately, with the 21 UTC (1300 LT) frame from each flight date averaged together. For the 850 hPa level and below, we instead analyse the higher resolution WRF model because we expect the pressure fields at these altitudes to be influenced heavily by the complex terrain. The 500 hPa NARR analysis is shown in Fig. 8, and Fig. 9 shows vertical profiles of observations from the aircraft, radiosonde, and radio acoustic sounding systems (RASS).
Winds up to 750 m a.g.l. are consistently north-westerly (aligned with the Central Valley axis) for all 6 days ( Fig. 9), which is typical of summertime conditions (Zhong et al. 2004). For the three July flights, enhanced north-westerly flows are observed between 750 m a.g.l. and veering to northerly up to 2000 m, while for the August flights, the up-valley winds weaken above 750 m (in both the scalar and vector mean sense) virtually stagnating near 1500 m and developing south-south-westerly flows aloft. Overall, there are weaker winds aloft and stronger valley winds in the warm anomaly of July, while in the cooler conditions of early August there are stronger synoptic winds aloft and weaker valley winds.
The broader synoptic scale in the July NARR data shows a large upper air high-pressure system at the 500 hPa level above Southern Nevada (Fig. 8) with an accompanying thermal low-pressure centre and higher air temperatures near the surface (Fig. 9) about 100 km to

A B
Fig. 8 500 hPa pressure surface height contours (m) and wind vectors from the North American Regional Reanalysis (NARR) for a July flights and b August flights. Each plot shows averages of the three snapshots at 1300 LT from each afternoon. Plots and data are provided by North American Regional Reanalysis (NARR) at NOAA/ Physical Sciences Laboratory the south-west of the upper-level high near the NV-CA border (data not shown). For the July flights, the orientation of the isobars indicates an enhanced up-valley component to the 850 hPa geostrophic winds (not shown), which are above the ABL and consistent with the observed profile above Visalia (Fig. 9). In the August flights, the 850 hPa geostrophic winds are cross-valley, and we see that the buffer-layer winds at VMA are significantly weaker than those of the July flights. This differing pattern of synoptic forcing may be what drives the stronger entrainment zone winds in July (mean = 2.40 m s −1 ) versus August (mean = 1.37 m s −1 ), as the synoptic forcing in July is aligned with the mountain-valley circulation driven up-valley winds. This distinction manifests as well in the inverse Richardson numbers from the July flights (mean = 1.00) and the August flights (mean = 0.34) seen in Fig. 6. This demonstrates a possible connection between entrainment strength above the SJV and the larger synoptic patterns, although far more observations would be needed to establish any causal relationships. The 500 hPa ridge pattern in July results from a lower-level thermal low that induces a geostrophic component that augments the background mountain-valley circulation. Zhong et al. (2004) show that the daytime winds above the SJV up to about 1500 m during summer are predominately out of the north-west and tend to back to westerlies at the upper levels especially in the afternoon as the integrated valley heating progresses (their Fig. 6). Furthermore, Faloona et al. (2020) show that above about 3 km, the climatological winds back further to south-westerly (their Figs. 3 and 9). Therefore, when a warm anomaly is present around the Nevada/California border near the surface, the lower-level depression will reinforce the usual up-valley flow strengthening the SJV mountain-valley circulation. Beaver and Palazoglu (2009) performed an analysis of synoptic patterns and ozone levels in SJV, and we note here that the 500 hPa synoptic setup for our July flights is similar to their "onshore high" pattern, and for our August flights it is similar to their "onshore high/trough" pattern. In the Beaver and Palazoglu (2009) study, high ozone was seen to be more likely for the onshore high set-up, and this is consistent with observations that show roughly 10 ppb greater ozone levels in the July CABOTS flights . It is very well A B Fig. 9 Average vertical profiles of wind vectors from the Radio Acoustic Sounding System (RASS) at the Visalia Municipal Airport embedded in a valley cross-cut showing topography. a 27-29 July, and b 4-6 August 2016. The colour of the windbarbs show the mean wind magnitude (i.e. average wind speed), and the vector mean is also shown. Averaged profiles of ozone from the flights are in blue and averaged profiles of potential temperature in red established that near surface ozone levels are strongly related to air temperature in polluted regions due to the fact that biogenic VOC emissions and many photochemical rates (e.g. the dissociation of peroxyacetyl nitrate) tend to be temperature dependent leading to higher ozone on warmer days. Relationships between daytime ozone and temperature reported in the literature for the SJV tend to fall between 1.5 and 2.0 ppb/K (Steiner et al. 2010;Pusede et al. 2014;Pan and Faloona 2022). Therefore, we note that the higher ozone during the warm anomaly in July is consistent with these temperature dependent enhancements, and in fact, the buffer layer ozone in both periods is not dramatically different (see Fig. 9). Thus, we find that the stronger entrainment induced by warmer temperatures serves as a countervailing mechanism on the near surface ozone levels, because lower ozone concentrations in the buffer-layer are being mixed into the boundary layer at a faster rate.

Conclusions
Our study successfully employed the scalar budget methodology to measure sensible and latent heat fluxes in the San Joaquin Valley during the summertime ozone season. The Bowen ratios, entrainment velocities, and boundary-layer heights are consistent with other literature reports in heavily irrigated semi-arid regions. Presently, the WRF model does not seem to accurately predict the ABL depths in the SJV, as shown in both our study and that of Alexander et al. (2022). The shallower than expected summertime boundary-layer heights in the valley can likely be accounted for by the robust mountain-valley induced subsidence, as well as the heavily irrigated land surface. While this led to large discrepancies between simulated and observed fluxes, the overall day-to-day correlations were encouraging. There is nonetheless a limitation in the correlation analyses with only six data points, but they are worth mentioning due to the large variability in day-to-day observations of this study. The scalar budgeting technique performed here can be used for future studies to evaluate and refine models in the SJV and other mountain-valley systems. Additionally, we have outlined a novel and relatively low-cost measurement technique for turbulence using aircraft wind and thermodynamic data, which can supplement the scalar budget technique to even better characterize the boundary layer environment without the need of costly gust probes.
In addition, our findings suggest that stronger than average winds above the ABL, as a consequence of differing synoptic patterns during warm spells, can enhance wind shear across the entrainment zone leading to more vertical entrainment mixing of ABL air. In particular, north-westerly geostrophic winds above the ABL are associated with low-level depressions and higher surface temperatures to the east of the SJV. The stronger entrainment associated with warmer temperatures can ventilate the ABL, in principle, counteracting the accelerated photochemical production of ozone. In this study, we aimed to provide a step towards linking the stories of synoptic and boundary layer meteorology in the context of the summertime polluted ABL of California's Central Valley. A better characterization of these interactions between different scales may help air quality forecasters in understanding the meteorological impacts on air pollution episodes.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.