Description and Implementation of a Mixed Layer Model (mxl, V1.0) for the Dynamics of the Atmospheric Boundary Layer in the Modular Earth Submodel System (messy)

We present a new submodel for the Modular Earth Submodel System (MESSy): the MiXed Layer (MXL) model for the diurnal dynamics of the convective boundary layer, including explicit representations of entrainment and surface fluxes. This submodel is embedded in a new MESSy base model (VERTICO), which represents a single atmospheric column. With the implementation of MXL in MESSy, MXL can be used in combination with other MESSy submodels that represent processes related to atmospheric chemistry. For instance, the coupling of MXL with more advanced modules for gas-phase chemistry (such as the Mainz Iso-prene Mechanism 2 (MIM2)), emissions, dry deposition and organic aerosol formation than in previous versions of the MXL code is possible. Since MXL is now integrated in the MESSy framework, it can take advantage of future developments of this framework, such as the inclusion of new process submodels. The coupling of MXL with submodels that represent other processes relevant to chemistry in the atmospheric boundary layer (ABL) yields a computationally inexpensive tool that is ideally suited for the analysis of field data, for evaluating new parametrizations for 3-D models, and for performing systematic sensitivity analyses. A case study for the DOMINO campaign in southern Spain is shown to demonstrate the use and performance of MXL/MESSy in reproducing and analysing field observations.


Introduction
In atmospheric chemistry, various types of models are used for studies on different spatio-temporal scales, ranging from box models representing a single point in space (e.g.Sander et al., 2011) to regional and global 3D models (e.g.Jöckel

Correspondence
to: R. H. H. Janssen (ruud.janssen@mpic.de) et al., 2006).Here we describe the implementation of a MiXed Layer model for the dynamics of the atmospheric boundary layer (ABL) and a generic 1D basemodel called 35 VERTICO (VERTIcal COlumn) in the Modular Earth Submodel System (MESSy, Jöckel et al., 2010).The implementation of the latter is necessary because of the structure of the MESSy framework: in this framework, a basemodel determines the basic configuration of the coupled model, which 40 can be a box, 1D or 3D model, and the individual processes are represented by basemodel independent process submodels that are coupled to the basemodel through an interface layer (Jöckel et al., 2010).
MXL represents the dynamics of the convective boundary 45 layer at diurnal time-scales with two vertically changing layers, one for the ABL and one for the free troposphere (FT).
In combination with MESSy submodels that represent the other processes relevant to atmospheric chemistry, it provides a missing link between box and 3D atmospheric chemistry 50 models, at a spatio-temporal scale typical for measurements of atmospheric chemistry.It is especially suited for studying the diurnal evolution of chemical species that have chemical lifetimes similar to the time scales of atmospheric mixing, so for which the effects of chemistry and dynamics should be 55 studied simultaneously (Fig. 1).Since the typical time-scale of mixing of a convective boundary layer is about 15 min., the evolution of species with a chemical time scale of minutes to several hours (e.g.O 3 , isoprene, NO, and NO 2 ) will be to a similar extent driven by both dynamics and chemistry.

60
The MXL model has been developed in the 1960's and 1970's (Lilly, 1968;Tennekes, 1973) to study the dynamics of clear and cloudy boundary layers and was first applied to atmospheric chemistry by Vilà-Guerau de Arellano et al. (2009), as MXLCH (MXL-CHemistry).Afterwards 65 MXLCH has been applied in studies of the diurnal evolution of gas-phase chemistry (Vilà-Guerau de Arellano et al., 2011;Ouwersloot et al., 2012;Van Stratum et al., 2012) and organic aerosol (Janssen et al., 2012(Janssen et al., , 2013) ) in mid-latitude, tropical and boreal areas.The model explicitly accounts for entrainment (ABL-FT exchange) and parametrizes turbulence in a simplified way.Although MXL uses simple parametrizations for turbulence and entrainment, the ABL dynamics as simulated by the model compare well with results from a turbulence resolving Large Eddy Simulation model (Pino et al., 2006) and with vertical profiles as observed from radio soundings (Ouwersloot et al., 2012).
With the implementation of MXL in MESSy (henceforth MXL/MESSy) as dynamical core and the subsequent coupling to the submodels that represent the other processes relevant to the evolution of chemical species in the troposphere, a versatile boundary layer chemistry model is created.MXL/MESSy can be used for different types of studies.First, it is well-suited to support the interpretation of field observations in an eulerian framework, e.g.observations from a tower.After obtaining a best fit with the observational data, a budget analysis can be performed, dividing the total tendency of species into the contributions of gas-phase chemistry, ABL dynamics, emission, deposition and/or gas/particle partitioning.The chemistry term can be further subdivided into the total production and loss terms, and the contributions of individual reactions to these terms.
Secondly, MXL/MESSy is well suited for performing systematic sensitivity analyses, in which the parameter space is explored and through which non-linearities can be studied.
Because it is computationally inexpensive, detailed and systematic sensitivity runs are possible.In that way, it can serve as a test bed for new parametrizations, especially in combination with a direct comparison to field data.It thus forms a link between theoretical/lab study results and 3D model applications.Finally, it can be used in theoretical/conceptual studies, for a qualitative evaluation of feedbacks and forcings (e.g.Vilà-Guerau de Arellano et al., 2012;Janssen et al., 2012).
In Sect. 2 we describe the MXL model and give its governing equations, in Sect. 3 we focus on the implementation of the MXL model in MESSy and in Sect. 4 we show an example application of MXL/MESSy using observations from the DOMINO campaign.

MXL submodel description
The version of the MiXed Layer model that we implement here has been developed in several stages: Lilly (1968) developed the mixed-layer equations and Tennekes (1973) the turbulence closure that we apply in this work.The model for the dynamics of the convective boundary layer that resulted from this work was later coupled to an atmospheric chemistry module by Vilà-Guerau de Arellano et al. ( 2009) and a land surface scheme by Van Heerwaarden et al. (2009, 2010).While the equations in this section have been published in one or more of the aforementioned references, we here present an overview of the specific and complete set of equations that are part of the implementation of MXL in MESSy.
A detailed derivation of the MXL governing equations can be found in Vilà-Guerau de Arellano et al. (2015). 125 The mixed layer theory states that under convective conditions, strong turbulent flow, driven by the surface heat fluxes, causes perfect mixing of quantities over the entire depth of the ABL (Tennekes, 1973;Vilà-Guerau de Arellano et al., 2009, 2015).Therefore, scalars and reactants in the convec-130 tive boundary layer are characterized by a well-mixed vertical profile over the whole depth of the ABL.The transition between the well-mixed ABL and the free troposphere is marked by an infinitesimally thin inversion layer.Typical profiles of potential temperature (θ), specific humidity (q) 135 and chemical species in the mixed-layer model are shown in Fig. 2. It shows that the ABL is represented by bulk values of scalars and chemical species mixing ratios, and capped by a infinitesimally thin layer over which these strongly change.Above this inversion layer lies the residual layer (the remain-140 der of the convective boundary layer of the previous day) during the early morning and the free troposphere later during the day, which are characterized in the model by a lapse rate for heat and moisture, and a concentration profile of chemical species which is constant with height.

145
In addition to the local surface heat fluxes, the boundary layer dynamics can be influenced by large-scale atmospheric flows which act as external forcings on the ABL development: for instance, Janssen et al. (2013) and Pietersen et al. (2014) studied cases for which the observed ABL dynam-150 ics could only be reproduced if the influence of advection, caused by meso-scale flows, and/or subsidence, caused by a high-pressure system, were taken into account.Therefore, advection and subsidence can be prescribed as forcings to MXL.

Governing equations for the heat budget
The main variable in the dynamics of the convective boundary layer is the potential temperature (θ), since it is used to quantify the convective turbulence.The evolution of the potential temperature of a dry convective boundary layer is driven by the heat input at the surface (surface heat flux), at boundary layer top (entrainment heat flux) and at the lateral boundaries (advection): Equation ( 1) is the result of a vertical integration of the 1-dimensional equation of the heat budget and we have assumed that the vertical profile of θ is in quasi-steady state (Lilly, 1968), which causes a linear vertical gradient of the 160 heat flux.In this equation, w θ s is the kinematic heat flux at the surface, which is related to the sensible heat flux (SH) as: w θ s = SH/(ρ • c p ), with ρ the density of air and c p the specific heat of air.The entrainment process, represented by w θ e , is defined as the process whereby air from the FT is mixed into the mixed layer, and it is therefore related to the θ jump at the inversion.The evolution of θ thus equals the input of heat into the ABL at the surface and due to entrainment over the ABL height h.Additionally, a heat advection term (adv θ ), which is a large-scale forcing on the boundary layer dynamics, can be prescribed to the model.
When calculating the entrainment flux, we assume that the transition from the ABL to the FT, the inversion, is represented by a sharp discontinuity, namely the zero-order jump closure (ZOJ; Tennekes, 1973).ZOJ closure defines this jump as ∆θ = θ FT − θ over an infinitely thin inversion layer.In this ZOJ approach, the entrainment flux is the product of the entrainment velocity w e (defined positive in the upward direction) and the potential temperature jump ∆θ at the inversion.First-order closure approaches exist, which include the explicit representation of the depth of the entrainment zone, but the ZOJ approach already gives satisfactory results (Pino et al., 2006).
The equation for the potential temperature entrainment flux reads: In our model, we calculate the subsidence velocity as: where ω represents the large scale vertical velocity that is a function of the horizontal wind divergence in s −1 .It can be thought of as the fraction with which the ABL is pushed down per second, due to large-scale subsiding air motions.Therefore, w l is per definition negative.By rewriting Eq. ( 2), we obtain an expression for the boundary layer growth ( ∂h ∂t ) as a function of the entrainment flux w θ e , the potential temperature jump (∆θ) and the subsidence velocity (w l ).It reads: This equation states that the mixed layer grows by entrainment of warm air from the free atmosphere (w e > 0) and that it is opposed by the vertical subsidence velocity (w l < 0) driven by high pressure situations.This growth is limited by the presence of a stably stratified layer defined by a potential temperature jump on top of the convective ABL.This jump at the inversion, represented by ∆θ, changes during the growth of the mixed-layer.Consequently, depending on a positive or negative tendency of ∆θ, the ABL growth increases or decreases.
As Eq. (2) shows, the entrainment flux for heat depends on the entrainment velocity and the potential temperature jump at the inversion.It is therefore necessary to obtain a prognostic equation for the potential temperature jump.This equation reads: where γ θ is the lapse rate (change with height) of θ in the FT.
The final set of equations, which describes the heat budget 200 in the diurnal atmospheric boundary layer, is therefore composed by Eqs. ( 1), ( 4) and (5).The three prognostic variables are θ , ∆θ and h.
In this set of equations, four variables still remain as unknowns.The surface heat flux w θ s is the result of landatmosphere interactions, but can be prescribed based on observations.The external variables w l and γ θ represent the influence of the free tropospheric conditions on the ABL development and their values are imposed on the mixed-layer model.Finally, the entrainment heat flux w θ e remains.Here, we assume an important closure and we relate the entrainment flux to the surface heat flux as: where the coefficient β can be imposed as a constant or as a parameter depending on the thermodynamic characteristics 205 at the inversion, for instance the presence of shear (Tennekes and Driedonks, 1981;Pino et al., 2003;Conzemius and Fedorovich, 2006).In our model, we assume a value of β equal to 0.2 (Pino et al., 2003), which physically means that the contribution of entrainment to the heat budget equals 20 % 210 of the contribution of the surface heat flux.The potential temperature is normally underestimated if this contribution is neglected.As we will show later on when introducing the moisture budget, Eq. ( 6) needs to be modified to make entrainment dependent on the buoyancy flux.

Governing equations for the moisture budget
By adding the moisture budget to the heat budget, we incorporate the effect of moisture on the buoyancy of air parcels and therewith complete the configuration of the thermodynamic variables in the ABL.The inclusion of the moisture effects on the dynamics of the ABL requires the introduction of two new equations.The first one (similar to Eq. 1) is the evolution of the mixed-layer specific humidity q : ∂ q ∂t = surface moisture flux w q s h − entrainment moisture flux w q e h + moisture advection adv q (7) where w q s represents the surface specific moisture flux, w q e the entrainment flux of moisture and adv q an optionally prescribed moisture advection term.The specific moisture flux is related to the latent heat flux (LE) following w q s = LE/(ρ • L v ), with L v the latent heat of vaporization.Similarly to Eq. (2), we represent the entrainment flux as: This equation relates the dynamics of the boundary layer growth, represented by the entrainment velocity, with the specific humidity jump at the interface between the ABL and the FT.Equation ( 8) requires an additional equation for the temporal evolution of the jump of q at the ABL-FT interface (∆q).It reads: where γ q is the lapse rate of q in the FT.At this point, we need to introduce a new variable, the virtual potential temperature, which accounts for both changes in the heat and moisture budgets and is used to quantify buoyancy.It is defined as: The virtual potential temperature thus accounts for the effect of water vapour on the density of air.Since moist air is less dense than dry air at the same conditions of temperature and pressure, θ v is always greater than the actual temperature, but only by a few degrees.The turbulent transport of this variable, the buoyancy flux, combines in one quantity the information of the potential temperature flux and the specific moisture flux.It reads: w θ v = w θ + 0.61 θ w q + q w θ + w θ q (11) This buoyancy flux expresses the production of turbulent kinetic energy in the ABL due to density differences.Turbulence driven by shear (mechanical turbulence) on the mixed-layer thermodynamic equations is not dealt with in our model.
The buoyancy flux directly enters in the ABL growth formulated in Eq. ( 2).Therefore, we rewrite Eq. ( 2) in the definitive form as: where ∆θ v is expressed in terms of the characteristics of the θand q-budgets: By introducing the buoyancy flux as the driver of the turbulent process in the determination of the boundary layer growth, we complete the main framework of our model formulation based on mixed-layer theory.

Governing equations for the horizontal wind budget
As the last part of the mixed-layer dynamics, we introduce the horizontal wind, or the momentum budget (Conzemius and Fedorovich, 2006).The wind speed is important for calculating surface-atmosphere exchange, and appears for instance in the calculations of the aerodynamic resistance which governs evapotranspiration (Sect.2.7.2) and dry deposition.

245
Similar expressions as for θ and q can be used for the two horizontal components of the wind speed.These equations also contain the Coriolis force that takes into account the rotation of the Earth.This gives another four equations: where u and v are the mixed-layer wind velocities in x and y direction, ∆u and ∆v are the jumps of these two variables, γ u and γ v are the lapse rates in the FT and f c is the Coriolis parameter.The free tropospheric wind velocities are assumed to be equal to the geostrophic wind at their height.Also for the wind budget, we assume zero-order closure: For applications such as in Eq. ( 42), u and v are combined to form the total horizontal wind speed: where w * is the free convection scaling velocity:

Governing equations for chemical species
The mixed-layer model allows us to investigate the influence of the boundary layer dynamics on reactive atmospheric compounds during daytime.It is important to stress that in mixed-layer theory, we assume that species are mixed instantaneously as soon as they are chemically produced, emitted or entrained from the FT.In other words, the MXL model acts as a reactive chamber with the additional advantage of accounting for boundary layer growth and ABL-FT exchange.
As for scalars, the inclusion of reactive species requires the introduction of two additional equations for each species.The expression for the evolution of the generic species C is similar to that for potential temperature (Eq. 1) and moisture (Eq.7), but includes a term for the chemical transformation (Vilà-Guerau de Arellano et al., 2009): By solving Eq. ( 24), we determine how C varies over time as a function of emission/deposition processes at the surface (represented by the term (w C ) s ), the dynamic effects (h and (w C ) e ) and the chemical transformation (S C ).Note that the surface fluxes of chemical species (emission or deposition) can be either prescribed (Sect.2.5) to MXL/MESSy or calculated by other MESSy submodels (Sect.3.1, Table 1).Advection of tracers is not considered, because horizontal gradients in species concentrations, which are necessary to calculate advection, are generally not constrained by observations.We are planning to implement the TNUDGE submodel (Kerkweg et al., 2006b) in the future to allow the relaxation of chemical species concentrations towards observations, which will also compensate for sources or sinks that are the result of advection of air masses with different concentrations of species.
The flux at the top of the boundary layer is represented in the same way as the entrainment flux for buoyancy, moisture and wind.For C it reads: By representing the exchange of C as the product of the entrainment velocity and the jump of C, we account for both the dynamics and chemistry, which together determine the exchange between the ABL and FT.Equation ( 25) requires an additional prognostic equation to solve the evolution of ∆C: where C FT is the concentration of C in the FT with ∂C FT /∂t driven by chemical production and loss only in MXL/MESSy.
By assuming instantaneous mixing, we assume that the intensity of segregation of chemical species due to possible 280 inefficient turbulent mixing is zero.It has been hypothesized that this segregation might be important for explaining model-measurements discrepancies for certain species, e.g.isoprene and the hydroxyl radical (Butler et al., 2008;Pugh et al., 2010).Later studies, however, showed that the effects 285 of imperfect mixing are only small, with a maximum reduction of the effective rate constant for the isoprene and OH by less than 15 % (Ouwersloot et al., 2011;Pugh et al., 2011), and not by as much as 50 %, which was necessary to explain the measurement-model discrepancy.Moreover, Large 290 Eddy Simulation results show well-mixed profiles throughout the convective boundary layer for species like isoprene and OH (Vilà-Guerau de Arellano et al., 2009Arellano et al., , 2011;;Ouwersloot et al., 2011) under homogeneous land surface conditions.A parametrization for the effect of the segregation of 295 species on the effective rate constants was developed by Vinuesa and Vilà-Guerau de Arellano ( 2003) and can in principle be implemented to account for this effect.

Simple emission functions
Besides online calculation of emissions, MXL/MESSy provides the possibility to prescribe emissions, following simplified diurnal emission patterns.These include constant, sinusoidal and cosine shaped emission fluxes, respectively: F C is the (maximum) daily emission flux (ppbms −1 ), where t is the time of day and t d the day length and For these functions, F C can be set by the user in the MXL 300 namelist (an example namelist is included in the Supplement).

Radiation model
In this section, we describe a simple radiation model that is implemented along with MXL/MESSy.The MESSy radi-305 ation submodel RAD is intended to be implemented in the future as a second option.
The net radiation at the surface R n is defined as the sum of the four radiation components: where S in and S out are the incoming and outgoing components of the short wave radation, and L in and L out the incoming and outgoing component of the long wave radiation.The incoming short wave radiation S in is calculated by: where S o is the constant solar irradiance at the top of the atmosphere, set to 1365 Wm −2 , T r the net sky transmissivity, that takes into account the influence of radiative path length and atmospheric absorption and scattering using: Through the solar elevation angle Ψ, both expressions depend on geographical location, day of the year and time of the day.Ψ is calculated using: (33) where t UTC is the universal time coordinate (UTC) and t d is the diurnal period of 24 h.The latitude φ (positive north of Equator) and longitude λ e (positive east of Greenwich) define the geographic location.Variable δ s is the solar declination, which is a function of the day number: where Φ r is the tilt of the Earth's axis relative to the elliptic, equal to 0.409 rad.The Julian day is represented by d and d r is 173, the day of the summer solstice.The number of days in a year d y is set to 365.
The outgoing shortwave radiation depends on the surface albedo: where α is the surface albedo.

315
The outgoing long wave radiation is calculated using the the Stefan-Boltzmann Law: where s is the surface emissivity, set to 1, σ SB the Stefan-Boltzmann constant equal to 5.67 × 10 −8 Wm −2 K −4 , and T s the surface temperature (K).
The incoming long wave radiation is computed using the same expression, but here it uses the temperature at the top of the surface layer T sl .This temperature is acquired by converting the potential temperature of the mixed-layer θ to absolute temperature using the height of surface layer top, which we define as 10 % of the boundary-layer height h (Jacobs and De Bruin , 1992;Van Heerwaarden et al., 2010): z sl = 0.1h.This gives the expression: where the atmospheric emissivity a is 0.8.

320
A land surface model (Van Heerwaarden et al., 2009, 2010;Van Heerwaarden, 2011) is included as an option in MXL/MESSy, which enables the interactive calculation of surface heat fluxes.With this inclusion, MXL/MESSy can be used to simultaneously and interactively calculate the ex-325 change of energy (sensible heat flux), water (latent heat flux), wind (momentum flux) and chemical species (emission and deposition) between the land surface and the ABL.In that way, a fully online coupled land surface-ABL chemistry model is obtained, in which all exchanges between the land 330 surface and the ABL at the diurnal time-scale are internal variables of the coupled system.This means that only forcings (drivers external to the system at the appropriate time scales) are prescibed to the model.These includes land surface forcings (e.g.leaf area index (LAI), roughness length, 335 soil moisture) and large-scale meteorological forcings on the ABL (e.g.incoming solar radiation, FT conditions).

The surface energy balance
The surface energy balance (SEB) forms the basis for the land surface model.It relates net radiation R n (Eq.30) to the surface heat fluxes: where SH is the sensible heat flux, LE is the latent heat flux and G is the ground heat flux.

340
The sensible and latent heat flux, which together supply the energy to the atmosphere that drives turbulent convection, are a function of land surface and atmospheric properties: where r a is the aerodynamic resistance, r s is the surface resistance, θ s and θ are the potential temperatures of the surface and the mixed-layer respectively, q sat (T s ) is the saturated specific humidity inside the canopy and q is the mixedlayer specific humidity.

345
The Penman-Monteith equation (Monteith, 1965) is then used to relate the evaporation flux to radiation, temperature, humidity, and aerodynamic and surface resistance: where dqsat dT , which is the slope of the saturated specific humidity curve with respect to temperature, and q sat are evaluated using the atmospheric mixed-layer temperature at the top of the atmospheric surface layer, defined as 10% of h.
The exchange of momentum and heat between the atmosphere with the surface is related to the strength of the turbulence.The aerodynamic resistance r a , that appears in Eqs. ( 39) and ( 40) is the inverse of this turbulent intensity and defined as: where C H is the drag coefficient for heat and U is the mixed-layer wind speed (Sect.2.3).It shows that stronger winds lead to lower aerodynamic resistance and enhanced turbulent exchange.
The turbulent exchange also depends on the drag, which is a function of the stability of the surface layer: a growing instability in the surface layer leads to stronger convection and mixing.The drag coefficient is calculated following: with: where κ is the von Kármán constant, z 0m and z 0h are the roughness lengths for momentum and heat, respectively, z sl is the depth of the atmospheric surface layer of 0.1h, L is the Monin-Obukhov length and Ψ M and Ψ H are the integrated stability functions for momentum and heat taken from Beljaars (1991).To find the values of L that are required in this function, the following implicit function is solved using a Newton-Raphson iteration method: where Ri B is the bulk Richardson number and θ v,s and θ v are the virtual potential temperatures of the surface and the mixed-layer atmosphere respectively.The calculations of both resistances r veg and r soil follow the method of Jarvis (1976) and employed similarly as in the ECMWF global forecasting model.The vegetation resistance is based on the following multiplicative equation.
where r s, min is the minimum surface resistance, LAI the leaf area index of the vegetated fraction, f 1 a correction function depending on incoming short wave radiation S in , f 2 a function depending on soil moisture w, f 3 a function depending on vapor pressure deficit (VPD) and f 4 a function depending on temperature T .Of these correction functions, the first three originate from the ECMWF model documentation and the fourth from Noilhan and Planton (1989) and are defined as: = min 1, 0.004S in + 0.05 0.81(0.004S in + 1) (46) where w wilt is the volumetric soil moisture at wilting point, w fc is the volumetric soil moisture at field capacity and g D is a correction factor for vapor pressure deficit that only plays 360 a role in high vegetation.
The soil resistance depends on the amount of soil moisture in the layer that has direct contact with the atmosphere.

Evapotranspiration calculation
The total evapotranspiration consist of three parts: soil evaporation, leaf transpiration and evaporation of liquid water from the leaf surface.The total evapotranspiration is therefore proportional to the vegetated fraction of the land surface: where LE veg the transpiration from vegetation, LE soil the evaporation from the soil and LE liq evaporation of liquid water.The fractions that are used are c veg which is the fraction of the total area that is covered by vegetation and c liq , which is the fraction of the vegetated area that contains liquid water.Since c liq is not constant in time, it is modeled following: where W max is the representative depth of a water layer that can lay on one leaf and W l the actual water depth.The evolution of W l is governed by the following equation: where ρ w is the density of water.
A soil model is available for the calculation of soil temperature and moisture changes that occur on diurnal time scales, and that therefore affect the surface heat fluxes.It is a forcerestore soil model, based on the model formulation of Noilhan and Planton (1989) with the soil temperature formulation from Duynkerke (1991).The soil model consists of two layers, of which only the thin upper layer follows a diurnal cycle.The soil temperature in the upper layer (T soil1 ) is calculated using the following equation, where the first term is the force term and the second the restore term: where C T is the surface soil/vegetation heat capacity, G the soil heat flux already introduced in the SEB, τ the time constant of one day and T soil2 the temperature of the deeper soil layer.The soil heat flux is calculated using where Λ is the thermal conductivity of the skin layer.The heat capacity C T is calculated following Clapp and Hornberger (1978) using: where C T, sat is the soil heat capacity at saturation, w 2 the volumetric water content of the deeper soil layer and b an empirical constant that originates from data fitting.The evolution of the volumetric water content of the top soil layer w 1 is calculated using: where w 1,eq is the water content of the top soil layer in equilibrium, d 1 is a normalization depth and where C 1 and C 2 are two coefficients related to the Clapp and Hornberger parameterization (Clapp and Hornberger, 1978), that are calculated using: (58) where C 1,sat and C 2,ref are constants taken from Clapp and Hornberger (1978).The water content of the top soil layer in equilibrium is: with a and p two more fitted constants from Clapp and Hornberger (1978).Finally, we introduce parametrizations for the momentum fluxes at the surface (Stull , 1988) and the friction velocity (u * ), which are required to calculate the horizontal wind budget in the mixed layer (Sect.2.3).
where C M is the drag coefficient for momentum, defined as: 3 Implementation in the MESSy structure MESSy is an interface to couple submodels of earth system processes, which offers great flexibility to choose between different geophysical and -chemical processes.In the 375 first version of MESSy, only the General Circulation Model ECHAM5 could be used (Jöckel et al., 2005).The second round of development also enabled the use of basemodels with different dimensions (Jöckel et al., 2010).For the current implementation, we used MESSy version 2.50.

380
Here, we give a brief overview of the MESSy structure, more details can be found in Jöckel et al. (2005Jöckel et al. ( , 2010)).In MESSy, each FORTRAN95 module belongs to one of the following layers: -Base Model Layer (BML): defines the domain of the 385 model, which can be box, 1D or 3D.This can be a complex atmospheric model, for instance a General Circulation Model like ECHAM5 (Jöckel et al., 2005), but in our case it consists of two boxes stacked on top of each other (Fig. 3).

390
-Base Model Interface Layer (BMIL): this layer manages the calls to specific submodels, data input and output, and data transfer between the submodels and the base model.Global variables are stored in structures called 'channel objects'.

395
-SubModel Interface Layer (SMIL): the submodelspecific interface layer collects relevant data from the BMIL, transfers it to the SMCL and sends the calculated results back to the BMIL.The SMIL contains the calls of the respective submodel routines for the initializa-400 tion, time integration, and finalizing phase of the model.
-SubModel Core Layer (SMCL): this layer contains the code for the base model-independent implementation of the physical and chemical processes or a diagnostic tool.
For the implementation of MXL, a generic 1D basemodel is created in MESSy, called VERTICO.VERTICO contains calls to modules for time and tracer management, the time loop which integrates the model equations and in VERTICO tracer concentrations are updated each timestep, combining the tracer tendencies from each active submodel.It is de facto a 3D basemodel in which the horizontal resolution has been reduced to a single grid box, to facilitate the submodel coupling within the MESSy framework.This also facilitates the possible development of a column model that includes more vertical levels in both boundary layer and free troposphere.
The current implementation of VERTICO consists of two domains (see Fig. 3): the lower one represents the wellmixed boundary layer during daytime, as represented by the MXL equations, and the upper box contains a simplified description of the free troposphere on top of the ABL.The FT is represented as a well-mixed box in which only chemistry and aerosol partitioning takes place.We have fixed the top of the FT to 10 km, so changes in the ABL height only slightly influence the volume of the FT box.FT values of temperature, moisture and pressure that are needed for the calculations of FT reaction rates are taken at the inversion.The prognostic equations for the ABL dynamics and the land surface scheme are integrated by an Euler forward solver.

Coupling with other MESSy submodels
Through the MESSy framework, three types of submodels are coupled: the first type are the generic submodels, which constitute the model infrastructure that is not directly related to actual physical processes, such as time, tracer and data management.The second type are the process submodels, which represent the individual processes that contribute to the evolution of chemical species in the ABL (see Eq. 24 and Fig. 4).Finally, there are diagnostic submodels, which are used when additional post-processing of simulation output is needed.Table 1 shows all generic, process and diagnostic submodels that are currently coupled to form MXL/MESSy, and their tasks.In principle, thanks to such an implementation, any MESSy submodel can be used in MXL/MESSy.
The modular nature of the MESSy interface allows full flexibility: processes can easily be switched on and off through switches in a namelist.Emission and deposition of species takes place only in the lower box, while chemistry and gas/particle partitioning take place in both ABL and FT.The chemical production and loss of a species is calculated in the MECCA submodel (Sander et al., 2011), which uses the Kinetic PreProcessor (KPP; Sandu and Sander, 2006) for the numerical integration of the chemical reaction mechanism.MECCA uses photolysis rates calculated by JVAL (Sander et al., 2014).
For the lower boundary conditions, there are two possibilities.On one hand, MXL/MESSy offers the possibility to use interactive emission (via the ONEMIS, Kerkweg et al., 2006b andMEGAN Guenther et al., 2006 submodels), dry deposition (DDEP Kerkweg et al., 2006a) and land surface parametrisations (Sect.2.7).In these submodels, land 460 surface-ABL exchange is calculated as function of land surface and ABL characteristics, like stomatal resistance and air temperature.On the other hand, it is possible to prescribe emissions and surface heat fluxes following simplified functions for the users who like to keep full control 465 over the boundary conditions of the model (Sect.2.5).In that way, MXL/MESSy can for instance be used to evaluate the sensitivity of the chemistry in the ABL to uncertainties in emission estimates.Further, OFFEMIS (Kerkweg et al., 2006b) allows for the extraction of emission time series from 470 an emission database, which is read by IMPORT.Additionally, the organic aerosol submodel ORACLE (Tsimpidi et al., 2014), allows for the representation of the organic aerosol composition and evolution.

DOMINO case study 475
To evaluate MXL/MESSy, we revisited the case study of Van Stratum et al. (2012), which is based on observations from the DOMINO campaign in the south of Spain.Van Stratum et al. (2012) selected one day from this campaign, the 23th of November 2008, as a case study for dis-480 entangling ABL dynamics and chemistry.On this particular day the influence of synoptic scale flows on the ABL dynamics was relatively small, there were no clouds, and wind speeds and background pollution levels were low.This assures that local processes (ABL dynamics, chemistry, emis-485 sion and deposition) were the main drivers of the observed chemistry.Moreover, this day was during an intensive observation period, assuring that both ABL dynamics and chemistry are well-characterized by observations.MXL/MESSy, with initial conditions as in Table A1, rep-490 resents the dynamics of the boundary layer well, as compared to observations from a tower and radio soundings (Fig. 5).
Since the equations and the initial and boundary conditions are equal to those of Van Stratum et al. (2012), the dynamics of θ , q and h are identical to those shown in their Fig. 3.

495
As a result of the positive heat fluxes (Fig. 5d), the initial potential temperature inversion is broken after 09:00 LT, and the boundary layer starts to grow rapidly.During this period of strong ABL growth, air from the FT is entrained, which causes (1) a strong increase in θ , since warm air 500 from the FT is entrained and (2) a decrease in q despite the positive latent heat flux (see Table A1), since dry air from the FT is entrained into the ABL.After 13:00 LT, the ABL growth slows down, and the effect of entrainment on scalars and chemical species becomes smaller.The correct repre-505 sentation of the boundary layer dynamics ensures that we also simulate the effect of ABL dynamics on the evolution of chemical species well (through h and w e in Eqs.24 and 25).
R. H. H. Janssen and A. Pozzer: Boundary layer dynamics with MXL/MESSy Figure 6 shows the diurnal evolution of several gas-phase chemical species, from measurements and from three model runs: the MXL/MESSy code with MIM2 chemistry (Taraborrelli et al., 2009), with initial conditions as in Van Stratum et al. ( 2012) (SameIC), the MXL/MESSy code with MIM2 chemistry (Taraborrelli et al., 2009), with initial conditions that gave the best fit with the measurements (BestFit), the MXLCH code as in Van Stratum et al. (2012).
The MXL/MESSy SameIC and MXLCH simulations where performed with identical initial conditions as listed in Table B1.In the MXL/MESSy SameIC run we overestimate NO, NO 2 and isoprene and underestimate O 3 , HO 2 and H 2 O 2 .However, these results are influenced by the surface fluxes (emission and deposition) and the mixing ratios in the FT of several species.Since these were not measured, this gives us some degrees of freedom to adjust them, within realistic bounds, to obtain an improved fit of MXL/MESSy with the observations.Figure 6 also shows the MXL/MESSy BestFit results (with initial conditions as in Table B1).It shows that MXL/MESSy with the MIM2 chemical mechanism is able to reproduce the diurnal dynamics of the main gas-phase chemical species well.While it reproduces the timing and peak values of OH and HO 2 well, the mixing ratios of both radicals are underestimated in the early morning and late afternoon.Also H 2 O 2 is underestimated from 10.00 LT to the late afternoon, which is related to the underestimation of HO 2 .
For this case study under conditions of low NOx and isoprene concentrations, MXL/MESSy with MIM2 chemistry and MXLCH with its reduced chemical scheme perform equally well, given the degrees of freedom available to tune the boundary conditions for the chemical species for each of these chemical schemes.
To give more insight in the diurnal evolution of chemical species, it is useful to analyse the contributions of the different processes to the total tendency (Eq.24).The modular nature of MXL/MESSy facilitates the evaluation of the contributions of the individual processes, thanks to the diagnostic tools present in the MESSy framework.
In Fig. 7, an example is shown for O 3 , based on the optimally tuned (BestFit) simulation.The upper panel shows how the total O 3 tendency is built up from the contributions of chemistry, entrainment and dry deposition (although the latter is set to zero), and the lower panel shows the chemistry term and how gas-phase and photolysis reactions contribute to the loss and production of ozone in the ABL.The contribution of entrainment to the O 3 budget is strongest during the morning when the ABL grows rapidly and a large quantity of ozone-rich air is mixed in into the ABL from the FT.
Except for the early morning and late afternoon, the sum of the chemical production and destruction of O 3 is positive, and the net O 3 production peaks around noon when photochemistry is strongest.Note that while we chose to show total production and loss terms for O 3 here, it is possible to 565 split up these terms into the contributions of the individual reactions.

Summary
We have implemented the MXL model as a new submodel in MESSy, in order to represent the diurnal dynamics of the 570 atmospheric boundary layer and their effect on atmospheric chemistry.The comprehensiveness of MXL/MESSy in representing the processes relevant to atmospheric chemistry in the convective boundary layer, while keeping computational requirements low, makes it an ideal tool for applica-575 tions in atmospheric chemistry that ask for systematic sensitivity analyses.These include the interpretation of observations from field campaigns and the evaluation of new process parametrizations under ambient conditions, as well as theoretical studies on the coupled land surface-boundary layer-580 atmospheric chemistry system.Expansion of the model with additional relevant MESSy submodels and conversion into a multi-layer column model is planned for the future: e.g.including the surface layer, entrainment zone, account for imperfect mixing in ABL (e.g.due to clouds, aerosol layers) 585 and the stable nocturnal boundary layer.Butler, T. M., Taraborrelli, D., Brühl, C., Fischer, H., Harder, H., Martinez, M., Williams, J., Lawrence, M. G., and        Jkg −1 latent heat of vaporization p − Clapp-Hornberger retention curve parameter q gkg −1 specific humidity q gkg −1 mixed-layer specific humidity q F T gkg −1 FT specific humidity q sat (T ) gkg −1 saturated specific humidity at temperature T q sat (T s ) gkg  the chemical production and loss with the chemical loss split up in a gas-phase destruction and a photolysis term. 215

235 4 .
closure assumption relating the surface buoyancy flux to the entrainment buoyancy flux w θ v e = −β w θ v s .
6 R. H. H. Janssen and A. Pozzer: Boundary layer dynamics with MXL/MESSy

Fig. 1 .Fig. 2 .Fig. 3 .Fig. 4 .
Fig.1.Spatio-temporal scales in atmospheric chemistry, typical species associated with these scales and the type of model that should be used to study their behaviour in the atmosphere.

Fig. 6 .
Fig. 6.Observed and modeled mixed-layer mixing ratios of a) O3, b) NO2, c) NO, d) isoprene, e) HNO3, f) OH, g) H2O2 and h) HO2 and i) NO2 photolysis rate for the DOMINO case, comparing MXL/MESSy and MIM2 chemistry with MXLCH and reduced chemistry.Results for both the SameIC and the BestFit runs with MXL/MESSy are shown.

Fig. 7 .
Fig. 7. O3 budget for the DOMINO case showing the contributions of (a) the processes entrainment, chemistry and dry deposition and (b) the chemical production and loss with the chemical loss split up in a gas-phase destruction and a photolysis term.

Table A2 .
Initial mixing ratio in ABL and FT, and surface emission fluxes of the reactants for the different runs in MXL/MESSy and MXLCH.Species in the reaction mechanism that are not included in this table have an initial concentration of zero and no surface emissions.For O2 and N2 we have imposed the values 2 × 10 8 and 8 × 10 8 ppb, respectively.

Table B1 .
List of acronyms used in this work

Table B2 :
List of variables and constants used in this work ms −1 ms −1 entrainment momentum flux in x-direction u w s ms −1 ms −1 surface momentum flux in x-direction u * gkg −1 ms −1 entrainment specific moisture flux w q s gkg −1 ms −1 surface specific moisture flux w θ v e