The role of subsidence in a weakly unstable marine boundary layer: a case study

. The diurnal evolution of a cloud free, marine boundary layer is studied by means of experimental measurements and numerical simulations. Experimental data belong to an investigation of the mixing height over inner Danish waters. The mixed-layer height measured over the sea is generally nearly constant, and does not exhibit the diurnal cy-cle characteristic of boundary layers over land. A case study, during summer, showing an anomalous development of the mixed layer under unstable and nearly neutral atmospheric conditions, is selected in the campaign. Subsidence is identiﬁed as the main physical mechanism causing the sudden decrease in the mixing layer height. This is quantiﬁed by comparing radiosounding proﬁles with data from numerical simulations of a mesoscale model, and a large-eddy simulation model. Subsidence not only affects the mixing layer height, but also the turbulent ﬂuctuations within it. By ana-lyzing wind and scalar spectra, the role of subsidence is further investigated and a more complete interpretation of the experimental results emerges.


Introduction
Measurements of large-scale divergence in the atmospheric boundary layer (ABL) are difficult and often contaminated by error (Lenschow et al. , 2007). Large-scale divergence in ABL is governed by subsidence, which depends mainly on synoptic-scale conditions. Although subsidence velocity rarely exceeds a few cm s −1 , it may significantly influ-ence mass conservation, material advection and mixing layer growth (Stull, 1988). Thus, considering that consequences of subsidence can be relevant, it is crucial to model the physical process correctly and estimate it quantitatively by indirect methods and/or numerical simulations, which can be used to integrate the knowledge coming from the experimental observations. Being associated with synoptic-scale variation, ABL subsidence velocity is treated as a mean field, unaffected by turbulence or by rapidly varying fluctuations. In a nutshell, subsidence parametrization corresponds to estimating a negative vertical velocity, generally assumed to be constant over ABL space and timescales (Stull, 1988).
Due to the lack of accurate divergence data from meteorological measurements, different approaches are adopted. In some ABL studies, subsidence velocity is -for simplicity -neglected or considered to be negligible (Batcharova and Gryning, 1991;Margulis and Entekhabi, 2004); while in other studies it is explicitly considered (see, e.g., Batcharova and Gryning, 1994;Yi et al., 2001;Bellon and Stevens, 2012). When this is the case, a common parametrization is to assume horizontal divergence constant with height. By mass continuity, this implies that subsidence velocity is proportional to the height z (Stull, 1988;Sempreviva and Gryning, 2000;Stevens et al., 2001;Letzel and Raasch, 2002;Mirocha and Kosović, 2010) , Recently, lidar measurements have revealed their potential to study boundary layer height variation and evolution (Eichinger et al., 2005;Di Liberto et al., 2012). These studies often rely on prognostic equations of the boundary layer height evolution such as the one derived in Batcharova and Gryning (1994), where large-scale subsidence velocity is needed as an input parameter. A possible choice (see Eichinger et al., 2005) is to use a relation such as where w subs (z i ) is the subsidence velocity at the boundary layer height z i , w RL is the negative vertical velocity at the top of the residual layer and z RL is the height at the top of the residual layer, thus estimating subsidence velocity from the residual layer of the day before. In Flagg (2005), a comprehensive discussion of different subsidence parametrizations is done. It emerges that, particularly when trying to model multiple-day evolution of the ABL, a constant value of the subsidence parameter can not account for change in synoptic regime or other local effects, creating a potentially inaccurate parametrization.
The development of a unified modelization of large-scale subsidence has been hindered by the difficulty of having accurate measurements of low-magnitude vertical velocities at synoptic and subsynoptic scales (Muschinski et al., 1999). Parametrizations often include its effects together with those of, e.g., radiation and turbulence (Carlson and Stull, 1986), or large-scale advection.
While subsidence contributes to reducing the boundary layer height, entrainment acts to increase it by mixing stably stratified air from above into the unstable boundary layer. More generally, the relative weight of the top fluxes, due to entrainment and subsidence, to surface ones may lead to different regimes, thus stressing the importance of a detailed description of all phenomena possibly present in the evolution of the boundary layer.
Here we want to disentangle the role of subsidence only, by considering a case study of a cloud free, marine boundary layer under weakly unstable conditions. We show that a numerical approach coupling mesoscale and large-eddy simulation (LES) modeling is appropriate for quantifying the effect of subsidence on the mixed layer inversion growth. Subsidence is identified as the key factor responsible for the observed collapse of the mixed layer. Moreover, by comparing the output of a LES run with subsidence to that of a control simulation without subsidence, we are able to quantify turbulent fluctuation evolution, otherwise unaccessible.
The paper is organized as follows. Section 2 shortly describes the experimental site and the apparatus, together with the case study. Sections 3 and 4 report on the numerical simulations with a mesoscale model and with a large-eddy simulation one, respectively. The rationale for using both is to have a quantitative control both on mean profiles and on small-scale turbulent fluctuations. Results are presented  Sempreviva and Gryning (2000). in Sect. 5, while conclusions and perspectives are discussed in the last section.

The experiment
A meteorological measuring station on the island of Anholt in the Kattegat Sea (lat = 56.7 • N, lon = 11.57 • E), between Denmark and Sweden (see Fig. 1), was operational from September 1990 to October 1992, as a part of the -90 Hav-90 marine research program funded by the Danish National Agency of Environmental Protection.
The goal was twofold: (i) investigating the climatology of the mixed layer height and the structure of the turbulence in the marine boundary layer (MBL) over inner Danish waters; (ii) quantifying the pollutants transport from the mainland and typical deposition rates into the sea. To monitor Nonlin. Processes Geophys., 21, 489-501, 2014 www.nonlin-processes-geophys.net/21/489/2014/ turbulent fluctuations associated with marine conditions (corresponding to wind blowing from the sector between 240 and 360 degrees), a 22 m-high meteorological mast was placed as close as possible to the shoreline, i.e., approximately at 10 m, on the western part of the island. The mast was equipped with instrumentation for standard measurements of wind speed U and direction DIR, temperature T , specific humidity q, pressure P , and solar radiation R. Pressure and solar radiation are measured at the surface; while high-frequency time series (20 Hz) of wind speed components, temperature and humidity were performed at the height of 22 m. Wind and temperature time series were recorded by a sonic anemometer (Kaijo-Denki DAT/TR-6lB); humidity time series were recorded by a fast response humidiometer (OPHIR Corporation, Lakewood, CO). All parameters were averaged over 10 and 30 min lapses. Radiosondes, type RS-80 by Vaisala, were on average released three times a day providing vertical profiles of wind direction and speed, temperature, humidity and pressure. The vertical profiles were recorded with a frequency of 0.5 Hz. With a radiosonde ascent velocity of 2.5 m s −1 , the frequency corresponded to a vertical resolution of approximately 5 m. In Sempreviva and Gryning (2000), a statistical study of the growth of the mixing height over two years was presented, based on the data recorded by the mast at Anholt Island. One of the key results is that the mixing layer growth mostly depends on the temperature gradient at the air-sea interface, i.e., the temperature difference between sea and the air mass just above it. In Fig. 2, we report data measured by the mast during four consecutive days in the summer period (15-18 June 1992). We note that the wind is constantly blowing from the west, i.e., from the sea. It is also to be noted that heat flux at the surface shows very little variation. As reported in Sempreviva and Gryning (2000), between 15 and 16 June, a mixed layer starts to grow after midnight and continues until the afternoon of 16 June, when it sinks. Correspondingly, a lower inversion develops. A similar behavior of the lower inversion is found in the measurements between 17 and 18 June. The interpretation of the observations is that the passage of cold and dry air masses from the west sector, from about 18:00 UTC, 15 June (their Fig. 5), associated with an increase in the temperature difference between the air and the sea surface, could be responsible of the observed abrupt collapse of the mixed layer height on 16 June.
Since this is not at the core of our investigation, we just mention that a second inversion is often detected over the marine boundary layer. There is no actual agreement about its origin and different phenomena have been proposed as a possible cause: the presence of a residual inversion from the previous mixing layer or a convective layer over the island where the measurements were taken (Sempreviva and Gryning, 2000); the presence of a boundary layer over land advected over the sea or the development of strato-cumulus clouds in weak frontal zones connected to low-pressure systems (Johansson et al., 2005). Here we focus on the MBL evolution during 16 June, to disentangle and quantify the role of subsidence in both mean fields and small-scale turbulent fluctuations further. With this aim, we first describe the mesoscale atmospheric condition obtained from a numerical simulation lasting 60 h (14 June, 12:00 UTC-17 June, 00:00 UTC) obtained with the WRF model, supporting the presence of large-scale subsidence and giving a quantitative measure of the subsidence velocity. We then refine our analysis by means of a large-eddy simulation of the boundary layer evolution at Anholt Island during the morning of 16 June, and lasting 9 h approximately.

The mesoscale conditions via a WRF numerical simulation
The WRF-ARW model, version 3.0, has been implemented to simulate the meso-and large-scale features of the case study. Initial and boundary conditions are taken from the ERA-Interim reanalysis (T255 spectral resolution approximately corresponding to 0.75 • ) (Untch et al., 2006). The model run starts at 12:00 UTC, 14 June and lasts for 60 h. The number of vertical levels is 40, extending up to 20 km, but more closely spaced in the atmospheric boundary layer. Two two-way nested domains, with horizontal resolutions respectively of 16 and 4 km, are employed. The number of grid points in the two domains are, respectively, 109 × 109 in the outer grid, and 161 × 161 in the inner grid (Fig. 3). The domains are centered in the location of the measurement site. The model configuration is the same implemented and tested in Miglietta and Regano (2008), and Moscatello et al. (2008), which includes the following parametrization schemes: Yonsei University PBL non-local scheme (Hong et  al., 2006), Thompson microphysics (Thompson et al., 2006), Kain-Fritsch convection scheme (only in the coarser grid) (Kain, 2004), Monin-Obukhov surface layer, 5-layer thermal diffusion for soil (Skamarock et al., 2005), Rapid Radiation Transfer Model for longwave radiation (Mlawer et al., 1997), and Dudhia scheme for shortwave radiation (Dudhia, 1989).
The simulation shows that, starting from around 18:00 UTC, 15 June, the circulation changes significantly, as the westerly wind component and the wind speed increase, and the flow progressively becomes westerly and then northerly (in agreement with Fig. 2), after a trough crosses the domain and a ridge reinforces over the British islands. The simulation also shows a cold front rapidly moving from north to south across the inner domain, thus responsible for cold air advection in the region; moreover, it suggests the presence of few low clouds, explaining the relative minimum in radiation during the morning, and why the peak in radiation is generally smaller compared to the previous and the next days, as shown in Fig. 5 of Sempreviva and Gryning (2000). However in Sempreviva and Gryning (2000) the presence of clouds for the day here considered is not reported, and we assume that they do not have an influence on the evolution of the boundary layer of 16 June.
The front is followed by an anticyclonic circulation, associated with a significant reduction of humidity. In the in-ner domain, at the station location, the simulated 2 m relative humidity decreases by more than 40 % in 15 h. In the same time interval, i.e., from 06:00 to 21:00 UTC, 16 June, the WRF simulated 1000 hPa temperature increases by about 3 K. The low level warming produces a progressive decrease in the temperature difference between the sea and the air, that explains the observed weak and slightly decreasing turbulent fluxes during the day, as reported in Sempreviva and Gryning (2000) (Fig. 5, lower right panel). The situation of 16 June is hence characterized by a variation in the synoptic conditions, due to incoming of the high pressure responsible for wind rotation.
The region of subsidence nearly corresponds to the area affected by the ridge, thus several hundred km along the main axis and a few hundred across. The vertical velocity decreases with time from values of about 0 to about w subs −0.07 m s −1 , in correspondence with the transit of the ridge. In order to estimate the uncertainty associated with the WRF model simulation, another experiment has been performed by changing the boundary layer scheme (by adopting the local Mellor-Yamada-Janjic closure; Janjic, 2001), the land-surface model (by adopting the so-called community Noah model; see http://gcmd.nasa.gov/records/NOAA_ NOAH.html) -which are the parameterization schemes most relevant for the present study, and the initial conditions (starting date: 12:00 UTC of 15 June instead of 12:00 UTC of 14 June). Simulations show that the subsidence velocity is weakly affected by such changes, as a minimum intensity of ∼ −0.09 m s −1 is extracted from the new experiment. By comparing the WRF model estimate with the vertical velocity shown in the large-scale analysis (e.g., the NCEP/NCAR reanalysis map at 12:00 UTC, 16 June), a similar vertical velocity of w subs = −0.1 m s −1 can be derived. However, the fields provided by the WRF are more accurate, since the large-scale analysis does not take into account the mesoscale effect of the orography of Norway, which, in the presence of a northerly wind, as in the present study, can modify the wind field at low and medium levels in a non negligible way. Summarising, the arrival of the ridge suggests that conditions of subsidence affect the area in the second part of the day, with negative vertical velocity of the order of 0.1 m s −1 .

Detailed evolution of the Marine Boundary Layer: a LES study
Turbulent motions, whose length scale can be much smaller than the horizontal grid spacing employed in mesoscale models, cannot be solved explicitly in mesoscale models, but they can only be parametrized. The impact of these subgrid-scale motions on grid-scale variables is relevant, particularly in the low levels, where they may significantly alter the atmospheric status through mixing. Especially in situations with strong spatial inhomogeneities (e.g., at the land-sea transition zone, where the structure of the ABL flow is more complex due to the abrupt changes in the surface roughness or thermal forcing) and rapid temporal variations, mesoscale models are not able yet to simulate the structure of PBL in all its complexity (De Tomasi et al., 2011), with significant discrepancies among different parametrization schemes. Generally, sub-grid fluxes are parametrized using two categories of closure schemes (Shin and Hong, 2011). The firstorder closure schemes do not include any additional prognostic equation to express the effects of turbulence. In addition to the simple local diffusion, they also consider non-local turbulent mixing in the ABL, which incorporates the contribution of the large-scale eddies to the total flux in terms of a correction to the local gradient of the prognostic variables (e.g., Hong et al., 2006). In the other category of schemes, an additional prognostic equation for the the Turbulent Kinetic Energy (TKE) is considered (Janjic, 2001). Thus, they are classified as TKE closure (one-and-a-half order closure) schemes. The different nature of the two categories of boundary layer schemes (local versus nonlocal turbulent diffusion), may affect the mesoscale flows well as the vertical thermal gradient of the atmosphere . Their advantages and disadvantages were examined in some recent studies (e.g., Shin and Hong, 2011;Rögnvaldsson et al., 2011), exploiting the different parametrization schemes options provided with the WRF model.
Also, for mesoscale simulation with horizontal scales of O(1 km), large eddies begin to blend with the parametrized mixing from the PBL scheme (Stensrud, 2007). As a result, the ability of the actual mesoscale models to reproduce atmospheric phenomena on such scales accurately can be questionable: we are close to the no man's land separating classical PBL schemes from large-eddy simulations (Weisman et al., 2008). For this reason, a numerical model at a finer scale is needed to simulate the marine boundary layer evolution in our study properly.

The LES model
A large-eddy simulation model (Moeng, 1984) is applied to compare numerical predictions with experimental data better. In large-eddy simulations, Eulerian fields are decomposed into their resolved and subgrid components, indicated with an overbar and a prime, respectively. The former are associated with space-time fluctuations whose evolution is directly described by the equation of motions; the latter take place at space-time scales smaller and faster than some cutoff scales and are modeled in terms of a turbulent closure. For instance, for the ith component of the velocity field it holds u i (x, t) = u i (x, t) + u i (x, t).
In the case of atmospheric flows, the governing equations are the incompressible Navier-Stokes equations, with Boussinesq approximation for the velocity field, and the advection-diffusion equations for the scalar fields (potential temperature, θ, and specific humidity, q). The LES model equations, obtained by low pass filtering of the physical equations, are Here the indexes i and j are running over x, y, z, and repeated indexes are retained summed, δ ij and ijk are the Kronecker delta and the Levi-Civita symbol, respectively. Note that x is the stream-wise direction along the geostrophic wind, and y is the span-wise direction, transverse to it. The other variables represent: g, the acceleration due to gravity, directed along z; θ v = θ 0 (1 + 0.61q 0 ), with θ 0 and q 0 the initial surface values of potential temperature and specific humidity, a reference virtual potential temperature; f c the Coriolis parameter; U gj the j component of the geostrophic wind. τ d ij is the deviatoric part of the subgrid scale strain tensor τ ij , which is defined according to τ ij = u i u j − u i u j = u i u j + u iūj + u i u j . The isotropic component of the strain is included in the pressure term: P = p/ρ 0 + τ kk /3 with p the physical pressure and ρ 0 the density of air.
The SGS stress for the scalar θ (or q) is defined as τ The buoyancy term, gθ /θ v (1 + 0.61q), couples the temperature and the humidity fields to the momentum in the Navier-Stokes equations. The closure of the equations is done by modeling the subgrid scale (SGS) terms through the resolved field. Note that condensation is not allowed, hence the described ABL is cloud free.
In the present work, we adopt the dynamic model of Germano et al. (1991). The main advantage with respect to Smagorinsky type of closures (Smagorinsky, 1963;Lévêque et al., 2007) is that, once fixed the cut-off scales, there are no tunable parameters in the SGS scheme. The use of Germano scheme requires the introduction of an additional test filter. Details on the Large-Eddy simulation model and on the SGS closure can be found, respectively, in Moeng (1984), Mazzitelli and Lanotte (2012), and Lanotte and Mazzitelli (2013).
The effect of subsidence is included by adding the largescale term Fφ on the right-hand side of the governing equations www.nonlin-processes-geophys.net/21/489/2014/ Nonlin. Processes Geophys., 21, 489-501, 2014 where w subs (z) is the subsidence velocity, andφ is to be replaced withū x andū y in the momentum equations (3), and withθ andq, in the Eqs. (5) and (6), respectively. Figure 4 shows the profile of the subsidence velocity that we adopted in the LES runs: it is a polynomial curve, which is maximal at the top of the boundary layer, and which goes to zero at the surface and above the inversion (see Appendix B of Siebesma et al., 2003).
To study the effect on subsidence we performed two series of LES, in the same domain, by changing the spatial resolution. The simulated domain is L x × L y × L z = (5 × 5 × 2.2) km 3 . We run the LES model with N x × N y × N z = 64 3 with mesh spacing ∼ 78 m×78 m×35 m, and 128×128×192 grid points, with mesh spacing ∼ 39 m × 39 m × 11 m.
We recall just few details of the LES numerical integration. Momentum and scalar fields equation are discretized on a regular grid in the horizontal planes, where periodic boundary conditions are applied and hence pseudo-spectral methods are used. Dealiasing is performed on horizontal directions applying the 2/3 rule to the nonlinear terms in the equations of motion and to the SGS model terms. A finite-centered difference scheme is adopted along the inhomogeneous vertical direction. Time integration is based on a third-order Runge-Kutta algorithm. A two-dimensional sharp spectral cutoff kernel is applied for both the grid and the test filters in the homogeneous directions. The width of the grid filter is¯ = (¯ x¯ y¯ z ) 1/3 , where, taking into account the dealiasing procedure,¯ i = (3L i )/(2N i ), i = x, y, and¯ z = L z /N z . The width of the test filter is˜ i with i = x, y is about the same of the grid filter. No explicit test filtering is applied along the vertical direction.
Results shown in the sequel are from the run at higher resolution. We verified that the main characteristics of the ABL are unchanged by varying the resolution, which can be easily understood since doubling the resolution does not alter the equilibrium response of the large-eddy simulations as a function of large-scale parameters.
We start the LES runs at t 0 = 08:48 UTC, 16 June 1992. Initial conditions for the temperature and humidity are plotted in Fig. 5: the profiles approximate the experimental ones. The velocity field is initialized with a barotropic geostrophic wind profile, approximating the one obtained by radiosounding (not shown). Scalar equations are forced by the surface-fluxes, w θ S and w q S , that are obtained from the available experimental measurements with 10 min frequency. The surface sensible and latent heat fluxes, see Fig. 5 of Sempreviva and Gryning (2000), do not display the diurnal variation typical of land boundary layer, but they stay positive and almost constant in time during the simulation period: typical values are 0.027 K m s −1 and 0.05 g kg −1 m s −1 , respectively, with the sensible heat flux exhibiting a slight decrease. Hence, we have a weakly convective marine boundary layer.  Other input parameters and simulation variables are summarized in Table 1. Large-eddy simulations are carried out with and without the subsidence term (7) in the equations of motion (3), (5), and (6). According to the experimental observations and the output of the WRF runs, subsidence is included from time t = 12:00 UTC till the end of the run. We fix its maximum intensityw max = 0.07 m s −1 , in agreement with the estimation of the WRF model. Note that the subsidence has a constant profile throughout the runs. Changes in the evolution of large-scale subsidence over time -that are in principle possible -are not taken into account in the present work. Results from the LES runs are compared with radiosoundings, when available. Table 1. Simulation parameters. The symbols indicate: U g the imposed geostrophic wind, fa the free atmosphere lapse rate, z i the boundary layer height, z i /L, with L the Monin-Obukhov length, the stability parameter, ω * the convective velocity and τ * = z i /ω * the convective timescale. Variables z i , z i /L, ω * and τ * are temporal averages from 10:30 to 12:00 UTC, before subsidence is introduced in the simulation.

Mean profiles
We start by plotting the results obtained from the large-eddy simulations for the boundary layer height z i , shown in Fig. 6. The boundary layer depth z i is estimated as the height at which the sensible heat flux is minimal. In the subsidencefree run the height of the boundary layer remains nearly constant in time, in agreement with a sensible heat flux slightly decreasing during the day. Differently, in the presence of subsidence a rapid and intense decrease in z i is observed. Note that, within error bars, the experimental measurements show a trend similar to the one numerically estimated. We recall that the experimental estimates of the BL height reported in Fig. 6 were obtained in Sempreviva and Gryning (2000) ( Fig. 6 of the paper), and roughly correspond to the upper inversion visible from Fig. 7. Moreover, since the WRF runs revealed the passage of a cold front, we can exclude that the observed boundary-layer height evolution might be due to the advection of warm air in the lower troposphere. In Figs. 7 and 8, we compare the plots of temperature and humidity profiles obtained from the LES at 15:30 UTC with the radiosoundings and with profiles obtained from the WRF at about 12:00 and 15:00 UTC (insets). In the LES, subsidence starts at 12:00 UTC, and clearly needs some time before being effective over the whole boundary layer. Hence it is reasonable to compare the radiosoundings with the LES temperature and humidity profiles at 15:30 only.
Few comments can be done. The main feature is the abrupt change in the observed vertical profiles at about 750 m, which appear smoother in the WRF simulations. Up to the first inversion, the LES profiles closely resemble those from the radiosoundings. In particular, the position of the first inversion is correctly reproduced, while there is some discrepancy in the profiles between the first and second inversion. We ascribe this discrepancy to the following facts. The first is that in our runs subsidence is maximal at the mixed layer edge (Siebesma et al., 2003): another possible choice is to define it with an exponential decrease as it is done in Bellon and Stevens (2012). This is unimportant for what concerns  mixing layer properties that we investigate here, while it can change the troposphere status over the ML. Also, we keep subsidence constant in time, while there could have been a slow evolution in the mesoscale conditions. Finally, as common, our ABL is barotropic; however, as experimental observations show, in addition to subsidence, 16 June is characterized by some wind variability. Baroclinicity, added even in the simplest form of an external, time-dependent geostrophic forcing (see, e.g., Zilitinkevich and Esau, 2003;Rizza et al., 2013), could improve our results for a weakly unstable ABL. We comment that in Rizza et al. (2013), the use of geostrophic wind profiles from the WRF improved the prognostic capability of LES in reproducing the wind field pattern in the boundary layer. The Monin-Obukhov length, the friction velocity and the surface fluxes were significantly modified by the inclusion of a baroclinic term in LES equations, while its effect on vertical profiles of temperature and humidity was negligible. It is reasonable to suppose that, also in the present case, characterized by a rapid evolution of largescale patterns, the inclusion of a baroclinic term might affect the simulation results. We leave the investigation along this direction for future work.
It is important to note that the disagreement between the WRF model profiles and observed soundings is larger than for LES, showing that LES represent a useful tool for representing the evolution of the atmospheric boundary layers better. On the other hand, the observed departure of the WRF profiles from the experimental ones can be mainly attributed to the initial and boundary conditions, which are based on ERA-INTERIM re-analysis, whose horizontal grid spacing is very coarse, being about 80 km.
Let us now examine the evolution of the boundary layer and the influence of subsidence on the statistics of the potential temperature and of the specific humidity. As previously remarked, the surface fluxes are positive and almost constant during the day in exam. At the top of the boundary layer, entrainment and subsidence effects compete and it is important to account for their importance.
Starting from Fig. 6, we can measure the entrainment velocity as w e = dz i dt , from the evolution of the boundary layer height in the LES run without subsidence. This gives w e = 0.0052±0.0005 m s −1 . Alternatively, as in Lanotte and Mazzitelli (2013), where the influence of different entrain-ment fluxes on scalar statistics in convective boundary layers was studied, we estimate the entrainment velocity at the top of the boundary layer by means of the equation (Stull, 1976): where θ 0 is a reference potential temperature, d 1 is the difference between the boundary layer height z i and the height of zero heat flux, ( EZ θ ) is the average (over the horizontal directions) temperature difference over the entrainment zone, ( EZ U ) is the difference for the magnitude of the wind, w * is the convective velocity, and u * the surface friction velocity.
In the LES runs, we calculated the entrainment velocity at 12:00 UTC when the subsidence is turned on: the result with this second method is w e (z i ) = 0.0050 ± 0.0005 m s −1 . So on the basis of these estimates, the entrainment velocity results to be an order of magnitude smaller than subsidence velocity: we can hence conclude that the role of entrainment is negligible with respect to that of subsidence in affecting turbulent fluxes at the top of the mixed layer. To account for this, we assume for simplicity a simple slab or bulk model for the MBL (see Stull, 1988): where body source terms have been neglected, and w θ bottom and w θ top are the bottom and top fluxes, respectively. In the above equation, the average θ is taken over the homogeneous directions and over the depth of the mixed layer z i . At the top of the boundary layer, the combined effects of the entrainment and the subsidence influence the turbulent top fluxes: as we have quantified, in our casestudy, the former is negligible in comparison to the latter. Since, because of subsidence, the mixed layer depth decreases, from Eq. (9) we have that the magnitude of temperature time derivative increases, i.e., net MBL warming. In the reasoning we have kept all fluxes positive and slowly varying with respect to z i . Numerical and experimental observations, in Fig. 7, indeed confirm that, in the presence of a slightly decreasing, but positive sensible heat flux at the surface, the MBL warms up because of the subsidence. Moreover, the observed warming is higher than the one obtained in the control LES run, without subsidence.
Concerning the specific humidity, a similar equation can be applied. Since the jump in the specific humidity profile across the entrainment layer is generally negative in diurnal conditions, and because the magnitude of the dry air entrainment flux in fair weather might exceed the surface fluxes, we could expect that a reduced mixed layer produces a net MBL drying. By looking at Fig. 8, we note that the MBL drying takes place in the morning; while, similar to what happens for the potential temperature, the effect of the subsidence is such to lead to a net MBL moistening (stronger than in the control case). Note that entrainment latent heat flux (estimated as 0.021 m s −1 g kg −1 ), that could lead to MBL drying, is half of the surface one, being the entrainment velocity very small in the case study. Hence even in the presence of a slight decrease in sensible surface heat flux, we have a net MBL moistening due to the negligible contribution of entrainment.
To summarize, the large-scale descending current is responsible of a net mixed-layer warming and moistening observable in the second part of 16 June, both in the radiosoundings and in the LES (note that the advection associated with the frontal system has an opposite effect, producing a cooling of the low troposphere). Moreover, subsidence acts to compress the mixed layer as a quasi adiabatic lid, which also leads to the increase in scalar fluctuations, as documented in the next section.

Turbulent statistics
One of the advantages of geophysical modeling by means of large-eddy simulations is that they give access to turbulent quantities scale-by-scale, from the large scale of the motions down to the cut-off scale. Here, we are interested in characterising the effect of subsidence onto the mixed-layer turbulent fluctuations, since these quantities are likely to enjoy a higher degree of universality. Before looking at the outcomes of the numerical simulations for the turbulent spectra of temperature, humidity and vertical velocity, we can try to have a physical intuition in terms of ABL similarity theory.
In the presence of subsidence, and assuming that surface fluxes have a slower time evolution with respect to the time variation of the boundary layer height z i (t), we can expect that the convective velocity scale w * decreases, while the convective temperature scale θ * increases. Hence if we look at variances at fixed values of z/z i , in the presence of subsidence we expect to observe a smaller value for the vertical velocity variance w 2 , and a higher value of the temperature (or specific humidity) variance θ 2 , with respect to the same situation without turbulence.
In Figs. 9 and 10, we plot the one dimensional spectra of vertical velocity component (transverse spectrum), relative humidity and potential temperature, respectively, calculated in the horizontal plane at wavenumber between n h = n 2 x + n 2 y and n h + dn h . Statistical convergence is obtained by analysing 16 equispaced snapshots of the Eulerian fields, between 12:30 and 15:30 UTC, and by averaging over a slab approximately 100 m thick. These spectra are obtained by looking at fluctuations at fixed value z/z i ∼ 0.7, since mixed layer physics is investigated. First, we observe that for the vertical component of the velocity, a reduction in the intensity of the turbulent fluctuations is indeed observed, but limited at the large scale of motion. The crossover between the average spectrum of the case study with subsidence, and that obtained without subsidence, takes place at a scale Fig. 9. Vertical velocity one-dimensional spectra for wavenumbers n h in the horizontal plane, measured from the LES runs with (thick line) and without (thin line) subsidence. Spectra are measured in both cases in the mixed layer, at the fixed height z/z i = 0.7. The dashed line gives the Kolmogorov slope for the inertial range of scale. Fig. 10. The same as Fig. 9, but for the potential temperature and for the humidity fields. The dashed line gives the Kolmogorov slope for the inertial range of scale. L 0 ∼ 800 m, approximately equal to the mixed layer depth in the presence of subsidence. This points to the fact that, in the horizontal plane, subsidence affects the most energetic eddies of scale larger than z i , while leaves unchanged turbulent fluctuations at horizontal scales smaller than or equal to the boundary layer height. This is confirmed by visually inspecting in Fig. 11, a vertical cut of the vertical velocity, measured in the ABL evolution before turning on the subsidence (at 12:00 UTC), and at the end of the run (at about 15:30 UTC), where it is clearly seen that largest vertical structures are dumped.
On the other hand, in agreement with our expectation for the case of temperature (and similarly for the humidity), we observe a net and global increase in the scalars' fluctuation intensity. Note that the integral scale variance of scalar turbulent fluctuations in convective boundary layers can be five to six times larger than that of the vertical velocity (Lenschow and Stankov, 1986).
In Mirocha and Kosović (2010), the effects of subsidence on the stream-wise velocity and temperature spectra were discussed. In particular, it was found that that the main effect of subsidence is to cause a shift in spectral power to higher frequencies, which is more visible on the velocity than in the temperature signals. Here, we observe that the effect of subsidence is to increase scalar fluctuations at any spatial scale as a result of having a shallower convective region, at fixed  Finally, to analyze the turbulence structure further, in Figs. 12 and 13, we plot the evolution in the LES with subsidence and in the control run of the sensible heat flux, and of the vertical velocity variance flux, respectively. They confirm previous observations, namely that turbulent exchanges of scalar fluctuations are enhanced by the action of subsidence, while turbulent transport of vertical velocity variance is reduced.

Conclusions
High-pressure regimes are most likely associated with largescale divergence and subsidence. As we have shown with the case study of a cloud-free marine boundary layer, subsidence can be responsible for shrinking the mixed layer depth by 600 m, from the height of about 1250 m recorded in the morning to that of about 500-600 m in the afternoon. By running a numerical experiment with the WRF model, we have characterized the mesoscale situation during the case study. This is affected by a cold front rapidly moving from north to south across the region of interest. Mesoscale model runs quantified a subsidence velocity following the front of the order of 0.1 m s −1 . Then, by performing two series of large-eddy simulations, we characterized mean and fluctuating field evolution. In particular, we find that the (i) the slowly decreasing sensible heat flux at the sea surface can not be responsible for the boundary layer evolution during the case study, as shown by the control run; (ii) by means of a polynomial profile for the subsidence velocity, we are able to reproduce the mean field evolution of the scalars and hence the observed collapse of the boundary layer height, associated with global air warming; (iii) when looking at turbulent fluctuations, quantified in terms of the second-order moments, we find that subsidence modifies their amplitude and spatial organization. In particular, subsidence dumps vertical motions at scales larger than the ML height, while keeping turbulent fluctuations at smaller scales unchanged. This can be relevant when estimating turbulent kinetic energy budgets scale-by-scale. Scalars exhibit increased turbulent fluctuations at all scales, and a simple similarity argument can be formulated to explain the observation. This is an interesting result of the work, since it points to the fact that subsidence can largely affect scalar turbulent fluctuations in the mixed layer, and not only mean profiles. More generally, our work confirms the importance of having an accurate estimate of subsidence velocity when modeling the atmospheric boundary layer. It is clear however that major improvements will be achieved only in the presence of a large data set of reliable vertical velocity measurements at synoptic and sub-synoptic scales, under different meteorological conditions. This implies performing studies of the statistical effects caused by varying subsidence forms and amplitudes, in a way similar to what has been previously done, e.g., in Sorbjan (1996), where the potential temperature lapse rate variation was examined, or as in Zilitinkevich et al. (2006Zilitinkevich et al. ( , 2007, where the role of the Brunt-Väisälä frequency and the Coriolis parameter on the height of boundary layer conditions was systematically analyzed under different stability conditions. Here, we focused on a specific situation, leaving a methodical investigation of the effects of subsidence to a further publication. In particular, beyond amplitude variation, it will be important to test the time and spatial variability of the subsidence velocity, beyond simple representations such as the one adopted in this study.