Near-surface temperature inversion growth rate during the onset of the stable boundary layer

This study aims tofind the typical growth rate of the temperature inversionduring theonset of the stable boundary layer around sunset. The sunset transition is a very challengingperiod for numericalweather prediction, since neither accepted theories for the convective boundary layer nor those for the stable boundary layer appear to be applicable. To gainmore insight in this period, a systematic investigation of the temperature inversion growth rate is conducted. A statistical procedure is used to analyze almost 16 years of observations from the Cabauw observational tower, supported by observations from twoadditional sites (DomeCandKarlsruhe). The results show that, on average, the growth rate of the temperature inversion (normalized by the maximum inversion during the night) weakly declines with increasing wind speed. The observed growth rate is quantitatively consistent among the sites, and it appears insensitive to various other parameters. The results were also insensitive to the afternoon decay rate of the net radiation except when this decay rate was very weak. These observations are compared to numerical solutions of three models with increasing complexity: a bulk model, an idealized single-column model (SCM), and an operational-level SCM. It appears only the latter could reproduce qualitative features of the observations using a first-order closure. Moreover, replacing this closure with a prognostic TKE scheme substantially improved the quantitative performance. This suggests that idealized models assuming instantaneous equilibrium flux-profile relations may not aid in understanding this period, since history effects may qualitatively affect the dynamics.


Introduction
During the sunset transition the incoming shortwave radiation from the sun gradually lessens. As a result, the convective boundary layer (CBL) weakens and a temperature inversion starts to form near the surface, which is the onset of the stable boundary layer (SBL). Typically, the temperature inversion grows within a few hours and then remains relatively constant until sunrise. It is the early growth that we are interested in for this study. In particular, we use both observations and numerical models to investigate what parameters determine how fast the SBL evolves during this period.
During the late afternoon, turbulence in the boundary layer is typically decaying (Stull 2000). Much research effort has been aimed at understanding this period on the full range from fundamental studies (van Heerwaarden and Mellado 2016), to practical models (Nadeau et al. 2011), and observations (Blay-Carreras et al. 2014;Jensen et al. 2016). The same is true for the established SBL (Armenio and Sarkar 2002;Sorbjan 2010;Ansorge and Mellado 2016;Mahrt 2014;van Hooijdonk et al. 2017).
These previous studies of the SBL have mostly relied on Monin-Obukhov similarity theory (MOST; Monin 1970), which, once the SBL is established, provides a reasonable description for a broad range of conditions (Sorbjan 2010;Grachev et al. 2013). The decaying turbulence itself, before the onset of the SBL, may also exhibit self-similar behavior, at least for strongly convective situations with weak synoptic forcing (e.g., Nadeau et al. 2011;van Heerwaarden and Mellado 2016). Thus, our understanding of both the (quasi) steady CBL and the (quasi) steady SBL has progressed significantly, despite their own challenges (Holtslag et al. 2013;Lothon et al. 2014).
This understanding and resulting parameterizations may be of limited use during the sunset transition itself (e.g., Sun et al. 2003), since the onset of the SBL typically occurs a few hours before the net radiation becomes negative (van der Linden et al. 2017). Consequently, the decay of convective turbulence and the SBL coexist and interact during the sunset transition (Nadeau et al. 2011;Blay-Carreras et al. 2014). This period is therefore very challenging for numerical weather prediction models . Moreover, large-eddy simulation is complicated owing to the changing resolution and domain requirements during the transition (Basu et al. 2008) and the importance of other processes, such as radiative transfer (Edwards 2009). Therefore, it is very important to gain observationally based insight into the period just after the onset of the SBL. As an initial step, we investigate the evolution of the temperature inversion for a broad range of conditions.
Most studies that have focused on the near-surface temperature (inversion) in the SBL have related the instantaneous temperature inversion to other flow variables (André and Mahrt 1982;Sorbjan 2006;Vignon et al. 2017), on detailed investigation of individual cases (Sun et al. 2003), or on mechanistic understanding of the ''steady state'' (van de Wiel et al. 2012;van Hooijdonk et al. 2015;van de Wiel et al. 2017). However, the temporal behavior starting at the onset of the SBL has received relatively little attention.
A few studies have explicitly aimed at the temporal evolution of the temperature during this period. Earlier work often assumed a square root dependence of the inversion strength with time (e.g., Stull 1983). Whiteman et al. (2004) studied the onset of the SBL within sinkholes, as measured by the evolution of the temperature inversion. Since the wind was always weak inside the sinkholes, the evolution could be modeled considering radiative processes only, and hence the sky-view factor was identified as a key parameter for the temporal evolution inside the sinkholes. A follow-up study by De Wekker and Whiteman (2006) extended to plains and basins, using an exponentially decaying function to fit the time series of the absolute temperature to determine a ''typical'' time scale during clear and calm nights.
Another study by Pattantyús-Ábrahám and Jánosi (2004) also studied the absolute temperature evolution during 2 years at two distinct sites. They specifically limited their analysis to cases where an exponential fit provided excellent agreement with the time series, which again mostly corresponded to clear and calm nights with wind speed less than 2 m s 21 . Despite the large number of (predominantly moisture-and radiation-related) quantities they investigated, no clear relation of the time scale was found with any of these. Edwards (2009) used an idealized single-column model (SCM) to study the early evolution of the SBL at weak winds. His main focus was the importance of radiative processes on the near-surface heat budget, while the time scales of the boundary layer evolution did not receive attention. Later, Edwards (2011) imposed the tendency of the surface temperature in a dimensionless model to study the evolution of the sensible heat flux. Nieuwstadt and Tennekes (1981) used a similar imposed tendency to model the evolution of the boundary layer height.
We employ a similar approach to De Wekker and Whiteman (2006) in the sense that a typical time scale is determined for each night. In our case this time scale is based on the growth of the temperature inversion rather than the decay of the absolute temperature. Therefore, we use a logistic growth function, instead of an exponential function, to fit the time series of the temperature inversion (see section 2). Realizing that our approach is heuristic, the typical time scale is considered to be representative of the rate at which the SBL evolves.
The main analysis focuses on a large number of observations from the Cabauw Experimental Site for Atmospheric Research (CESAR) observatory (the Netherlands). These are supported by observations from the Karlsruhe station (Germany) and the Dome C observatory (Antarctica). The Cabauw and Dome C sites are essentially flat, and the Karlsruhe site is located in a small forest clearing in a wide valley. To aid the interpretation, we use three types of numerical simulations: a bulk model for the SBL, which is based on van de Wiel et al. (2017); a single-column model for the SBL (A. M. Holdsworth and A. H. Monahan 2017, unpublished manuscript); and the Regional Atmospheric Climate Model (RACMO)-SBL, the single-column version of an operational-level regional climate model of the Royal Dutch Meteorological Institute (Baas et al. 2017, manuscript submitted to Bound.-Layer Meteor.).
This analysis is used to answer the following questions: To what physical quantities is the rate at which the SBL evolves during the hours around sunset most sensitive? How may this (in)sensitivity be explained? And how similar are the results for three contrasting sites? A major part of the analysis is dedicated specifically to the wind speed dependence, since this quantity has been identified as a key parameter for the characterization of the SBL van Hooijdonk et al. 2015;Monahan et al. 2015;Acevedo et al. 2016;van de Wiel et al. 2017).
This paper is organized as follows. In section 2 the important meteorological quantities in this study, the observations sites, and the numerical methods are discussed. Section 3 contains the results, which are then discussed in section 4. The paper is concluded in section 5.

a. Parameter definition
We investigate the temporal evolution of the nearsurface temperature inversion Du 5 u(z 2 ) 2 u(z 1 ), where u(z) is the potential temperature at height z (and z 2 . z 1 ). For each night we obtain a characteristic temperature evolution rate of Du, which we denote › t Du.
In this section, we introduce how › t Du is determined, followed by several key parameters that may affect › t Du.
First, time is defined relative to the time the net radiation Q N becomes negative; that is, we define sunset at t 5 0 when Q N 5 0 [see also van Hooijdonk et al. (2015)]. Typically this occurred several hours before the meteorological sunset, and 1-2 h after the onset of the SBL (cf. van der Linden et al. 2017).
To define › t Du, a logistic growth function was used to fit the time series of Du ( Fig. 1a; details in the appendix). Previously, an exponential function was used by De Wekker and Whiteman (2006) to fit the time series of the cumulative (vertically integrated) cooling. For the temperature inversion we found the logistic growth more suitable, since it captures both the initial and the ''final'' stage of inversion growth with reasonable accuracy owing to the inflection point. Moreover, the logistic growth was more robust to variations in the interval on which the fit was based.
The logistic growth function is defined as where T, k, and t 0 are parameters to be estimated: T represents the estimated maximum of the temperature inversion, k is the steepness of the curve, and t 0 is the time at which the inflection point occurs. We define the typical evolution of Du as which is the derivative of f at the inflection point at t 5 t 0 . Additionally, a temperature scale Du max is defined as the maximum temperature inversion (based on 10-min averages) during each night. Using 30-min averages instead did not have a large impact (typically &3% decrease). To reduce the effect of synoptic changes during the night, the detection of the maximum was limited to t , 5 h and to 0000 LST at Cabauw and Karlsruhe, respectively. For most nights this had no effect on the value of Du max , since the maximum inversion typically occurs early in a night at these locations, probably because pressure acceleration may slightly reduce the inversion strength at later times. At the Dome C the evolution is more gradual and less prone to synoptic variability (Vignon et al. 2017). Therefore, no time restriction was imposed for this site. The temperature scale is used to normalize › t Du, which results in the inverse time scale Du 21 max › t Du. Alternatively T could be used as a temperature scale. However, we found that measuring Du max directly was more robust.
The quality of the fit greatly varied between nights owing to irregular time series. Therefore, we defined a measure for the quality of the fit, denoted by h 2 , which is the mean-squared residual with respect to the observations normalized by Du max . Only nights in which h 2 , 0:02 were used for the analysis (see Table 1). This value is an arbitrary trade-off between the fit quality and the quantity of the nights used for the analysis. Regular checks showed that the results were largely insensitive to the precise value of the threshold, with one caveat. During nights with clouds, the time series of the observations were typically more irregular, which affected the fit quality. Therefore, the selection on h 2 resulted in a sampling bias, since a relatively large fraction of the cloudy cases were affected compared to clear-sky cases (Fig. 1b). Similarly, the selection on h 2 introduced a minor bias toward lower relative humidity, larger (negative) net radiation, and a slightly weaker wind speed. However, the remaining set of nights in the present study still spans a much broader range of conditions than the stringent selection of Pattantyús-Ábrahám and Jánosi (2004).
Since wind speed has been identified as a key parameter for the SBL van de Wiel et al. 2012), we characterize each night by the mean wind speed around (the net-radiation based) sunset U sunset . At this time the wind speed has a strong vertical correlation, and thus the wind speed at a single altitude may be used as a measure for the mechanical forcing. In each case the highest level on which Du was based was used to obtain U sunset . As an averaging period we chose the interval t 2 [22, 1] h. We investigated the sensitivity to this choice by arbitrarily choosing different averaging periods within the time range t 2 [24, 1] h. Between these choices the resulting values of U sunset were strongly correlated (not shown), indicating insensitivity to the particular choice.
The established SBL may be characterized further using the net radiative emission Q N (defined positive toward the surface). Several hours after sunset the net radiation minimizes (maximizes in absolute sense) and remains relatively constant. This minimum value Q N,min was used as a characteristic measure. Since Q N evolves during the sunset transition itself, we also define › t Q N as the decay rate of Q N . Using a similar procedure as for › t Du, › t Q N was determined based on a linear fit to the time series of Q N (details are in the  For each quantity listed in the first column, the minimum (min) and maximum (max) values, the mean m, and standard deviation s are listed as an indication of the observational spread of each quantity for each site. The final row contains the number of nights that had sufficient data and fit quality to be used in the analysis and two fractions (f 1 , f 2 ); f 1 indicates how many nights could be used as fraction of the total number of nights, and f 2 indicates a similar fraction relative to the number of nights that had sufficient quality of the data (e.g., preselected on missing data points). appendix). A linear decrease provided a reasonable representation of the time series between several hours before sunset until shortly after sunset. Moreover, we did not find significant improvement when a more elaborate fit function (e.g., a sine function) was used. The slope of the linear fit provides an estimate of › t Q N . The ratio Q 21 N,min › t Q N defines an inverse time scale characterizing how fast Q N evolves toward its steady state.
Because of the differences in the relative importance of radiative, turbulent, and soil heat fluxes in the two regimes, it is possible that the character of the growth of the inversion may differ fundamentally in the very stable boundary layer (VSBL) from that in the weakly stable boundary layer (WSBL). As such, we will investigate the growth of the inversion separately in these two regimes. To separate the two regimes in observations, we make use of a hidden Markov model (HMM) analysis [as in Monahan et al. (2015), who also provided an introductory example for this method]. This approach models the observed variable x i , i 5 1, . . . , N, as being dependent on an unobserved Markov chain z i taking a set of K discrete values (the so-called hidden states). At each time, the HMM associates the observation with one of a number (in our case, 2) of probability distributions, resulting in a time series of HMM states. More technically, conditioned on z i residing in state k, the distribution of x i is described by a specified probability density function p(x i , l k ) (where l k , k 5 1, . . . , K, is a state-dependent set of parameters). For a specified number K of hidden states, the HMM algorithm finds maximum likelihood estimates both of the parameters l k and the hidden state trajectory z i . Having obtained this discrete hidden-state trajectory, conditional probability distributions of any observations concurrent with those used in the HMM can be estimated.
A benefit of the HMM approach is that the separation of states is determined by the intrinsic structure of the data, rather than relying on subjective thresholds (which may be impossible to define in highdimensional datasets). In this study, we follow the approach of Monahan et al. (2015), distinguishing the VSBL and the WSBL using a two-state HMM applied to three-dimensional data from Cabauw (the wind speed at 10 and 200 m and the potential temperature difference between 2 and 200 m, all at 10-min resolution).

b. Observation sites
Data from three observational sites were used for this study. The main focus lies on the CESAR observatory, while the Karlsruhe and Dome C stations are used to verify consistency of results among sites. To provide insight in the climate of each site, several characteristic values of the quantities introduced in the previous section are listed in Table 1.
The CESAR observatory near Cabauw, the Netherlands (51.9718N, 4.9278E), and its surroundings are extensively described in van Ulden and Wieringa (1996). Further details on the measurement equipment and the tower configuration may be found at http://www.cesarobservatory.nl. The surface is mostly covered with short grass and the surrounding land is relatively flat. We use observations collected between January 2001 and October 2016. The wind speed and direction (cup anemometers) and the (potential) temperature (KNMI Pt500 Element) are measured at 10, 20, 40, 80, 140, and 200 m. Additionally, the temperature is measured at 1.5 (all years) and 0.1 m (since 2013). The relative humidity is derived from the air and dewpoint temperature (Vaisala HMP243) at 1.5 m. All measurements are based on 10-min averages. The net radiation Q N is composed of the longwave and shortwave incoming and outgoing radiation components measured at 1 m. Measurements of the cloud cover have been obtained using a nubiscope (Wauben et al. 2010 Larger hills are approximately 50 km to the east and north. The tower is placed in a small forest clearing, a distance of 10-25 m in each direction from the tree line. The canopy height of the forest is 30 m. Further details on the site and the instrumentation may be found in Kalthoff and Vogel (1992) and at http://imkbemu. physik.uni-karlsruhe.de/;fzkmast/. To the northeast the urbanized area of the Karlsruhe Institut für Technologie is located, with buildings higher than 30 m at about 200 m from the tower. Since the forest affects the radiation measurements during sunset and sunrise (M. Kohler 2017, personal communication), these data were not used in our analysis. The sunset was therefore estimated based on the temperature inversion. For our analysis at this site the exact time of sunset is not important.
The Dome C observatory is located on Antarctica (75. 068S, 123.2008E). The station is located on a continental-scale dome (3233 m AGL), and the local slope is very weak (,1%). We use 30-min averages of the temperature (Campbell HMP155 thermohygrometers in mechanically ventilated shields) and the wind measurements (Young 05106 aerovanes) from a 45-m tower at z 5 [18,25,33,41] m. Additionally, the surface skin temperature is estimated from longwave radiation measurements. Detailed information on the instrumentation and processing are given by Genthon et al. (2013) and Vignon et al. (2016) (see also http://www.institut-polaire. fr/ipev-en/infrastructures-2/stations/concordia/). The analysis is limited to nights with a clear diurnal cycle in the months November-February during the period November 2011-December 2016 (i.e., the Antarctic summer). During the Antarctic summer the sun never sets; that is, the shortwave incoming radiation is always nonzero. Nonetheless, the net radiation does become negative during the ''night'' owing to variations in the zenith angle. It should be noted that the diurnal cycle of Q N is much weaker than at Cabauw, owing to the lower zenith angle and the higher surface albedo.

c. Numerical models
This section provides an overview of the three numerical models that were used to aid the interpretation of the observations. These models have been described elsewhere, and our main interest is to what extent each model reproduces elements of the observations. As such, only a brief summary of each model is given here, and we do not discuss the merits and drawbacks of specific model choices. For detailed descriptions of the models, the reader is referred to the cited papers.
The bulk model is used with the aim to obtain rudimentary understanding of the SBL dynamics. The model is based on the idealized model by van de Wiel et al. (2017), and it uses a single evolution equation for the temperature inversion: where Du is the time-dependent temperature difference between the surface and a constant air temperature at reference height z ref , c y is the heat capacity of the surface per square meter, and a proportionality constant A 5 rc p c D U, in which r is the density of air, c p is the specific heat capacity of air at constant pressure, and U is the wind speed. The neutral drag coefficient is defined as c D 5 k 2 / ln 2 (z ref /z 0 ), with k the von Kármán constant and z 0 the surface roughness length. The net radiation is a prescribed function of time Q N (t) 5 2qt when t , jQ N,min /qj and Q N (t) 5 Q N,min when t $ jQ N,min /qj, with q a proportionality constant controlling the decay rate of Q N . Both q and Q N,min are considered known.
The stability function is defined as Here, a is a fit parameter of the log-linear stability functions (Businger et al. 1971), and R b 5 (g/u 0 )(z ref Du)/U 2 is the bulk Richardson number, with g the acceleration by gravity and u 0 a reference temperature. This model is a strong idealization of reality, especially since MOST may not be valid during the sunset transition (Sun et al. 2003). Nonetheless it is instructive to determine whether or not the model reproduces aspects of the observations. Equation (3) is numerically integrated for various combinations of z ref and U.
The second model is an idealized SCM for the SBL (SBL-SCM). It solves the Reynolds-averaged equations for the two components of the horizontal wind vector profile (including the Coriolis force), the potential temperature profile, and the ground temperature as set up by Blackadar (1979). The model height is 5000 m and it uses a stretched grid with 100 vertical levels and the highest resolution near the surface. The turbulent diffusivities are modeled using the Prandtl mixing-length hypothesis and the standard Businger-Dyer relation to account for stable and slightly unstable conditions (Businger et al. 1971). The soil model is based on Blackadar (1976), and the radiation scheme is based on the effective atmospheric emissivity (Staley and Jurica 1972). Details may be found in A. M. Holdsworth and A. H. Monahan (2017, unpublished manuscript). Apart from an idealized treatment in the radiation scheme, the model neglects moisture effects.
RACMO is used with the aim to match the observations as closely as possible. ). The main difference with the IFS model is that the first-order closure model was replaced with a prognostic turbulent kinetic energy (TKE) closure model, such that TKE is not forced to be in equilibrium (details in Lenderink and Holtslag 2004;Baas et al. 2017, manuscript submitted to Bound.-Layer Meteor.). The model has 90 vertical levels up to a height of 20 km. The grid spacing near the surface is approximately 6 m, with the lowest model level at 3 m. To investigate the effect of using a TKE scheme, the model was also run using the default IFS first-order closure for the years 2014-15. Further details on the computational details and the model physics may be found in ECMWF (2007). The model was initialized at 1200 UTC using the daily forecast output of the three-dimensional RACMO, version 2.1 (van Meijgaard et al. 2008). The total model run spanned 48 h, but only the second 24-h period was used for the analysis.

Results
a. Inversion growth rate versus wind speed Figure 2 shows the joint frequency distributions of Du 21 max › t Du (based on the potential temperature difference between 40 and 1.5 m), and its individual components Du max and › t Du, with the sunset wind speed U sunset (also at 40 m) at Cabauw. Consistent with previous studies (e.g., Sun et al. 2012;van Hooijdonk et al. 2015;Monahan et al. 2015;van de Wiel et al. 2017) Du max suddenly increases when U sunset is below 5-7 m s 21 . The red and blue curves indicate the joint frequency distributions using subsets of the observations based on the hidden Markov model. Broadly speaking, the VSBL (red) occurs when the wind is weak, while the WSBL (blue) occurs when the wind is strong. Note that these observations contain both cloudy-and clear-sky cases.
It is not surprising that the absolute growth rate › t Du is closely related to Du max (Fig. 2b), since Du max is, in essence, the time integral of the growth rate. To investigate the qualitative differences, the inverse time scale is used. When this normalization was applied, no strong distinction between the regimes was evident (Fig. 2c). This fact is investigated further below.
To gain a more systematic insight into the wind dependence of the inverse time scale, Fig. 3a shows the median value of Du 21 max › t Du as a function of U sunset , where for each case U sunset was measured at z 2 (the highest level on which Du was based). The median is based on the subset of nights for which the sunset wind speed was in the interval U sunset 6 d, with d 5 0:5 m s 21 for Cabauw and Karlsruhe and d 5 1:0 m s 21 for Dome C. Varying d did not have a significant effect on the results, provided that sufficient nights were available per subset (not shown). Since at Dome C the observations spanned fewer years, and since only Antarctic summer nights could be used, a larger value of d was chosen for this site. At the same time, the temporal structure at Dome C is generally smoother than at the other locations, and › t Du could typically be determined with greater accuracy than at the other two sites. We also verified that using the mean value, instead of the median, did not have a significant effect on the results (not shown). Figure 3a shows how Du 21 max › t Du depends on U sunset . It shows an overall weak declining trend, with an apparent plateau at intermediate values of U sunset . There is no indication that the development of nights that become very stable (red) is qualitatively different from nights that  Figure 3b shows that the results from Karlsruhe and Dome C are remarkably consistent (i.e., in the range 0.2-0.4) with the observations at Cabauw, considering the differences among the sites. At Karlsruhe the range of U sunset was limited, such that we could not investigate if a similar decrease occurred at higher U sunset . The width of the distributions (the error bars) also complicate the interpretation of the observed trends. The observations at Dome C typically show lower values of Du 21 max › t Du, which may be explained by the weaker diurnal cycle (Argentini et al. 2014); that is, the magnitude of the cycle of Q N during a 24-h period is lower. This results in lower values of › t Q N (Table 1), which is further investigated below. Nonetheless the trend appears to be weakly downward at this location as well. The error bars in Fig. 3 also show that Dome C exhibits much less variability. This results in narrower distributions, further demonstrating the suitability of this site for idealized, and conceptual, research. Figure 4 indicates that the results in Fig. 3 are not sensitive to the measurement height of Du and U sunset . A closer look suggests that the plateau as observed at Cabauw for Dz 5 40-1.5 m is not a generic feature. For example, when the lower level is taken at z 5 0.1 m, the decreasing trend appears more gradual. However, these differences fall well within the observed error bars.
At Karlsruhe it is not clear if a trend over the relatively restricted range of U sunset values is present at any measurement height. At Dome C the downward trend appears to be present, and it appears to be stronger below 20 m. A possible explanation is that the boundary layer is typically much shallower at this location than at the two midlatitude locations. As a result, the top of the tower may be outside of the boundary layer (Vignon et al. 2017).

b. Width of the distributions
The size of the error bars in Fig. 3b suggests that other parameters than U sunset play an important role in determining Du 21 max › t Du as well. It also shows that the error bars at Dome C are smaller, which may be explained by the infrequent variations of moisture and cloud cover at this site (Lanconelli et al. 2011;Vignon et al. 2016). Moreover, measurements were used from a single season only.
The common factor of moisture, clouds, and season is that they are coupled (or at least correlated) to the evolution of Q N and to each other. Therefore, it is difficult to investigate directly how each quantity is related to Du 21 max › t Du. This is in contrast to U sunset , which is a relatively independent parameter, in the sense that U sunset is primarily a function the synoptic pressure gradient, while other parameters play a secondary role.
First, we investigate if high or low values of Du 21 max › t Du are systematically related to various Q N -related quantities. For this purpose we use the percentile division as in Fig. 3b: the 25% slowest sunset transitions (based on Du 21 max › t Du) of each subset are collected in a single group FIG. 3. (a) The median of Du 21 max › t Du as function of U sunset at Cabauw. The median was calculated for each subset based on a range of U sunset . The horizontal axis represents the center value of each subset. Colors represent all cases (black), all VSBL cases (red), and all WSBL cases (blue). (b) Comparison of the median inversion growth rate as function of U sunset for three different sites: Cabauw (black), Karlsruhe (blue), and Dome C (green). The respective vertical ranges over which the temperature inversions were measured are 40-1.5, 60-2, and 18-0 m. The thin lines in both panels represent the number of nights (310 3 ) available for each subset with the same color (right axes). The error bars illustrate the typical width of the distributions, which is measured by the 25th and 75th percentiles (i.e., 50% of the observations fall within the error bars).
(and similarly for each subsequent 25%). Figure 5a shows the ensemble-averaged time series of each group for the cloud cover, the relative humidity at 1.5 m, and the net radiation. This shows that, on average, nights with the lowest values of Du 21 max › t Du are associated with higher levels of the cloud cover and the relative humidity. This may explain the lower (absolute) values of Q N , both during the day and the night. However, Fig. 5b shows that winter (December-February) nights are disproportionally represented in the slowest 25%, which may also explain the weaker trend of Q N .
A similar analysis was performed using several other quantities. We found no systematic relation with the wind direction, the absolute temperature, the friction velocity, the wind speed at 200 m, or the sensible heat flux (all before sunset). Therefore, we do not present results related to those quantities.
Next, we focus directly on how the evolution of the net radiation is related to Du 21 max › t Du. A similar approach is taken as for U sunset , which essentially corresponds to treating › t Q N (section 2) as an external parameter without explicit concern for the processes (such as clouds) that determined › t Q N in a particular night. We define subsets based on nights with values for › t Q N that are in the range › t Q N 6 d, where d 5 2.5 W m 22 h 21 for Cabauw and d 5 0.5 W m 22 h 21 for Dome C. Figure 6a shows the median Du 21 max › t Du as a function of › t Q N . Consistent with Fig. 5, we observe a decrease when › t Q N is small (in absolute sense) at Cabauw.
At Dome C the diurnal cycle of Q N is weak, such that only a very narrow range of › t Q N occurs. Consequently, no systematic comparison with Cabauw may be made. In fact, the close agreement is somewhat surprising. At Dome C the sky is often clear, and › t Q N was determined with great accuracy. Conversely, at Cabauw low values of j› t Q N j were typically observed during overcast conditions, which often resulted in poor estimates of › t Q N , owing to the irregular time series of Q N . Although Fig. 6a is suggestive of a direct relationship of Q N and its evolution with Du 21 max › t Du, the robustness of this relationship remains unclear. Figure 6b shows essentially the same result based on Q 21 N,min › t Q N , although it is less clear if a downward trend at low Q 21 N,min › t Q N is present owing to the size of the error bars. One may hypothesize that when the net radiation takes very long to reach its final state (Q 21 N,min › t Q N / 0), also the evolution of the temperature inversion should be very slow (Du 21 max › t Du / 0). Following this reasoning, possible support for a trend in Fig. 6b is that for low values of Q 21 N,min › t Q N , the system evolves so slowly that it is close to local equilibrium and Du 21 max › t Du is directly related to Q 21 N,min › t Q N . Conversely, when Q 21 N,min › t Q N / ' (i.e., the sunset transition occurs instantaneously), the system is not in local equilibrium with the evolution of Q N ; hence, Q 21 N,min › t Q N is no longer a relevant time scale for Du 21 max › t Du, and the curve levels off.

c. Numerical models
The observations suggest that Du 21 max › t Du takes a relatively small range of values (0.2-0.4 h 21 ) and it decreases slightly with increasing U sunset . It is not clear what physical processes control this relation, as U sunset is not an externally imposed control parameter, but it is a boundary layer quantity dynamically coupled to the various other quantities we are considering. Therefore, we use three types of numerical models, such that a more controlled test may be performed. Information on which type(s) of model(s) reproduce(s) the trend of the observations could give an indication of what dynamics are responsible.
Given the idealized nature of the bulk model, its use for studying the effect of the wind speed provides qualitative insight only. Nonetheless, it would be highly advantageous if such a model is sufficient to explain the observations, since it provides the most direct mechanistic insight. As a first step, one could approach the SBL growth as a simple diffusion problem, which is governed by a time scale z 2 ref /K, with K 5 ku * z ref the turbulent diffusivity (chosen equal for heat and momentum). The effect of the wind speed is twofold. First, a larger wind speed would increase the turbulent diffusivity and, therefore, reduce the time scale. Conversely, increasing U sunset would also increase the characteristic length scale, such as the Obukhov length, the boundary layer height The bulk model is used to investigate both effects. Figure 7a shows that, at a fixed z ref , the evolution of the normalized inversion strength in all cases follows more or less the same slope, which is dictated by the magnitude of Q N /c y [Eq. (3)]. This indicates that the turbulent heat flux and the net radiation are always close to balance, which is due to the simplicity of the model. The small differences in slope between the curves are caused by variations in U 40 , which represents U sunset at 40 m. When U sunset was increased, while keeping z ref constant, the initial growth of Du 21 max › t Du was indeed larger. To test the effect of z ref , U sunset should not be kept constant, since then increasing z ref would effectively decrease the turbulent mixing. Therefore, we vary z ref while modifying U sunset such that the neutral friction velocity u * N 5 c 1/2 D U sunset is kept constant. Figure 7b shows opposite results from Fig. 7a; that is, the growth rate of Du 21 max Du was reduced at higher z ref and U. In summary, since both effects work in opposite direction, the bulk model remains inconclusive on the compound effect on the evolution of Du 21 max Du. In the SBL-SCM both the wind speed and a vertical length scale are dynamic quantities that result from the imposed geostrophic wind speed, the surface properties, and the initial conditions. If the results would be consistent with the observations, this would imply that the compound effect of U sunset on the turbulent mixing and the vertical length scale is sufficient to explain the observed trend in Fig. 3. However, Fig. 8a shows that an increased geostrophic wind (and thus U sunset ) accelerates the growth of Du 21 max Du. Since the rudimentary soil model is not tuned to any of the observational sites, we investigate the sensitivity of the evolution of Du/Du max to the soil model. Figure 8a shows that the results are comparable between dry clay and fresh snow. In fact, a wide variety of other soil types gave very similar results initially in the sense that the normalized inversion grew faster with increasing wind speed (not shown). As such, we cannot attribute the poor representation of the sunset transition in this model to the particular choice of soil type.
The qualitative disagreement with the observations may be explained partially by the invalidity of the equilibrium flux-gradient relations during the sunset transition. Moreover, the model is initialized right before the onset of the SBL, such that the SBL growth is very sensitive to the initial conditions (Fig. 8b). In fact, model simulations are more sensitive to the initial conditions than to the soil type. Thus, although the model provides useful insight in the relation between the model parameters and dynamics after the sunset transition, it appears less suitable to model the transition itself.
Since the idealized models do not provide qualitative similarities with the observations, we also investigate a more realistic model. The RACMO-SCM with a TKE scheme (as in Baas et al. 2017, manuscript submitted to Bound.-Layer Meteor.) does not rely on a first-order parameterization of turbulence. Additionally, this model is more elaborate than the SBL-SCM, in the sense that it models a full diurnal cycle with a simplified representation of horizontal advection and moisture schemes and more advanced radiation and soil schemes. Therefore, we may expect that the model provides a more realistic representation of the SBL transition. Figure 9a shows the median of Du 21 max › t Du (similar to Fig. 4). It shows that reasonable agreement with the observations was obtained, especially at higher wind speeds. Modeled time series are typically smoother than observational time series. Therefore, a slightly smaller value of Du max in the model is not surprising and the overestimation of the inversion growth rate by the model could be due to this fact. However, the variation among various model levels is in contrast with observations (cf. Fig. 4) and suggests that the model is still prone to other biases. Moreover, the increase at very low wind speeds should be approached with caution, since these types of models typically do not perform well in this range. Figure 9b shows that there also is a systematic overestimation by the model, when the sensitivity to the evolution of the net radiation was investigated. The shapes of the curves are nonetheless qualitatively similar.
The reasonable agreement between the RACMO-SCM and the observations could be explained by the more realistic prognostic representation of the TKE in this model. To test this hypothesis, a new set of runs was performed (spanning 2 years of observations) using the first-order IFS closure for turbulence of the ECMWF. Figure 9a shows that this severely reduces the model performance. The reduction of the performance is almost completely due to the poor prediction of › t Du, while reasonable agreement between the two schemes was found for Du max (not shown). Nonetheless, it is clear that between the SBL-SCM and the RACMO-SCM with the IFS scheme a substantial qualitative improvement was already achieved. The use of a prognostic TKE scheme provided another significant improvement step (Fig. 9a).
In summary, the bulk model and the SBL-SCM appeared to be unable to provide qualitative agreement with the observations. This suggests that the physical processes of particular importance are either absent or poorly represented in these models. A first qualitative important improvement is observed when the RACMO-SCM with the IFS scheme is used, which we attribute to the fact that a full diurnal cycle was modeled, and thus that a realistic and coherent state of the boundary layer was known at the time of the onset of the SBL. Further improvement by using the TKE scheme suggests that the history of the boundary layer is important in a more general sense. This seems to be consistent with the notion that turbulence is not in local equilibrium (e.g., Blay-Carreras et al. 2014).

Discussion
In this paper, we investigated the relation between the growth of the SBL across the evening transition and other physical quantities. For this approach it was necessary to use a statistical fitting procedure, which neglects the large variety in how the temperature inversion may develop in individual cases.
Previous studies found strong relations between the wind shear and the temperature inversion van Hooijdonk et al. 2015), which could be explained using a strongly idealized model (van de Wiel et al. 2012. Furthermore, evidence was found of at least two distinct regimes in the SBL (Grachev et al. 2005;Banta et al. 2007;van Hooijdonk et al. 2015;Monahan et al. 2015), of which the most stable regime seems incompatible with MOST-type scaling laws (Sorbjan 2010;Mahrt 2014). Considering the distinct behavior between the regimes later in the night, it would not have been surprising if signals of distinct development were detected at an earlier stage. However, once normalized by the maximum temperature inversion, it appears there is little evidence of distinct regimes (Fig. 3). A possible explanation is that the initial growth is dominated by the same processes (e.g., continuous turbulence or remnants of convective motions). At a later stage turbulence may be weak and intermittent, such that other processes (e.g., the ground heat flux) may be more significant . In fact, the wind shear appeared to have only a very mild effect on Du 21 max › t Du. We hypothesize that the weak effect of the wind may be due to counteracting processes: a stronger wind shear leads to more mixing but also to a deeper boundary layer (i.e., more available warm air). Other possibilities may be a relation to the weaker decay of convective turbulence in case wind shear is strong (Pino et al. 2006). Pattantyús-Ábrahám and Jánosi (2004) suggested that the tropospheric water content was the most important factor in the evolution of the nocturnal temperature. Our observations indicate that very low (absolute) values of › t Q N (associated with clouds, moisture, and season) may lead to slower growth of Du 21 max › t Du. Additionally, the results showed that moisture-related dynamics (e.g., clouds) may partly account for the spread in Du 21 max › t Du. It is unfortunate that Pattantyús-Ábrahám and Jánosi (2004) did not investigate the time scale in relation to Q N , since Q N (and its evolution), moisture and clouds, and the seasonal cycle are coupled and/or correlated. As a result, disentangling effects of individual parameters is complicated.
To further complicate the problem, we found indications that the daytime evolution (e.g., how the CBL decays) may be an important factor that affects the SBL growth (Fig. 9). This would explain why no, or only weak, relations were found with many contemporaneously measured variables (e.g., wind speed, temperature, humidity). Since these quantities contain little information on the SBL growth, it seems unlikely that a simple physical model assuming instantaneous equilibrium of turbulent fluxes to the resolved state could qualitatively reproduce the early SBL growth, and how this growth depends on the various parameters.
Consistent with the present results, Edwards (2009) observed clear temporal differences depending on whether or not the full diurnal cycle was modeled. A plausible explanation is that when a full diurnal cycle is being modeled, the flow constitutes a physically realizable state at the time of the onset of the SBL. The archetypal neutral boundary layer (a logarithmic wind profile and no temperature gradients) may not be a realistic state of the boundary layer at any time around sunset. Consequently, they cannot reproduce important aspects of the transient SBL ( Figs. 7 and 8). This suggests that at least the presunset history (i.e., the afternoon decay of convective turbulence) should be taken into account in models that aim to study the onset of the SBL.
The results of Jensen et al. (2016) also suggest that the presunset state is important for the dynamics at sunset. They investigated detailed flow characteristics during this presunset period, and they found that the (third order) budget terms of the sensible heat flux contain clues about the nature of the countergradient fluxes around sunset. However, it remains unclear if these budget terms qualitatively impact the growth rate of the temperature inversion, since the RACMO-SCM is able to reproduce the inversion growth with reasonable accuracy, even though it does not resolve third-order terms.
The studies of De Wekker and Whiteman (2006) and Pattantyús-Ábrahám and Jánosi (2004) shared similar approaches to determine a typical evolution of the SBL. These approaches were distinct from this study, but they are related in the sense that, in first-order approximation, a fixed functional form was used to fit the temporal evolution of a temperature-based quantity. For a set of 18 individual nights from 9 different locations, De Wekker and Whiteman (2006) used an exponential function to the normalized cumulative cooling during the whole night. The time scale that they found ranged from 3 to 7.5 h, corresponding to a rate of ;0.13-0.33 h 21 . The related approach by Pattantyús-Ábrahám and Jánosi (2004) yielded values for the inverse time scale of 0.1-0.5 h 21 , compared to the scale of 0.2-0.4 h 21 that we found at three distinct sites. As such, there appears to be some generality to the evolution rate of the SBL among studies, locations, measurement heights (Fig. 4), and a broad range of conditions. So despite the complexity of the dynamics, a simple ad hoc model for the evolution of the SBL (or at least the temperature inversion) may be sufficient to be useful in studies that do not require detailed insight in the boundary layer itself. The consistency among sites could be further verified by performing a similar analysis on a large number of sites [e.g., as in Monahan et al. (2011)]. Of course, such an analysis could not achieve the same level of detail per site, but it could reveal a systematic dependence on other parameters, such as soil characteristics.
Since at Dome C clouds are mostly absent and conditions are dry in general, the diurnal cycle of Q N is mostly determined by the zenith angle. Therefore, an alternative approach would be to conduct a more detailed investigation of this site. This might not only be beneficial for understanding the evolution of the SBL, but also for the decay of the CBL. With respect to the latter, Sorbjan (1997) used large-eddy simulation to find that turbulence decreases faster when Q N decreases more abruptly (i.e., smaller › t Q N ; note that › t Q N is negative), which is consistent with the bulk model and may explain the slower dynamics when › t Q N is less negative (i.e., at Dome C, or when there are clouds).
Finally, we found that an operation-level SCM results in reasonable agreement in terms of statistics. Unfortunately these models are typically complex and may behave unexpectedly if dynamical features (such as clouds) are turned off, which makes them less versatile than more idealized models. It may be possible, however, to investigate the effect of small perturbations of preexisting cases to gain insight in the effect on the growth of the temperature inversion.
Cabauw observations. The following funding agencies are acknowledged for their contribution: IvH is supported by the Netherlands Organisation for Scientific Research (NWO, ALW Grant 832010110); AHM, AH, and CA are supported by the Natural Sciences and Engineering Research Council Canada; and BvdW and PB are supported by an ERC Consolidator grant (648666). Finally, the University of Victoria is thanked for hosting IvH during several weeks for this research.

Fitting Procedure
The fitting procedure was aimed at estimating a ''typical'' inversion growth rate for each night, measured by › t Du. Considering the irregularity of the time series, different choices of fitting procedure may lead to a better fit for some nights and a poorer fit for others. The procedure described below appeared to be reasonably robust, in the sense that the quality of the fit resembled the irregularity within a particular night.
The procedure used Du (measured over a certain vertical range) and the net radiation. Furthermore the time is defined relative to sunset; that is, we defined t 5 0 when Q N 5 0. There are several main steps in the procedure: First, the time was located when the stable boundary layer (SBL) starts-that is, when Du crosses zero just before sunset. This time is defined as t min . Next, the maximum inversion Du max was determined, and t max was defined as the time when Du 5 0:9Du max for the first time after sunset. We found that using 90% of the maximum leads to less sensitivity to spikes.
Third, the minimum value of Du was subtracted from the time series, since Eq. (1) starts at 0, while Du is slightly negative during the day.
The lsqcurvefit function of MATLAB was used to fit Eq. (1) to the shifted Du-time series of each day between t 5 t min 1 a and t 5 t max 1 b, where a 5 [0, 22, 21, 22, 21, 21, 0] and b 5 [0, 0, 0:5, 1, 0, 1, 1]. This means that seven distinct fits were made for each day. No fits were made when fewer than five observations were available. Figure 1 in section 2c shows two examples of a fit for two different nights. For each fit the typical evolution › t Du 5 kT/4 was determined.
Finally, for each night, the extremes were removed, i.e., the two highest and two lowest estimates of › t Du were omitted. The remaining three estimates were then averaged to give a final estimate of › t Du, which was used for the analysis.
To determine › t Q N a similar procedure was used. First, we define Q N,max the maximum of Q N in the time interval 24 , t , 0 h and Q N,min the minimum in the interval 0 , t , 5 h. Again we define t min and t max . The time when Q N becomes less than 0:9Q N,max was used as t min and the time when Q N becomes less (more negative) than 0:9Q N,min was used as t max . A linear function f (t) 5 a 1 t 1 a 0 was used to fit the time series between t 5 t min 1 a and t 5 t max 1 b. The mean value of the middle three values of a 1 was then used as an estimate of › t Q N .