Estimations of the Mexicali Valley (Mexico) Mixing Height

: We report estimations of the Mexicali Valley (Mexico) mixing height for three seasons. Surface and upper air meteorological measurements were carried out nearby Cerro Prieto geothermal power plant during July 2010 (summer), January 2012 (winter), and October 2016 (autumn). Four di ﬀ erent methods were applied to estimate the convective boundary layer (CBL) height from radiosonde (RS) proﬁles: the parcel method, the gradients method, the least-squares variational approach based on the slab model of the CBL structure, and a covariance method. For nocturnal conditions, we used diagnostic models based on friction velocity and Monin–Obukhov length. Under unstable conditions, we obtained (on average) mixing heights of 497 m at 06:00 LST, 1242 m at 12:00 LST, 1404 m at 15:00 LST, and 482 at 18:00 LST during summer; 754 m at 12:00 LST during winter; and 1195 m at 12:00 LST and 13:00 m between the 14:00 and 15:00 LST during the autumn. The results allowed adjusting a semiempirical model to evaluate mixing height from turbulent sensible heat ﬂux and friction velocity data. Our results provide practical tools that could facilitate the application of regulatory dispersion models to assess air quality in the region.


Introduction
Pollution in the Mexicali Valley (Mexico) is not a simple phenomenon. It involves factors that play essential roles in the region, both in the economic development and the environment, such as landfills, cement and smelting industries, agricultural burning, and power plants nearby [1][2][3][4]. The Cerro Prieto geothermal power plant, which is the primary source of energy in the state of Baja California (BC, Mexico), is frequently considered as the principal source of pollution in the region. However, the environmental impact of the geothermal plant is controversial, with divergent opinions from the Comisión Federal de Electricidad (CFE), scholars, and social sectors [5][6][7][8].
Air quality impact assessment is an essential activity for determining the relative contribution to ground-level pollutant concentrations of specific current or future source emissions at receptor sites. Its main tools are air quality monitoring and modeling techniques. Air pollution models also play a pertinent role in risk analysis, emergency planning, and source apportionment studies.

Materials and Methods
We performed this study with the support of surface and upper-air meteorological data obtained from short-term experimental campaigns carried out nearby of Cerro Prieto power plant in summer, autumn, and winter. In the first part of this section, we described the experimental campaigns and data obtained, whereas in the second part, we presented the methods we used to estimate the boundary layer heights.

Meteorological Campaigns
We carried out three short-term experimental campaigns in an open rural site located at 32°24.0015' N and 115°14.0065' W, with an altitude of 10 meters (above sea level), within the facilities of the Cerro Prieto geothermal power plant, 700 m to the west of the geometrical center of the emission sources distribution, approximately. We made surface meteorological measurements throughout July 2010 (summer), January 2012 (winter), and October 2016 (autumn), and we launched atmospheric radiosondes during the periods of 17-28 July 2010, 25-29 January 2012, and 11-15 October 2016.

Surface Meteorological Measurements
We used a 3D ultrasonic anemometer (USA-1, METEK GmbH, Elmshorn) to measure the wind velocity components and temperature at z = 12 m above ground with a sampling rate of 10 Hz. No significant obstacles were around the monitoring site. From this high-frequency data, following the methods described in [14], we calculated 1-h averages of the turbulent parameters, such as friction velocity ( * u ), sensible heat flux (H0), Monin-Obukhov length (L), and turbulent kinetic energy (ε), among others. Using 1-min sets of data, we carried out two rotations on the ordinary velocity components (vx, vy, vz) to obtain the mean streamwise velocity   , 0, 0 u and the turbulent fluctuations Following the climatology [12,13], the mean seasonal temperatures in the region are 13 • C in winter, 19 • C in spring, 28 • C in summer, and 19 • C in autumn. The windier part of the year extends from March to June, with average wind speeds of more than 3.3 m/s, and the calmer period lasts from July to February, with average wind speeds around 2.8 m/s. The cloudy period in Mexicali extends from November to March. The day length varies significantly over the year. On average, the sunrise occurred at 06:30 LST (sunset: 17:15 LST) in winter, at 05:00 LST (sunset: 18:25 LST) in spring, at 05:00 LST (sunset: 18:30 LST) in summer, and at 06:00 LST (sunset: 17:00 LST) in autumn. The average daily incident shortwave solar energy experiences extreme seasonal variation over the year. The brighter period of the year extends from middle April to the end of July, and the darker period runs from beginning November to beginning February. The Mexicali local standard time (LST) is UTC-8 (Coordinated Universal Time 8).

Materials and Methods
We performed this study with the support of surface and upper-air meteorological data obtained from short-term experimental campaigns carried out nearby of Cerro Prieto power plant in summer, autumn, and winter. In the first part of this section, we described the experimental campaigns and data obtained, whereas in the second part, we presented the methods we used to estimate the boundary layer heights.

Meteorological Campaigns
We carried out three short-term experimental campaigns in an open rural site located at 32 • 24.

. Surface Meteorological Measurements
We used a 3D ultrasonic anemometer (USA-1, METEK GmbH, Elmshorn) to measure the wind velocity components and temperature at z = 12 m above ground with a sampling rate of 10 Hz. No significant obstacles were around the monitoring site. From this high-frequency data, following the methods described in [14], we calculated 1-h averages of the turbulent parameters, such as friction velocity (u * ), sensible heat flux (H 0 ), Monin-Obukhov length (L), and turbulent kinetic energy (ε), among others. Using 1-min sets of data, we carried out two rotations on the ordinary velocity components (v x , v y , v z ) to obtain the mean streamwise velocity (u, 0, 0) and the turbulent fluctuations ( u, v, w) and T of the velocity components and (sonic) temperature (T). Then, we used the following equations to estimate the turbulent parameters from the covariance between the turbulent fluctuations [14]: u * = u w 2 + v w 2 1 4 (1) Here, ρ and c p are the density and specific heat of the air, k is the von Karman constant, and g is the gravity acceleration.
We calculated the wind speed (WSP) and wind direction (WDR) from the wind velocity components. Using conventional instruments, we also measured atmospheric pressure (p), temperature (T), and relative humidity (RH) with a 1 Hz sampling rate. We calculated the 1-h averages of all variables.

Atmospheric Soundings
From the atmospheric soundings, we obtained the vertical profiles of pressure, temperature, relative humidity, and dew point temperature. From the radiosonde (RS) data profiles, we calculated the virtual potential temperature (VPT) and specific humidity (SH) profiles. We used the Vaisala's MW21 DigiCORA III upper air sounding system with radiosondes RS92-SGP [15], for the sounding campaigns of July 2010 and January 2012 and a GRAW ground station GS-B sounding system with radiosondes DFM09 [16], for the October 2016 campaign. The numbers of launched radiosondes were 25 in the period of 17-28 July 2010; 13, during 25-29 January 2012; and 14 during 11-15 October 2016. These launchings were performed mainly at the sunrise, noon, afternoon, sunset, and midnight.

Estimation Methods
We provided in this section a brief description of the methods we used to estimate the atmospheric boundary layer (ABL) depth for both the diurnal (convective) and nocturnal (stable) conditions from the RS profiles and the surface meteorological measurements.
We used the surface measurements of Monin-Obukhov length and friction velocity to estimate the stability parameters ζ = z/L and µ = u * / f L, where z is the height where L is measured and f the Coriolis parameter. We also estimated the near-surface temperature lapse rate of the RS temperature profiles. With these three parameters, we organized the RS according to the atmospheric stability. Near the transition hours (sunrise and sunset), we used only the near-surface temperature lapse rates of the T-profiles for reckoning the stability conditions of the soundings. Qualitatively, the near-surface atmospheric conditions were labeled as unstable, neutral, and stable conditions when ζ < 0, ζ = 0, and ζ > 0, respectively. In terms of the other stability parameter, following the proposal of Arya [17], we labeled the atmospheric conditions according to Table 2. Table 2. Stability classes, as suggested by Arya [17]. Sonic at z = 12 m. To avoid mistaking free tropospheric features for the top of the ABL, we restricted the RS profiles to the lowest 4000 meters above ground level (magl). We smoothed all the RS profiles to reduce the effect of extreme values (outliers). Specifically, if P = (z k , T k )|k = 1, 2, . . . , N is the set of data points for the temperature profile obtained by RS, the smoothed value at a given z k was the sum of the values of temperature at z k-1 multiplied by 0.25, at the z k of interest multiplied by 0.5, and at z k+1 multiplied by 0.25:

Convective Boundary Layer
We applied four different methods to estimate the convective boundary layer (CBL) height from the radiosonde data. Two of them were traditional approaches often encountered in the ABL literature: the parcel method (PM) [10,18,19], which requires only the temperature profiles, and the gradient method (GM), which is applied to the VPT [20][21][22] and SH [23] vertical profiles. The other two methods were the least-squares variational approach (LSVA) [24] that carries out the best fitting of slab model profiles to the VPT and SH vertical profiles and a covariance method (CM) that identifies the mixing height with the level of the minimum covariance between the VPT and SH profiles.

The Parcel Method
The PM is a procedure widely used to estimate the mixing height from RS data, with no need for wind profile. The original method (Holzworth) determines the maximum daily value of mixing height from the intersection between the dry adiabatic profile, which starts at the maximum surface temperature, and the RS temperature profile observed at sunrise (generally, from 05:00 through 07:00 LST, approximately) [10,18,19]. It is ordinary an interpretation of the estimated height as the equilibrium level of a hypothetical rising parcel of air representing a thermal. This method strongly depends on the surface temperature, and high uncertainty in the estimated height may result in situations without a pronounced inversion at the top of the CBL. A fixed meteorological station measures the surface temperature.
In this work, we did not attempt estimating the maximum daily value of mixing height [18], which sometimes finds application in air quality modeling, or the minimum daily value [19], which provides a worst-case scenario for estimating surface concentrations of air pollutant emissions. Instead, we used the parcel method as described by Seidel et al. [25], and we evaluated the mixing Atmosphere 2020, 11, 505 6 of 26 height at the time of the atmospheric sounding. No intersection between the dry adiabatic profile and the RS temperature profile will exist if the lower part of the T-profile exhibits a subadiabatic lapse rate; therefore, no direct estimation of mixing height can be provided by the PM in this case.

Gradient Methods
Several methods are available to determine the convective mixing height from the gradient profiles of the RS data. In this work, for estimating the mixing height under unstable conditions, we applied gradient methods to the vertical profiles of virtual potential temperature and specific humidity. The identification procedures are as follows: (a) Virtual potential temperature: The CBL height is estimated to be the level of the maximum vertical gradient of virtual potential temperature, which is characteristic of a transition from a convectively less stable region below to a more stable region above. [20,23,26].
(b) Specific humidity: The CBL height is estimated to be the height where the minimum of the vertical gradient of the specific humidity profile occurs [23].
The gradient approaches are most effective when the CBL is capped with a well-defined inversion layer with the VPT increasing and the SH decreasing sharply above the top. However, the local gradient approach could lead to an ambiguous height determination for scenarios where a sharp inversion layer is absent, or where a cloud layer tops the ABL, or where there are multiple layers in the troposphere with strong gradients.

Slab Model of CBL and the Least-Squares Fitting Approach to the Mixing Height
The basic structure of the CBL comprises a few stacked atmospheric layers: the surface layer, where the vertical profiles of the principal meteorological quantities follow the similarity laws; the mixing layer, which is in the middle, characterized by intense turbulence of thermal origin (convection); and the entrainment layer, marked by a temperature inversion that suppresses the turbulence and buffers the free atmosphere above [20][21][22]. This layered structure of CBL is frequently evident in the vertical profiles of VPT and SH. This feature especially holds for VPT, but with minor restrictions also for SH [27].
The slab model assumes that, under convective conditions (z/L < 0), the entrainment layer has an infinitesimal thickness, and that the vertical profiles of VPT and SH display simple shape profiles ( Figure 2), which we can use to estimate the convective mixing height [20].
Atmosphere 2019, 9, x FOR PEER REVIEW 6 of 27 (a) Virtual potential temperature: The CBL height is estimated to be the level of the maximum vertical gradient of virtual potential temperature, which is characteristic of a transition from a convectively less stable region below to a more stable region above. [20,23,26].
(b) Specific humidity: The CBL height is estimated to be the height where the minimum of the vertical gradient of the specific humidity profile occurs [23].
The gradient approaches are most effective when the CBL is capped with a well-defined inversion layer with the VPT increasing and the SH decreasing sharply above the top. However, the local gradient approach could lead to an ambiguous height determination for scenarios where a sharp inversion layer is absent, or where a cloud layer tops the ABL, or where there are multiple layers in the troposphere with strong gradients.

Slab Model of CBL and the Least-Squares Fitting Approach to the Mixing Height
The basic structure of the CBL comprises a few stacked atmospheric layers: the surface layer, where the vertical profiles of the principal meteorological quantities follow the similarity laws; the mixing layer, which is in the middle, characterized by intense turbulence of thermal origin (convection); and the entrainment layer, marked by a temperature inversion that suppresses the turbulence and buffers the free atmosphere above [20][21][22]. This layered structure of CBL is frequently evident in the vertical profiles of VPT and SH. This feature especially holds for VPT, but with minor restrictions also for SH [27].
The slab model assumes that, under convective conditions (z/L < 0), the entrainment layer has an infinitesimal thickness, and that the vertical profiles of VPT and SH display simple shape profiles ( Figure 2), which we can use to estimate the convective mixing height [20]. We applied the least-squares variational approach (LSVA) proposed in [24] to estimate the mixing height from the VPT profile under convective conditions. We also used an extension of this method with the SH profile. The basis of LSVA is on the slab model of the CBL structure [20], but with the additional assumption that the thickness hs of the surface layer is infinitesimal (negligible).
Radiosonde ascents provide the atmospheric vertical profiles of temperature (T), pressure (p), and dew point temperature (Tdew), and these basic profiles allow us to calculate the specific humidity q and the virtual potential temperature  using standard relations [20,22].
the set of data points obtained by RS for the VPT profile ( ). z  Then, according to LSVA, if no residual layer is present, the profile ( ) z  exhibits only three sections [20,24]. The mixing layer, where the VPT is assumed to be constant and equal to  from 0 z  (earth We applied the least-squares variational approach (LSVA) proposed in [24] to estimate the mixing height from the VPT profile under convective conditions. We also used an extension of this method with the SH profile. The basis of LSVA is on the slab model of the CBL structure [20], but with the additional assumption that the thickness h s of the surface layer is infinitesimal (negligible).
Radiosonde ascents provide the atmospheric vertical profiles of temperature (T), pressure (p), and dew point temperature (T dew ), and these basic profiles allow us to calculate the specific humidity q and the virtual potential temperature θ using standard relations [20,22].
Let P = (z i , θ i )|i = 1, 2, . . . , N the set of data points obtained by RS for the VPT profile θ(z). Then, according to LSVA, if no residual layer is present, the profile θ(z) exhibits only three sections [20,24]. The mixing layer, where the VPT is assumed to be constant and equal to θ m from z = 0 (earth surface) to z = h m (mixing height); the entrainment zone, whose thickness is assumed infinitesimal, located at z = h m , where the VPT has a discontinuity, jumping from θ = θ m to θ = θ 0 ; and the free atmosphere, for z > h m , where VPT increases linearly with a slope γ.
The following equation describes the model profile (solid line in Figure 3): Here, H(x) is the Heaviside function, which is 0 for negative argument values and one otherwise.
Atmosphere 2019, 9, x FOR PEER REVIEW 7 of 27 Here, H(x) is the Heaviside function, which is 0 for negative argument values and one otherwise.
under the next five restrictions, define the best choice of the VPT model profile. The limit b of the integral in Equation (7) is a cutting height subjectively imposed on the balloon ascent. It is an adjusting parameter worked (by a skilled user) to select the portion of the virtual potential temperature profile free of residual layers. In principle, the following set of nonlinear equations provide the solution of the problem [24]: However, due to the experimental nature of the VTP profile ( ) z  , we used a numerical From the standpoint of the LSVA, the estimating procedure of mixing height consists of finding the values of the parameters h m , θ m , θ 0 , γ , which define the best fitting of the model profile to the set P of the data points obtained from RS for the VPT profile θ(z). It is a parametric estimation problem that the LSVA solves through nonlinear least-squares [24].
The values of h m , θ m , θ 0 , γ that minimize the functional under the next five restrictions, define the best choice of the VPT model profile. The limit b of the integral in Equation (7) is a cutting height subjectively imposed on the balloon ascent. It is an adjusting parameter worked (by a skilled user) to select the portion of the virtual potential temperature profile free of residual layers. In principle, the following set of nonlinear equations provide the solution of the problem [24]: Atmosphere 2020, 11, 505 8 of 26 However, due to the experimental nature of the VTP profile θ(z), we used a numerical procedure to obtain the minimizing set h m , θ m , θ 0 , γ . With the Equations (9)-(12), we reduced the functional (7) to another one that depends only on the parameters h m and b. Given b, we varied h m over the interval 0 < h m < b until obtaining the minimum of F. The minimizing set h m , θ m , θ 0 , γ defines the best fitting of the model VTP profile to the experimental data. For a given b, the minimizing h m represents the best choice for the mixing height.
This approach takes into account all the VTP profile up to the cutting height b. Then, if the VPT profile shows the existence of multiple layers in the troposphere with strong gradients, different values of b can result in different estimations of mixing height. In this work, we also used a direct extension of the LSVA to estimate the mixing height from the SH profile.

Covariance Method
One of the main characteristics of the mixing layer is that an intense vertical mixing tends to leave the conserved properties such as virtual potential temperature and specific humidity nearly constant with height [27]. Otherwise, in the entrainment layer (the transition zone between the mixing layer and the stably stratified, quasi-nonturbulent free atmosphere above), the typical features are a significant positive temperature lapse rate, a sharp decrease of specific humidity, and sometimes considerable vertical wind shear [27]. The sharp gradients of VPT and SH (positive, the first and negative, the second) in the entrainment layer compel to a joint covariance of VPT and SH with an ample negative value in the capping inversion of the CBL. The level it occurs identifies the mixing height. We refer to this approach as the covariance method (CM).

Stable Boundary Layer
In contrast to the daytime convective boundary layer, under stable nocturnal conditions, the situation is more complicated and less straightforward to determine the height of the stable boundary layer (SBL). Turbulence is gradually suppressed during stable conditions and may show an intermittent behavior [28]. Some other processes, such as longwave radiation flux divergence [29], gravity waves, wave breaking [30], baroclinicity [31], and low-level jet [32], may affect the structure of the nocturnal boundary layer (NBL). Therefore, no straightforward interpretation exists of temperature and wind speed profiles, and both the measurement and prediction of the nighttime mixing height have proven to be far more complicated and less reliable. When direct measurements (using instruments like SODAR, RASS, and LIDAR) of mixing height are not available, as it is the case for most routine meteorological observations, one estimates the height of the SBL through parametrical relations that express it as a function of some boundary layer parameters [33,34].

Parametrizations of the Turbulent SBL
Many parameterizations for the height of the turbulent SBL have been suggested in the literature. Both diagnostic and prognostic relationships have been proposed, but no agreement exists on which type is the most suitable [10,27]. Several diagnostic models for estimating the height of the turbulent SBL from turbulence parameters such as friction velocity and Monin-Obukhov length have been proposed [17,[35][36][37][38][39][40][41][42]. The main models to estimate the SBL height were discussed and summarized by Seibert et al. [27], and Lena and Desiato [43] performed an intercomparison of several indirect appraisal methods to estimate of the stable mixing height for urban air pollution modeling. They also compared the results of the models against mixing heights derived from wind (SODAR) and temperature (RASS) profiles measured in the Milan urban area during spring and summer 1996.
In Table 3, we presented three of the most popular diagnostic equations (based on scaling arguments) [10]. These equations correspond to the classical parametrization models of Zilitinkevich [35], Nieuwstadt [38], and Arya [17]. Most of the verification studies done in the past do not seem to favor the application of more elaborated parameterizations [10]. Table 3. Stable mixing height diagnostic algorithms. The equations give stable mixing height (h m ) as a function of friction velocity (u * ) Monin-Obukhov length (L), and the Coriolis parameter (f ).

Model
Authors ID [38,42] N1981 Zilitinkevich [35]. Arya [17]: a = 0.74. Mahrt, Andre, and Heald [39]: a = 0.6. Nieuwstadt [41]: The Nieuwstadt parametrization, N1981, was derived for the height of the stationary boundary layer during stable lapse rate conditions. It satisfies the conventional limits for neutral conditions and large values of stability [38]. The parametrization Z1972 was first suggested by Zilitinkevich in 1972 [35] as a model of the nocturnal atmospheric boundary layer, and several authors confirmed it from 1974 to 1978 [44][45][46][47]. The parametrization A1981 was suggested by Benkley and Schulman [36] as a simple relation for routine operational use. Both simple relations were studied by Arya [17], who found functional correlations of h m with (Lu * / f ) 1 2 and (u * / f ). Arya concluded that Z1972 is adequate for moderately stable conditions, whereas A1981, is satisfactory for all but extremely stable (µ > 100) nighttime conditions. Both diagnostic models, Z1972 and A1981, can be obtained from the N1981 as asymptotic limits of h m /L 1 (very stable conditions) and h m /L → 0 (neutral conditions), respectively [27]. A1981 provides the simplest parametrization for the nocturnal mixing height, and it is suggested frequently for practical applications, but with the prescription of a convenient lower threshold value [43].
Venkatram [37] and Nieuwstadt [41] proposed other simpler diagnostic algorithms where the SBL height is proportional to the 3/2-power of friction velocity and the 10 m-wind-speed. The US-EPA processors AERMET [48] and CALMET [49] use these relations; however, they show underestimation or overestimation in low and high wind conditions, respectively [43].
Other diagnostic models have been proposed using the time scale 1/N BV , defined by the Brunt-Vaisala frequency, instead of 1/f or where the diagnostic formulae considering the bulk structure of the SBL are based on the assumption that turbulence production must vanish at the SBL top and the Richardson number must, therefore, exceed its critical value [10]. However, these and other available methods (see [10] for a review) involve information not available from radiosondes, such as cloud base height, micrometeorological covariance statistics, or turbulence flux profiles.
In this work, we used Equations (N1981, Z1972, A1981) to estimate the nocturnal mixing height from surface measurements of friction velocity and Monin-Obukhov length. In practice, we considered Z1972 and A1981 with the proportionality constants equal to the average of the respective constants recommended by the different authors. In the case of the monitoring site in the Mexicali Valley (with latitude around 32 • 24' N), the Coriolis parameter is f = 7.79 × 10 −5 s −1 . Then, Z1972 became to and Equation A1981 became to These equations, which will be referred to as Z1972S and A1981S, give h m in meters if u * is in m/s, and L is in meters. Friction velocity in A1981S is limited to 0.05 < u * < 1.2 m/s [50].

Surface-Based Inversions
Longwave radiation flux divergence may affect the structure of the NBL [29]. At nighttime, under clear skies, the net radiative-flux divergence is a potentially crucial cooling phenomenon that affects the boundary layer structure. When the sun goes down, the ground loses heat very quickly, and this cools the air above the ground. This cooling of the ground continues over more extended periods during the wintertime nights, which are much longer than the summertime nights. Weak winds prevent air mixing near the surface, and clear skies increase the rate of cooling at the Earth's surface. Stable conditions inhibit vertical and horizontal mixing near the ground and, consequently, favor the development of a strong surface-based temperature inversion, where the near-surface temperature lapse rate is positive. A surface-based T inversion is a clear indicator of a stable boundary layer, whose top can define an ABL height. If a surface-based inversion (SBI) is found in a sounding, the previous methods are not evaluated, as they assume a different NBL structure [25]. The top of the SBI can be defined by the height above which temperature first decreased with elevation [25,51].

Estimation of Mixing Height from the Sensible Heat Flux
A straightforward semiempirical model was proposed in 2001 by Yi, Davis, Berger, and Bakwin (YDBB) [52] to estimate the convective mixing height as a linear function of the square root of the energy accumulation due to the turbulent sensible heat flux H 0 : Here, a and b are constants to be empirically determined by fitting to measurements. The YDBB model can be used to predict mean mixing height during that period of the day when the mixing layer grows due to surface heating.
The lower limit t 0 in the integral of Equation (15) denotes the instance when H 0 changes sign from negative to positive (at the sunrise). When t = t 0 , we have h m (t 0 ) = a. Then, the constant a represents the value of mixing height at the sunrise. The upper limit t is the time of interest, for which, we want to evaluate the mixing height, t 0 < t ≤ t m + 2 h. t m is the time when the sensible heat flux reaches its maximum value [52]. From the standpoint of YDBB, the mixing height gets the most significant value when t ≈ t m + 2 h, approximately. However, in [52], the authors underlined that it might occur (in northern Wisconsin) that the weather is generally transparent or partly cloudy on days when maximum mixing height is reached 4 h after the maximum H 0 . Moreover, with clear skies in the summertime, the mixing layer can continue to grow for 5 or 6 h after the maximum H 0 .

The Time Evolution of the Mixing Height
The time behavior of mixing height follows a diurnal cycle. It also depends on the meteorology parameters, the turbulent surface fluxes, and the physiographic characteristics of the site. The diurnal cycle of heating of the earth's surface by the direct influence of the insolation creates thermal turbulence and causes the recurrent mixing layer growth with a period of 24 h. However, the underlying diurnal behavior of mixing height is known not to be purely harmonic (i.e., it does not occur in a simple sinusoidal way) [20]. For convective conditions, the YDBB model can be used to obtain an estimation of the temporal evolution of the mixing height during its growing phase from the surface values of turbulent sensible heat flux. For nighttime conditions, on the other hand, estimations of mixing height can be obtained from the surface values of friction velocity and Monin-Obukhov length using one of the models presented in Table 3.
Then, we can use the YDBB and the stable parametrizations, complementarily, to establish a procedure that allows the hourly estimations of the diurnal mixing heights from the surface values of sensible heat flux, and nocturnal mixing heights from friction velocity and Monin-Obukhov length: Here, π(t) is the rectangle function, which is equal to 1 if t 0 < t ≤ t m + 2 h and 0 otherwise. The function h s is one convenient parametrization (N1972, Z1981, A1981, or any other) that allows estimating the nocturnal value of h m . Here, π(t) is the rectangle function, which is equal to 1 if t0 < t ≤ tm + 2 h and 0 otherwise. The function hs is one convenient parametrization (N1972, Z1981, A1981, or any other) that allows estimating the nocturnal value of hm.

Surface Meteorological Data
Figures 4 and 5 show the hourly trends of the minimal, mean, and maximal values of temperature, relative humidity, wind speed, sensible heat flux, friction velocity, and the stability parameters ( = ⁄ and = * ⁄ ). For each sounding campaign (17-28 July 2010, 25-29 January 2012, and 11-15 October 2016), we determined hour by hour the minimum, mean, and maximum values of the surface meteorological variables measured during its period. Figure 5 also shows the hourly trends of sensible heat flux and friction velocity for the corresponding month.
In Figure 4, we observed that the larger values of temperature and relative humidity prevailed in summer, while in winter, the smaller ones. On average, the minimum, mean and maximum temperatures were 30, 33, and 35 °C in summer, 16, 18, and 19 °C in winter, and 26, 27, and 29 °C in autumn. For relative humidity, the values were 39%, 48%, and 55% in summer, 17%, 30%, and 44% in winter, and 24%, 38%, and 58% in autumn. These results are consistent with the Mexicali Valley climatology [13]. For wind speed, the averages of the minimum, mean, and maximum values were 2.2, 3.9, and 5.5 m/s for summer, 1.7, 4.3, and 6.8 m/s for winter, and 1.3, 2.1, and 2.9 m/s for autumn, with the larger values in winter and the smaller ones in autumn. Regarding the Mexicali climatology [13], the mean wind speeds obtained for the summer and autumn campaigns fell within the 25th to 75th percentile band, but for the winter campaign fell in the 10th to 90th percentile band.  winter, and 06:00-16:30 LST in autumn. During these periods, the surface heating produced thermal turbulence (convection). These results with H0 > 0 are consistent with the observations of ζ < 0 and μ < 0, corroborating the unstable atmospheric conditions. On average, H0 reached its maximum at noon, independently of the seasonal period, with values of 239, 123, and 96 W/m 2 in the summer, winter, and autumn RS-campaigns, respectively. At midday, the surface upward turbulent transport of thermal energy was more abundant in summer than in winter and autumn campaigns. It is striking that the maximum value of H0 at noon was higher in winter than in autumn, but this is reflecting only that on the day January 28, the average value of H0 registered from 10:00 LST to 13:00 LST was 142 W/m 2 against the average of 88 W/m 2 registered during the rest of the campaign. Moreover, we can observe that the hourly trends of H0 averaged over the month revealed maximums of 205 W/m 2 in summer, 82 W/m 2 in winter, and 106 W/m 2 in autumn at noon, which are consistent with the expected seasonal differences. Accordingly, we expect that the convective mixing heights reveal the seasonal differences correctly.
The mean values of friction velocity also were more significant in summer than winter and autumn. On average, the values were 0.32, 0.17, and 0.13 m/s in the summer, winter, and autumn In Figure 4, we observed that the larger values of temperature and relative humidity prevailed in summer, while in winter, the smaller ones. On average, the minimum, mean and maximum temperatures were 30, 33, and 35 • C in summer, 16, 18, and 19 • C in winter, and 26, 27, and 29 • C in autumn. For relative humidity, the values were 39%, 48%, and 55% in summer, 17%, 30%, and 44% in winter, and 24%, 38%, and 58% in autumn. These results are consistent with the Mexicali Valley climatology [13]. For wind speed, the averages of the minimum, mean, and maximum values were 2.2, 3.9, and 5.5 m/s for summer, 1.7, 4.3, and 6.8 m/s for winter, and 1.3, 2.1, and 2.9 m/s for autumn, with the larger values in winter and the smaller ones in autumn. Regarding the Mexicali climatology [13], the mean wind speeds obtained for the summer and autumn campaigns fell within the 25th to 75th percentile band, but for the winter campaign fell in the 10th to 90th percentile band.
In Figure 5, the plots in rows 3 and 4 show that, with very few exceptions, the stability parameters ζ = z/L (at z = 12 m) and µ = u * / f L were negative from 05:30 to 17:00 LST in summer, from 08:00 to 16:15 LST in winter, and from 06:30 to 16:15 LST in autumn. Therefore, during these periods, unstable atmospheric conditions prevailed in the surface with possible near-neutral conditions at the ends [53]. Moreover, the plots of the first row of Figure 5 show that the sensible heat flux was positive on average during the periods 05:00-17:45 LST in summer, 07:30-16:00 LST in winter, and 06:00-16:30 LST in autumn. During these periods, the surface heating produced thermal turbulence (convection). These results with H 0 > 0 are consistent with the observations of ζ < 0 and µ < 0, corroborating the unstable atmospheric conditions. On average, H 0 reached its maximum at noon, independently of the seasonal period, with values of 239, 123, and 96 W/m 2 in the summer, winter, and autumn RS-campaigns, respectively. At midday, the surface upward turbulent transport of thermal energy was more abundant in summer than in winter and autumn campaigns.
It is striking that the maximum value of H 0 at noon was higher in winter than in autumn, but this is reflecting only that on the day January 28, the average value of H 0 registered from 10:00 LST to 13:00 LST was 142 W/m 2 against the average of 88 W/m 2 registered during the rest of the campaign. Moreover, we can observe that the hourly trends of H 0 averaged over the month revealed maximums of 205 W/m 2 in summer, 82 W/m 2 in winter, and 106 W/m 2 in autumn at noon, which are consistent with the expected seasonal differences. Accordingly, we expect that the convective mixing heights reveal the seasonal differences correctly.
The mean values of friction velocity also were more significant in summer than winter and autumn. On average, the values were 0.32, 0.17, and 0.13 m/s in the summer, winter, and autumn measurement periods, respectively. Then, the downward turbulent transport of momentum was more significant in the summer than in the winter and autumn RS-campaigns. We observed the smaller values of friction velocity in the nighttime, especially during the dawn hours. On the monthly average, the temporal evolutions of the mean friction velocity exhibit trends with peaks of 0.4514, 0.2255, and 0.2326 m/s around the 16:00 LST, 13:00 LST, and 11:00 LST during the summer, winter, and autumn, respectively.

Vertical Profiles of Temperature
We carried out 52 RS under different stability conditions during the three campaigns (25 in summer 2010, 13 in winter 2012, and 14 in autumn 2016). In Table 1, we highlighted the sunlight conditions. In Figure 6, we presented the average temperature profiles obtained from the RS carried out under sunlight conditions. The only case we carried out diurnal RS for all seasons (summer, autumn, and winter) was that of the 12:00 LST soundings (Figure 6b,e,f). Because of the variability of the sunrise moment at the Mexicali latitudes, we obtained diurnal RS early in the morning (06:00 LST) only in the summer 2010 campaign.
For each RS, from the surface measurement campaigns, we determined the mean values of Monin-Obukhov length, friction velocity, sensible heat flux, and the stability parameters (ζ = z/L and µ = u * / f L) over the sounding duration. We also estimated the near-surface temperature lapse rates (Γ) of the T-profiles. In Tables 4-6, we presented the values of the main parameters for all the RS that we considered during unstable conditions. The RS revealed quasi-adiabatic and superadiabatic lapse rates, Γ −0.0099 K/m. Otherwise, in Tables 7-9, we gave the values of the same parameters for the RS that we considered under stable conditions. These RS revealed subadiabatic lapse rates, Γ > −0.0099 K/m. Atmosphere 2019, 9, x FOR PEER REVIEW 14 of 27

Estimation of the Convective Mixing Height with the LSVA, GM, CM, and PM
In Tables 4-6 The RS presented in Tables 4-6 were carried out during sunlight hours. All of them, excepting the soundings carried out at 18:00 LST during the summer 2010 campaign, were done under atmospheric surface conditions with ζ < 0, µ < 0, and H 0 > 0, which indicate unstable atmospheric conditions near the ground. The soundings of the 18:00 LST were performed under surface conditions with ζ > 0, µ > 0, and H 0 < 0, which indicate stable conditions. However, as we can observe, the near-surface temperature lapse rates of these RS were superadiabatic, indicating that they found unstable conditions in the first 100-200 m above ground. The radiosoundings of the 18:00 LST were carried out during the transition period of the sunset, which explains the small positive values of the stability parameters and the small negative values of H 0 .
The atmospheric conditions for the 06:00 LST soundings carried out during the summer 2010 campaign were unstable according to the stability criteria (ζ < 0, µ < 0), but the near-surface lapse rates of these temperature profiles were quasi-adiabatic or subadiabatic, indicating that they found stable conditions in the first 100-200 m above ground. In this case, however, the diagnostic parametrizations based on the friction velocity and Monin-Obukhov length cannot be applied since the formulae assume L > 0. Therefore, we used the LSVA, GM, CM, and PM to obtain the mixing height from the RS profiles. In order to apply the PM, however, we built the dry adiabatic profile with the surface temperature increased (subjectively) in 3 K (otherwise, no intersection would be found between the temperature profiles of the RS and dry adiabatic).
The graphs in Figures 7-9 illustrate the application of the LSVA, GM, CM, and PM to the VPT, SH, and T profiles obtained from the RS carried out at 12:00 LST during the sounding campaigns. In these Figures, the red lines indicate the VTP and SH profiles, the green lines represent their gradient profiles, and the blue lines represent the VPT and SH slab model profiles given by the LSVA. The graphs in the third column show the covariance (magenta lines) between the VPT and SH profiles. The last column illustrates the application of the PM for the 12:00 LST temperature profiles (solid lines). The dotted line denotes the dry adiabatic profiles.  In Figure 8, the VPT profile of the RS identified as 20120129-1200 reveals the presence of three inversions aloft at 730, 1210, and 2110 m, with VPT gradients of 0.022, 0.032, and 0.055 K/m, respectively. Probably, the third one is an inversion by subsidence, which are common in the Imperial County and the Mexicali Valley between November and June and relatively absent between July and October [54,55]. The presence of the semipermanent cyclonic cell of the Pacific Ocean can cause the fall of air masses. When these air masses descend, they undergo heating by compression, reaching temperatures that are higher than those of the layers that are closer to the ground surface. The intensity of these highly stable atmospheric conditions makes them difficult to disrupt and can persist for one or more days. These inversions may occur at heights between 2000 and 6000 m. The second inversion resembles a reminiscence of the day before. Therefore, we accepted the mean value of 706 m from the LSVA/VPT, GM/VPT, and PM as the valid one for 12:00 LST of 29 January 2012 (Table 4). In Figure 8, the VPT profile of the RS identified as 20120129-1200 reveals the presence of three inversions aloft at 730, 1210, and 2110 m, with VPT gradients of 0.022, 0.032, and 0.055 K/m, respectively. Probably, the third one is an inversion by subsidence, which are common in the Imperial County and the Mexicali Valley between November and June and relatively absent between July and October [54,55]. The presence of the semipermanent cyclonic cell of the Pacific Ocean can cause the fall of air masses. When these air masses descend, they undergo heating by compression, reaching temperatures that are higher than those of the layers that are closer to the ground surface. The intensity of these highly stable atmospheric conditions makes them difficult to disrupt and can persist for one or more days. These inversions may occur at heights between 2000 and 6000 m. The second inversion resembles a reminiscence of the day before. Therefore, we accepted the mean value of 706 m from the LSVA/VPT, GM/VPT, and PM as the valid one for 12:00 LST of 29 January 2012 (Table 4). Comparing the plots of the VPT and SH profiles in Figures 7-9, we observed that the SH profiles were less regular than the VPT profiles. Then, the idealized slab model (Figure 2) provided a better representation of the actual VPT profiles than the actual SH profiles. Nevertheless, the application of the LSVA to the SH profiles gave reasonable mixing height values, similar to those of the other methods, excepting the 12:00 LST case of 29 January 2012.
In Figure 8, the VPT profile of the RS identified as 20120129-1200 reveals the presence of three inversions aloft at 730, 1210, and 2110 m, with VPT gradients of 0.022, 0.032, and 0.055 K/m, respectively. Probably, the third one is an inversion by subsidence, which are common in the Imperial County and the Mexicali Valley between November and June and relatively absent between July and October [54,55]. The presence of the semipermanent cyclonic cell of the Pacific Ocean can cause the fall of air masses. When these air masses descend, they undergo heating by compression, reaching temperatures that are higher than those of the layers that are closer to the ground surface. The intensity of these highly stable atmospheric conditions makes them difficult to disrupt and can persist for one or more days. These inversions may occur at heights between 2000 and 6000 m. The second inversion resembles a reminiscence of the day before. Therefore, we accepted the mean value of 706 m from the LSVA/VPT, GM/VPT, and PM as the valid one for 12:00 LST of 29 January 2012 (Table 4).
Tables 4-6, and the plots in Figure 10, reveal similar results among the convective mixing height estimations obtained with the LSVA/VPT, GM/VPT, LSVA/SH, GM/SH, CM, and PM, for the same atmospheric sounding.
All the methods we used to estimate the convective mixing height gave similar results when the CBL was capped with a well-defined inversion layer with the VPT increasing and the SH decreasing sharply above the top, indicating a marked transition from a convectively less stable region below to a more stable region above. The local gradient approach and the LSVA lead to ambiguous estimations with scenarios where a sharp inversion layer is absent or where multiple layers in the troposphere with strong gradients were present. In general, the VPT profile under unstable conditions fits better the slab model profile than the SH profile, then the LSVA/VPT provided more reliable results than LSVA/SH. The reliability of the CM depends on the regularity of both profiles, VPT and SH.
From Tables 4-6, we observed that, on average, the convective mixing height in the Mexicali Valley ranged from 400 to 1400 m. These results are in agreement with the mixing height values (ranging from 300 to 1500 m) that Bei et al. [56] Tables 4-6, and the plots in Figure 10, reveal similar results among the convective mixing height estimations obtained with the LSVA/VPT, GM/VPT, LSVA/SH, GM/SH, CM, and PM, for the same atmospheric sounding.
All the methods we used to estimate the convective mixing height gave similar results when the CBL was capped with a well-defined inversion layer with the VPT increasing and the SH decreasing sharply above the top, indicating a marked transition from a convectively less stable region below to a more stable region above. The local gradient approach and the LSVA lead to ambiguous estimations with scenarios where a sharp inversion layer is absent or where multiple layers in the troposphere with strong gradients were present. In general, the VPT profile under unstable conditions fits better the slab model profile than the SH profile, then the LSVA/VPT provided more reliable results than LSVA/SH. The reliability of the CM depends on the regularity of both profiles, VPT and SH.
From Tables 4-6, we observed that, on average, the convective mixing height in the Mexicali Valley ranged from 400 to 1400 m. These results are in agreement with the mixing height values (ranging from 300 to 1500 m) that Bei et al. [56]

Summer 2010 Campaign
During the summer 2010 campaign, we carried out nocturnal atmospheric soundings at 00:00 LST. Therefore, the atmospheric conditions for these RS were expected to be stable. On average, during this set of RS, the observed values of the sensible heat flux and the stability parameters were H0 = −9 W/m 2 , ζ = 0.525, and μ = 59 (Table 7). These observations indicate that the 00:00 LST soundings were done under very stable atmospheric conditions. The near-surface lapse rates of the temperature profiles were subadiabatic ( = −0.0027 K/m, on average). We used N1981, Z1972S, and A1981S to estimate the Mexicali Valley midnight stable mixing height from the friction velocity and the Monin-

Summer 2010 Campaign
During the summer 2010 campaign, we carried out nocturnal atmospheric soundings at 00:00 LST. Therefore, the atmospheric conditions for these RS were expected to be stable. On average, during this set of RS, the observed values of the sensible heat flux and the stability parameters were H 0 = −9 W/m 2 , ζ = 0.525, and µ = 59 (Table 7). These observations indicate that the 00:00 LST soundings were done under very stable atmospheric conditions. The near-surface lapse rates of the temperature profiles were subadiabatic (Γ = −0.0027 K/m, on average). We used N1981, Z1972S, and A1981S to estimate the Mexicali Valley midnight stable mixing height from the friction velocity and the Monin-Obukhov length. Table 7 and Figure 11 summarize the results. Table 7. The summer stable mixing height of the Mexicali Valley estimated with the surface mean values of friction velocity and Monin-Obukhov length according to the N1981, Z1972S, and A1981S. Here, we also included the mean surface values of the stability parameters and sensible heat flux, and the near-surface temperature lapse rate of the radiosonde (RS).

Summer 2010
Mixing Height (m)   Figure 11. Estimations of the Mexicali Valley mixing height under the 00:00 LST stable atmospheric conditions during the summer 2010 campaign. We observe that the A1981S parametrization, Equation (17), overestimated the mixing height in comparison with the other two models during very stable conditions (μ > 50).
In Table 7 and Figure 11, we noticed that N1981, Z1972S, and A1981S provided similar estimations for the 00:00 LST stable mixing height. The smaller estimations were obtained with N1981, which was derived for stable lapse rate conditions with a wide range of stabilities, from near neutral to very stable [38]. However, as A1981S was derived in the limit of neutral or near-neutral conditions, it overestimated the mixing height in comparison with the other two models for the cases with very stable conditions (50 < μ < 100), but revealed a small difference for Z1972S in the case (20100718-00:00) with a moderately stable condition (10 < μ < 50). On average, the nocturnal (00:00 LST) mixing height value was 160 m.

Winter 2012 and Autumn 2016 Campaigns
For the winter 2012 and autumn 2016 campaigns, all the nocturnal RS were done under stable conditions. It is evidenced by the average surface values of the stability parameters and the nearsurface temperature lapse rates of the RS profiles. Tables 8 and 9 summarize the results.
During the winter 2012 sounding campaign, we carried out nocturnal soundings at 00:00 LST, 06:00 LST, 18:00 LST, and 21:00 LST. For this part, during the autumn 2016 campaign, we carried out soundings at 06:00 LST and 18:00 LST, again under nocturnal conditions. In Tables 8 and 9, we summarized the nocturnal stability conditions in the Mexicali Valley during the periods of these experimental campaigns. For all the RS, we observed that the surface values of the stability Figure 11. Estimations of the Mexicali Valley mixing height under the 00:00 LST stable atmospheric conditions during the summer 2010 campaign. We observe that the A1981S parametrization, Equation (14), overestimated the mixing height in comparison with the other two models during very stable conditions (µ > 50).
In Table 7 and Figure 11, we noticed that N1981, Z1972S, and A1981S provided similar estimations for the 00:00 LST stable mixing height. The smaller estimations were obtained with N1981, which was derived for stable lapse rate conditions with a wide range of stabilities, from near neutral to very stable [38]. However, as A1981S was derived in the limit of neutral or near-neutral conditions, it overestimated the mixing height in comparison with the other two models for the cases with very stable conditions (50 < µ < 100), but revealed a small difference for Z1972S in the case (20100718-00:00) with a moderately stable condition (10 < µ < 50). On average, the nocturnal (00:00 LST) mixing height value was 160 m.

Winter 2012 and Autumn 2016 Campaigns
For the winter 2012 and autumn 2016 campaigns, all the nocturnal RS were done under stable conditions. It is evidenced by the average surface values of the stability parameters and the near-surface temperature lapse rates of the RS profiles. Tables 8 and 9 summarize the results.
During the winter 2012 sounding campaign, we carried out nocturnal soundings at 00:00 LST, 06:00 LST, 18:00 LST, and 21:00 LST. For this part, during the autumn 2016 campaign, we carried out soundings at 06:00 LST and 18:00 LST, again under nocturnal conditions. In Tables 8 and 9, we summarized the nocturnal stability conditions in the Mexicali Valley during the periods of these experimental campaigns. For all the RS, we observed that the surface values of the stability parameters were positive and large, indicating that the atmospheric conditions on the surface were from very stable to extremely stable. We also observed that the near-surface temperature lapse rates of the RS profiles were all positive, revealing the presence of surface-based inversions (SBI).
An SBI is a clear indicator of a stable boundary layer, whose top can define an ABL height [25]. If an SBI is found in a sounding, none of the previous diagnostic parametrizations can be applied because they assume a different ABL structure [25]. Turbulence is suppressed, and no mixing height exists.
We evaluated the SBI as the height where the temperature changes from increasing to decreasing with the elevation. Ideally, at this height, the vertical gradient of the temperature profile becomes zero: (∂T/∂z) SBI = 0. Tables 8 and 9 and Figures 12 and 13 show our estimations of the SBI. On average, the SBI values were from 53 m at 18:00 LST to 183 m at 06:00 LST during the winter 2012 campaign, whereas in the autumn 2016 campaign, the values were from 132 m at 18:00 LST to 370 m at 06:00 LST. As we can observe in Figure 14 and the first row of Table 6, these inversions seem to break-off 1 h after sunrise (0 8:00 LST). An SBI is a clear indicator of a stable boundary layer, whose top can define an ABL height [25]. If an SBI is found in a sounding, none of the previous diagnostic parametrizations can be applied because they assume a different ABL structure [25]. Turbulence is suppressed, and no mixing height exists.
We evaluated the SBI as the height where the temperature changes from increasing to decreasing with the elevation. Ideally, at this height, the vertical gradient of the temperature profile becomes zero: ( ) = 0. ⁄ Tables 8 and 9 and Figures 12 and 13 show our estimations of the SBI. On average, the SBI values were from 53 m at 18:00 LST to 183 m at 06:00 LST during the winter 2012 campaign, whereas in the autumn 2016 campaign, the values were from 132 m at 18:00 LST to 370 m at 06:00 LST. As we can observe in Figure 14 and the first row of Table 6, these inversions seem to break-off 1 h after sunrise (~ 08:00 LST).        Table 6, the radiosonde (RS) profile 20161015-08:00 revealed a subadiabatic near-surface temperature lapse rate of  = -0.0046 K/m.

Estimation of the Mixing Layer Growth with the YDBB Model
The YDBB model [52] allows calculating the mixing height during its growing phase in terms of the time series of sensible heat flux H0(t). We determined the sensible heat flux as the covariance between the turbulent fluctuations of temperature and the vertical wind component measured with a 3D ultrasonic anemometer at z = 12 m, using Equation (2).
With the H0 data collected during the surface meteorological campaigns and the estimations of the convective and stable mixing heights, the parameters a and b of the YDBB model were determined. In the case of the summer 2010 campaign (the longer one), the results for the convective conditions included 14 estimations of the mixing height for the 06:00 LST (6 RS), 12:00 LST (6 RS), and 15:00 LST (2 RS). For the winter 2012 campaign, we obtained only three estimations of the convective mixing height for the 12:00 LST. For the autumn 2016 campaign, we had estimations of     Table 6, the radiosonde (RS) profile 20161015-08:00 revealed a subadiabatic near-surface temperature lapse rate of  = -0.0046 K/m.

Estimation of the Mixing Layer Growth with the YDBB Model
The YDBB model [52] allows calculating the mixing height during its growing phase in terms of the time series of sensible heat flux H0(t). We determined the sensible heat flux as the covariance between the turbulent fluctuations of temperature and the vertical wind component measured with a 3D ultrasonic anemometer at z = 12 m, using Equation (2).
With the H0 data collected during the surface meteorological campaigns and the estimations of the convective and stable mixing heights, the parameters a and b of the YDBB model were  Table 6, the radiosonde (RS) profile 20161015-08:00 revealed a subadiabatic near-surface temperature lapse rate of Γ = −0.0046 K/m.

Estimation of the Mixing Layer Growth with the YDBB Model
The YDBB model [52] allows calculating the mixing height during its growing phase in terms of the time series of sensible heat flux H 0 (t). We determined the sensible heat flux as the covariance between the turbulent fluctuations of temperature and the vertical wind component measured with a 3D ultrasonic anemometer at z = 12 m, using Equation (2).
With the H 0 data collected during the surface meteorological campaigns and the estimations of the convective and stable mixing heights, the parameters a and b of the YDBB model were determined.
In the case of the summer 2010 campaign (the longer one), the results for the convective conditions included 14 estimations of the mixing height for the 06:00 LST (6 RS), 12:00 LST (6 RS), and 15:00 LST (2 RS). For the winter 2012 campaign, we obtained only three estimations of the convective mixing height for the 12:00 LST. For the autumn 2016 campaign, we had estimations of convective mixing height for the 08:00 LST (1 RS), 12:00 LST (1 RS), 14:00 LST (3 RS), and 15:00 LST (2 RS). No RS soundings were available for sunrise during the winter and autumn campaigns. Therefore, we used the friction velocity surface measurements to estimate the mixing height at 07:00 LST (each day of these campaigns) with the diagnostic parametrization A1981S. In Table 10, we presented the results obtained for the YDBB model parameters, including the determination coefficient R 2 . In Figure 15, we showed the results of the application of the YDBB model (Equation (16)) to calculate the mixing height for all hours during the sounding campaigns.

Conclusions
The main goal of this work was to show that simple and straightforward methods can be applied to obtain reasonable estimations of the hourly mixing height that regulatory dispersion models require to carry out air pollution modeling studies. The use of sophisticated equipment (SODAR, RASS, and LIDAR) is more and more frequent in experimental research campaigns to study the structure of the atmospheric boundary layer; however, for air quality assessment based on short-term pollutants dispersion modeling, usually one estimates the mixing height from a more conventional approach, such as atmospheric sounding. Unfortunately, radiosondes involve some interference with aircraft flight plans, and long campaigns are not always allowed, nor with the necessary

Conclusions
The main goal of this work was to show that simple and straightforward methods can be applied to obtain reasonable estimations of the hourly mixing height that regulatory dispersion models require to carry out air pollution modeling studies. The use of sophisticated equipment (SODAR, RASS, and LIDAR) is more and more frequent in experimental research campaigns to study the structure of the atmospheric boundary layer; however, for air quality assessment based on short-term pollutants dispersion modeling, usually one estimates the mixing height from a more conventional approach, such as atmospheric sounding. Unfortunately, radiosondes involve some interference with aircraft flight plans, and long campaigns are not always allowed, nor with the necessary characteristics. Therefore, it often occurs that one has to carry out the air pollution modeling with a limited set of mixing height observations, which are mainly representative of daylight conditions.
In this paper, we reported the estimation of mixing height in the Mexicali Valley from RS temperature profiles and surface data of sensible heat flux, friction velocity, and Monin-Obukhov length, which we measured during the periods of 17-28 July 2010, 25-29 January 2012, and 11-15 October 2016. For unstable atmospheric conditions, we used four different approaches: the parcel and gradient methods [27], a least-squares variational approach [24], and a covariance method. For the turbulent SBL height, we used the parametrizations of Nieuwstadt [38], Zilitinkevich [35], and Arya [17], which allow estimations of mixing height from friction velocity and Monin-Obukhov length.
These complementary methods provided a set of mixing height values, which we used to implement a semiempirical model (YDBB) for the diagnostic of the mixing height growing as a function of sensible heat flux [52]. These models allowed estimating the hourly time evolution of mixing height during the periods of interest. A comparison of the YDBB estimations against the results of the PM, GM, and LSVA approaches showed reasonable agreements. In practice, the regulatory models (such as AERMOD and US-EPA) require hour by hour, the value of mixing height as input to perform the estimation of the hourly surface concentrations of the pollutants. However, an atmospheric sounding campaign will provide only a small number of RS a day during the campaign. In general, it would be too expensive and technically prohibitive to carry out one RS each hour of the day, even only during the sunlight hours. Then, the YDBB model may be considered as a practical tool to build up the hourly series of mixing height from the available RS. We also found good agreement with the mixing height values (ranging from 300 to 1500 m) that Bei et al. [56] determined in 2010 at the Parque Morelos (PQM, Tijuana Municipal System Theme Park in the center of Tijuana).
On the other hand, the nocturnal RS we carried out during the winter and autumn campaigns revealed the development of surface-based temperature inversions from sunset to 1 h before sunrise of the next day, with their breaking-off 1 h after sunrise, approximately.