The CHIMERE chemistry-transport model v2023r1

A new version of the CHIMERE model is presented. This version contains both computational and physico-chemical changes. The computational changes make it easy to choose the variables to be extracted as a result, including values of maximum sub-hourly concentrations. Performance tests show that the model is 1.5 to 2 times faster than the previous version for the same setup. Processes such as turbulence, transport schemes and dry deposition have been modified and updated. Optimization was also performed for the management of emissions such as anthropogenic and mineral dust. The impact of fires on wind speed, soil properties and leaf area index (LAI) was added. Pollen emissions, transport and deposition were added for birch, ragweed, olive and grass. The model is validated with a simulation covering Europe with a 60 km × 60 km resolution and the entire year of 2019. Results are compared to various measurements, and statistical scores show that the model provides better results than the previous versions.

1 Introduction

Chemistry-transport modelling is useful for the analysis and forecast of atmospheric pollution events. Regional deterministic models have proven their value for studies about climate change, emissions regulations, day-to-day pollution evolution and health effects, among many other topics. These models have to both be accurate and have a low computational time.

The CHIMERE chemistry-transport model has been under continuous development since 1999. Several versions have been released as open-source software since the first one. The model is now used for analysis and forecast in many institutes, research laboratories and private companies. For the last 10 years, releases were distributed at a frequency of a new major one every 3 years. Each version of the model is distributed with other ancillary codes. These ancillary codes are developed by the same team and are intended to help users prepare the input data required for a simulation. These codes are detailed in the documentation. The model is distributed with other codes dedicated to helping users prepare for forcings such as weather and broadcasts. The model is also distributed with a test case to show the user the necessary input data and what the model's output looks like.

The first version with full documentation of all the processes taken into account in the model was version 2014. It included the possibility of reading emission fluxes from fires, calculated using a specific model in pre-processing; reading daily satellite data of fire radiative power and burnt areas; and converting this information into chemical species hourly fluxes using additional information, such as vegetation type and related fuel load for each. It also includes the Statewide Air Pollution Research Center (SAPRC) chemical mechanism (Carter, 2010), an adaptive time step and the new datasets for the chemical boundary conditions (Menut et al., 2013). Version 2017 (Mailler et al., 2017) added the possibility of using a hemispheric domain (and the way to manage the grid and the transports schemes), the new datasets to calculate mineral dust emissions in any global locations, the addition of the Fast-JX module to calculate photolysis rates and diagnose the lidar profiles, and the addition of a new resuspension scheme. Version 2020 provides coupling with the Weather Research and Forecasting (WRF) meteorological model through the Ocean Atmosphere Sea Ice Soil code using the Model Coupling Toolkit (OASIS3-MCT), the external coupler described in Craig et al. (2017). In addition, many chemical processes were added in the model: the volatility basis set aerosol schemes, the dimethyl sulfate (DMS) emissions and chemistry, the H²O aerosol scheme, the mineral dust mineralogy, the sub-grid-scale variability in anthropogenic emissions, and an updated strategy for vertical advection of pollutants (Menut et al., 2021a).

The novelties in CHIMERE v2023r1 relative to previous model versions concern the physical and chemical processes as well as the technical implementation of the model. In this article, we summarize all new developments and we show examples of their use. The main changes regarding the technical implementation are (i) new file format and databases for anthropogenic and fire emissions, (ii) addition of the XIOS2 (XML-IO-Server version 2) library, (iii) modularization of subroutines, (iv) possibility of outputting a partial aerosol optical depth (AOD) for each aerosol species, and (v) a new Fortran 90 implementation of the ISORROPIA (meaning equilibrium in Greek) model. Regarding the physico-chemical processes, the main changes are the (i) new calculation for the turbulence and dry deposition, (ii) new advection schemes, (iii) implementation of pollen emission, (iv) parameterization of the impact of fires on mineral dust emissions and leaf area index, and (v) externalization of the modules included in the SSH-aerosol package.

In this article, we present in detail all these novelties. The model structure is detailed, and the computing time is illustrated. Test cases are performed and comparisons to observations are presented in order to quantify the model ability to reproduce real pollution events. Finally, a general conclusion about this new version is presented in Sect. 9, including some perspectives for future development.

2 Model numerical structure

One of the most important changes for this model version is a complete re-organization of the code structure. It could appear as if it is not impacting the users in their daily way of using the model, but, in fact, it changes a lot the simulations. The simulations are faster, and the implementation of the XIOS (XML-IO-Server) library enables the user to better manage the outputs as well as their size. All these changes are described in this section.

2.1 Modularization

The code was rewritten to have programs organized into Fortran 90 modules. This does not change the order of the processes' calculations or the results of the simulation. It is just an internal choice to have more easily optimizable code for future developments.

Figure 1 Organization chart of CHIMERE v2023r1. The first yellow boxes present the initialization phase: static namelist and initial values are selected only on one processor (inichimere). Second, on all processors, the main fields are initialized for the first hour: boundary conditions, meteorology and concentrations of all chemical species (iniworker). Finally, the main temporal loop is performed (worker). It includes an hourly time step (in offline mode), then a physical time step, and then a chemical time step. Depending on the mode (offline or online), some exchanges are performed between the WRF model and the CHIMERE model.

Figure 1 presents the organization chart of this model version. The main program is called chimere. This program mainly manages three different blocks. First, the inichimere routine is called and is dedicated to the initialization of each module: reading input parameters (namelists or aerosol size distribution and list of chemical species, among others) and computing fixed parameters and array allocations. Depending on the processes requested by the user, modules are initialized or not; this way we prevent unnecessary file reading or memory burden. In the case of online coupling, an additional namelist for OASIS is read. Second, another initialization part has to be done in the iniworker routine. This is performed only for the first simulated hour in order to obtain values at zero time for all input variables: meteorology, turbulence, anthropogenic and fire emissions (depending on input databases), and boundary conditions. Third, the loop over time starts and is managed for each processor in the worker routine. This loop is performed at the physical time step, a sub-hourly division depending on the domain cell sizes. Depending on whether the simulation is offline (no feedback between chemistry and weather) or online (in blue), meteorological fields are read (hourly; in red) or received from WRF (physical time step; in blue). In the case of it being online, the concentrations are sent to WRF to take into account direct and indirect effects for the next physical time step. If offline, if the time step coincides with an entire hour, the hourly boundary conditions and anthropogenic and fire emissions are updated. For these hourly fields, a linear time interpolation is performed to retrieve information at the corresponding physical time step. Next, all processes related to the physical time step are performed: emissions, resuspension, aerosol physics, mixing, photolysis, deposition and chemical rates. Once we have all the necessary information, a loop is performed at the chemical time step. Two chemical solvers are available: twostep and splitting, as already explained in Menut et al. (2021a). For example, in twostep, the first step is performed, and then an additional loop is performed with Gauss–Seidel iterations for the second step. Finally, the newly implemented XIOS routine is called to write the results on the disc.

The model has thus three main steps: initialization of modules, initialization of first table values at zero time and time integration. The time integration contains three levels: (i) a sub-hourly physical time step, mainly depending on the grid size for meteorological variable accuracy; (ii) a chemical time step depending on the Courant–Friedrichs–Lewy (CFL) condition (Courant et al., 1952); and (iii) the Gauss–Seidel iterations (at least two are recommended).

2.2 XIOS I/O management library

Version 2 of XIOS (XIOS2) has been implemented in CHIMERE. XIOS (XML-IO-Server) is an open-source C++ library dedicated to I/O management that is usable in conjunction with numerical models written in Fortran/C/C++ (XIOS, 2023). It is a powerful tool for both output history file management and spatial and temporal post-processing of data at runtime. XIOS was developed with several goals in mind: (i) simplifying the I/O management inside the source code by reducing the number of calls needed to send a variable to be written; (ii) offering flexibility to the user without having to make modifications inside the source code; and (iii) making writing and computing happen simultaneously using asynchronous calls, meaning that the writing will not slow down the computation. XIOS manages to achieve these goals by setting up servers and clients; one or several CPUs can be dedicated to XIOS servers, which manage I/O processes exclusively. These servers receive asynchronous data from what are called clients (here the chimere ranks) and write them into one or multiple files; they can also manage any additional temporal or spatial post-processing requested by the user. XIOS has been implemented and evaluated in several models. Among many others, Institut Pierre-Simon Laplace (IPSL) Earth system models use XIOS exclusively (Boucher et al., 2020).

File and variable management is done by XML files. For standard use of the model, all necessary files, variables and parameters have been prepared in standard XML files, and the use of XIOS is transparent for the user. However, one of the perks of using XIOS is the flexibility it offers the users. This means that the user can create almost any combination of output files they want by adding or removing species or dimensions. Indeed, the XML files can be modified or new ones can be created at any time by the user without the need to recompile the model or, largely, even code in the model. Here, we list some common tasks that XIOS can perform:

• Hourly minima, maxima or averages of a particular species can be processed by XIOS and written to the disc without adding more computation time. For example, maxima can be useful for calculating air quality threshold exceedances.

• Species can be written only for a specific part of the domain (like a sounding or a subdomain) or a specific layer (the surface, for instance).

• Arithmetic operations can be performed to sum up variables in a chosen dimension.

• The user is completely free to define an output file only with the variables that are important for their specific usage.

Examples of such personalized cases have been provided in the code. It is important to keep in mind that customizing an XML file to create personalized XIOS output files is not limited to the examples provided by the development team, and each user can create output files that correspond to their specific needs.

Figure 2 Runtime comparison between several I/O options implemented in the model. Each boxplot contains thirty 1 d simulations, and the black line on each boxplot shows the average for each configuration. Five options are shown – no_output: model does not write any output files; XIOS_slim: an example of a personalized XIOS output; XIOS_full: the standard output files of the model written with XIOS; nc4: the standard output files of the model written with the standard I/O module in the model; and v2020r3: the previous distributed model version. The additional time (in %) is shown above each boxplot. For the no_output case, the average time in minutes is shown.

The performance of XIOS and the standard I/O system using netCDF 4 in parallel in the model (called nc4 hereafter) has been compared extensively. Figure 2 shows runtime comparisons of thirty 1 d simulations in different configurations. Two configurations of XIOS are shown here: the standard configurations that replicates the nc4 output files to the letter (called XIOS_full) and an example of a personalized configuration which only writes a list of species necessary for the usage in mind (called XIOS_slim). These are compared to the model runtime when it does not write any output, which can be considered the pure calculation time of the model (including reading of input data). The figure shows that using XIOS offers an important reduction in the runtime of the model, adding only 8.8 % and 15.6 % to the runtime for the XIOS_slim and the XIOS_full configurations, respectively. Compared to the added 69.4 % computation time in the case of the standard I/O module, XIOS offers a good performance increase in the model. Finer tuning of the performance is also possible by adapting the number of I/O servers or the size of buffers. Finally, this new model version is compared to the previous one, the distributed v2020r3 version: with a time increase of +105.4 %, the new model version is 2 times faster than the previous one.

2.3 Update of the user parameter namelist

The CHIMERE parameter file chimere.par was updated in order to remove some unused parameters and to add novelties. The complete list of parameters is presented in the model documentation (https://www.lmd.polytechnique.fr/chimere/docs/CHIMEREdoc_v2023.pdf, last access: 9 July 2024) for all possible flags. In this section, only the new flags corresponding to new functionalities are explained.

1. forecast. The forecast is re-activated in this model version. It is mainly the way to manage the WRF files when preparing the meteorology.

2. nproc_chimere, nproc_wrf and nproc_xios. These options are for the management of the simulation as a function of the available processors, for CHIMERE, WRF and XIOS, respectively. Obviously, the total of the three must be less than or equal to the total number of the available processors on the user's computer. For CHIMERE and WRF, the number of processors corresponds to the number of subdomains for the parallel computation. If you are offline, the number of subdomains is calculated for CHIMERE only. For example, in online mode, for a cluster core of 64 processors (classical configuration), it is recommended to select 15 for WRF, 45 for CHIMERE and 2 for XIOS.

3. ifire2dust, ifire2lai and ifire2drydep. These options concern the effect of the fires on other processes. This takes into account the possible changes in soil erodibility and surface wind speed and the leaf area index (LAI) for the biogenic emissions and for the dry deposition of gaseous species following Menut et al. (2023).

4. ilai. This option is used to select the LAI database between ilai = MOD (MODIS; Yuan et al., 2011) and ilai = COP (Copernicus; Fuster et al., 2020).

5. ip_ragweed. It checks whether the ragweed pollen calculation is activated (1) or not (0). Three schemes are available for this pollen type: (1) for the Sofiev et al. (2013) scheme, (2) for the Efstathiou et al. (2011) scheme, and (3) for the Menut et al. (2021b) scheme. Calculations are explained in Sect. 6.2.

6. ip_birch, ip_grass, ip_alder and ip_olive. These are dedicated to activating (1) or not activating (0) the pollen emissions for birch, grass, alder and olive, respectively.

7. ipolresus. This option activates the resuspension only for pollen.

8. pollendir. This option indicates the directory where all pollen databases are stored. For this model version, and with no other available dataset, pollen forcing data are available only over western Europe.

9. gtrc and ptrc. These activate the gaseous and aerosol tracers, respectively. It is necessary to complete two ASCII files (see Sect. 3.2): TRACERS_GAS.data and TRACERS_PART.data.

10. dgrb. This is the directory where the global meteorological files in GRIB format are stored (ECMWF or NCEP).

11. geodata. This is the directory of surface input data for the WRF Pre-Processing System (WPS). In this version, WPS (from WRF) is distributed since the domain file used by CHIMERE is the geog file of WRF. This avoids some interpolations of meteorological fields and then reduces errors in the meteorology and speeds up the model.

12. dmin and dmax. These options enable the definition of the limits of the aerosol size distribution. dmin and dmax are the size of the first interval of the first bin and the last interval of the last bin, respectively. These limits are the same for all aerosols. Values of 0.01 and 40 µm are recommended.

13. iforcut. If this flag is equal to 1, the calculation of the size distribution is forced to have retained the diameter of 2.5 and 10 µm as the cutoff.

14. timeverb. This flag manages some text printed on the screen, providing the time step in each module of the model. It is mainly useful to developers to see if new calculations have an important numerical cost.

15. AODPerSpecies. The aerosol optical depth is the outputted aerosol species per species.

The chemical module SSH-aerosol (Sartelet et al., 2020) is now compatible with CHIMERE to represent aerosol dynamics. SSH-aerosol is distributed independently of the CHIMERE model, and it will evolve in the next few years. It means that in order to be able to use this specific option, the user needs to download the SSH-aerosol module and compile it independently of all CHIMERE programs. For this, a specific script is available in the CHIMERE main folder to download SSH-aerosol. In the CHIMERE namelist, the use of SSH-aerosol is activated when values of soatyp are h2o or h2or. A specific namelist for SSH-aerosol is proposed in the namelist_ssh folder, and it can be adapted by the user according to the level of complexity suitable for computing aerosol dynamics. Note that, in this version, the SSH-aerosol package is currently only compatible with the operator splitting option and not with the twostep option. For consistency between the secondary organic aerosol (SOA) compounds in gas and particle phases, the MELCHIOR 2 chemical mechanism is recommended.

3 New calculations for the transport

3.1 Changes in the advection formulation

The advection schemes have been corrected and updated, and new schemes were added as a user's choice. The user can select different schemes for the horizontal and the vertical transport, the two being calculated separately.

In the vertical direction, the Després and Lagoutière (1999) scheme has been implemented. As described in Lachatre et al. (2020), two options are now available for the representation of vertical velocity: diagnose it either from the horizontal wind components, as it was historically done in CHIMERE, or using the vertical wind speed provided by the meteorological model. Lachatre et al. (2022) have shown that reducing numerical diffusion using the Després and Lagoutière (1999) advection scheme and using the vertical wind speed, w, from the meteorological model (is_diagwinw = .false. in chimere.par) also modifies the chemical processes in thin plumes; due to the nonlinearity of tropospheric chemistry, diluting a plume over a larger volume changes not only its distribution, but also the kinetics and relative importance of the chemical processes.

For horizontal advection, five different schemes are available: (1) the simple first-order advection scheme of Godunov and Bohachevsky (1959). It is the fastest scheme and may be used in the case of low resolution or forecast when results are expected soon. However, it is also diffusive and is not recommended in the case of thin plume studies. (2) The Van Leer scheme is fast and accurate (Van Leer, 1977). (3) The parabolic piecewise method (PPM) scheme (Colella and Woodward, 1984) has been revised and is now faster than the previous implementation. (4) The Walcek (2000) advection scheme has been implemented, and, finally, (5) the PPM+W advection scheme, a combination of the PPM and Walcek schemes, has been implemented (Mailler et al., 2023c). In an academic framework, Mailler et al. (2023c) have shown that the PPM+W performs best among these schemes in terms of accuracy, while the Van Leer and Walcek schemes are computationally cheaper. The PPM+W scheme has a computational cost similar to PPM.

3.2 Changes in tracer management

The introduction of “tracer” emissions inside model simulations including other (realistic) emissions or not is a possibility offered by the CHIMERE model. In CHIMERE, the word tracer does not necessarily imply that the species is inert chemically. The emitted species can be chemically either inert (as in Lachatre et al., 2020) or active (as in Lachatre et al., 2022) as soon as the emission occurs over a finite number of horizontal points (but possibly spread in the z direction). The user can specify the vertical repartition of tracer emissions and their evolution in time. This framework is well suited not only to include emissions due to volcanic eruptions or industrial accidents, for example (e.g. Lachatre et al., 2022), but also to study the behaviour and motion of air masses to and from key areas such as cities or mountain valleys (Lapere et al., 2023; Deroubaix et al., 2019). Lachatre et al. (2022) show an example of a real species (SO₂) emitted as a tracer over a specific point, along with its normal surfacic anthropogenic emissions from emiSURF.

In the new version of the CHIMERE model, gas and/or particulate tracers are set in two distinct input files (TRACERS_GAS.data and TRACERS_PART.data). The users decide to activate gas tracers or not by acting on two distinct parameters in the chimere.par file:

• gtrc. It is to be set to 1 (resp. 0) to activate (resp. deactivate) gas tracers.

• ptrc. It is to be set to 1 (resp. 0) to activate (resp. deactivate) particulate tracers.

Compared to earlier model versions, the new model version permits more flexibility and user control for particulate tracer emissions: for each tracer emission event, up to three different emission modes can be specified (from the median diameter and the multiplicative spread of the log-normal distribution), independently from the modes defined for anthropogenic emissions. In the TRACERS_PART.data, for each tracer emission, the user can define the following parameters:

• longitude and latitude of the emission point

• onset date and termination date of emission in the YYYYMMDDHH format (in v2023r1, the emission duration must be an entire number of hours and at least 1 h)

• emitted tracer mass

• minimum and maximum altitude of emission (the tracer mass will be distributed evenly between these two altitudes)

• f₁, f₂ and f₃ – fraction of tracer mass in the first, second and third size mode, respectively (the user must check that $f_{\mathrm{1}}+f_{\mathrm{2}}+f_{\mathrm{3}}=\mathrm{1}$)

• D₁, D₂ and D₃ – mass median diameter of the first, second and third size mode, respectively, for the emission.

• σ₁, σ₂ and σ₃ – geometric standard deviation of the first, second and third size mode, respectively, for the emission (σ_i>1)

The proportion of tracer mass between diameters D_m and D_M>D_m is

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle \frac{m(D_{\mathrm{m}}\to D_{\mathrm{M}})}{m}}& =\sum _{i=\mathrm{1}}^{\mathrm{3}}{f}_i\left[\mathrm{erf}\left({\displaystyle \frac{\ln D_{\mathrm{M}}-\ln D_i}{\sqrt{\mathrm{2}}\ln \mathit{\sigma }}}\right)\right.\\ & \left.-\mathrm{erf}\left({\displaystyle \frac{\ln D_{\mathrm{m}}-\ln D_i}{\sqrt{\mathrm{2}}\ln \mathit{\sigma }}}\right)\right].\end{array}\end{array}$$

These three emitted modes are then distributed into the nbins size bins for chemistry and transport, as described in Menut et al. (2013).

4 New land use

The pre-processing for the CHIMERE-specific nine classes of land use was removed in this version. We made the choice to directly use the land use generated by the WRF Pre-Processing System (WPS) and its output geog file. This permits it (i) to be faster to prepare in pre-processing; (ii) to be consistent with the meteorology when the two models, WRF and CHIMERE, are coupled; and (iii) to be consistent with the sub-grid-scale approach of the mineral dust emissions. Changing the land use classes imposed an update of the surface information required for turbulence, resuspension, albedo calculations and dry deposition. The surface information required for a simulation is presented in Table 1.

Lana et al. (2011)Sindelarova et al. (2014); Yuan et al. (2011)Prigent et al. (2012); Beegum et al. (2016); Journet et al. (2014)

Table 1 Files and variables required to initiate a CHIMERE simulation.

The use of a new tabulated dynamical roughness is done using the values proposed by the WRF model, which are already fitted for the United States Geological Survey (USGS) land use. The seasonality is considered using the WRF green vegetation fraction (GVF) as a proxy, with the following equation:

$$\begin{array}{}\text{(2)}& z_{\mathrm{0}\mathrm{m}}=(\mathrm{1}-\text{GVF})\times z_{\mathrm{0}\mathrm{m},\min }+\text{GVF}\times z_{\mathrm{0}\mathrm{m},\max }.\end{array}$$

In Table A1, the last z_0  m value, representative of water surface, is here set to zero but later calculated using the Charnock scheme (Charnock, 1955) (see Sect. 5.2 for the detailed calculation). Following the same process, the albedo values are tabulated and are presented in Table A2. It impacts the calculation of the resuspension: with the new land use, the resuspension is applied for all cell percentages related to urban and agricultural areas following USGS land use classification. It also impacts the dry deposition, and a new correspondence matrix was defined between the European Monitoring and Evaluation Programme (EMEP) (Emberson et al., 2001) and the USGS land use classes.

5 New calculation for the turbulence

All calculations dedicated to turbulence were grouped into a single Fortran module. Some parameterizations were updated at the same time. Compared to CHIMERE v2020r1, this new module includes the following changes:

• The roughness length over the sea, z₀, is not tabulated but calculated using the Charnock equation (Charnock, 1955).

• The parameterization of the surface sensible heat flux, Q₀, and the aerodynamic resistance, r_a, were upgraded to a more recent scheme. For r_a, its calculation is now on the CHIMERE vertical levels and not on the meteorological vertical grid.

• The calculation of turbulent parameters is now on a loop to ensure consistency between all variables (it includes Monin–Obukhov length calculation).

• It is now possible to diagnose the 10 m wind speed $|U{|}_{\mathrm{10}\mathrm{m}}$, the one from WRF that is dependent on the boundary layer height scheme and not always satisfying.

5.1 Variables to read or diagnose

In the CHIMERE model, many boundary-layer-related variables can be diagnosed. However, if the meteorological driver model provides some of them, they could be read directly and then used to diagnose the missing ones. The variables that can be read are selected in the chimere.par file and are the following:

• idiagblh. These variables use boundary layer height of the weather driver (idiagblh = 0) or diagnose it (idiagblh = 1) following Troen and Mahrt (1986).

• idiagflux. These variables use the surface turbulent characteristics as the friction velocity u_* and the surface sensible heat flux Q₀ of the weather driver (idiagflux = 0) or diagnose it (idiagflux = 1) following Troen and Mahrt (1986).

Depending on this choice, the list of variables that can be read or diagnosed and those that can only be diagnosed is presented in Table 2.

Table 2 List of variables required to describe the turbulence in and above the boundary layer. R stands for read and D for diagnostic.

For all calculations, several constants are defined in the code. They can be modified by the user, but it is not really recommended. These constant values are listed in Table 3.

Table 3 Constant values used for the calculation of additional parameters in the routine calc_turb.

5.2 Diagnostic of meteorological variables

Before starting the calculation of the turbulent parameters, it is necessary to estimate some other variables. A specific order is required, with some variables being estimated using other variables.

First, the surface pressure, P_surf, is simply diagnosed using the pressure values at the two first model levels:

$$\begin{array}{}\text{(3)}& P_{\mathrm{surf}}=P_{\mathrm{1}}+{\displaystyle \frac{z_{\mathrm{1}}(P_{\mathrm{1}}-P_{\mathrm{2}})}{z_{\mathrm{2}}-z_{\mathrm{1}}}}.\end{array}$$

The same calculation is applied to estimate the surface values of specific humidity, cloud liquid content, ice content and rain. The surface values of wind components are set to zero. The total amount of water in the atmosphere is

$$\begin{array}{}\text{(4)}& \text{wa}\left(z\right)=\text{cliq}\left(z\right)+\text{cice}\left(z\right)+\text{rain}\left(z\right),\end{array}$$

where cliq, cice and rain are the cloud, ice and rainwater mixing ratio (kg kg⁻¹). The potential temperature θ is calculated as

$$\begin{array}{}\text{(5)}& \mathit{\theta }\left(z\right)=T\left(z\right).{\left({\displaystyle \frac{P_{\mathrm{0}}}{P\left(z\right)}}\right)}^{R_a/C_p},\end{array}$$

with T being the air temperature (K), P₀=100 000 Pa the reference pressure, R_a=287.04 J K⁻¹ kg⁻¹ the specific gas constant for dry air and C_p=1005 J K⁻¹ kg⁻¹ the specific heat capacity of air. Note that the ratio $R_{\mathrm{a}}/C_p$ is also noted and used as γ=0.2857. The virtual potential temperature, θ_v, is the equivalent of θ but accounting for dry and humid air as

$$\begin{array}{}\text{(6)}& {\mathit{\theta }}_v\left(z\right)=(\mathrm{1}+\mathrm{0.61}q(z\left)\right)T\left(z\right){\left({\displaystyle \frac{P_{\mathrm{0}}}{P\left(z\right)}}\right)}^{R/C_p}.\end{array}$$

The wind speed module, $\left|U\right|$ (m s⁻¹), is expressed as

$$\begin{array}{}\text{(7)}& \left|U\right|\left(z\right)=\sqrt{u^{\mathrm{2}}+v^{\mathrm{2}}}.\end{array}$$

Air density is calculated (in molecules m⁻³) as

$$\begin{array}{}\text{(8)}& \mathit{\rho }={\displaystyle \frac{P}{R_{\mathrm{a}}\cdot T}}{\displaystyle \frac{N_A}{{\mathrm{10}}^{-\mathrm{3}}{M}_{\mathrm{air}}}},\end{array}$$

with $N_A=\mathrm{6.022}\times {\mathrm{10}}^{\mathrm{23}}$ mol⁻¹ and M_air=29 g mol⁻¹ the molar mass of air.

For the calculation of relative humidity (RH), the partial pressure of water vapour, e_* (Pa), is first estimated following Tetens' formula (Tetens, 1930) as

$$\begin{array}{}\text{(9)}& e_{*}=\mathrm{611}\cdot \exp \left[\mathrm{17.27}\left({\displaystyle \frac{T-T_{\mathrm{0}}}{T-\mathrm{35.86}}}\right)\right],\end{array}$$

with T₀=273.15 K. The dew point (in g g⁻¹) is deduced as

$$\begin{array}{}\text{(10)}& \text{de}=\mathrm{0.622}{\displaystyle \frac{e_{*}}{P-e_{*}}}.\end{array}$$

Finally, the relative humidity is calculated as

$$\begin{array}{}\text{(11)}& \text{RH}={\displaystyle \frac{q}{\text{de}}},\end{array}$$

with q being the specific humidity. The relative humidity is protected against extreme values and bounded to the rhmax constant value defined in Table 3.

One key point for turbulence parameterizations is to diagnose the near-surface virtual potential temperature, here noted as θ_v,0. To access this value, we extrapolate the vertical temperature profile using the 2 m temperature, T_2 m, and the temperature at the first model vertical level, T₁, to obtain T₀. The following expression is used:

$$\begin{array}{}\text{(12)}& T_{\mathrm{0}}=T_{\mathrm{2}\mathrm{m}}-{\displaystyle \frac{z_{\mathrm{2}\mathrm{m}}(T_{\mathrm{1}}-T_{\mathrm{2}\mathrm{m}})}{z_{\mathrm{1}}-z_{\mathrm{2}\mathrm{m}}}}.\end{array}$$

Over sea, it is common to use the Charnock equation to estimate the momentum roughness length, which is dependent on the sea surface and the wave height:

$$\begin{array}{}\text{(13)}& z_{\mathrm{0}}={\mathit{\alpha }}_{\mathrm{c}}{\displaystyle \frac{u_{*}^{\mathrm{2}}}{g}},\end{array}$$

where α_c is a constant (α_c≈0.015). Pena and Gryning (2008) note that this constant could range between 0.008 and 0.06 depending on the study. An iterative approach is necessary to quantify the roughness length depending on the friction velocity and vice versa of the two. Fortunately, Hersbach (2011) proposes an efficient fit to estimate z₀ without knowledge of u_*:

$$\begin{array}{}\text{(14)}& z_{\mathrm{0}}={\displaystyle \frac{z}{\exp \left(b_n\right)-\mathrm{1}}},\end{array}$$

with $b_n={\left[(b_n^{\mathit{\nu }}{)}^p+(b_n^{\mathit{\alpha }}{)}^p\right]}^{\mathrm{1}/p}$ and $p=-\mathrm{12}$. The two parameters expressed that z₀ depends on kinematic viscosity, ν, for light wind and on the Charnock relation for strong wind. They are estimated with

$$\begin{array}{}\text{(15)}& \begin{array}{rl}{b}_n^{\mathit{\nu }}& =-\mathrm{1.47}+\mathrm{0.93}\ln \left(R\right),\\ R& ={\displaystyle \frac{z}{{\mathit{\alpha }}_M\mathit{\nu }}}\left(k{u}_n\right),\\ b_n^{\mathit{\alpha }}& =+\mathrm{2.65}-\mathrm{1.44}\ln \left(A\right)-\mathrm{0.015}\left(\log \right(A){)}^{\mathrm{2}},\\ A& ={\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{ch}}}{gz}}(k{u}_n{)}^{\mathrm{2}},\end{array}\end{array}$$

with k=0.41 being the von Kármán constant, $\mathit{\nu }=\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹ the kinematic viscosity, u_n the wind speed at height z (here z=10 m), α_M=0.11 and α_ch=0.018 the Charnock constant. However, this constant varies following studies and ways to determine it as, for example, explained in Feng et al. (2021). This roughness length calculated for three different Charnock constants is displayed in Fig. 3.

Figure 3 Roughness length over water surfaces calculated for three different Charnock constants.

5.3 Optical depth above model top

In order to have the optical thickness due to clouds for the full tropospheric column, an optical depth is estimated above the CHIMERE domain top using the meteorological model data. It is generally possible since the CHIMERE model currently has a model top around 200 hPa (user's choice), while the meteorological model has a model top around 50 hPa. For the atmospheric column in CHIMERE, the Fast-JX module calculates the cloud optical depth (COD) directly. The top level of CHIMERE, ptop_chimere, is first estimated in the meteorological model grid. The following equation is then used:

$$\begin{array}{}\text{(16)}& {\text{COD}}_{\mathrm{above}}=\sum _{{\mathrm{ptop}}_{\mathrm{chimere}}}^{{\mathrm{ptop}}_{\mathrm{meteo}}}{\displaystyle \frac{\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}.(\stackrel{\mathrm{‾}}{\text{RH}}-\mathrm{0.95}).\mathrm{\Delta }z}{\mathrm{1}-\mathrm{0.95}}},\end{array}$$

where 0.95 is the value of high cloud attenuation and $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}$ is the high cloud optical depth for RH = 1 (in m⁻¹). ptop_meteo is the pressure of the top level of the meteorological model. Finally, the attenuation due to cloud above the model top is estimated using the following empirical formula:

$$\begin{array}{}\text{(17)}& A_{\mathrm{top}}=\exp (-\mathrm{0.11}\cdot {\text{COD}}_{\mathrm{above}}^{\mathrm{2}/\mathrm{3}}).\end{array}$$

5.4 Urban correction

A simple urban correction is available in the model. This is very simple and only implemented to retrieve some observed trends: lower surface wind speed and higher surface sensible heat flux. This option may be activated in the chimere.par parameter file by changing the urbancorr flag. This activation induces a slight change for the wind speed in the surface layer (from ground to surface layer height) and Q₀, the surface sensible heat flux, only over urbanized areas as

$$\begin{array}{}\text{(18)}& {\displaystyle }{U}_{\mathrm{10}\mathrm{m},\mathrm{urb}}=\text{factU}\times U_{\mathrm{10}\mathrm{m}},\text{(19)}& {\displaystyle }{Q}_{\mathrm{0},\mathrm{urb}}=\text{addflux}+Q_{\mathrm{0}}.\end{array}$$

The factor for the wind speed is equal to the proportion of urban surface in the cell multiplied by a factor defined in the module calc_turb.F90 as a parameter, factU. A standard default value is 0.8. The change for Q₀ is an additional flux called addQ0 and equal to 30 W m⁻². These values may be changed directly in the model code. Note that this correction is applied to the wind profile, (u,v), before the diagnostic of turbulent fluxes and after the diagnostic of Q₀. The values correspond to a concatenation of wind and temperature observations found in numerous publications on urban meteorology, such as Oke (1982), Calbo et al. (1998) and Arnfield (2003).

5.5 Diagnostic of turbulent variables

It is possible to directly use the turbulent variables calculated by the meteorological forcing model. For that, two conditions have to be fulfilled: (i) the flag idiagflux = 1 in chimere.par and (ii) the corresponding variables, friction velocity, u_*, and surface sensible heat flux, Q₀, must be available in the meteorological model outputs. If these two conditions are not reached together, these variables have to be diagnosed, with each variable being dependent on the others. Due to these dependencies, it is an iterative process.

5.5.1 The Richardson number

First, a calculation of the surface layer turbulent variables must be done. It is necessary to calculate the stability of the atmosphere. Following the scheme of Jiménez et al. (2012), the bulk Richardson number is calculated as follows:

$$\begin{array}{}\text{(20)}& {\mathit{Ri}}_{\mathrm{b}}={\displaystyle \frac{gz}{{\mathit{\theta }}_{\mathrm{a}}}}{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{va}}-{\mathit{\theta }}_{\mathrm{vg}}}{U^{\mathrm{2}}}},\end{array}$$

with θ_vg and θ_va (K) being the virtual potential temperature at ground and first model levels, respectively. The first vertical model level is at height z. g is the acceleration of gravity. U is the mean wind speed at altitude z. To avoid values of Ri_b that are too large and unrealistic, a minimum value of U=0.1 m s⁻¹ is forced.

The Monin–Obukhov length, L, is iteratively estimated using the first-guess value of Ri_b and the following equation:

$$\begin{array}{}\text{(21)}& {\mathit{Ri}}_{\mathrm{b}}={\displaystyle \frac{z}{L}}{\displaystyle \frac{{\mathrm{\Psi }}_T}{{\displaystyle }{\left[\ln \left({\displaystyle \frac{z+z_{\mathrm{0}\mathrm{m}}}{z_{\mathrm{0}}}}\right)-{\mathit{\psi }}_{\mathrm{m}}\left({\displaystyle \frac{z+z_{\mathrm{0}\mathrm{m}}}{L}}\right)+{\mathit{\psi }}_{\mathrm{m}}\left({\displaystyle \frac{z_{\mathrm{0}\mathrm{m}}}{L}}\right)\right]}^{\mathrm{2}}}},\end{array}$$

with Ψ_T being equal to

$$\begin{array}{}\text{(22)}& {\mathrm{\Psi }}_T=\ln \left({\displaystyle \frac{z+z_{\mathrm{0}\mathrm{h}}}{z_{\mathrm{0}\mathrm{h}}}}\right)-{\mathit{\psi }}_{\mathrm{h}}\left({\displaystyle \frac{z+z_{\mathrm{0}\mathrm{h}}}{L}}\right)+{\mathit{\psi }}_{\mathrm{h}}\left({\displaystyle \frac{z_{\mathrm{0}\mathrm{h}}}{L}}\right),\end{array}$$

where z_0  m is the momentum roughness length and is tabulated. For the heat roughness length, several formulations and uses exist. For example, Beljaars and Holtslag (1991) use $z_{\mathrm{0}\mathrm{h}}=z_{\mathrm{0}\mathrm{m}}/\mathrm{10}$. In Jiménez et al. (2012), all equations use only z₀ without distinction between momentum and heat. We thus apply the relation z_0h=z_0m here because the use of a lower z_0h gives too low of a surface sensible heat flux. The stability functions ${\mathit{\psi }}_{\mathrm{m}}(z/L)$ (for momentum) and ${\mathit{\psi }}_{\mathrm{h}}(z/L)$ (for heat) are calculated using the modified formulation of Jiménez et al. (2012).

1. If Ri_b>0 (nocturnal stable boundary layer), then

   $$\begin{array}{}\text{(23)}& \left\{\begin{array}{l}{\mathit{\psi }}_{\mathrm{m}}=-a\ln \left[\frac{z}{L}+{\left[\mathrm{1}+{\left(\frac{z}{L}\right)}^b\right]}^{\mathrm{1}/b}\right],\\ {\mathit{\psi }}_{\mathrm{h}}=-c\ln \left[\frac{z}{L}+{\left[\mathrm{1}+{\left(\frac{z}{L}\right)}^d\right]}^{\mathrm{1}/d}\right],\end{array}\right.\end{array}$$

   with a=6.1, b=2.5, c=5.3 and d=1.1.

2. If Ri_b≤0 (free unstable conditions), then

   $$\begin{array}{}\text{(24)}& {\mathit{\psi }}_{\mathrm{h},\mathrm{m}}={\displaystyle \frac{{\displaystyle }{\mathit{\psi }}_{\mathrm{K},\mathrm{h},\mathrm{m}}\left({\displaystyle \frac{z}{L}}\right)+{\left({\displaystyle \frac{z}{L}}\right)}^{\mathrm{2}}{\mathit{\psi }}_{\mathrm{C},\mathrm{h},\mathrm{m}}\left({\displaystyle \frac{z}{L}}\right)}{{\displaystyle }\mathrm{1}+{\left({\displaystyle \frac{z}{L}}\right)}^{\mathrm{2}}}},\end{array}$$

   where ${\mathit{\psi }}_{\mathrm{C},\mathrm{h},\mathrm{m}}$ represents the convective contribution as

   $$\begin{array}{}\text{(25)}& {\mathit{\psi }}_{\mathrm{C},\mathrm{h},\mathrm{m}}={\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}\ln \left({\displaystyle \frac{y^{\mathrm{2}}+y+\mathrm{1}}{\mathrm{3}}}\right)-\sqrt{\mathrm{3}}\arctan \left({\displaystyle \frac{\mathrm{2}y+\mathrm{1}}{\sqrt{\mathrm{3}}}}\right)+{\displaystyle \frac{\mathit{\pi }}{\mathrm{3}}},\end{array}$$

   with $y={\left[\mathrm{1}-{\mathit{\alpha }}_{\mathrm{h},\mathrm{m}}\left(z/L\right)\right]}^{\mathrm{1}/\mathrm{3}}$. α_m=10 and α_h=34, after Grachev et al. (2000). ${\mathit{\psi }}_{\mathrm{K},\mathrm{h},\mathrm{m}}$ represents the Kansas-type contribution as initially proposed by Paulson (1970).

   $$\begin{array}{}\text{(26)}& \left\{\begin{array}{l}{\mathit{\psi }}_{\mathrm{m}}={\displaystyle }\mathrm{2}\ln \left({\displaystyle \frac{\mathrm{1}+x}{\mathrm{2}}}\right)+\ln \left({\displaystyle \frac{\mathrm{1}+x^{\mathrm{2}}}{\mathrm{2}}}\right)-\mathrm{2}\arctan \left(x\right)+{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}},\\ {\mathit{\psi }}_{\mathrm{h}}={\displaystyle }\mathrm{2}\ln \left[{\displaystyle \frac{\mathrm{1}+\sqrt{\mathrm{1}-\mathit{\gamma }{\displaystyle }\left({\displaystyle \frac{z}{L}}\right)}}{\mathrm{2}}}\right],\end{array}\right.\end{array}$$

   with $x={\left[\mathrm{1}-\mathit{\gamma }(z/L)\right]}^{\mathrm{1}/\mathrm{4}}$ and γ=16.

5.5.2 The friction velocity (u_*)

The friction velocity, u_*, and the potential temperature scale, θ_*, are calculated as follows (Beljaars and Holtslag, 1991):

$$\begin{array}{}\text{(27)}& {\displaystyle }{u}_{*}={\displaystyle \frac{kU}{\ln \left({\displaystyle \frac{z+z_{\mathrm{0}\mathrm{m}}}{z_{\mathrm{0}\mathrm{m}}}}\right)-{\mathit{\psi }}_{\mathrm{m}}\left({\displaystyle \frac{z+z_{\mathrm{0}\mathrm{m}}}{L}}\right)+{\mathit{\psi }}_{\mathrm{m}}\left({\displaystyle \frac{z_{\mathrm{0}\mathrm{m}}}{L}}\right)}},\text{(28)}& {\displaystyle }{\mathit{\theta }}_{*}={\displaystyle \frac{k\left({\mathit{\theta }}_{\mathrm{a}}-{\mathit{\theta }}_g\right)}{\ln \left({\displaystyle \frac{z+z_{\mathrm{0}\mathrm{h}}}{z_{\mathrm{0}\mathrm{h}}}}\right)-{\mathit{\psi }}_{\mathrm{h}}\left({\displaystyle \frac{z+z_{\mathrm{0}\mathrm{h}}}{L}}\right)+{\mathit{\psi }}_{\mathrm{h}}\left({\displaystyle \frac{z_{\mathrm{0}\mathrm{h}}}{L}}\right)}}.\end{array}$$

Note that we also added the correction term here, with $z+z_{\mathrm{0}\mathrm{m},\mathrm{h}}$ in place of z. For the first iteration, we can calculate u_* and θ_* under neutral conditions, considering ψ_m and ψ_h are equal to unity. The Monin–Obukhov length, L, is then updated as

$$\begin{array}{}\text{(29)}& L={\displaystyle \frac{{\mathit{\theta }}_{\mathrm{a}}{u}_{*}^{\mathrm{2}}}{kg{\mathit{\theta }}_{*}}}.\end{array}$$

The iterative algorithm may be summarized as follows:

• 1.

  Ri_b is calculated using the meteorological variables in Eq. (20).

• 2.

  L is estimated using Eq. (21) and considering neutral conditions to have ${\mathit{\psi }}_{\mathrm{m}}={\mathit{\psi }}_{\mathrm{h}}=\mathrm{0}$.

• 3.

  A loop is then performed with the

  • (a)

    calculation of ψ_m and ψ_h with the updated $(z/L)$ value,

  • (b)

    calculation of u_* and θ_* using the updated ψ_m and ψ_h, and

  • (c)

    calculation of L using the updated u_* and θ_*.

The iterations stop when two successive values of L are lower than a pre-defined threshold. Here, we use ΔL=0.1. After several tests, the convergence is generally reached after four iterations.

The 10 m wind speed is then diagnosed by

$$\begin{array}{}\text{(30)}& u_{\mathrm{10}\mathrm{m}}=u_{\mathrm{a}}{\displaystyle \frac{\ln \left({\displaystyle \frac{\mathrm{10}+z_{\mathrm{0}\mathrm{m}}}{z_{\mathrm{0}\mathrm{m}}}}\right)-{\mathit{\psi }}_{\mathrm{m}}\left({\displaystyle \frac{\mathrm{10}+z_{\mathrm{0}\mathrm{m}}}{L}}\right)+{\mathit{\psi }}_{\mathrm{m}}\left({\displaystyle \frac{z_{\mathrm{0}\mathrm{m}}}{L}}\right)}{\log \left({\displaystyle \frac{z+z_{\mathrm{0}\mathrm{m}}}{z_{\mathrm{0}\mathrm{m}}}}\right)-{\mathit{\psi }}_{\mathrm{m}}\left({\displaystyle \frac{z+z_{\mathrm{0}\mathrm{m}}}{L}}\right)+{\mathit{\psi }}_{\mathrm{m}}\left({\displaystyle \frac{z_{\mathrm{0}\mathrm{m}}}{L}}\right)}}.\end{array}$$

5.5.3 The surface heat fluxes

The surface heat fluxes are not always provided in forcing models. But we need this quantity to evaluate w_* and $\stackrel{\mathrm{‾}}{h}$. For its quantification, there are two ways. First, we have the information about the sensible H (W m⁻²) and latent LE (W m⁻²) turbulent fluxes (the “eddy fluxes”). We can thus estimate the fluxes close to the surface as

$$\begin{array}{}\text{(31)}& \left\{\begin{array}{l}{\stackrel{\mathrm{‾}}{w^{\prime }{T}^{\prime }}}_{\mathrm{0}}=\stackrel{\mathrm{‾}}{(H+H_{\mathrm{urb}})}.{\displaystyle \frac{R}{C_p}}.{\displaystyle \frac{T_{\mathrm{2}\mathrm{m}}}{P_{\mathrm{surf}}}},\\ {\stackrel{\mathrm{‾}}{w^{\prime }{q}^{\prime }}}_{\mathrm{0}}=\stackrel{\mathrm{‾}}{\mathit{LE}}.{\displaystyle \frac{R}{L_{\mathrm{v}}}}.{\displaystyle \frac{T_{\mathrm{2}\mathrm{m}}}{P_{\mathrm{surf}}}},\end{array}\right.\end{array}$$

with R being 287.04 J K⁻¹ kg⁻¹, C_p the specific heat capacity (C_p=1005 J K⁻¹ kg⁻¹) and L_v the latent heat of vaporization ($L_{\mathrm{v}}=\mathrm{2.45}\times {\mathrm{10}}^{\mathrm{6}}$ J kg⁻¹). H_urb is an additional surface sensible heat flux that represents an urban contribution and is chosen to be H_urb=0 by default. A value of H_urb=30 W m⁻² is often used in regional or urban modelling. The total (heat + latent) surface flux, Q₀, is then estimated as

$$\begin{array}{}\text{(32)}& Q_{\mathrm{0}}={\stackrel{\mathrm{‾}}{w^{\prime }{T}^{\prime }}}_{\mathrm{0}}\times (\mathrm{1}+\mathrm{0.61}q)+{\stackrel{\mathrm{‾}}{w^{\prime }{q}^{\prime }}}_{\mathrm{0}}\times \mathrm{0.61}\times {\mathit{\theta }}_{\mathrm{0}}.\end{array}$$

Second, if we have no information, we use the previous parameterized values of u_* and θ_*. Q₀ is estimated to be

$$\begin{array}{}\text{(33)}& Q_{\mathrm{0}}=-u_{*}{\mathit{\theta }}_{*}.\end{array}$$

5.5.4 The boundary layer height $\stackrel{\mathrm{‾}}{h}$

If the boundary layer height ($\stackrel{\mathrm{‾}}{h}$) is available in the meteorological model and if it is the user's choice, the variable is just read and not calculated. If the user prefers to have a diagnostic made in CHIMERE, the calculation is done as follows. To select the use of the meteorological model, the following two conditions must be reached: (i) the boundary layer height, $\stackrel{\mathrm{‾}}{h}$, must be available in the meteorological model outputs, and (ii) the flag idiagblh must be 0 in chimere.par.

In the case of diagnostics, the following calculation is performed. Formulations are different depending on the atmospheric static stability. When stable, i.e. when L>0, $\stackrel{\mathrm{‾}}{h}$ is estimated to be the altitude when the Richardson number reaches a critical number, here chosen to be Ri_c=0.5, following Troen and Mahrt (1986). Under unstable situations (i.e. convective), $\stackrel{\mathrm{‾}}{h}$ is estimated using a convectively based boundary layer height calculation. This is based on a simplified and diagnostic version of the approach of Cheinet and Teixeira (2003), which consists of the resolution of the (dry) thermal plume equation with diffusion. The in-plume vertical velocity and buoyancy equations are solved, and the boundary layer top is taken to be the height where calculated vertical velocity vanishes. Thermals are initiated with a non-zero vertical velocity and potential temperature departure, depending on the turbulence similarity parameters in the surface layer. The two formulations are used, and the maximum of the two diagnosed $\stackrel{\mathrm{‾}}{h}$ values is retained.

5.5.5 The free convection velocity scale (w_*) and the aerodynamical resistance (r_a)

Now that we have the boundary layer height, $\stackrel{\mathrm{‾}}{h}$, it is possible to more properly estimate the free convection velocity scale, w_*:

$$\begin{array}{}\text{(34)}& w_{*}={\left({\displaystyle \frac{g\cdot Q_{\mathrm{0}}.\stackrel{\mathrm{‾}}{h}}{{\mathit{\theta }}_{\mathrm{v},\mathrm{a}}}}\right)}^{\mathrm{1}/\mathrm{3}},\end{array}$$

where $\stackrel{\mathrm{‾}}{h}$ is the convective boundary layer height and $\stackrel{\mathrm{‾}}{{\mathit{\theta }}_{\mathrm{v},\mathrm{a}}}$ the mean virtual potential temperature at the first model vertical level. For the aerodynamical resistance, we use the Bruin and Holtslag (1982) equation:

$$\begin{array}{}\text{(35)}& r_{\mathrm{a}}={\displaystyle \frac{\mathrm{4.72}{\left[\ln \left(\frac{z}{z_{\mathrm{0}\mathrm{m}}}\right)\right]}^{\mathrm{2}}}{\mathrm{1}+\mathrm{0.54}U}},\end{array}$$

with U being the wind speed at height z (in our case, the first vertical model level). Note that Bruin and Holtslag (1982) use a fit with the momentum roughness length z_0 m, while the usual expression of r_a uses the heat roughness length. This fit replaces the regular expression of the aerodynamical resistance (Beljaars and Viterbo, 1994) by

$$\begin{array}{}\text{(36)}& {\displaystyle \frac{\mathrm{1}}{r_{\mathrm{a}}}}={\displaystyle \frac{k{u}_{*}}{\ln \left(\frac{z+z_{\mathrm{0}\mathrm{h}}}{z_{\mathrm{0}\mathrm{h}}}\right)-{\mathrm{\Psi }}_{\mathrm{h}}\left(\frac{z+z_{\mathrm{0}\mathrm{h}}}{L}\right)}}.\end{array}$$

5.5.6 Momentum vertical turbulent diffusivity (K_z)

The momentum turbulent diffusivity, K_z (m² s⁻¹), has different formulations depending on the period of the day (stable or unstable conditions) and the altitude (surface layer, mixed layer or free troposphere).

In the surface layer, when conditions are stable, the turbulence is the most difficult to diagnose, with the variances of meteorological variables being very small. Troen and Mahrt (1986) remind us that when the surface layer is considered to be in equilibrium, then the parameterization of Louis (1979) is used. In the mixed layer, K_z is calculated following the parameterization of Troen and Mahrt (1986) without the counter-gradient term by

$$\begin{array}{}\text{(37)}& K_z=k{w}_{\mathrm{s}}\stackrel{\mathrm{‾}}{h}{\displaystyle \frac{z}{\stackrel{\mathrm{‾}}{h}}}{\left(\mathrm{1}-{\displaystyle \frac{z}{\stackrel{\mathrm{‾}}{h}}}\right)}^{\mathrm{2}},\end{array}$$

where k=0.41 is the von Kármán constant, z the altitude, $\stackrel{\mathrm{‾}}{h}$ the boundary layer height and w_s a vertical velocity scale. Note that the formulation of Eq. (1) in Troen and Mahrt (1986) is false. The calculation of w_s depends on the atmospheric stability. For the stable case, the following applies:

$$\begin{array}{}\text{(38)}& w_{\mathrm{s}}={\displaystyle \frac{u_{*}}{(\mathrm{1}+\mathrm{4.7}z/L)}}.\end{array}$$

On the other hand, for the unstable case, the following applies:

$$\begin{array}{}\text{(39)}& w_{\mathrm{s}}=(u_{*}^{\mathrm{3}}+\mathrm{7}\mathit{ϵ}k{w}_{*}^{\mathrm{3}}{)}^{\mathrm{1}/\mathrm{3}}.\end{array}$$

The constant ϵ is determined to be 0.1 in Troen and Mahrt (1986). For the model, this value is defined as $e=min(\mathrm{0.1},z/\stackrel{\mathrm{‾}}{h})$ to ensure a robust calculation in the first few vertical model levels. To avoid numerical problems, values of K_z are bound using pre-defined constant values, as explained in Table 3.

Above the boundary layer and in the presence of clouds, an additional diffusivity is added. This diffusivity is estimated if the relative humidity is larger than 0.9. For each model vertical level, a new vertical diffusivity, K_z, is estimated. The formulation is inspired by the Betts et al. (1995) study and the implementation of the turbulence scheme in the NCEP operational Medium-Range Forecast (MRF) model. The concept is based on the Louis (1979) formulation, defining stability function.

The vertical diffusion coefficient, $K_{z,m,h}$ (as a function of altitude z and for momentum m and heat h), is expressed here as

$$\begin{array}{}\text{(40)}& K_{z,m,h}=l_{\mathrm{mix}}^{\mathrm{2}}{f}_{m,h}\left({\mathit{Ri}}_{\mathrm{b}}\right)\left|{\displaystyle \frac{\partial U}{\partial z}}\right|,\end{array}$$

where Ri is the local gradient Richardson number, U the wind speed module, z the altitude, f_m,h the stability function for momentum and heat, and l the mixing length diagnosed as

$$\begin{array}{}\text{(41)}& {\displaystyle \frac{\mathrm{1}}{l}}={\displaystyle \frac{\mathrm{1}}{kz}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{c}}}},\end{array}$$

with k being the von Kármán constant and λ_c the asymptotic length scale, with λ_c=150 m (but 250 m in Betts et al., 1995). The stability functions, f_m and f_m,h, are estimated to be

$$\begin{array}{}\text{(42)}& f_{\mathrm{m}}=\left\{\begin{array}{ll}\mathrm{1}+\mathrm{8}\left({\displaystyle \frac{-Ri}{\mathrm{1}+C_{m,h}\sqrt{\mathit{Ri}}}}\right)& \text{for }\mathit{Ri}<\mathrm{0}\text{ (unstable)},\\ {\displaystyle \frac{\mathrm{1}}{(\mathrm{1}+\mathrm{5}\mathit{Ri}{)}^{\mathrm{2}}}}& \text{for }\mathit{Ri}\ge \mathrm{0}\text{ (stable)},\end{array}\right.\end{array}$$

with C_m=1.746 and C_h=1.286.

6 Changes in emissions

Emissions are the main chemical forcings of the model. In parallel with the development of the CHIMERE model itself, the various emissions taken into account are also updated between two versions. This section therefore presents emissions newly taken into account as well as updated and improved ones.

6.1 Changes in anthropogenic emissions

The management of the anthropogenic emissions has changed in this version of the model. The goal is to have fewer files to manage, which is a sore point of supercomputers. The new format corresponds to one file per day of the week, so seven files per month. Indeed, we have temporal data to differentiate the months of the year, the days of the week in a month and the 24 h of a day.

All chemical emitted species are grouped together in each file (before it was one file per species). For each species, the resolution is hourly and the fluxes are projected on the horizontal grid of the model. Using the provided emiSURF program, the user can also choose to differentiate their emission fluxes by activity sector according to the GNFR nomenclature. This differentiation will then be necessary in the case of simulations using the sub-grid variability option.

Finally, the user is encouraged to use the emiSURF model (described in Mailler et al., 2017) developed by the CHIMERE model team. This model is able to ingest multiple input databases, as described in Table 4. In its latest version (v2023r1), emiSURF takes into account the most recent time factor data (Guevara et al., 2021; Crippa et al., 2020).

Table 4 The different anthropogenic emission inventories processed by the emiSURF model. NMVOC: non-methane organic volatile compound, BC: black carbon, OC: organic carbon.

6.2 Implementation of pollen emissions

Pollen emissions were taken into account in a parallel version of CHIMERE operated by INERIS for the Copernicus daily forecast over Europe. The code was cleaned up and is now implemented in the official CHIMERE version. It means that from now on, pollen emission, transport and deposition will always be available in CHIMERE. This development comes in addition to the simulation of gases and aerosols: pollens are inert and there are no interactions with other pollutants. Model scores are therefore not modified if pollens are calculated. Additional details about this calculation in CHIMERE are presented in Sofiev et al. (2015) (for birch) and Menut et al. (2021b) (for ragweed).

6.2.1 The pollen grain specifics

Pollen grains are very different from the usual atmospheric aerosols. Regarding their transport, Sofiev et al. (2006) showed that they can be transported together with air masses and may be mixed due to small turbulent eddies. In this model version, five different pollen species were implemented: birch, olive, ragweed, grass and alder. The pollen density depends on each species and varies between 800 and 1200 kg m⁻³, which yields the sedimentation velocity of 1.2–1.3 cm s⁻¹. The main dry deposition process for pollens is gravitational settling, and they can be washed out by rains. For pollen, only the “big” mode of the aerosol is designed. The pollen species are considered non-chemically reactive. Optical properties are prescribed for short-wave and long-wave wavelengths.

6.2.2 Main principle of emissions for birch, grass and olive

To calculate pollen emissions, several surface databases are needed. In this model version, the databases are limited to western Europe. Three parameters are needed for each species: the surface fraction, the heat sum threshold and a scale factor. For the grass pollen, the data come from the SILAM model (Sofiev et al., 2006). Some additional information is provided: the mean diameter, D_p=32 µm; the density, 800 kg m⁻³; and the height above surface, h=1 m. Note also that there is no time dependency in this database. It is probably to improve later since the ragweed locations increase every year in Europe (Schaffner et al., 2020).

The pollen emission modelling is challenging because it requires fair knowledge of plant distribution and phenology and adequate assumptions regarding the sensitivity to meteorological factors, such as humidity, temperature, wind speed and turbulence. Unlike industrial pollutants or mineral aerosols, pollen emissions depend on not only the instantaneous meteorological conditions, but also the conditions during the pollen maturation within the plants before the pollination starts: it means that an accurate simulation is required to start several months before the targeted period. For all species, emissions schemes have the same basis:

$$\begin{array}{}\text{(43)}& E(x,y,t)=D(x,y)\times P(x,y,t)\times \mathit{\varphi }(x,y,t)\times R(x,y,t),\end{array}$$

where D(x,y) is the ragweed density distribution in number of individual plants per square metre. $P(x,y,t)$ is the annual production in grains per individual plant. $\mathit{\varphi }(x,y,t)$ is the phenology factor (in s⁻¹), considering its yearly integrated value is unity. This factor represents the knowledge of the start and end dates of the pollen season as well as the shape of these potential emissions. $R(x,y,t)$ is the daily or sub-daily weather-dependent release of pollen grains in the atmosphere which depends on the hourly (or daily) meteorological variables. $R(x,y,t)$ is unitless. These different terms correspond to two different temporal scales: D(x,y), $P(x,y,t)$ and $\mathit{\varphi }(x,y,t)$ represent “annual” information and $R(x,y,t)$ represents “short-term” information for which we want to evaluate correlation with the meteorological variables.

To calculate these terms, we used the scheme of Prank et al. (2013) (referred to as P2013 hereafter). It is the scheme used in SILAM (Sofiev et al., 2013). Originally developed for birch pollen, it was adapted to ragweed by Prank et al. (2013). It is also used for olive, grass and alder. This SILAM version is directly implemented in CHIMERE and used without any changes. This scheme is based on a double-heat-sum concept (Linkosalo et al., 2010). The temperature accumulates from the 60th day of the year (29 February or 1 March), and when the accumulated heat sum reaches the first threshold, the flowering season starts. The end of the flowering season occurs when there is no more pollen left in the catkins (the open-pocket principle). The model not only predicts the start and the end of the flowering season, but also relates the amount of pollen ready to be released to the stage of bud development. The actual pollen release during the flowering season is determined by meteorological factors (Sofiev et al., 2013). The model also takes into account the fact that different trees within the same model grid cell get ready to emit pollen at slightly different times, with a few days' difference. As a result, the model calculates the pollen emission flux as a function of time. The emission flux is determined by the spatial distribution of birch trees, by the heat sum threshold distributions, and by meteorological conditions during the flowering season.

In terms of seasonal emission flux modulation, three time periods are distinguished:

• Before the beginning of the flowering season, the emission rate is zero and the heat sum, H, is accumulated starting from a certain date, D_s, until it reaches some start-of-flowering threshold, H_start.

• During the flowering season, as H continues to accumulate, the pollen emission rate, R(t), grows.

• The pollen season is over as soon as there is no more pollen in the catkins (the open-pocket principle).

Summing up, the pollen emission is determined by the following factors:

1. The heat sum exceeds a certain threshold value, H_start, in a grid point (H(t)>H_start).

2. The temperature exceeds the cutoff value, T_co (T>T_co).

3. The meteorological parameters inhibit emissions, unless they have their optimal values (humidity and rain) or stimulate emissions via mechanical action (wind speed and convective turbulence).

4. The available pollen is in catkins (N_tot≠0).

5. Trees are flowering (F_in and F_out≠0).

The resulting emission flux, E, is a product of the terms described above:

$$\begin{array}{}\text{(44)}& \begin{array}{rl}E(t,i,j)& =N_{\mathrm{tot}}\times \mathit{\varphi }(i,j)\times {\displaystyle \frac{T-T_{\mathrm{co}}}{\mathrm{\Delta }H}}\times f_{\mathrm{wind}}(U,w^{*})\\ & \times F_{\mathrm{in}}\left({\displaystyle \frac{H}{H_{\mathrm{start}}}}\right)\times F_{\mathrm{out}}\left({\displaystyle \frac{R\left(t\right)}{N_{\mathrm{tot}}}}\right)\\ & \times f_{\mathrm{th}}(q,q_{\mathrm{low}},q_{\mathrm{high}})\times f_{\mathrm{th}}(P,P_{\mathrm{low}},P_{\mathrm{high}}),\end{array}\end{array}$$

where

• N_tot (grains m⁻² yr⁻¹) is the total number of pollen grains released from 1 m² during the whole year,

• ϕ(i,j) is the pollen plant fraction in grid cell (i,j),

• $\frac{T-T_{\mathrm{co}}}{\mathrm{\Delta }H}$ is a relative emission intensity as a linear function of temperature,

• F_in is the emission flux fraction depending on the relative heat sum,

• F_out is the emission flux fraction depending on pollen quantity emitted since the start of the flowering season,

• f_wind is the wind-speed-dependent correction,

• f_th is the coefficient modulating emissions by humidity and precipitation (between 0 and 1).

Once the flowering season starts, the actual pollen release is determined by meteorological conditions that can either biologically inhibit emissions with respect to an optimal level (humidity and precipitation) or stimulate emissions due to mechanical forces (wind and convective turbulence) on the daily basis. The key variables are the wind speed, U; the convective velocity scale, w^*; the relative humidity, q; precipitation rate, P; and daily evolving temperature, T.

Precipitation and humidity impacts on the model are governed by the following equation, where x_low may be P_low or RH_low and x_high may be P_high or RH_high. Depending on the pollen species, the values are fixed in Table 5.

$$\begin{array}{}\text{(45)}& f_{\mathrm{th}}=\left\{\begin{array}{ll}\mathrm{1}& \mathrm{if}x⩽x_{\mathrm{low}},\\ {\displaystyle \frac{x_{\mathrm{high}}-x}{x_{\mathrm{high}}-x_{\mathrm{low}}}}& \mathrm{if}{x}_{\mathrm{low}}<x⩽x_{\mathrm{high}},\\ \mathrm{0}& \mathrm{if}x>x_{\mathrm{high}}.\end{array}\right.\end{array}$$

Table 5 Fixed parameters for the pollen emissions. N_tot is the total number of grains per square metre and year, D_s is the flowering start in Julian days, T_co is the temperature cutoff value (in K) and Hc is the emissions height.

The impacts of wind speed, U, and convective turbulence, w^*, are described by

$$\begin{array}{}\text{(46)}& f_{\mathrm{wind}}=f_{\text{wind-max}}-\exp {\displaystyle \frac{-(U+w^{*})}{U_{\mathrm{satur}}}},\end{array}$$

where U_satur corresponds to the maximum value where the wind still strengthens emissions and f_wind-max describes the magnitude of wind impact. For birch, U_satur = 5 m s⁻¹ and f_wind-max = 1.5 (Table 5). All the trees of the same grid cell do not start flowering at the same time, with the transition between phenological stages being accomplished in a gradual manner. The gradual start of the flowering of the trees in a grid cell is described by the relative heat sum, $x=\frac{H\left(t\right)}{H_{\mathrm{start}}}$, with the corresponding term F_in(x) given by

$$F_{\mathrm{in}}\left(x\right)=\left\{\begin{array}{ll}\mathrm{0}& \text{if }x<\mathrm{1}-{\mathit{\delta }}_{\mathrm{in}},\\ \frac{x-\mathrm{1}+{\mathit{\delta }}_{\mathrm{in}}}{\mathrm{2}{\mathit{\delta }}_{\mathrm{in}}}& \text{if }\mathrm{1}-{\mathit{\delta }}_{\mathrm{in}}⩽x<\mathrm{1}+{\mathit{\delta }}_{\mathrm{in}},\\ \mathrm{1}& \text{if }x⩾\mathrm{1}+{\mathit{\delta }}_{\mathrm{in}},\end{array}\right.$$

where δ_in is the relative uncertainty in H_start. The parameter controlling the gradual flowering season termination is the number of grains remaining in catkins with respect to the initial total number of grains, N_tot, with $x=\frac{R\left(t\right)}{N_{\mathrm{tot}}}$. In the current model version, N_tot is fixed; for example, N_tot=10⁹ for birch. The corresponding term F_out(x) is given by

$$F_{\mathrm{out}}\left(x\right)=\left\{\begin{array}{ll}\mathrm{1}& \text{if }x<\mathrm{1}-{\mathit{\delta }}_{\mathrm{out}},\\ \frac{x+\mathrm{1}-{\mathit{\delta }}_{\mathrm{out}}}{\mathrm{2}{\mathit{\delta }}_{\mathrm{out}}}& \text{if }\mathrm{1}-{\mathit{\delta }}_{\mathrm{out}}⩽x<\mathrm{1}+{\mathit{\delta }}_{\mathrm{out}},\\ \mathrm{0}& \text{if }x⩾\mathrm{1}+{\mathit{\delta }}_{\mathrm{out}},\end{array}\right.$$

where δ_out is the relative uncertainty in N_tot.

To simulate the pollen emission flux, the tree spatial distribution maps for each species are required. For birch, it was compiled by Sofiev et al. (2006) from a combination of national forest inventories of western Europe and satellite images of broadleaf forests when the former information was not available. The forest inventories reveal irregularities, which are probably related to the specifics of national methodologies for forest surveys. The use of satellite images allows one to smooth out the data between the countries. Some interpolations were made based on assumptions like if between 62° N and 65° N, the fraction of birch is expected to be equal to that found in Finland (80 %). In France, birch is mainly present in the Massif Central (middle mountain region), the west (Aquitaine Basin), and the north-east (the Vosges). As noted by Sofiev et al. (2013), there are noticeable regional variations in the temperature sum threshold that are necessary to trigger the flowering season. Hence, the temperature sum thresholds were fitted into phenological and aerobiological observations independently inside 33 sub-regions of Europe. This yielded a temperature sum threshold map for the flowering season start.

6.2.3 The heat sum concept

The heat sum (HS) concept follows the concept of Chapman et al. (2014), which is based on biological days. Some input information is necessary to update HS at each model time step: the current 2 m temperature, the current day, its length in hours and the time step. The HS is updated only if the current day is after the first day, which in our case is always equal to 79 (i.e. 20 March). Two ramp functions are estimated, one depending on the 2 m temperature and the second one on the day length. All data used are the same as in SILAM, as described in Prank et al. (2013), using parameters defined in Deen et al. (1998). The HS is then equal to

$$\begin{array}{}\text{(47)}& \text{HS}=\text{HS}+r_T\times l_{\mathrm{h}}\times {\displaystyle \frac{\mathrm{\Delta }T}{\mathrm{86}\mathrm{400}}},\end{array}$$

with

$$\begin{array}{}\text{(48)}& r_t=\left\{\begin{array}{lr}{\displaystyle \frac{T-T_{\mathrm{low}}}{T_{\mathrm{opt}}-T_{\mathrm{low}}}}& \mathrm{if}{T}_{\mathrm{low}}<T<T_{\mathrm{opt}},\\ {\displaystyle \frac{T_{\mathrm{high}}-T}{T_{\mathrm{high}}-T_{\mathrm{opt}}}}& \mathrm{if}{T}_{\mathrm{opt}}<T<T_{\mathrm{high}},\end{array}\right.\end{array}$$

where the l_h variable is equal to 1 except if the day length (d_l) is greater than the photo period parameter, which is 14.5 here, and for the following values:

$$\begin{array}{}\text{(49)}& l_{\mathrm{h}}=\left\{\begin{array}{lr}\exp \left(\right(d_l-N_{\mathrm{h},\mathrm{ph}})\times \log (\mathrm{0.95}\left)\right)& \mathrm{if}\mathrm{11.5}<\mathrm{HS}<\mathrm{16},\\ \exp \left(\right(d_l-N_{\mathrm{h},\mathrm{ph}})\times \log (\mathrm{0.94}\left)\right)& \mathrm{if}\mathrm{16}<\mathrm{HS}<\mathrm{20.5},\end{array}\right.\end{array}$$

where the threshold temperature values are fixed. Note that this HS may be set to zero again depending on specific meteorological conditions as follows:

• T_2 m is lower than T_thr (here, T_thr=273.15 K).

• Daily mean T_2 m is lower than dayT_thr (here, dayT_thr=280.65 K). Note that in this model version, the daily mean 2 m temperature is the running average for the last 24 h.

• HS is lower than StartHSThr. This value is fixed to StartHSThr = 25.0 here following Deen et al. (1998).

With this scheme, there are no emissions during the night. The calculation of sunrise and sunset is necessary. At the end, note that emissions are injected in the first model level only in CHIMERE, whereas they are injected in the whole boundary layer in SILAM.

6.2.4 The ragweed emission

The ragweed pollen emission is a specific case. It is possible to use the P2013 scheme, but two others are also implemented. The first one is the scheme developed by Efstathiou et al. (2011) that corresponds to a modified version of Helbig et al. (2004). The second one is the scheme of Menut et al. (2021b), which corresponds to a modified version of Efstathiou et al. (2011).

The approach of Efstathiou et al. (2011), hereafter called E2011, is a mix between the schemes of Helbig et al. (2004) and Sofiev et al. (2006), adapting these formulations to ragweed. In the following, the specific notations of the publications are used to have a reference. The terminology is different from P2013 and Helbig et al. (2004), but the principle is the same. The pollen emission flux, E (grains m⁻² s⁻¹), is calculated as

$$\begin{array}{}\text{(50)}& E=c_e\times c^{*}\times K_e\times u_{*}.\end{array}$$

The flux depending on the surface grid cell and time is $E(x,y,t)$. The terms of the equation are c_e, a plant-specific factor, and c^*, the grain production factor. The release factor, R_e, is represented by $K_e\times u_{*}$. K_e is a time-varying factor that depends on weather and u_* is the friction velocity. The plant-specific factor, c_e, is calculated as

$$\begin{array}{}\text{(51)}& c_e=\mathrm{4}{c}_b\left({\displaystyle \frac{d}{S}}-{\displaystyle \frac{d^{\mathrm{2}}}{S^{\mathrm{2}}}}\right),\end{array}$$

where c_b is approximated to 10⁻⁴. S is the pollen season duration in days, with S=60. d is the Julian day, which varies during the simulation. In practice, the starting and ending dates of the ragweed period are fixed to 210 and 270 Julian days here, respectively. The grain production factor, c^*, is calculated as

$$\begin{array}{}\text{(52)}& c^{*}={\displaystyle \frac{q_{\mathrm{p}}}{\text{LAI}\times h_{\mathrm{c}}}},\end{array}$$

where q_p is the total production of grains per year and is 10⁹ grains m⁻². It corresponds to the maximum number of emitted grains. Day after day, q_p is reduced by the number of the already emitted grains the day before. LAI is the leaf area index. LAI is a map in E2011, but in this CHIMERE version, we are using the value of LAI = 3. h_c=1 m is the canopy height. The K_e variable is estimated by

$$\begin{array}{}\text{(53)}& K_e=K_{\mathrm{h}}\times K_{\mathrm{w}}\times K_{\mathrm{r}},\end{array}$$

with the three meteorological limiting factors: K_h for humidity, K_w for 10 m wind speed and K_r for precipitation. For relative humidity, the limiting factor is expressed as Eq. (45), with q_low=50 % and q_high=80 %. For precipitation, we use the fit in E2011:

$$\begin{array}{}\text{(54)}& K_{\mathrm{r}}=\mathrm{1.0}-{\displaystyle \frac{p}{\mathrm{0.5}}},\end{array}$$

with p being the precipitation rate (in mm h⁻¹). For the wind speed, K_w is based on the P2013 wind speed correction, previously described by Eq. (46).

The Menut et al. (2021b) scheme, hereafter called M2021, has the same formulation as the E2011 scheme, except that the release term is reformulated. The emission flux is expressed as

$$\begin{array}{}\text{(55)}& E=c_e\times c^{*}\times R_{\mathrm{TS}},\end{array}$$

where c_e and c^* are the same functions as in the E2011 scheme. R_TS is the modified instantaneous release factor. Based on the correlation results in M2021, it appears that the main driving factors are those related to thermodynamical processes – namely the 2 m temperature, T_2 m; vertical velocity scale, w_*; and short-wave radiation, SW_d. The pollen emissions may be moderated by the precipitation rate, P_r, and 2 m specific humidity, q_2 m.

The differences between birch and ragweed emissions could be explained by the plant typology itself: birch is a tree, with the pollen source up to 10 m above the ground. At this level, the wind may be considered a dominant process for the emission of grains. Ragweed rarely exceeds 1 to 2 m above the ground, where the wind speed is moderate. In this case, the dominant factor could be the temperature, considering the grains are emitted at the highest temperature when they are sufficiently dry (Holmes and Bassett, 1963). The precipitation rate is a limiting factor but not the most important one: even if it rains during the night, the grains can dry out and can be pulled off the plant in the morning. R_TS is thus estimated as

$$\begin{array}{}\text{(56)}& R_{\mathrm{TS}}={\displaystyle \frac{{\displaystyle \frac{T_{\mathrm{2}\mathrm{m}}}{T_{\mathrm{2}\mathrm{m},\mathrm{0}}}}\times {\displaystyle \frac{w_{*}}{w_{*,\mathrm{0}}}}\times {\displaystyle \frac{{\text{SW}}_d}{{\text{SW}}_{d,\mathrm{0}}}}}{r_{q_{\mathrm{2}\mathrm{m}}}+r_{\mathrm{Pr}}}},\end{array}$$

where the values of T_2 m, w_* and SW_d correspond to the mean daily values. These values are normalized in order to keep the release term non-dimensional. The normalization factors are $T_{\mathrm{2}\mathrm{m},\mathrm{0}}=\mathrm{10}$ °C, $w_{*,\mathrm{0}}=\mathrm{1}$ m s⁻¹ and ${\text{SW}}_{d,\mathrm{0}}=\mathrm{200}$ W m⁻².

In order to moderate these fluxes when meteorological conditions are not favourable, resistance terms are added. These resistances are mainly due to the 2 m specific humidity, q_2 m, and the precipitation rate, P_r. Each resistance is expressed as a sigmoid function ranging between 0 and 1, depending on the minimum and maximum value of the x parameter. The resistance has to reflect the fact that these parameters inhibit ragweed pollen emissions.

$$\begin{array}{}\text{(57)}& r_x=\mathrm{1}+\exp \left[{\displaystyle \frac{-b_{\mathrm{f}}(i_{\max }-i_{\min })}{\mathrm{2}}}\left({\displaystyle \frac{x}{x_{\max }-x_{\min }}}-\mathrm{1}\right)\right],\end{array}$$

with b_f being a constant, chosen here to be b_f=10, that determines the curve of the sigmoid function. i_min and i_max represent the range of the sigmoid and are chosen here to be i_min=0 and i_max=1 in order to use a normalized function for each resistance. The critical issue here is to choose the minimum and maximum value for each x meteorological parameter. These boundaries have to reflect the best possible range of variations in meteorological variables for all locations over Europe and for the whole year. The maximum values must be moderate enough in order to provide a realistic resistance: a too low of a maximum value would give a resistance of 1 too often, while a too high of a maximum value would give resistances that are too low. Based on all meteorological values used in this study, the boundaries for the 2 m specific humidity are q_2 m(min)=0 and $q_{\mathrm{2}\mathrm{m}}\left(\text{max}\right)={\mathrm{5.10}}^{-\mathrm{3}}$ g g⁻¹ and for the precipitation rate are P_r(min)=0 and P_r(max)=1.5 mm h⁻¹.

6.3 Changes in mineral dust emission

Two main novelties were developed in this model version. First, the way to calculate the fluxes is now following a sub-grid-scale variability, blending the percentages of soil in each percentage of land use and the Weibull distribution for the surface wind speed. It enables the calculation of a flux whatever the horizontal resolution is and as if it was at high resolution. Second, a link was added between vegetation fire emissions and mineral dust emissions (user's choice). We consider that vegetation fires are changing the soil and land use properties then create more erodible surfaces. In addition and during the fire, a vertical injection profile is applied to dust emission, considering that the fires induce a local convection.

6.3.1 Sub-grid-scale variability in the soil and land use

The principle of the calculation of emission flux was reshaped to be more efficient. The principle of sub-grid-scale variability is used. The main goal is to remove, as far as possible, the problem associated with the horizontal resolution when using emission schemes based on the threshold calculation principle. Thus, emission fluxes are calculated at the highest possible resolution and then summed up in the grid cell. For that, the input databases are pre-processed to have sub-grid tables, as presented in Fig. 4. In this figure, the big square represents a CHIMERE grid cell (here, for example, with a size of 8 km, with each sub-square at a resolution of 1 km). For each grid cell, we count the number of different land use classes and the number of different soil types. It defines a table as being a matrix and containing the relative part of each land use class in each soil type. For each value, we also attach the value of the saltation roughness length, z₀, and the MODIS erodibility, also with a 1 km resolution. After calculation of all values, the table is sorted as a function of the land use for each soil. Later, in CHIMERE, the sorted table enables us to optimize the calculation speed by exiting the loop when the land use relative part becomes zero.

Figure 4 The three datasets (kilometric resolution) are read together to build sub-grid-scale variability tables. For each cell, the relative percentage of land use for each soil type is calculated (ratiol) as well as the corresponding aeolian roughness length (z0chim).

The interesting part is that we can calculate mineral dust flux using threshold parameters such as aeolian roughness length as if we had a 1 km × 1 km horizontal resolution. With the sorted table, the calculation is fast. The weakness is that with this sub-grid-scale approach, we lost the spatial position of the high resolution: we know there is n % of land use, but we do not know where it is in the grid cell. Note it is not very important because we use the mean wind speed for the emissions calculation. For each sub-grid part, the emission flux calculation is done by iterating the Weibull distribution.

6.3.2 Impact of the fires on dust

Vegetation fire emission is parameterized using the emifiresCAMS program, which reads the CAMS daily data. For each day, the data provide information on the chemical species flux and on fire radiative power. Based on the hypothesis (Menut et al., 2022), the burnt area is deduced from this fire radiative power (FRP) value. The burnt area surface is cumulated from a starting date and then considered surface, having lost its vegetation. In this model version, the surface recovery is not taken into account, but it is considered that this recovery can take years, whereas we only simulate a few months during and just after the fires. The information needed is only valid using the CAMS vegetation fire emissions because the specific information of the burnt area surface was specifically developed in the CHIMERE emifiresCAMS program.

The fires have several impacts. First, short-term impacts are modelled and active only during the fire: (i) the local pyroconvection creates a surface pressure gradient and then increases the surface wind speed and (ii) the pyroconvection induces a local thermal, quickly mixing the surface emission in the whole atmospheric column. It has the effect of increasing the surface emissions and then transporting these emissions vertically, when the usual emission model only distributes fluxes into the first vertical model level. Second, long-term effects are also parameterized: by burning the surface, vegetation fires create additional barren surfaces, with more erodible properties for the mineral dust emissions.

7 Changes in deposition and chemistry

7.1 The settling velocity

The calculation of the settling velocity for aerosols has been updated using the AerSett v1 formulation (Mailler et al., 2023a). As described in Mailler et al. (2023b), the terminal fall speed of aerosols is calculated as

$$\begin{array}{ll}\text{(58)}& {\displaystyle }{\stackrel{\mathrm{̃}}{v}}_{\mathrm{\infty }}^{\mathrm{Stokes}}& {\displaystyle }=C_{\mathrm{c}}{\displaystyle \frac{D^{\mathrm{2}}\left({\mathit{\rho }}_p-{\mathit{\rho }}_{\mathrm{a}}\right)g}{\mathrm{18}\mathit{\mu }}},\text{(59)}& {\displaystyle }\stackrel{\mathrm{̃}}{R}& {\displaystyle }={\displaystyle \frac{{\mathit{\rho }}_{a}D{\stackrel{\mathrm{̃}}{v}}_{\mathrm{\infty }}^{\mathrm{Stokes}}}{\mathrm{2}\mathit{\mu }}},{\displaystyle }{v}_{\mathrm{\infty }}& {\displaystyle }={\stackrel{\mathrm{̃}}{v}}_{\mathrm{\infty }}^{\mathrm{Stokes}}\\ \text{(60)}& {\displaystyle }& {\displaystyle }\times \left[\mathrm{1}-{\left(\mathrm{1}+{\left({\displaystyle \frac{\stackrel{\mathrm{̃}}{R}}{\mathrm{2.440}}}\right)}^{-\mathrm{0.4335}}\right)}^{-\mathrm{1.905}}\right].\end{array}$$

In Eq. (58), the slip correction factor, C_c, is obtained following the classical Davies (1945) formulation, where $\mathit{Kn}=\frac{\mathrm{2}\mathit{\lambda }}{D}$ is the Knudsen number, D the particle diameter and λ the mean free path of molecules in the air:

$$\begin{array}{}\text{(61)}& C_{\mathrm{c}}=\mathrm{1}+\mathit{Kn}\left(\mathrm{1.257}+\mathrm{0.4}\exp \left(-{\displaystyle \frac{\mathrm{1.1}}{\mathit{Kn}}}\right)\right).\end{array}$$

As discussed in Mailler et al. (2023b), the value of v_∞ obtained from Eq. (60) always stays within 2 % of the “theoretical” v_∞ value obtained from the Clift and Gauvin (1971) parameterization but is 4 times as fast as a brute-force iterative calculation of v_∞. For small particles ($\stackrel{\mathrm{̃}}{R}<\mathrm{0.116}$), Eq. (60) is replaced by just $v_{\mathrm{\infty }}={\stackrel{\mathrm{̃}}{v}}_{\mathrm{\infty }}^{\mathrm{Stokes}}$, changing the result by less than 1 %. This typically applies to all particles with D<10 µm, while for large particles, Eq. (60) takes into account the Clift and Gauvin (1971) large-particle drag correction factor.

7.2 Impact of fires on dry deposition

A new process was added in the model. The dry deposition fluxes depend on the leaf area index for the calculation of the surface resistance, r_c (Simpson et al., 2012). More precisely, the bulk canopy resistance is expressed as

$$\begin{array}{}\text{(62)}& G_{\mathrm{c}}=\text{LAI}\times g_{\mathrm{sto}}+G_{\mathrm{ns}},\end{array}$$

where LAI is the leaf area index (in m² m⁻²), g_sto the stomatal conductance and G_ns the bulk non-stomatal conductance.

In the case of using fire emissions (user's choice), we consider a new impact on this model version: the fires destroy the vegetation and modify the soil properties (a soil that is drier and more erodible). It has an impact on not only the mineral dust emissions (as described in Sect. 6.3.2), but also the vegetation and then the LAI. The effect of the fires is cumulated in time for that moment from the beginning of each year. For the LAI, its decrease due to the fires is expressed as

$$\begin{array}{}\text{(63)}& \text{LAI}=\text{LAI}\times {\displaystyle \frac{S_{\mathrm{total}}-S_{\mathrm{burnt}}}{S_{\mathrm{total}}}},\end{array}$$

where S_total is the total surface of the grid cell when a fire is active or was previously diagnosed. S_burnt is the time-cumulated burnt surface.

7.3 Update to the ISORROPIA module

The ISORROPIA module has been used for a long time in the model in order to simulate the partitioning of inorganic aerosols (Nenes et al., 1998). This module was originally written in Fortran 77 and is used in both its forward and reverse modes in CHIMERE. In order to solve possible instability issues as well as increase the readability of the module, we opted to use a Fortran 90 version of ISORROPIA.

This was partially done by adapting the ISORROPIA module used in the GEOS-Chem (The International GEOS-Chem User Community, 2021) model for CHIMERE. Since GEOS-Chem only uses the forward mode of ISORROPIA, we opted to add the reverse mode into the module. The entirely translated version was then tested both outside (in a zero-dimensional configuration) and inside CHIMERE. It was concluded that the two modules provide the same results, while the Fortran 90 version of the module is more stable than the Fortran 77 version. This stable version is now part of the CHIMERE model.

8 Model validation

8.1 Simulation and validation setup

The model goes through continuous quality control validations. Six simulations are designed: first, with two model versions, (i) with the previous model version, v2020r3, and (ii) with this new version, v2023r1, and second, for each version, with three different configurations: (i) in offline mode and forced by ECMWF meteorological fields; (ii) in offline mode and forced by WRF 4.3; and (iii) in online mode, also forced by WRF 4.3. These three configurations are chosen since they show the results of the model with two major meteorological forcings as well as the effect of coupling with an online meteorological model on the simulated concentrations.

The simulations are performed over a domain that covers most of Europe and northern Africa. This domain is usually used for validation purposes since it takes into account a region containing all the possible emission sources in the model, i.e. anthropogenic, biogenic, biomass burning, sea salt and dust sources. The horizontal resolution of the model is 60 km × 60 km. The vertical resolution of the model has 15 levels, starting from the surface and going up to 300 hPa. The simulations are performed for the entire year of 2019 in order to take into account seasonal changes. For simulations using WRF, this model is forced by the CFSv2 meteorological inputs (Saha et al., 2011); for simulations using ECMWF, the meteorology from the IFS/ECMWF fields, with a 0.25×0.25° global resolution and 3-hourly time frequency, is used. Anthropogenic emissions come from the CAMS global emission inventory (v5.3; Soulie et al., 2024). Mineral dust, salt, DMS, anthropogenic, biogenic and NO_x emissions from lightning have all been calculated. The CAMS global reanalysis 6-hourly fields are used as initial/boundary conditions. The CHIMERE model is run with 10 aerosol size bins starting from 10 nm and going up to 40 µm.

For the calculation of statistical scores, the EvalTools module is used. It is a Python package developed by INERIS and Météo-France and freely available on a website (version 1.0.6r; https://opensource.umr-cnrm.fr/projects/evaltools, last access: 9 July 2024). Being also used by INERIS in the framework of the daily operational forecast using CHIMERE over France (PREV'AIR) and Europe (CAMS Copernicus), the advantage of using this tool is to be completely homogeneous with the validation made for the operational forecast.

Measurements used for these comparisons are collected from multiple sources. The model is validated on a regular basis with a vast amount of data for a large list of species downloaded from several sources. Starting with the meteorological data, 2 m temperature, wind speed, precipitation and relative humidity are compared. These measurements come from E-OBS, the British Atmospheric Data Centre network (BADC; http://data.ceda.ac.uk/badc, last access: 6 December 2023) and EBAS (https://ebas-data.nilu.no/, last access: 6 December 2023). The E-OBS dataset (prepared by the ECA&D project) provides daily data on a regular 0.1×0.25 grid covering the European area. For all variables provided by E-OBS, data for the European Environment Agency (EEA) ozone (O₃) stations are extracted from the gridded format and compared to the simulations.

For atmospheric pollutants, a long list of species is compared to the simulations on a regular basis, but only certain comparisons are presented in this paper. Ground-level measurements of O₃, PM₁₀ and PM_2.5 are downloaded from the European Environment Agency (EEA; https://www.eea.europa. eu, last access: 6 December 2023), EBAS and OpenAQ (https://openaq.org/, last access: 6 December 2023). The observations for aerosol optical depth (AOD) downloaded from the Aerosol RObotic NETwork global remote sensing network (AERONET; https://aeronet.gsfc.nasa.gov/, last access: 6 December 2023) are also included.

The comparisons have been performed for all types of stations. It should be noted that all the data mentioned here are at an hourly resolution, apart from the E-OBS data. Only stations with at least 90 % data availability threshold for the entire year of 2019 have been used. An exception has been given to the AOD data, where the data availability during the year has been fixed at 25 %. The number of stations used for each of these sources is given in Table 6.

We have also included validation for vertical profiles of pressure, temperature and O₃ using the data acquired from the World Ozone and Ultraviolet Radiation Data Centre (WOUDC; https://woudc.org/home.php, last access: 6 December 2023) resulting from ozonesonde measurements. There are 700 balloon flights inside the simulation domain that have sufficient data for the year from nine locations. Only stations with over 25 flights during the year were used for the comparisons.

Table 6 Number of stations for all data sources per species.

8.2 Statistical scores

Statistical scores are presented for O₃, PM_2.5 and PM₁₀ surface concentrations and aerosol optical depth (AOD). The scores are calculated for six model configurations: the previous model version, v2020r3, and the new model version, v2023r1, with three different meteorological forcings – ECMWF in offline mode, called offecm, and WRF in offline and online modes, respectively, called offwrf and onwrf.

Results for ozone surface concentrations are presented in Fig. 5. In the top part of the figure, the hourly mean bias is presented for the six simulations, with the three v2020r3 configurations on the left and the three v2023r1 configurations on the right. For each, results are calculated for the three datasets (EBAS, EEA and OpenAQ). With the old version, the bias is large and negative when compared to EBAS but close to zero when compared to EEA and OpenAQ. With the new version, the bias becomes low for EBAS and positive and much larger for EEA and OpenAQ. It means that independently of the surface dataset and independently of the offline or online mode, the new version simulates more hourly ozone than the previous one. This difference may reach +12.19 and +10.4 µg m⁻³ when v2023r1_onwrf simulations are compared to EEA and OpenAQ, respectively. It means that compared to the previous version, there is a positive shift of $\approx +\mathrm{9}$/+10 µg m⁻³ for the surface ozone concentrations.

For the correlation (bottom of Fig. 5), the values show an increase between the old model version and the new one. For example, for onwrf, correlation increases from 0.55 to 0.65 with EBAS, from 0.61 to 0.70 with EEA and from 0.63 to 0.71 with OpenAQ. This increase is noted for the three model configurations, with or without direct/indirect effects. On average, the new model version provides correlation with an increase of +0.07 for ozone.

Figure 5 Statistical scores for ozone.

Figure 6 presents the same kind of validation but for PM₁₀ surface concentrations. The mean bias is reduced for all simulations and changes from largely negative to slightly positive. For example and for the online coupling, the bias of −8.08 and −9.89 µg m⁻³ for the comparison to EEA and OpenAQ changes to −2.18 and −2.71 µg m⁻³ with the v2023r1 model version. On the other hand and by comparison to EBAS data, the bias changes from 0.35 to 4.6 µg m⁻³. On average, the surface concentrations of PM₁₀ increase by ≈5 µg m⁻³ for all model configurations of v2023r1 compared to v2020r3. As for ozone, the correlation systematically increases with the new model version. If we consider the validation by comparison to EEA data, we can see an increase from 0.46 to 0.50 when the meteorological forcing is offecm. The increase is much larger when the meteorological forcing is the WRF model, with, for example, an increase from 0.26 to 0.44 for onwrf between the old and the new model version. On average, the gain in correlation is roughly +0.05 when the forcing is offecm and +0.2 when the forcing is offwrf or onwrf.

Figure 6 Statistical scores for PM₁₀.

The same kind of validation is performed for aerosol optical depth, and results are presented in Fig. 7. The only dataset here is AERONET, providing a large geographical cover of the AOD. Regarding the bias, it is slightly reduced between the old and the new model version, but the differences between the old and the new version are very low. The bias was also smaller when using offwrf and onwrf in place of offecm. This result remains the same with the new model version. The correlation changes with the new version. With offecm, one can note an increase from 0.55 to 0.59. With offwrf and onwrf, the increase is much larger, with values from 0.33 and 0.32 for v2020r3 to 0.64 and 0.64 for v2023r1. Depending on the meteorological forcing, the net gain in AOD correlation ranges from +0.05 to +0.3. The use of the new version is always a gain.

Figure 7 Statistical scores for AOD.

Another way to validate the simulations with the new model version is to compare model outputs to vertical soundings. Results are synthesized in Fig. 8 for ozone. All soundings recorded during the simulation period (1 year) and over the whole domain are aggregated to calculate bias and correlation. Values are presented at the vertical model levels. The values for observations are thus presented as the mean average (the black bar and dot) and the standard deviation as a black dot in a coloured bar. At each level, the correlation is presented as coloured text: the colour corresponds to the model configuration with the best correlation. It enables us to easily visualize where model configurations outperform the average. For this simulation, the best scores are obtained with v2023r1 forced offline, with WRF close to the surface and then with the offline ECMWF version for the rest of the tropospheric column (except two model levels, where the best score is obtained with v2020r3). Globally, the best scores are with v2023r1, in the configuration forced by ECMWF. It is also noticeable that the standard deviation values at each model level are very close between observations and the model. It changes only close to the top, where the model variability becomes larger than the observations. At these levels, the boundary conditions play a role in this variability.

Figure 8 Comparison between observations and model results as a mean averaged vertical profile of ozone (in ppb). The data are averaged over the whole modelled period and for all available vertical soundings recorded over Europe.

9 Conclusions

The new CHIMERE version called v2023r1 is presented in this article. It contains both new processes and new numerical tools. More precisely for the researchers, the pollen emissions are now distributed and can be used for any local studies in Europe, the turbulence and the transport were updated with more recent parameterizations, and the impact of the fires on the surface was added. For the users performing long simulations or forecast, the integration of the XIOS model enables us to have a faster model, with a computational cost that is 40 % lower than the previous version. It is also possible to easily manage the output, including the possibility of writing sub-hourly fields, such as the maximum concentration during an hour.

The perspectives for the model are to continue to modularize the code and to optimize the calculation. The future of such a model being to go toward high-performance computing (HPC) architecture, it is important to follow the fast evolution in numeric methodology, including optimization. The use of XIOS for reading input data will also be implemented to allow some on-the-fly interpolation in the case of complex grids.

Appendix A: Tabulated roughness length and albedo

Table A1 USGS land use types and minimum and maximum values of momentum roughness length z₀.

Table A2 USGS land use types and tabulated albedos for bare and snow surfaces.