HERMESv3, a stand-alone multi-scale atmospheric emission modelling framework – Part 2: The bottom–up module

We describe the bottom–up module of the High-Elective Resolution Modelling Emission System version 3 (HERMESv3), a Python-based and multi-scale modelling tool intended for the processing and computation of atmospheric emissions for air quality modelling. HERMESv3 is composed of two separate modules: the global_regional module and the bottom_up module. In a companion paper (Part 1, Guevara et al., 2019a) we presented the global_regional module. The bottom_up module described in this contribution is an emission model that estimates anthropogenic emissions at high spatial- (e.g. road link level,) and temporal- (hourly) resolution using state-of-the-art calculation methods that combine local activity and emission factors along with meteorological data. The model computes bottom–up emissions from point sources, road transport, residential and commercial combustion, other mobile sources, and agricultural activities. The computed pollutants include the main criteria pollutants (i.e. NO_x, CO, NMVOCs (non-methane volatile organic compounds), SO_x, NH₃, PM₁₀ and PM_2.5) and greenhouse gases (i.e. CO₂ and CH₄, only related to combustion processes). Specific emission estimation methodologies are provided for each source and are mostly based on (but not limited to) the calculation methodologies reported by the European EMEP/EEA air pollutant emission inventory guidebook. Meteorologically dependent functions are also included to take into account the dynamical component of the emission processes. The model also provides several functionalities for automatically manipulating and performing spatial operations on georeferenced objects (shapefiles and raster files). The model is designed so that it can be applicable to any European country or region where the required input data are available. As in the case of the global_regional module, emissions can be estimated on several user-defined grids, mapped to multiple chemical mechanisms and adapted to the input requirements of different atmospheric chemistry models (CMAQ, WRF-Chem and MONARCH) as well as a street-level dispersion model (R-LINE). Specific emission outputs generated by the model are presented and discussed to illustrate its capabilities.

1 Introduction

The development of reliable anthropogenic emission inventories is crucial to understanding air pollution sources and designing effective emission abatement measures (Day et al., 2019). Emission inventories are also key inputs for air quality research and forecasting applications and represent one of the largest source of uncertainty in the air quality modelling chain (Russell and Dennis, 2000). For their use in atmospheric chemistry models, emissions need to be spatially distributed over a gridded domain, temporally resolved with (typically) hourly intervals and mapped to the species defined in the gas-phase and aerosol chemical mechanisms of the atmospheric chemistry model.

The last decades have seen large efforts to develop emission inventories for global and regional scales using a variety of input datasets and approaches. These inventories have become key in supporting scientific research and policymaking (e.g. the Air Quality Modelling Evaluation International Initiative, AQMEII; Pouliot et al., 2015). At a global scale, some of the most frequently used inventories are the Air Pollutants and Greenhouse Gases Emission Database for Global Atmospheric Research (EDGAR; Crippa et al., 2018) and the dataset derived from the Evaluating the Climate and Air Quality Impacts of Short-Lived Pollutants project (ECLIPSEv5.a; Klimont et al., 2017). There are also widely known regional emission inventories such as the European Monitoring and Evaluation Programme (EMEP) (Mareckova et al., 2019) and the TNO-MACC_III (Kuenen et al., 2014) or the Regional Emissions inventory in Asia (REAS), which covers China, Japan and other Asian countries (Kurokawa et al., 2013). More recently, and as part of the European Copernicus Atmosphere Monitoring Service (CAMS), new updated global and regional emission datasets covering both anthropogenic and natural sources have been developed (Granier et al., 2019).

These inventories provide estimates of emissions either globally or regionally for a variety of sectors, pollutants and years in a consistent way. However, they are usually limited to high-resolution modelling applications; assessing urban air quality or the local impact of emission reduction measures are two good examples (e.g. Timmermans et al., 2013). Such a limitation is due to the insufficient level of detail of the data used to estimate the emissions, typically national statistics, the uncertainties associated with their spatial and temporal distribution and the lack of flexibility for computing specific scenarios (e.g change in road speed limits). Emissions are first estimated at the annual and national level and the spatial proxies assigned to each pollutant source (e.g. population, city lights, land uses) are usually empirical and may not be representative of the real-world spatial emission patterns (Andres et al., 2016; Geng et al., 2017). Hourly emissions are computed through the application of temporal profiles (i.e. monthly, weekly and diurnal) to the original annual inventories, which are usually static (i.e. spatially constant) and do not account for the dynamical component of the emission processes (e.g. volatilization of ammonia as a function of meteorological parameters; Backes et al., 2016). Moreover, existing global and regional emission inventories are reported following sector aggregations (e.g. Gridded Nomenclature For Reporting, GNFR) that can hinder the application of detailed speciation profiles. For instance, the EMEP inventory, which compiles emissions from the parties of the Convention on Long-range and Trans-boundary Air Pollution (CLRTAP), reports all road transport emissions under one single category (i.e. GNFR F_RoadTransport) without discriminating by vehicle type (e.g. passenger cars, heavy duty vehicles), fuel type (i.e. petrol and diesel vehicles) or EURO (European emission standards) category. These factors become crucial when assigning the fraction of total nitrogen oxides (NO_x) emitted directly as nitrogen dioxide (NO₂) (e.g. Carslaw and Rhys-Tyler, 2013).

Dedicated emission models combining detailed databases with novel calculation methods that represent the main factors influencing the emission processes (e.g. meteorology, soil properties) can overcome these limitations. Some recent examples are the French VOLT'AIR model (Hamaoui-Laguel et al., 2014), which simulates ammonia (NH₃) volatilization fluxes after the application of fertilizer taking into account agro-environmental factors (e.g. meteorology, agricultural practices, soil properties), the Brazilian VEIN model (Ibarra-Espinosa et al., 2018), which provides bottom–up exhaust and evaporative vehicular emissions at street and hourly levels using different sets of emission factors (e.g. COmputer Programme to calculate Emissions from Road Transport version 5 (COPERT), EPA), and the Norwegian MetVed model (Grythe et al., 2019), which estimates residential wood combustion hourly emissions on a 250 m grid resolution considering several influencing factors such as outdoor temperature, type of available heating technologies and the number, type and size of dwellings.

Despite presenting highly accurate and detailed modelling methodologies, there are still some shortcomings associated with these modelling tools, mainly in terms of their usability for air quality modelling. On the one hand, each model covers only a specific pollutant sector, which means that they may have to be manually combined with other existing inventories. On the other hand, the models are usually designed to provide emissions for specific frameworks, which results in limitations in terms of the number of atmospheric chemistry models compatible with the emission outputs, map projections and type of working domains supported, as well as regions or countries in which the model can be applied. Finally, these tools are not always distributed under open-access licenses, which limits their usage within the scientific community.

This paper is the second part of the description of an open-source, Python-based, parallel, stand-alone and multi-scale atmospheric emission modelling framework. The High-Elective Resolution Modelling Emission System version 3 (HERMESv3) estimates atmospheric emissions for use in multiple air quality models (i.e. CMAQ, Appel et al., 2017; WRF-Chem, Grell et al., 2005; and MONARCH, Badia et al., 2017) as well as map projections and model grids (i.e. regular and rotated latitude–longitude, Lambert conformal conic, Mercator). The system reflects the learning from previous versions of HERMES developed by the Earth Sciences Department of the Barcelona Supercomputing Center (BSC) during the last decade (i.e. Baldasano et al., 2008; Ferreira et al., 2013; Guevara et al., 2013, 2017).

HERMESv3 is composed of two independent modules named global_regional and bottom_up. The global_regional module (HERMESv3_GR) is a customizable emission processing system that combines existing gridded inventories with user-defined vertical, temporal and speciation profiles for the generation of global and regional air quality model-ready emission files. A complete description of HERMESv3_GR can be found in Guevara et al. (2019a).

The bottom_up module, (described in this paper and further referred to as HERMESv3_BU) is an emission model that computes high spatial (e.g. road link, point source) and temporal (i.e. hourly) resolution anthropogenic emissions using state-of-the-art calculation methods that combine local activity and emission factors along with meteorological data. The model covers the estimation of bottom–up emissions from multiple sources, including power and manufacturing industries, road transport (i.e. exhaust and non-exhaust sources), residential and commercial combustion, other mobile sources (i.e. agricultural machinery, landing and take-off cycles at airports, shipping activities in ports, and recreational boats), and agricultural activities (manure management, fertilizer application and crop operations). The computed pollutants include the main criteria pollutants (i.e. NO_x; CO; NMVOCs; SO_x; NH₃; PM₁₀ and PM_2.5) and greenhouse gases (i.e. CO₂ and CH₄, only related to combustion processes). HERMESv3_BU provides specific estimation methodologies and emission factors for each source, which are mostly based on (but not limited to) the calculation methodologies reported by the European EMEP/EEA air pollutant emission inventory guidebook (EMEP/EEA, 2016). Users are allowed to load their own emission factors or apply specific tuning factors to the default dataset. With respect to the input activity data (e.g. geographic location of the industrial facilities and corresponding activity factors and temporal profiles), the user is responsible for providing the required information on the region of interest following the formats established in HERMESv3_BU. The application of the model is not restricted to a specific country; it can be run for any European region if the corresponding input data are provided. HERMESv3_BU includes a variety of global and regional state-of-the-art datasets to increase the usability of the tool and minimize the amount of input information that needs to be provided by the user. In the same line of thinking, HERMESv3_BU includes functionalities similar to geographic information systems (GISs) for automatically manipulating and performing spatial operations on geometric objects (e.g. remap spatial data from one spatial domain to another). The model counts as well with a flexible speciation mapping functionality, which allows the user to speciate the original pollutants to any desired chemical mechanism. Besides the aforementioned mesoscale atmospheric chemistry models that are compatible with HERMESv3, the road transport emission outputs of HERMESv3_BU can also be used within the R-LINE research-grade dispersion model (Snyder et al., 2013).

2 General description

2.1 Overview

A schematic of the model structure and execution workflow is shown in Fig. 1. The characteristics of the working domain, execution dates, emission sectors and pollutants to be included in the calculation process, input data paths, atmospheric chemistry output format, and number of calculation and writing processors have to be specified by the user in the general configuration file (see Sect. 2.2). Once all this information is compiled, HERMESv3_BU starts the general initialization process, which includes (i) the creation of the working grid where emissions will be calculated (grid function), (ii) the creation of the polygon feature that will be used to perform polygon clipping operations of the input data when needed (clip function) and (iii) the distribution of the computational resources among the different emission sectors (sector manager function). The data generated with the grid and clip functions, which are specific to each working domain, are stored as auxiliary files by default after their creation so that they can be reused in subsequent executions.

Figure 1 Schematic representation of the general structure of HERMESv3_BU. All the parameters described by the user in the general configuration file (e.g. working domain, execution dates, input data paths) are used to execute the initialization process (i.e. creation of working grid, clip polygon and auxiliary files and distribution of the computational resources). Once the general initialization is finished, HERMESv3_BU starts the execution of the emission core of the model. For each pollutant sector, and considering the input data provided by the user, HERMESv3_BU performs a (i) sector initialization, (ii) sector calculation, (iii) spatial mapping, (iv) temporal distribution and (v) speciation. Once the emission estimation of all sectors is finished, the writing core proceeds to gather all the data and writing it according to the atmospheric chemistry output format selected by the user.

The input data used by HERMESv3_BU include sector-specific and general data. The first set of information is divided into individual folders (<sector>), each one containing the activity and emission factor files of each sector, and the <profiles> folder, which contains the temporal, vertical and speciation profiles associated with each sector. Regarding the general data, this section includes meteorological files, which are used by several sectors, as well as global and regional datasets such as population maps. The different types of files used in the model are described in detail in Sect. 3.

The emission core of the model is composed of individual and independent submodules that calculate the emissions of each sector following specific estimation methodologies (see Sect. 3.1 to 3.6). Most of the submodules share a common procedure, which consist of the following steps: (i) sector initialization, which creates sector-specific auxiliary files, (ii) sector calculation, (iii) spatial mapping, (iv) temporal distribution and (v) speciation.

Once the execution of all the sectors is finished, HERMESv3_BU starts the data writing process. This function consists in first gathering the emissions estimated by each submodule (provided as 4D matrices with information on emissions across space, time and vertical levels) and secondly writing the merged data in an output NetCDF file according to the conventions of the atmospheric chemistry model of interest (see Sect. 2.5).

2.2 General configuration file

The general configuration options are passed to HERMESv3_BU via a configuration file, which is divided into seven different sections.

1. General: this section defines the main paths of the model (i.e. input, output, general data, auxiliary files), the name of the output emission file, time step configuration parameters (i.e. start and end dates and number of hourly time steps) and the option of removing the existing auxiliary files at the beginning of the execution.

2. Domain and output format: this section defines the characteristics of the working grid (e.g. spatial coverage, horizontal resolution and vertical description) as well as the structure and naming convention of the output NetCDF emission file. Currently, HERMESv3_BU supports four map projections (i.e. regular lat–long, rotated lat–long, Lambert conformal conic and Mercator) and three atmospheric chemistry model file conventions (i.e. CMAQ, WRF-Chem and MONARCH) (see Sect. 2.5).

3. Clipping: this section defines the polygon feature that will be used to perform the clipping operation during the general initialization process (see Sect. 2.4).

4. Sector management: this section defines the number of computational processors that will be assigned to each pollutant sector during the emission calculation process. Sectors can be individually deactivated setting their corresponding numbers to 0. This section also defines the number of processors that will be assigned for the writing process (see Sect. 4).

5. General shapefiles and raster: this section defines the path to the general shapefiles (i.e. administrative boundaries) and rasters (i.e. population, land use, livestock and soil property maps) used in the model.

6. Pollutant sector data: this section contains individual subsections for each pollutant sector, in which the user defines (i) the list of pollutants to be calculated, (ii) the data paths that point to the specific-sector input files (see Table B1) used for the emission calculation process (i.e. users can freely define their own file data storage convention) and (iii) an optional subset of pollutant categories to be considered for the calculation process and only available for certain sectors (i.e. list of vehicle categories in the road transport sector, list of crop categories in the fertilizer application sector, list of animal categories in the livestock sector, list of airport codes in the aircraft sector, list of port codes in the shipping sector, list of fuels in the residential and commercial sector). This last option can be very useful when studying the contribution of certain pollutant categories to total sectoral emissions (e.g. diesel vehicles in road transport) or when performing source attribution modelling studies (e.g. Pay et al., 2019).

7. Meteorology: this section define the paths to the gridded meteorological files used as input.

2.3 Input files and cross-referencing

HERMESv3_BU uses four types of input file formats.

• Esri shapefiles (points, polygons and polylines): used to provide spatial georeferenced information, including road transport networks, collections of point source facilities, infrastructure boundaries (i.e. airport, port) and administrative boundaries.

• Geotiff raster files: used to provide spatially gridded information, including land use information, population and livestock distributions, and soil properties (i.e. pH and cation exchange capacity).

• CSV files: used to store non-georeferenced activity data and emission factors as well as sets of temporal, vertical and speciation profiles.

• NetCDF files: used to provide modelled meteorological data (e.g. temperature, wind speed). HERMESv3_BU can currently use gridded meteorological data provided by MONARCH and ERA5 (C3S, 2017).

Depending on the pollutant sector one or more types of data are combined during the emission calculation process. For each sector, cross-references and spatial relationships between files are used for matching all the different input data. Figure 2 shows an example of the files used in the point source sector (see Sect. 3.1) and how all the information is linked. Most of the point source's input information (e.g. activity and emission factors, stack parameters, geographical coordinates) used by the model is provided in a multipoint shapefile, each row containing the information on a specific facility. The shapefile includes temporal and speciation profile IDs (e.g. “MXXX” for the monthly profiles, “XXX” being a three-digit numeric code that starts at “001”), which are cross-referenced with temporal and speciation CSV files where the numeric profiles are stored. On the other hand, the geographical coordinates of each facility are used to identify the closest grid cell of the meteorological NetCDF file and subsequently associate with them the required meteorological information.

Figure 2 Illustration of how the different input files used by the point source sector are linked and cross-referenced. The multipoint shapefile (point_source_facilities.shp), contains temporal and speciation profile IDs for each facility (e.g. “MXXX” for the monthly profiles, “XXX” being a three-digit numeric code that starts at “001”), which are cross-referenced with the corresponding temporal and speciation CSV files (month.csv, week.csv, hour.csv, speciation_cb05_aero5.csv). The shapefile also contains the geographical coordinates of each facility (long, lat), which are used to identify the closest grid cell of the meteorological NetCDF file (meteo_file.nc) and derive the surface temperature data.

The input data required to run HERMESv3_BU can be classified into three main categories.

• User-dependent data files: files that contain local information for the domain of study (e.g. energy consumption statistics, cultivated crop areas) and that need to be provided by the user.

• Built-in data files: store information that is not tied to a specific domain (e.g. emission factors, temporal profiles) and that is provided by default with the model. Users can modify the files provided by default if needed or add new ones (e.g. speciation profiles for a new chemical mechanism) but it is not mandatory in order to correctly run the model.

• External data files: open-source files reported by third parties (e.g. Joint Research Centre, Copernicus Land Monitoring Service) that contain global or regional information (e.g. population density, land use map) and that allow minimizing the amount of local information that needs to be provided by the user

All the model input data (user-defined, built-in and external) are needed to correctly execute the emission core of the system, as illustrated in Fig. 1. A classification of all the input files according to the three categories described above can be found in Table B1. Section 3 provides more information on the input datasets used for each pollutant sector. A complete description of all the input files and corresponding information fields can be found in the wiki of the model (see Sect. 6).

2.4 Spatial operations

HERMESv3_BU includes multiple functionalities for manipulating and performing spatial operations on geometric objects that are spatially referenced and have associated attributes (i.e. shapefiles and raster files) without requiring the users to have a GIS software. This allows the model to automatically manipulate georeferenced input data files, as well as to create spatial surrogates that have lately been used to map the estimated emissions onto the grid cells of the desired working domain. The following operations are implemented in the model:

• read, write or create – reads, writes and creates vector-based spatial data including Esri shapefiles and Geotiff raster files. These operations also allow changing the map projection of the original files.

• clip – overlays a polygon on a target feature and extracts from it only the data that lie within the area outlined by the clip polygon. By default, the clip polygon is defined as the outline of the user-defined working domain, but users can optionally use an existing shapefile (e.g. administrative boundary) or define a costume polygon providing a set of latitude–longitude coordinates. The clipped data become a new feature.

• unary union – returns a representation of the union of the given geometric objects. This function is used to create the outline of the user-defined gridded working domain.

• data conversion – converts raster files to a polygon feature class.

• spatial intersection – computes a geometric intersection of two features. The operation returns only those geometries that are contained by both targeted features. Unlike the clip operation, in the spatial intersection the associated attribute values from the input feature classes are copied to the output feature class.

• spatial difference – as opposed to the spatial intersection, in this case the operation returns only those geometries that are not contained by both targeted features.

• spatial join – joins attributes from one feature to another based on the spatial relationship. The target features and the joined attributes from the join features are written to the output feature class. Unlike the spatial intersection, the spatial join does not modify the geometry of the target feature.

• nearest point – Calculate the nearest point in a pair of geometries. This operation is used to assign to each emission source (e.g. point source, road link) the closest meteorological data reported in the NetCDF input files.

As an illustration, Fig. 3 shows the steps performed by HERMESv3_BU to generate the gridded fuel consumption data used by the residential and commercial combustion emission submodule (see Sect. 3.4). In the example, Spanish natural gas and wood consumption data obtained at the province level (IDAE, 2018; MITECO, 2018) are mapped onto a 4 km by 4 km regional Lambert conformal conic grid covering the Iberian Peninsula. For creating the gridded data, the model uses the population maps reported by the Global Human Settlement Layer (GHSL) project (JRC and CIESIN, 2015; Florczyk et al., 2019). The GHSL provides global Geotiff raster files at a resolution of 1 km by 1 km on the distribution and density of population, expressed as the number of people per cell (Schiavina et al., 2019), and on the classification of human settlements on the basis of the built-up area and population density, expressed as high- and low-density clusters (i.e. large and small urban areas, here remapped under a single category expressed as urban areas) and rural areas (Pesaresi et al., 2019). In the example, a clip of the original GHSL population density raster is performed using a shapefile of the administrative borders of Spain (Fig. 3a). The resulting clipped raster is converted to a polygon feature (Fig. 3b, zoom over the region of Madrid), to which new information is appended performing two spatial joins: one with a shapefile of the Spanish Nomenclature of Territorial Units for Statistics level 3 (NUTS3) administrative boundaries to append the province code to each source grid cell and another one with the GHSL settlement classification layer (that has also been previously converted from raster to shapefile) to append the population type information (Fig. 3c). Once each grid cell of the polygon has information on the population, the NUTS3 code and type of settlement, HERMESv3_BU spatially distributes the annual fuel consumption input data, which are provided by the user in a CSV file. For that, the following expression is applied (Eq. 1):

$$\begin{array}{}\text{(1)}& {\mathrm{FC}}_f\left(\stackrel{\mathrm{‾}}{x}\right)={\mathrm{FC}}_{f,n}\cdot {\displaystyle \frac{{\mathrm{Pop}}_{n,t}\left(\stackrel{\mathrm{‾}}{x}\right)}{\sum _{x=\mathrm{1}}^N{\mathrm{Pop}}_{n,t}\left(\stackrel{\mathrm{‾}}{x}\right)}},\end{array}$$

where ${\mathrm{FC}}_f(\stackrel{\mathrm{‾}}{x}$) is the annual fuel consumption (GJ yr⁻¹) of fuel f on the source grid cell $\stackrel{\mathrm{‾}}{x}$; FC_f,n is the total annual fuel consumption (GJ yr⁻¹) of fuel f in the NUTS3 n; and $\mathrm{Pop}(\stackrel{\mathrm{‾}}{x}{)}_n$ is the amount of population (inhabitant per cell) of type t (urban, rural) from NUTS3 n on the source grid cell $\stackrel{\mathrm{‾}}{x}$. In the example provided, urban and rural population are considered for the distribution of natural gas consumption (Fig. 3d), and only rural population for the distribution of wood consumption (Fig. 3e).

Figure 3 Examples of the spatial operations performed by HERMESv3_BU during the initialization of the residential and commercial combustion emission sector, including (a) a clip of the original GHSL population density raster (population per pixel) using a shapefile of the administrative borders of Spain, (b) a conversion of the clipped raster to a polygon feature (population per cell) (zoom over the area of Madrid), (c) the spatial joins performed to append the NUTS3 administrative boundary codes (ES300, Madrid; ES424, Guadalajara and ES425, Toledo) and the GHSL settlement categories (urban, rural) to each source grid cell, and the spatial intersections applied to remap the fuel consumption data (natural gas in urban and rural areas and wood in rural areas) (GJ per cell per year) from the source (d, f) domain to the user-defined (e, g) destination domain (4 km by 4 km Lambert conformal conic grid).

In the final step, the resulting polygon features are spatially intersected with the 4 km by 4 km gridded domain in order to remap the fuel consumption data from the source domain to the destination domain (Fig. 3f and g). The remapping is performed taking into account the ratio of the area of the region of intersection between the source and destination grid cells ($A\left(\stackrel{\mathrm{‾}}{x}x\right))$ to the total area of the source grid cell ($A\left(\stackrel{\mathrm{‾}}{x}\right))$, as expressed in Eq. (2):

$$\begin{array}{}\text{(2)}& {\mathrm{FC}}_f\left(x\right)={\mathrm{FC}}_f\left(\stackrel{\mathrm{‾}}{x}\right)\cdot {\displaystyle \frac{A(\stackrel{\mathrm{‾}}{x},x)}{A\left(\stackrel{\mathrm{‾}}{x}\right)}}.\end{array}$$

Similar operations are applied for the spatial manipulation of other georeferenced datasets such as land use categories, livestock maps or digitalized traffic networks.

2.5 Air quality model-ready files

HERMESv3_BU creates NetCDF emission output files in a format that is compatible with the conventions used by a number of air quality models, including CMAQ, WRF-Chem and MONARCH. For each model, an independent writing function has been implemented (e.g. writing_cmaq.py) to perform the required conversion of units and inclusion of mandatory global attributes. This modular approach allows us to easily extend the writing capabilities of the model to other atmospheric chemistry model conventions. Alternatively, the user can also estimate the emissions in a so-called DEFAULT format, which stores the emissions (g h⁻¹) in a NetCDF file that follows the Climate and Forecast (CF1.6) Metadata Conventions. In the case of road transport emissions, a dedicated writing function was designed so that the computed link-level emissions can be used by the R-LINE Gaussian dispersion model. This functionality allows HERMESv3_BU to be used for modelling air pollution at the urban (street level) scale (Benavides et al., 2019).

3 Emission sectors

Table 1 summarizes the major characteristics of each pollutant sector considered in HERMESv3_BU, including source type, categories and processes considered, pollutants involved and temporal, and vertical distribution. The following subsections provide a detailed description of the emission estimation methodologies implemented within HERMESv3_BU for each pollutant sector. Unless otherwise stated, all the equations reported in the following subsections are derived from the emission estimation expressions reported by EMEP/EEA (2016). Original expressions have been reformulated to produce high-resolution emissions (i.e. gridded or source-specific and hourly) instead of total national annual emissions. Some illustrative examples of the outputs that can be generated with the tool are also presented and compared against other existing emission datasets. In all the cases, the presented results were estimated for Spain or a Spanish region or city. The reason for this is mainly because the access to the required local, regional and national data is easy for the authors in this country. Compiling data with the same level of detail for other countries is beyond the scope of this work. Nevertheless, and as mentioned before, HERMESv3_BU is designed so that it can be applicable to other European countries or regions where similar input data are available.

Table 1 Summary of the main characteristics of each pollutant sector included in HERMESv3_BU.

¹ EMEP/EEA (2016) (chap. 1.A.3.b.i–iv, Tier 3 approach). ² EMEP/EEA (2016) (chap. 1.A.3.d, Tier 3 approach). ³ EMEP/EEA (2016) (chap. 1.A.3.a, Tier 3 approach).
⁴ EMEP/EEA (2016) (chap. 1.A.5.b, Tier 3 approach).

3.1 Point sources

This submodule estimates hourly emissions from process and combustion activities occurring in energy and manufacturing industrial point sources Eq. (3):

$$\begin{array}{}\text{(3)}& E_{p,i}\left(h\right)={\mathrm{AF}}_p\cdot {\mathrm{EF}}_{p,i}\cdot \mathrm{FM}(m{)}_p\cdot \mathrm{FW}(d{)}_p\cdot \mathrm{FH}(h{)}_p,\end{array}$$

where E_p,i(h) is the hourly emissions of pollutant i at point source p and hour h (g h⁻¹); AF_p is the annual activity factor (energy or material produced, fuel consumed) associated with point source p (GWh yr⁻¹ or GJ yr⁻¹ or grams of product yr⁻¹); EF_p,i is the emission factor linked to point source p and pollutant i (g GWh⁻¹ or g GJ⁻¹ or gram per gram of product); FM(m)_p is the monthly factor associated with month m and point source p (0 to 1); FW(d)_p is the weekly factor associated with day d and point source p (0 to 1); and FH(h)_p is the hourly factor associated with hour h and point source p (0 to 1).

As previously mentioned, most of the input data are provided in a georeferenced multipoint shapefile, each row containing the information for each specific facility (see Sect. 2.4). Emission factors are derived from facility-level emission reports when available, as recommended by the Tier 3 approach of EMEP/EEA (2016) (chap. 1.A.1 and 1.A.2). Alternatively, Tier 2 technology and fuel-dependent emission factors provided by the European guidelines are proposed. For each point source, emissions are horizontally allocated to the nearest grid cell of the destination working domain. Regarding the vertical allocation, HERMESv3_BU explicitly uses plume rise calculations to determine for each hour and each point source the effective emission heights. For this, the plume rise formulas as described by Gordon et al. (2018) are implemented. The algorithm takes into account stack and meteorological parameters, including stack height, stack diameter, exit temperature at the stack outlet, stack emission exhaust velocity, air temperature at stack height, wind speed at stack height, surface temperature, boundary-layer height, friction velocity and Obukhov length. Emissions are uniformly allocated across all the vertical layers that are included between the top and the bottom of the calculated plume. The plume rise function can be deactivated in HERMESv3_BU. In that case, the model will allocate the emissions to the layer closest to the stack height.

Alternatively to Eq. (3), HERMESv3_BU can directly ingest measured hourly emissions if available. For this, the user needs to provide a separate CSV file that contains the point source's measured emission fluxes per hour of the day (g h⁻¹) and to define the AF_p parameter in the shapefile as “−1”. This functionality becomes very relevant when assessing the impact of point source's plumes for specific days and under specific meteorological conditions (e.g. Baldasano et al., 2014).

Figure 4 shows an example of hourly and vertically distributed SO₂ emissions (kg h⁻¹) estimated by HERMESv3_BU for the As Pontes coal-fired power plant (Spain) during the months of January and July 2015. As Pontes is the largest power plant in Spain (1468.5 MW) and its exhaust stack (356 m) is the largest in the country and the second largest in Europe. The emission fluxes are directly derived from measurements reported by the Spanish Research Centre for Energy, Environment and Technology (CIEMAT, personal communication, 2018). The meteorological parameters for the plume rise calculations are obtained from the MONARCH model. It can be seen that there are significant differences between the vertical profiles obtained for January (winter) and July (summer), the emissions being injected to lower altitudes in the first case. The average plume thickness and plume top in January are 219.3 and 685.5 m, respectively, while in July the values are 259.4 and 745.7 m (18.3 % and 8.8 % larger, respectively). This is mainly due to meteorological differences between July and January in terms of air temperature at the stack height (+6.5 ^∘C) and boundary-layer height (−430.5 m). The results are in line with other plume rise calculations performed in other facilities located in similar climate zones (Bieser et al., 2011).

Figure 4 (a) Location of the As Pontes Spanish coal-fired power plant (sources of the ArcGIS Online basemap: Esri, DigitalGlobe, GeoEye, i-cubed, USDA FSA, USGS, AEX, Getmapping, Aerogrid, IGN, IGP, swisstopo and the GIS user community); (b) modelled average plume thickness (Δh) (m), air temperature at the stack height (T) (^∘C), boundary-layer height (planetary boundary layer, PBL) (m) and wind speed at the stack height (W) (m s⁻¹) for January and July (monthly averages); results of the hourly and vertically distributed SO₂ emissions (kg h⁻¹) estimated by HERMESv3_BU for the months of (c) January and (d) July 2015. The map in (a) was created using ArcGIS^® software by Esri. ArcGIS^® and ArcMap™ are the intellectual property of Esri and are used herein under license. Copyright © Esri. All rights reserved. For more information about Esri^® software, please visit https://www.esri.com (last access: February 2020).

3.2 Road transport

3.2.1 Hot and cold exhaust

Hot-exhaust emissions are estimated following Eq. (4):

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\mathrm{Ehot}}_{l,i}\left(h\right)=\sum _{v=\mathrm{1}}^n{\mathrm{AADT}}_{v,l}\cdot L_l\cdot \mathrm{EFhot}\left(V\right(h{)}_l{)}_{v,i}\\ & \cdot \mathrm{Mcorr}\left(V\right(h{)}_l{)}_v\cdot \mathrm{FM}(m{)}_l\cdot \mathrm{FW}(d{)}_l\cdot \mathrm{FH}(h{)}_l,\end{array}\end{array}$$

where Ehot_l,i(h) is the hourly hot-exhaust emissions of pollutant i at road link l and hour h (g h⁻¹); AADT_v,l is the annual average daily traffic for vehicle category v at road link l (no. vehicles per day); L_l is the length of the road link l (km); EF(V(h)_l)_v,i is the hot-exhaust emission factor linked to vehicle category v and pollutant i (grams per kilometre per no. of vehicles) as a function of the hourly mean vehicle travelling speed (V(h)_l at road link l and hour h (km h⁻¹); Mcorr(V(h)_l)_v is the mileage correction factor associated with vehicle category v, also estimated as a function of the hourly mean travelling speed; FM(m)_l is the monthly factor associated with month m and road link l (0 to 12); FW(d)_l is the weekly factor associated with day d and link l (0 to 1); and FH(h)_l is the hourly factor associated with hour h and link l (0 to 1). The number of vehicle categories is n.

Most of the activity input data (e.g. average daily traffic flow, mean vehicle speed) are provided in a multiline shapefile, each row containing the information on a specific road link. The shapefile includes vehicle fleet composition, temporal and speciation profile IDs, which are cross-referenced with the corresponding CSV files where all the numeric profiles are stored (similarly to the example shown for point sources; see Sect. 2.3). These profiles are used to distribute the total traffic flow among the different vehicle categories, temporally disaggregate the traffic flow and average speed at the hourly level and speciate the estimated emissions, respectively.

Both the emission and mileage correction factors implemented in HERMESv3_BU are the ones reported by the Tier 3 methodology of EMEP/EEA (2016) (chap. 1.A.3.b.i–iv), which correspond to the values reported by the European COPERT 5 (https://copert.emisia.com/, last access: February 2020). A total of 491 vehicle categories are considered, discriminated by vehicle type (i.e. mopeds, motorcycles, passenger cars, light duty vehicles, heavy duty vehicles, buses), fuel type (i.e. diesel, gasoline, liquefied petroleum gas (LPG), hybrid, electric), EURO category, engine power and gross weight class.

In the case of cold-start emissions, the calculation expression used is the following one (Eq. 5):

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& {\mathrm{Ecold}}_{l,i}\left(h\right)=\sum _{v=\mathrm{1}}^n{\mathrm{Ehot}}_{l,i,v}\left(h\right)\cdot {\mathit{\beta }}_{i,v}\left(T\right(h{)}_l)\\ & \cdot \left({\mathrm{Qcold}}_{l,i,v}\right(V(h{)}_l,T(h{)}_l)-\mathrm{1}),\end{array}\end{array}$$

where Ecold_l,i(h) is the hourly cold-exhaust emissions of pollutant i at road link l and hour h (g h⁻¹); ${\mathrm{Ehot}}_{l,i,v}\left(h\right)$ is the hourly hot-exhaust emissions of pollutant i for vehicle category v at road link l and hour h (g h⁻¹); β_i,v( T(h)_l ) is the fraction of mileage driven with a cold-engine pollutant i and vehicle category k (0 to 1) as a function of the hourly outdoor temperature T(h)_l  for hour h and road link l (^∘C); and ${\mathrm{Qcold}}_{l,i,v}\left(V\right(h{)}_l$, T(h)_l) is the cold or hot emission quotient for pollutant i and vehicles category v [≥1] as a function of the hourly outdoor temperature T(h)_l (^∘C) and hourly mean vehicle travelling speed (V(h)_l for hour h and road link l (km h⁻¹). The number of vehicle categories is n.

As in the case of the hot-exhaust emissions, the β_i,v(T(h)_l ) and ${\mathrm{Qcold}}_{l,i,v}\left(V\right(h{)}_l,T(h{)}_l$) parameters are estimated following the expressions and constants reported by the EMEP/EEA (2016) Tier 3 methodology (chap. 1.A.3.b.i–iv).

Besides the COPERT 5 constants used to calculate vehicle- and pollutant-specific hot and cold emission factors, HERMESv3_BU also includes scaling factor parameters (defined as 1 by default) that the user can modify to tune the original emission factors. This functionality can be useful for adjusting the default factors based on the insights reported by measurements performed under real-world driving conditions (e.g. underestimation of COPERT NH₃ cold-start emissions according to Suarez-Bertoa et al., 2017). The mileage correction factors reported by COPERT 5, which only apply to gasoline vehicles, are expanded to diesel vehicles (i.e. deterioration of tailpipe NO_x emissions of 22 % and 10 % on EURO 2 and 3 diesel passenger cars), following the results reported by Chen and Borken-Kleefeld (2016).

The estimated link-level vehicle emissions are mapped onto the user-defined gridded working domain by applying a spatial intersection. Once the intersection is performed, emissions are automatically gathered at the grid cell level and the total sum is computed. Figure 5 shows an example of the hourly PM_2.5 road transport emissions estimated for an area of Barcelona city (09:00 UTC), both at the road link level (kg km ⁻¹ h⁻¹, Fig. 5a) and grid cell level (1 km×1 km) (kg h ⁻¹, Fig. 5b). The total annual NO_x and PM₁₀ road transport emissions were estimated for the city of Barcelona using HERMESv3_BU and the results were compared against the latest available local emission inventory developed by the Barcelona City Council (AB, 2015) (Fig. 5c). Information on the traffic flow data was obtained from the local automatic traffic counting network (Barcelona city council, mobility and transport department, personal communication, 2017) and the TomTom historical average speed profiles product (https://www.tomtom.com, last access: February 2020), whereas vehicle fleet composition profiles were derived from a remote-sensing campaign (RACC, 2017). It is observed that HERMESv3_BU results are 60 % and 39 % higher than the ones reported by Barcelona City Council (AB). This is due to a combination of several factors, including (i) the different years of reference (2017 for HERMESv3_BU and 2013 for AB), (ii) the inclusion of Barcelona's port-area-associated road transport in HERMESv3_BU (large amount of heavy duty vehicles that contribute more than 250 and 15 t yr⁻¹ of NO_x and PM₁₀, respectively), (iii) the use of COPERT 5 real-world adjusted NO_x emission factors for EURO 5 and 6 diesel vehicles in HERMESv3_BU (AB inventory is based on COPERT 4, which does not consider the dieselgate effect), and (iv) the consideration of deterioration factors on old diesel passenger cars in HERMESv3_BU.

Figure 5 Hourly PM_2.5 road transport emissions estimated for an area of Barcelona city (09:00 UTC) at (a) the road link level (kg km⁻¹ h⁻¹) and (b) grid cell level (1 km×1 km) (kg h⁻¹) (sources of the ArcGIS Online basemap: Esri, DigitalGlobe, GeoEye, i-cubed, USDA FSA, USGS, AEX, Getmapping, Aerogrid, IGN, IGP, swisstopo and the GIS user community). (c) Barcelona city total annual NO_x and PM₁₀ road transport emissions (t yr⁻¹) estimated by HERMESv3_BU and reported by the Barcelona City Council (AB, 2015). Maps in (a) and (b) were created using ArcGIS^® software by Esri. ArcGIS^® and ArcMap™ are the intellectual property of Esri and are used herein under license. Copyright © Esri. All rights reserved. For more information about Esri^® software, please visit https://www.esri.com (last access: February 2020).

3.2.2 Non-exhaust (wear and resuspension)

HERMESv3_BU also estimates non-exhaust PM₁₀ and PM_2.5 traffic emissions, including road surface, tyre and brake wear, and resuspension. Emissions derived from processes of abrasion are estimated following Eq. (6):

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& {\mathrm{Ewear}}_{l,i}\left(h\right)=\sum _{v=\mathrm{1}}^n{\mathrm{AADT}}_{v,l}\cdot L_l\cdot {\mathrm{EFwear}}_{v,i}\cdot S\left(V\right(h{)}_l)\\ & \cdot \mathrm{FM}(m{)}_l\cdot \mathrm{FW}(d{)}_l\cdot \mathrm{FH}(h{)}_l,\end{array}\end{array}$$

where Ewear_l,i(h) is the hourly emissions of pollutant i at road link l and hour h (g h⁻¹); AADT_v,l is the annual average daily traffic for vehicle category v at road link l (n d⁻¹); L_l is the length of the road link l (km); EFwear_v,i is the emission factor linked to vehicle category v and pollutant i (g km⁻¹); S(V(h)_l) is the correction factor (0.902 to 1.39) estimated as a function of the hourly mean vehicle travelling speed (V(h)_l at road link l and hour h (km h⁻¹); FM(m)_l is the monthly factor associated with month m and road link l (0 to 12); FW(d)_l is the weekly factor associated with day d and link l (0 to 1); and FH(h)_l is the hourly factor associated with hour h and link l (0 to 1). Both the emission and correction factors are derived from the Tier 2 methodology proposed by EMEP/EEA (2016) (chap. 1.A.3.b.vi, Tables 3-4, 3-6 and 3-8). The number of vehicle categories is n.

In the case of resuspension, emissions are estimated as follows (Eq. 7):

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& {\mathrm{Eresus}}_{l,i}\left(h\right)=\sum _{v=\mathrm{1}}^n{\mathrm{AADT}}_{v,l}\cdot L_l\cdot {\mathrm{EFresus}}_{v,i}\\ & \cdot S\left(\mathrm{Hrain}\right(h{)}_l)\cdot \mathrm{FM}(m{)}_l\cdot \mathrm{FW}(d{)}_l\cdot \mathrm{FH}(h{)}_l.\end{array}\end{array}$$

All the parameters used in the expression are the same as the ones defined in Eq. (6) except for the resuspension emission factor linked to vehicle category v and pollutant i (g km⁻¹) (EFresus_v,i) and the correction factor S(Hrain(h)_l) (0 to 1). The resuspension emission factors proposed by default are vehicle type dependent (i.e. motorcycles, passenger cars, light duty vehicles, heavy duty vehicles) and derived from a measurement campaign performed in Barcelona (Amato et al., 2012a). The correction factor is estimated as a function of the number of hours after a precipitation event at road link l and hour h (Hrain(h)_l), following the expression reported by Amato et al. (2012b) (Eq. 8):

$$\begin{array}{}\text{(8)}& S\left(\mathrm{Hrain}\right(h{)}_l)=(\mathrm{1}-e^{-r\cdot \mathrm{Hrain}(h{)}_l}).\end{array}$$

The formula, which is based on measurements undertaken in Barcelona (Spain) and Utrecht (the Netherlands), indicates that after a rainfall (when the mobility particles drop to values close to zero), the loading of mobile road dust mobility increases exponentially tending to reach again the maximum emission strength. The equation depends on a recovery rate (r) that varies according to the traffic characteristics and local climatic conditions. By default, HERMESv3_BU uses the recovery rate value derived from the Barcelona measurements, but the user can change it to other values if desired. The effect of precipitation on resuspension emissions is only applied when at least a 0.254 mm h⁻¹ rainfall occurs (US EPA, 2011).

3.2.3 Gasoline evaporation

NMVOC evaporative diurnal emissions are considered in HERMESv3_BU as follows (Eq. 9):

$$\begin{array}{}\text{(9)}& E_i(x,h)=\sum _{v=\mathrm{1}}^{n}N(x{)}_v\cdot \mathrm{EF}(T(x,d){)}_v\cdot \mathrm{FH}\left(T\right(x,h\left)\right),\end{array}$$

where E_i(x,h) is the hourly emissions of NMVOC at the destination grid cell x and hour h (g h⁻¹); N(x)_v is the number of registered vehicles of category v in the destination grid cell x (no. vehicles);  EF(T(x,d))_v is the emission factor for vehicle category v (grams per no. of vehicles) as a function of the daily mean outdoor temperature T(x,d) for day d and destination grid cell x (^∘C); and FH(T(x,h)) is the hourly factor associated with hour h as a function of the hourly mean outdoor temperature T(x,h) for hour h and destination grid cell x. The number of vehicle categories is n.

The gridded number of registered vehicles (N(x)_v) is obtained combining the GHSL gridded population map with information provided by the user on registered gasoline vehicles at NUTS level 3, following the spatial operations showed in Sect. 2.4. Contrary to the exhaust and wear emissions, evaporative emissions are considered as an area source and directly computed at the grid cell level.

Emission factors are derived from the Tier 2 method of EMEP/EEA (2016) (chap. 1.A.3.b.v, Tables 3-5 and 3-6). Neither running loss nor hot-soak emissions are currently considered in HERMESv3_BU. This is due to the fact that these emissions mainly occur in gasoline vehicles with carburettors, and the fraction of European passenger cars and light duty vehicles post-EURO 1 with this technology is almost zero.

3.3 Agriculture

3.3.1 Fertilizer application

Hourly and spatially disaggregated NH₃ emissions from agricultural fertilizers are estimated following the expression reported by Paulot et al. (2014) (Eq. 10):

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& E(x,h)={\sum }_{c=\mathrm{1}}^{n}A(x{)}_c\cdot C_c\cdot \mathrm{\Gamma }(x{)}_c\cdot \mathrm{EF}(x{)}_c\\ & \cdot \mathrm{FD}(x,d{)}_c\cdot \mathrm{FH}(h),\end{array}\end{array}$$

where E(x,h) is the hourly NH₃ emissions at destination grid cell x and hour h (g h⁻¹); A(x)_c is the annual cultivated area of crop c at destination grid cell x (ha yr⁻¹); C_c is the ratio of cultivated to fertilized area for crop c (0 to 1); Γ(x)_c is the fertilizer application rate for crop c at destination grid cell x (kg N ha⁻¹); EF(x)_c is the emission factor for crop c at destination grid cell x (g NH₃ kg N⁻¹); FD(x,d)_c is the daily factor for crop c at destination grid cell x and day d (0 to 1); and FH(h) is the hourly factor associated with hour h (0 to 1). The number of crop categories is n.

The distribution of the cultivated crop areas onto the destination grid cells (A_c(x)) is performed using the spatial operation capabilities of HERMESv3_BU (see Sect. 2.4). The model combines the land use Geotiff raster reported by the CORINE Land Cover (CLC) inventory 2018 version 18 at 250×250 m (CLMS, 2018) with cultivated crop area statistics at NUTS level 2 provided by the user. HERMESv3_BU performs a mapping between the different CLC and crop categories (Table A1) in order to spatially distribute the statistics across the space. One limitation of this approach is that the number of agricultural land use categories in CLC is limited and therefore certain crops are assigned to the same CLC category (e.g. maize, barley, wheat, oat and rye crop categories are all mapped to the “Permanently irrigated land” CLC category). Future work will include exploring the use of more detailed datasets such as the crop type map product included in the Sentinel-2 for Agriculture portfolio (http://www.esa-sen2agri.org, last access: February 2020), which provides maps of the main crop types at 10 m resolution based on Sentinel-2 and Landsat-8 imagery.

The emission factors (EF(x)_c) are calculated following the methodology proposed by Bouwman and Boumans (2002), which determines that the NH₃ volatilization is driven by soil pH and cation exchange capacity (CEC), the type of fertilizer used (e.g. urea, ammonium, ammonium sulfate, manure), type of crop (i.e. upland, flooded), and application mode (i.e. broadcast or injection). HERMESv3_BU takes the soil parameters from the International Soil Reference and Information Centre (ISRIC) World Soil Information database (Hengl et al., 2017), which reports global pH soil and CEC maps at 250 m resolution. Original maps are remapped onto the user-defined grid applying a spatial intersection operation. All crops are assumed to be upland except for rice, and the application mode is assumed to be broadcast in all cases, following Paulot et al. (2014). The input on the type of fertilizer used can be distinguished by crop type and NUTS level 2.

The daily factors (FD(x,d)_c) are estimated following the dynamical ammonia emission parameterization reported Gyldenkærne (2005) and Skjøth et al. (2004, 2011), which is dependent on the outdoor air temperature, wind speed and timing of the fertilizer application, the last parameter being described with a Gauss function (Eq. 11):

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& \mathrm{FD}(x,d{)}_c=e^{\mathrm{0.0223}\cdot T(x,d)+\mathrm{0.0419}\cdot \mathrm{WS}(x,d)}\\ & \cdot {\sum }_{a=\mathrm{1}}^{\mathrm{3}}{\displaystyle \frac{{\mathit{\beta }}_{a,c}}{{\mathit{\sigma }}_{c,a}\cdot \sqrt{\mathrm{2}\cdot \mathit{\pi }}}}\cdot e^{\left(\frac{(d-{\mathit{\tau }}_{c,a}{)}^{\mathrm{2}}}{-\mathrm{2}\cdot {\mathit{\sigma }}_{c,a}^{\mathrm{2}}}\right)},\end{array}\end{array}$$

where T(x,d) is the 2 m outdoor temperature at destination grid cell x and day d (^∘C); WS(x,d) is the 10 m wind speed at destination grid cell x and day d (m s⁻¹); β_a,c is the fraction of fertilizer applied to crop c at stage a (1: planting; 2: at growth; 3: after harvest); τ_c,a is the optimal application date for crop c at stage a (Julian day, $\mathrm{1}:\mathrm{365}/\mathrm{366}$); σ_c,a is the deviation around date τ_c,a (number of days); and d is the day of the year (Julian day, $\mathrm{1}:\mathrm{365}/\mathrm{366}$).

Concerning the hourly distribution of emissions, the fixed temporal profile for the agriculture sector reported by Denier van der Gon et al. (2011) is proposed by default.

Figure 6a shows the results for the Spanish total annual NH₃ fertilizer emissions calculated on a Lambert conformal conic grid of 4 km by 4 km resolution. Cultivated crop area statistics were obtained from MAPA (2017a), and information on the fertilizer application rate and type of fertilizer used by crop type were derived from both MAPA (2011, 2017b) and Mueller et al. (2012). The time series of daily NH₃ emissions for the region of Aragon and Catalonia is plotted in Fig. 6b. To calculate the daily distribution, the fraction of fertilizers applied to each crop and stage are obtained from Paulot et al. (2014), and the application dates and deviation values are derived from multiple sources, including Gyldenkærne et al. (2005), Sacks et al. (2010) and Skjøth et al. (2011). In the particular case of barley, rye, wheat, oats and maize, the growth application dates are determined using the accumulated growing degree days (GDDs) since planting (McMaster and Wilhelm, 1997). Meteorological parameters were derived from the ERA5 dataset (C3S, 2017). As shown in the example, HERMESv3_BU is able to discriminate the emission estimation by crop type, which allows quantifying the contribution of each source to the total NH₃. The annual emissions estimated for this Spanish region, which is considered to be one of the Europe's main NH₃ hotspots, are compared against the emission fluxes derived from Infrared Atmospheric Sounding Interferometer (IASI) satellite observations (Van Damme et al., 2018) (Fig. 6d). It is observed that results are in agreement, the emissions reported by HERMESv3_BU being just −6 % lower to the ones derived from the IASI instrument. It is important to note that for this comparison, livestock emissions estimated by HERMESv3_BU have also been considered (see Sect. 3.3.2).

Figure 6 Spanish annual NH₃ (a) fertilizer and (c) livestock emissions (t yr⁻¹) calculated on a Lambert conformal conic grid of 4 km by 4 km resolution and corresponding time series of daily NH₃ (t d⁻¹) for the regions of (b) Aragon and Catalonia per crop type and (d) Murcia per livestock category; (e) the total NH₃ emissions (fertilizers+livestock) estimated for these two hot spots (t yr⁻¹) are compared against IASI satellite-derived NH₃ emission fluxes (Van Damme et al., 2018).

3.3.2 Livestock

Hourly gridded NH₃ emissions derived from manure management activities are estimated according to the expression reported by Paulot et al. (2014) (Eq. 12):

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& E(x,h)=\sum _{a=\mathrm{1}}^{a=\mathrm{4}}{\sum }_{l=\mathrm{1}}^{l=n}D(x{)}_l\cdot {\mathrm{\Gamma }}_l\cdot {\mathit{\beta }}_l\cdot {\mathit{\gamma }}_{a,l}\cdot {\mathrm{EF}}_{a,l}\\ & \cdot \mathrm{FD}(x,d{)}_{a,l}\cdot \mathrm{FH}(h),\end{array}\end{array}$$

where E(x,h) is the hourly NH₃ emissions at destination grid cell x and hour h (g h⁻¹); D(x)_l is the animal density for livestock category l at destination grid cell x (no. of heads per cell); Γ_l is the nitrogen (N) excretion rate for livestock l (grams of N per head); β_l is the fraction of total ammoniacal nitrogen (TAN) content of the excreta from livestock l (0 to 1); γ_p,l is the fraction of total excreta associated with activity a (1: housing; 2: yarding; 3: storage; 4: grazing) for livestock l (0 to 1); EF_a,l is the emission factor for livestock l and activity a (g NH₃ g N⁻¹); $\mathrm{FD}(x,d{)}_{a,l}$ is the daily factor for livestock l and activity a at destination grid cell x and day d (0 to 1); and FH(h) is the hourly factor associated with hour h (0 to 1). The number of livestock categories is n.

The estimation methodology and emission factors follow the Tier 2 approach proposed by EMEP/EEA (2016) (chap. 3.B, Table 3.9), which uses a mass-flow approach based on the concept of a flow of TAN through the manure management system.

Regarding the animal density data (D(x)_l), HERMESv3_BU uses as a basis the gridded livestock population from the Gridded Livestock of the World version 3 (GLWv3; Gilbert et al., 2018), which provides raster maps of global population densities of cattle, buffaloes, horses, sheep, goats, pigs, chickens and ducks for 2010 at a spatial resolution of 0.083333^∘. HERMESv3_BU adjusts the original data to match province-level official records from most recent years. These official statistics, which need to be provided by the user, are also used to distribute each general livestock GLWv3 group (e.g. pigs) into specific categories (e.g. fattening pigs between 50 and 80 kg, boars, sows not yet covered). This disaggregation is relevant due to the different types of feeding received by each animal type, which subsequently affect the levels of N and TAN content in their excreta. The remapping and adjustment of the GLWv3 original data to the destination gridded domain is performed using the spatial operation tools described in Sect. 2.4. HERMESv3_BU currently considers a total of 36 livestock categories, which are grouped into five main groups: pigs (10), cattle (11), poultry (2), goats (6) and sheep (7). Other livestock groups (i.e. buffaloes, horses and ducks) are not currently included due to their low contribution to NH₃ emissions in Europe.

Daily factors ($\mathrm{FD}(x,d{)}_{a,l}$) are estimated following the dynamical emission parameterization reported in Gyldenkærne et al. (2005) and Skjøth et al. (2004, 2011). For housing operations, the daily factors are assumed to be dependent on the barn air temperature and ventilation rate, following Eq. (13):

$$\begin{array}{}\text{(13)}& \begin{array}{rl}\mathrm{FD}(x,d{)}_{\mathrm{housing},l}& =\left(T\right(x,d{)}_{\mathrm{housing},l}^{\mathrm{0.89}}\\ & +V(x,d{)}_{\mathrm{housing},l}^{\mathrm{0.26}}),\end{array}\end{array}$$

where $T(x,d{)}_{\mathrm{housing},l}$ is the barn temperature associated with the housing of livestock category l at destination grid cell x and day d (^∘C) and $V(x,d{)}_{\mathrm{housing},l}$ is the ventilation rate associated with the housing of livestock category l at destination grid cell x and day d (m s⁻¹). Both parameters are calculated as a function of the outdoor 2 m temperature and 10 m wind speed considering the parameterizations reported in Gyldenkærne et al. (2005), which take into account whether the livestock are kept in open or closed barns. HERMESv3_BU assumes that pigs and poultry are kept in closed barns, while cattle, sheep and goats are kept in open barns, following Backes et al. (2016). The daily factors for yarding and storage activities also follow Eq. (13), but using wind speed and air temperature. Finally, for grazing activities the temporal variability is linked to the availability of grass and to its growing period. Therefore, the daily distribution is estimated using Eq. (11) and considering only the growing stage of grass. Concerning the hourly distribution of emissions, the same profile proposed for fertilizer application is proposed.

Emissions from NMVOC, PM₁₀ and PM_2.5 are estimated by multiplying the animal density (D(x)_l) by the default Tier 1 annual emission factors reported per livestock category in EMEP/EEA (2016) (chap. 3.B, Tables 3.4 and 3.5). As noted in the guidelines, only housing emissions are considered due to the large uncertainties and lack of available information from other sources. A flat temporal distribution is assumed for both pollutants since emissions are related to feeding processes.

Figure 6c shows the results of the Spanish annual NH₃ livestock emissions calculated on a Lambert conformal conic grid of 4 km×4 km resolution. Animal number statistics at NUTS level 3 were obtained from MAPA (2017c), and information on the N excretion rate for each livestock category was derived from MAPA (2017b). The TAN content data were obtained from Antezana et al. (2016) for pig categories and EMEP/EEA (2016) (chap. 3.B, Table 3.9) for the remaining animals. The time series of daily NH₃ emissions for the region of Murcia is plotted in Fig. 6d. Meteorology is derived from ERA5. As shown in the example, HERMESv3_BU is able to discriminate the emission estimation by livestock group, which allows quantifying the contribution of each source to the total NH₃. The annual emissions estimated for this Spanish region, also considered a major European NH₃ hotspot, are again compared against the emission fluxes derived from IASI (Van Damme et al., 2018) (Fig. 6e). The HERMESv3_BU results include both livestock and fertilizers emissions. It is observed that results are in the same order of magnitude, HERMESv3_BU reporting 21 % less emissions.

3.3.3 Crop operations

Particulate matter emissions released during soil cultivation and crop harvesting activities are estimated following Eq. (14):

$$\begin{array}{}\text{(14)}& \begin{array}{rl}& E(x,h)=\sum _{o=\mathrm{1}}^m{\sum }_{c=\mathrm{1}}^n{A}_c\left(x\right)\cdot {\mathrm{EF}}_{c,o}\\ & \cdot \mathrm{FM}(m{)}_{c,o}\cdot \mathrm{FW}(d)\cdot \mathrm{FH}(h),\end{array}\end{array}$$

where E(x,h) is the hourly PM₁₀orPM_2.5 emissions at destination grid cell x and hour h (g h⁻¹); A_c(x) is the annual cultivated area of crop c at destination grid cell x (ha yr⁻¹); EF_c,o is the emission factor for crop c and operation o (g PM₁₀orPM_2.5 ha⁻¹); FM(m)_c,o is the monthly factor associated with crop c, operation o and month m (0 to 1); FW(d) is the weekly factor associated with day d (0 to 1); and FH(h) is the hourly factor associated with hour h (0 to 1). The number of crop categories is n.

The model uses the emission factors reported by the EMEP/EEA (2016) Tier 2 methodology (chap. 3.D, Tables 3.6 and 3.8). The methodology takes into account emissions happening in four different types of crops (i.e. wheat, rye, barley and oat). Gridded crop areas are estimated using the same approach described in the fertilizers sector (see Sect. 3.3.1). Considering the monthly distribution of emissions, specific weight factors are applied based on the soil cultivation and harvesting calendars reported by Sacks et al. (2010). The daily and hourly profiles used by default are the ones recommended by EMEP/EEA (2016) for temporally allocating agricultural machinery activities (chap. 1.A.4.c ii, Table 5.1).

3.4 Residential and commercial combustion

Emissions from residential and commercial small combustion plants are estimated as follows (Eq. 15):

$$\begin{array}{}\text{(15)}& \begin{array}{rl}& E_i(x,h)={\sum }_{f=\mathrm{1}}^n{\mathrm{FC}}_f\left(x\right)\cdot {\mathrm{EF}}_{f,i}\\ & \cdot \mathrm{FD}(x,d{)}_f\cdot \mathrm{FH}(h{)}_f,\end{array}\end{array}$$

where E_i(x,h) is the hourly emissions of pollutant i and hour h (g h⁻¹); FC_f(x) is the annual gridded fuel consumption data for fuel category f and destination grid cell x (GJ yr⁻¹); EF_f,i is the emission factor linked to the consumption of fuel f and pollutant i (g GJ⁻¹); FD(d,x)_f is the gridded daily factor associated with fuel type f, destination grid cell x and day d (0 to 1); and FH(h)_f is the hourly factor associated with fuel type f and hour h (0 to 1). The number of fuel categories is n.

The gridded consumption data are calculated combining the fuel statistic consumptions at NUTS level 3 provided by the user with the GHSL population data, as previously described in Sect. 2.4. HERMESv3_BU considers the following fuel types: natural gas, LPG, heating diesel oil, wood and coal. For all of them, emission factors are obtained from the EMEP/EEA (2016) Tier 2 approach (chap. 1.A.4, Tables 3-10, 3-18, 3-19, 3-21, 3-31 and 3-37). In the particular case of wood, the proposed emission factors are obtained as an average of the values reported for the different types of appliances (i.e. fireplace, conventional stove, conventional boiler, eco-labelled boiler and pellet stove; Tables 3-14, 3-17, 3-18, 3-24 and 3-25). The average emission factors are obtained considering the appliance shares reported by Denier van der Gon et al. (2015). It is also important to note that wood-related PM₁₀ and PM_2.5 emission factors take into account the condensable fraction of PM.

The temporal distribution of annual emissions is performed using gridded daily temporal profiles, which are derived according to the heating degree day (HDD) concept. The HDD is an indicator used as a proxy variable to reflect the daily energy demand for heating a building (Quayle and Diaz, 1980). The original expression is reformulated according to Mues et al. (2014) to consider those combustion processes that are not only related to space heating but also to other activities that remain constant throughout the year (i.e. water heating, cooking) (Eq. 16):

$$\begin{array}{}\text{(16)}& {\mathrm{FD}}_f(x,d)={\displaystyle \frac{\mathrm{HDD}(x,d)+f_f\cdot \stackrel{\mathrm{‾}}{\mathrm{HDD}}\left(x\right)}{(\mathrm{1}+f_f)\cdot \stackrel{\mathrm{‾}}{\mathrm{HDD}}\left(x\right)}},\end{array}$$

where HDD(x,d) is the heating degree-day factor for grid cell x and day of the year d (0 to 1); $\stackrel{\mathrm{‾}}{\mathrm{HDD}}\left(x\right)$ is the yearly average of the heating degree-day factor per grid cell x; and f_f is a constant offset that indicates the share of the fuel f that is used for activities not related to space heating (0 to 1). By default, f_f is considered to be 0 for wood and 0.2 for the other fuels, following the European household energy statistics reported by Eurostat (2018). HDD(x,d) and $\stackrel{\mathrm{‾}}{\mathrm{HDD}}\left(x\right)$ are estimated as shown in Eqs. (17) and (18) (Quayle and Diaz, 1980):

$$\begin{array}{}\text{(17)}& {\displaystyle }& {\displaystyle }\mathrm{HDD}(x,d)=\max (T_{\mathrm{b}}-T_{\mathrm{2}\mathrm{m}}(x,d),\mathrm{1})\text{(18)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{\mathrm{HDD}}\left(x\right)={\displaystyle \frac{{\sum }_{\mathrm{1}}^N\mathrm{HDD}(x,d)}{N}},\end{array}$$

where T_b is the threshold temperature above which a building needs no heating (i.e. heating appliances will be switched off) (^∘C), T_2 m(x,d) is the daily mean 2 m outdoor temperature for grid cell x and day d (^∘C), and N is the number of days of the simulation year (365 or 366). Following Spinoni et al. (2015), who developed gridded European degree-day climatologies, we assume T_b=15.5 ^∘C, a value also suggested by the UK Met Office. The HDD(x,d) value increases with increasing difference between the outdoor and base temperatures. Note that a minimum value of 1 is assumed instead of 0 to avoid numerical problems.

Two profiles are proposed in HERMESv3_BU for the hourly distribution of emissions. The first one applies only to residential wood burning emissions, and it is a combination of existing profiles derived from citizen interviews performed in Norway and Finland (Finstad et al., 2004; Gröndahl et al., 2010) as well as from long-term measurements of the wood burning fraction of black carbon in Athens (Athanasopoulou et al., 2017). The second profile (applicable to all other fuel types) is equivalent to the one proposed by Denier van der Gon et al. (2011) for residential combustion emissions. The wood-related profile presents an intense peak in the evening hours but not during the morning (Grythe et al., 2019). This fact is related to a common practice in Europe of using fireplaces and other types of wood-burning appliances mainly in the evening.

3.5 Other mobile sources

3.5.1 Shipping in port areas

Hourly emissions related to fuel combustion processes occurring in main and auxiliary maritime engines during manoeuvring and hotelling operations in port areas are estimated as follows (Eq. 19):

$$\begin{array}{}\text{(19)}& \begin{array}{rl}& E_{p,f,i}(x,h)={\sum }_{v=\mathrm{1}}^n\sum _{e=\mathrm{1}}^{\mathrm{2}}{N}_{p,v}\cdot S(x{)}_{p,f}\cdot t_{p,v,f}\cdot P_{v,e}\\ & \cdot {\mathrm{LF}}_{v,e,f}\cdot {\mathrm{EF}}_{v,e,f,i}\cdot \mathrm{FM}(m{)}_{v,p}\cdot \mathrm{FW}(d{)}_p\cdot \mathrm{FH}(h{)}_{v,p},\end{array}\end{array}$$

where $E_{p,f,i}(x,h)$ is the hourly emissions of pollutant i at port p, destination grid cell x and hour h during phase f (i.e. manoeuvring, hotelling) (g h⁻¹); $N_{p,v,f}\left(x\right)$ is the number of annual operations associated with vessel v (i.e. liquid bulk ship, dry bulk carrier, general cargo, ro-ro, cruiser, ferry, container, tug or others) at port p (operation yr⁻¹); S(x)_p,f is the spatial weight factor associated with port p and destination grid cell x during phase f (0 to 1); $t_{p,v,f}$ is the time spent by vessel v to complete phase f at port p (h); P_v,e is the average power of vessel v's engine e (1: main; 2: auxiliary) (kW); ${\mathrm{LF}}_{v,e,f}$ is the average load factor of vessel v's engine e of during phase f (0 to 1); ${\mathrm{EF}}_{v,e,f,i}$ is the emission factor for pollutant i of vessel v's engine e during phase f (g kWh⁻¹); FM(m)_v,p is the monthly factor associated with month m, vessel v and port p (0 to 1); FW(d)_p is the weekly factor associated with day of the week d and port p (0 to 1); and FH(h)_p is the hourly factor associated with hour h and port p (0 to 1). The number of vessel categories is n.

The estimation methodology and corresponding emission factors are derived from the EMEP/EEA (2016) Tier 3 approach (chap. 1.A.3.d, Table 3-10). Information on the vessel's technical characteristics are obtained from Trozzi (2010), including (i) the average power of a vessel's engines (as a function of the gross tonnage) and (ii) the engine type (i.e. slow-speed diesel, medium-speed diesel, high-speed diesel, gas turbine, steam turbine) and fuel class (i.e. marine diesel oil, marine gas oil and bunker fuel oil) assigned to each vessel category. The values of engine load factors and operation times per phase and vessel are obtained from Entec (2002). The remaining information, which is port-dependent, needs to be provided by the user.

The gridded weight factors (S(x)_p,f), which are used to spatially allocate the emissions, are calculated by the model through a spatial intersection between the destination domain and two shapefiles representing the areas of manoeuvring and hotelling of each port. Both shapefiles need to be provided by the user. Figure 7a shows an example of the total NO_x annual emissions (t yr⁻¹) estimated for the Port of Barcelona on a 1 km×1 km regular grid. Both the manoeuvring and hotelling shapefiles used for the spatial distribution were digitalized using as a basis an infrastructure map provided by the Barcelona's Port Authority (APB, personal communication, 2017). The hotelling layer (in red) consists of a multipolygon shapefile, each polygon representing one of the docks in which the ships operate. A specific weight was assigned to each dock as a function of its usage, allowing a more realistic distribution of the total emissions. On the other hand, manoeuvring emissions were distributed homogeneously between all the cells intersected by the corresponding shapefile (light blue layer). Annual emissions are compared against the results reported by the APB 2013 inventory (APB, 2016) (Fig. 7b). Results show that HERMESv3_BU estimates lower emissions both for NO_x and PM₁₀. The difference can be related to the different years of reference and subsequently the number of vessel operations considered, as well as to the fact that the APB inventory assumes a main engine load factor of 20 % during hotelling operations, while in HERMESv3_BU the load factor is null following the recommendations in Entec (2002).

Figure 7 (a) Annual NO_x emissions (t yr⁻¹) for the Port of Barcelona calculated on a regular 1 km by 1 km grid. Dark red and light blue layers represent the spatial proxies used to allocate hotelling and manoeuvring emissions, respectively. (b) Comparison between annual NO_x and PM₁₀ of Barcelona's port emissions (t yr⁻¹) estimated with HERMESv3_BU and reported by the Barcelona Port Authority inventory (APB, 2016). (c) Spatial distribution of Spanish marinas (with the associated number of docks) and recreational boats along the Costa Brava region (coast of Catalonia). (d) Comparison between annual NO_x, CO, NMVOC and PM_2.5 pleasure boat emissions (t yr⁻¹) estimated for Spain (HERMESv3_BU) and the Baltic Sea region (Johansson et al., 2020). (Sources of the ArcGIS Online basemap: Esri, DigitalGlobe, GeoEye, i-cubed, USDA FSA, USGS, AEX, Getmapping, Aerogrid, IGN, IGP, swisstopo and the GIS user community. Maps in (a) and (b) were created using ArcGIS^® software by Esri. ArcGIS^® and ArcMap™ are the intellectual property of Esri and are used herein under license. Copyright © Esri. All rights reserved. For more information about Esri^® software, please visit https://www.esri.com (last access: February 2020).)

3.5.2 Recreational boats

Emissions derived from pleasure boat activities can be estimated in HERMESv3_BU following Eq. (20):

$$\begin{array}{}\text{(20)}& \begin{array}{rl}& E_i(x,h)={\sum }_{b=\mathrm{1}}^n{N}_b\left(x\right)\cdot T_{\mathrm{b}}\cdot P_b\cdot {\mathrm{LF}}_b\cdot {\mathrm{EF}}_{b,i}\\ & \cdot \mathrm{FM}\left(m\right)\cdot \mathrm{FW}\left(d\right)\cdot \mathrm{FH}\left(h\right),\end{array}\end{array}$$

where E_i(x,h) is the hourly recreational boat emissions of pollutant i at destination grid cell x and hour h (g h⁻¹); N_b(x) is the number of pleasure boats associated with category b and destination grid cell x (no. boats); T_b is the number of hours that the pleasure boat of category b is used during 1 year (h); P_b is the engine nominal power associated with the pleasure boat of category b (kW); LF_b is the engine load factor associated with the pleasure boat of category b (0 to 1); EF_b,i is the emission factor linked to the pleasure boat of category b and pollutant i (g kWh⁻¹); FM(m) is the monthly factor associated with month m (0 to 1); FW(d) is the weekly factor associated with day of the week d (0 to 1); and FH(h) is the hourly factor associated with hour h (0 to 1). The number of recreational boat categories is n.

The parameters EF_b,i, LF_b, T_b and P_b are derived from the EMEP/EEA (2016) Tier 3 approach (chap. 1.A.5.b, Table 3-11). The gridded number of pleasure boats (N_b(x)) is computed by combining official statistics of registered recreational crafts with a raster file that simulates the distribution of recreational boat activities along the coast. Both datasets need to be provided by the user. Figure 7c shows the spatial distribution of Spanish recreational boats along the Costa Brava region (coast of Catalonia). The spatial proxy is based on the location of each Spanish marina and associated number of docks, which are also represented in the map (Fondear, 2019). A raster interpolation was performed to simulate the activities of recreational boats nearby the marinas, considering the number of docks as a weight and assuming that no operations are happening beyond the territorial waters (i.e. 12 nautical miles away from the coastline). A hot spot region is observed near Cap de Creus (headland located at the far northeast of Catalonia) due to the presence of the largest Spanish marina (Empuriabrava, with 5000 docks). The total number of recreational boats per type of boat were derived from ICOMINA (2016). Figure 7d shows the total annual emissions estimated with HERMESv3_BU for Spain. To the authors' knowledge, no other national emission inventory is currently available for this pollutant source. Emissions are contrasted against an inventory of pleasure boats estimated for the Baltic Sea (Johansson et al., 2020) to assess that the results are within the same range of magnitude. HERMESv3_BU reports higher NO_x (4.1 times) and lower CO and NMVOC (0.7 and 0.5 times, respectively) emissions. This is probably due to the different fleet characteristics of each region. While in Spain more than 40 % of the boats are related to large diesel motor sail boats, in the Baltic Sea region this category only accounts for less than 15 % of the total fleet.

3.5.3 Agricultural machinery

Emissions related to the use of agricultural equipment (i.e. two-wheel tractors, agricultural tractors and harvesters) are estimated following Eq. (21):

$$\begin{array}{}\text{(21)}& \begin{array}{rl}& E_i(x,h)={\sum }_{e=\mathrm{1}}^n{N}_e\left(x\right)\cdot T_e\cdot P_e\cdot {\mathrm{LF}}_e\cdot {\mathrm{EF}}_{e,i}\left(P_e\right)\\ & \cdot (\mathrm{1}+{\mathrm{DF}}_{e,i})\cdot {\mathrm{FM}}_e\left(m\right)\cdot \mathrm{FW}\left(d\right)\cdot \mathrm{FH}\left(h\right),\end{array}\end{array}$$

where E_i(x,h) is the hourly emissions of pollutant i at destination grid cell x and hour h (g h⁻¹); N_e(x) is the number of agricultural equipment associated with category e and destination grid cell x (no. equipment); T_e is the number of hours that the equipment of category e is used during 1 year (h); P_e is the engine nominal power associated with the equipment of category e (kW); LF_e is the engine load factor associated with the equipment of category e (0 to 1); EF_e,i(P_e) is the emission factor linked to the agricultural equipment of category e and pollutant i (g kWh⁻¹) as a function of the engine nominal power P_e; DF_e,i is the deterioration adjustment factor for the equipment of category e and pollutant i; FM_e(m) is the monthly factor associated with month m and equipment category e (0 to 1); FW(d) is the weekly factor associated with day d (0 to 1); and FH(h) is the hourly factor associated with hour h (0 to 1). The number of agricultural equipment categories is n.

HERMESv3_BU takes as an input the total number of agricultural equipment and corresponding engine nominal power registered at the NUTS 3 level. The spatial allocation of this data onto the destination domain is performed considering the CLC non-irrigated arable land category as a proxy and applying the spatial operations described in Sect. 2.4. Both the emission and deterioration adjustment factors are derived from the EMEP/EEA (2016) Tier 3 methodology (chap. 1.A.4, Tables 3-6 and 3-11). The load factor adjustments proposed by default are the ones reported by Winther and Nielsen (2006) (i.e. 0.4 for two-wheel tractors, 0.5 for agricultural tractors and 0.8 for harvesters). Emissions from two-wheel tractors and agricultural tractors are disaggregated at the monthly level considering the crop calendars associated with the cultivation operation. In the case of the harvester emissions, the monthly factors are associated with the harvesting period of the non-irrigated arable crops (Sacks et al., 2010).

3.5.4 Landing and take-off cycles at airports

Hourly aircraft landing and take-off (LTO) emissions occurring at airports are estimated according to Eq. (22):

$$\begin{array}{}\text{(22)}& \begin{array}{rl}& E_{a,f,i}(x,h)={\sum }_{p=\mathrm{1}}^{n}N(m{)}_{a,p,f}\cdot S(x{)}_{a,f}\cdot {\mathrm{EF}}_{p,f,i}\\ & \cdot \mathrm{FW}(d{)}_a\cdot \mathrm{FH}(h{)}_{a,f},\end{array}\end{array}$$

where $E_{a,f,i}(x,h)$ is the hourly emissions of pollutant i at airport a, destination grid cell x and hour h during phase f (i.e. approach, landing, taxi-in, post taxi-in, pre taxi-out, taxi-out, take-off and climb-out) of the LTO cycle (g h⁻¹); $N(m{)}_{a,p,f}$ is the number of monthly operations associated with aircraft p at airport a during phase f for month m (operations per month); S(x)_a,f is the spatial weight factor associated with airport a and destination grid cell x during phase f (0 to 1); ${\mathrm{EF}}_{p,f,i}$ is the emission factor for pollutant i associated with aircraft p and phase f (grams per operation); FW(d)_a is the weekly factor associated with day d and airport a (0 to 1); and FH(h)_a,f is the hourly factor associated with hour h phase f and airport a (0 to 1).

Depending on the LTO phase, different emission processes and subsequently emission factors are considered. For taxi-in, taxi-out, take-off, climb out and approach operations the emission factors correspond to the fuel combustion in the main engines (${\mathrm{EFmain}}_{p,f,i}$) (Eq. 23):

$$\begin{array}{}\text{(23)}& {\mathrm{EFmain}}_{p,f,i}=E_p\cdot t_{a,p,f}\cdot {\mathrm{EFengine}}_{p,f,i},\end{array}$$

where E_p is the number of engines associated with aircraft p (engine); $t_{a,p,f}$ is the time spent by aircraft p to complete phase f at airport a (s); and ${\mathrm{EFengine}}_{p,f,i}$ is the emission factor for pollutant i of the main engine associated with aircraft p and phase f (grams per second per engine per operation). The number of engines and emission factors associated with each aircraft category are derived from the EMEP/EEA (2016) Tier 3 approach (Annex 5 spreadsheets to the chap. 1.A.3.a). For taxi-in and taxi-out, times are also obtained from the same source, whereas for the other operations (take-off, climb out and approach) different values are assumed as a function of the type of aircraft (i.e. wide-body planes, narrow-body planes, business planes and light planes with piston engines) (Dellaert and Hulskotte, 2017).

For landing operations, particulate-matter emission factors (EFwear_p,i) related to the wear of the aircraft brakes and tyres are considered (Eq. 24) (Morris, 2006):

$$\begin{array}{}\text{(24)}& {\mathrm{EFwear}}_{p,i}={\mathrm{MTOW}}_p\cdot {\mathrm{EFwear}}_i,\end{array}$$

where MTOW_p is the maximum take-off weight associated with aircraft p (t) and EFwear_i is the emission factor for pollutant i (grams per tonne per operation), which is taken from Morris (2006).

Finally, emissions during the pre taxi-out and post taxi-in operations are linked to the fuel combustion in the auxiliary power units (APUs). The corresponding emission factors (EFaux_p,i) are estimated as follows (Eq. 25) (Watterson et al., 2004):

$$\begin{array}{}\text{(25)}& {\mathrm{EFaux}}_{p,i}=t\mathit{\_}{\mathrm{apu}}_{a,p,f}\cdot {\mathrm{EFapu}}_{p,i,f},\end{array}$$

where $t\mathit{\_}{\mathrm{apu}}_{a,p,f}$ is the APU running time associated with aircraft p, phase f and airport a (s) and EFapu_p,i is the emission factor for pollutant i of the APU engine associated with aircraft p and phase f (grams per second per operation). Data on the type of APU installed in each type of aircraft and corresponding emission factors are obtained from Watterson et al. (2004).

The spatial weight factors (S(x)_a,f) to spatially distribute the estimated emissions also vary according to the LTO cycle's phase. Taxi-in, post taxi-in, pre taxi-out and taxi-out emissions are assigned at the ground level; taking into account digitized airport areas. Emissions are mapped to the grid cells that spatially intersect with the airport polygons, which need to be provided by the user in a shapefile. Take-off and landing emissions are also allocated at the ground level, but they are horizontally distributed across the grid cells that intersect with the digitized runways (also provided by the user in a shapefile). The user has the option of defining specific weights to each runway as a function of its usage (for departures, arrivals or both). Finally, emissions from approach and climb-out operations are allocated on a 3D basis, considering the trajectory outlined between the origin or end of the runway and 1000 m of altitude. This operation is performed in HERMESv3_BU using as input a shapefile representing the air trajectories and information on the approach and take-off angles for each runway.

3.6 Speciation

This process disaggregates the calculated primary pollutants (i.e. NO_x, NMVOC, PM_2.5) into the more detailed species defined by the chemical mechanism of interest. Specific speciation CSV files, which contain a set of profiles with numerical factors for converting the primary pollutants (e.g. NO_x) to output model species (e.g. NO, NO₂ and HONO), are assigned to each pollutant sector. The number of speciation profiles considered varies according to the pollutant sector of interest. For instance, in the case of residential combustion the number of speciation profiles proposed is equal to the number of fuel types considered (e.g. natural gas, biomass), whereas in the case of road transport, specific profiles are assigned to each vehicle category. The assignment between profiles and pollutant categories is performed following the cross-referencing system described in Sect. 2.3.

The speciation factors are mass-based (i.e. grams of chemical species per gram of source pollutant) for NO_x and PM_2.5. For NO_x, emissions are divided into 90 % NO and 10 % NO₂ for all sectors (Houyoux et al., 2000) except for road transport, where vehicle-dependent speciation factors are derived from EMEP/EEA (2016) (chap. 1.A.3.b.i–iv, Table 3.87) and the work by Carslaw et al. (2016) and Rappenglueck et al. (2013). For PM_2.5, the speciation factors are derived from multiple sources including EMEP/EEA (2016) and the SPECIATE (Simon et al., 2010) and SPECIEUROPE (Pernigotti et al., 2016) databases. For the rest of the calculated primary pollutants (i.e. CO, NH₃, SO₂, PM₁₀) a default speciation factor of 1 is proposed for all sources. All computed gas-phase species are converted from mass to moles using a molecular weight CSV file that is included in the model database (built-in data file).

In the case of NMVOC, the speciation factors used by HERMESv3_BU are mol-based (mol of chemical species per gram of source pollutant) and are estimated following the expression proposed by Li et al. (2014) (Eq. 26):

$$\begin{array}{}\text{(26)}& {\mathrm{SF}}_{\stackrel{\mathrm{‾}}{e},p}=\sum _{j=\mathrm{1}}^n{\displaystyle \frac{X_{j,p}}{{\mathrm{MW}}_j}}\cdot C_{j,\stackrel{\mathrm{‾}}{e}},\end{array}$$

where X_j,p is the mass fraction of chemical compound j to total NMVOC emissions for speciation profile p, MW_j is the molecular weight of chemical compound j and $C_{j,\stackrel{\mathrm{‾}}{e}}$ is the mole-based conversion factor of chemical compound j to destination chemical species $\stackrel{\mathrm{‾}}{e}$. X_j,p values are obtained from multiple sources including EMEP/EEA (2016) and the NMVOC SPECIATE database (Simon et al., 2010), and MW_j and $C_{j,\stackrel{\mathrm{‾}}{e}}$ are obtained from the mechanism-dependent mapping tables proposed in Carter (2015). The number of individual chemical compounds considered in the total NMVOC is n.

For reference, the HERMESv3_BU test case input database currently includes speciation profiles for the Carbon Bond 05 (CB05) (Whitten et al., 2010) gas-phase mechanism and the fifth-generation aerosol module (AERO5) (Roselle et al., 2008). Following the formats of the associated input files, the user can create its own speciation profiles for other mechanisms of interest, also using alternative sources of information. As an illustration, Table 2 shows some examples of proposed CB05 and AERO5 speciation profiles for different pollutant sectors. As mentioned before, HERMESv3_BU allows using specific profiles for each of the source categories included in the different sectors. This feature enables the user to consider key factors influencing the splitting of primary pollutants into chemical mechanism species such as (i) fuel type (wood, natural gas) for NMVOC emissions in the residential sector (Simon et al., 2010), (ii) vehicle type (passenger car, motorcycle) and vehicle fuel (diesel, gasoline) for road transport NO_x (Carslaw et al., 2016; Rappenglueck et al., 2013), (iii) type of process (combustion, road wear) for road transport PM_2.5 emissions (EMEP/EEA, 2016) (chap. 1.A.3.b.i–iv, Table 3.88 and chap 1.A.3.b.vi, Table 3-11), and (iv) animal type (pigs, cattle) for livestock NMVOC emissions (EMEP/EEA, 2016) (chap. 3.B, Table A1.2).

Table 2 Example of speciation profiles included in HERMESv3_BU for speciating primary emissions to CB05 and AERO5 mechanisms. Different pollutant sectors and categories are shown to illustrate the degree of specificity allowed by the model. Speciation factors are mole-based (mol of chemical species per gram of source pollutant) for NMVOC and mass-based (grams of chemical species per gram of source pollutant) for NO_x, CO, SO_x, NH₃ and PM_2.5. Symbol “–“ denotes that no primary emissions are considered for that pollutant. All CB05 and AERO5 species are defined in Table A2 of Guevara et al. (2019a).

4 Technical implementation

HERMESv3_BU is coded using Python 3.7.X and requires numpy (> D1.16.0), NetCDF4 (> D1.3.1) under HDF5 in parallel mode, pandas (> D0.22.0), geopandas (> D0.4.0), pyproj (> D1.9.5.1), configargparse (> D0.11.0), cf_units (> D1.1.3), timezonefinder (> D2.1.0), mpi4py (> D3.0.0), pytest (> D3.6.1), shapely (D >1.6.4), scipy (D >0.14.1) and rasterio (D >1.0.21) Python libraries.

The emission core of each pollutant sector (i.e. calculation, temporal and spatial distribution and speciation) is parallelized using one of the following strategies:

• road link partition – this strategy is applied to the road transport sector. The original digitalized road transport network is split according to the number of processors to be used. Each processor computes the emissions of the road links subgroup that has been assigned with. The partition of the data is automatically balanced with respect to the number of processors used.

• point source partition – this strategy is used for the point source sector. The approach is equal to the one described above, but in this case the partition is applied to the point source facilities.

• polygon partition – this strategy is applied to the shipping and aviation sectors. Each processor is responsible for calculating the emissions of a subset of the total number of infrastructures (i.e. ports or airports) included in the working domain.

• grid partition – this strategy is applied to all the sectors that are computed at the destination grid cell level once the original input activity data have been mapped onto the working domain (i.e. agricultural sources, residential and commercial combustion, recreational boats and gasoline evaporative emissions). The destination working domain is divided into subgroups of row-major consecutive cells. Each subgroup or chunk contains only those grid cells that have activity data information. The number of divisions is balanced and equal to the number of processors to be used.

In all cases, emissions are computed independently in each processor and for each partition (i.e. group of road links, point sources, polygons or grid cells). As a result of the emission calculation process, each individual computational processor ends up containing the gridded hourly estimated emissions for the corresponding section of the working domain. Once this status is reached for all the processors (i.e. emissions have been estimated for all partitions and pollutant sectors), HERMESv3_BU starts the execution of the writing process. The model can perform this task in serial or in parallel mode, regardless of the processors used for the emission calculation process. If the user configures the model to run with more than one writing processor, HERMESv3_BU will decompose the destination working domain into horizontal sections (as many as processors selected), maintaining each row as indivisible. Each processor will first gather the estimated emissions of the corresponding writing section and then compute the total sum per grid cell, vertical layer and time step. During this process, unit conversions will also take place. Finally, each individual section will write the emission 4D array (hour, vertical layer, longitude, latitude) onto the corresponding memory section of the NetCDF output file. When only one writing processor is used, the same process is performed but without applying the aforementioned domain partition.

The number of processors used to calculate each pollutant sector and to perform the writing function are defined by the user in the general configuration file (see Sect. 2.2). The distribution of the defined computational resources is performed using a sector manager function under the framework of the Message-Passing Interface (MPI) protocol (mpi4py Python library). MPI creates a global communicator and the sector manager splits it into a group of sub-communicators (as many as pollutant sectors involved). This allows HERMESv3_BU to exclusively assign each group of processors to a specific sector and therefore isolate the corresponding calculation processes.

Figure 8 shows an illustrative example of the emission calculation and writing parallel implementations. The example considers four pollutant sectors: a point-type sector (e.g. point sources), a line-type sector (e.g. road transport), a polygon-type sector (e.g. aviation) and a grid-type sector (e.g. fertilizers). A total of 17 processors are defined for the emission calculation process (Fig. 8a), which are split as follows: seven for point sources, three for road transport, three for aviation and four for fertilizers. Each processor is assigned to an individual element (or group of elements) of the corresponding pollutant sector and performs the calculation of the gridded hourly emissions (Fig. 8b). Those cells that are not intersected with any emission source are left without information. Once all the emission calculation processors have finished, the writing process begins, in this case using a total of three processors. The working domain is accordingly split into three horizontal sections, the first two consisting of four completed rows and the last one of two (Fig. 8c). For each grid cell of each writing section, the corresponding processor gathers and sums the emissions estimated for all the different sources (Fig. 8d). The information stored in each writing processor is then written in the corresponding section of the output NetCDF file.

Figure 8 Illustrative example of how HERMESv3_BU parallelizes the calculation and writing processes and distributes the computational resources. (a) A total of 17 processors are selected to estimate emissions from point sources, road transport, aviation and an area source sector (e.g. residential combustion). (b) The calculation processors are distributed as follows: seven for the point sources (one processor per facility), three for road transport (one processor per road link), three for aviation (one processor per polygon) and four for residential combustion (one processor per subgroup of cells equally balanced, i.e. 10 cells per subgroup). Each processor estimates the emissions of the corresponding partitioned element and maps them onto the intersected grid cells (highlighted in grey). (c) During the writing process, HERMESv3_BU decompose the destination working domain into three horizontal sections (equivalent to the number of writing processors selected), maintaining each row as indivisible. (d) For each horizontal section, the model gathers the estimated gridded emissions of all the sources involved in that section and then computes the total sum per grid cell. Once this operation is finished, the resulting emissions are dumped into the corresponding section of the output file.

The supercomputer MareNostrum4, hosted by the BSC (http://www.bsc.es/marenostrum/marenostrum, last access: February 2020), was used to test the capability of HERMESv3_BU to scale up the emission calculation processes. HERMESv3_BU was executed using a number of cores from 1 to 256, doubling the number in each successive test. The proposed testing domain was a Lambert conformal conic domain of 4 km by 4 km with 397 rows, 397 columns and 15 vertical layers covering the Iberian Peninsula (Fig. 6a). Hourly emissions for Spain were estimated for 24 time steps and for each one of the pollutant sectors described in Sect. 3.

As shown in Fig. 9, for all sectors the computational time decreases as the number of cores is increased. Traffic is the pollutant sector that requires the largest amount of time to be completed, the total execution time decreasing from 48 585.5 to 450.3 s when changing from 1 to 256 processors. Similar time reduction rates are observed for other sectors when comparing the results obtained with 1 and 256 processors: livestock (from 17 530.2 to 53.2 s), point sources (from 5107.7 to 75.4 s), residential and commercial combustion (from 3230.7 to 27.0 s), and crop fertilizers (from 2814.7 to 39.1 s). For the rest of the sectors (i.e. agricultural machinery, crop operations, recreational boats, shipping in ports), execution times are already considerably low when only using one processor (between 19.5 and 326.5 s). Note that for aviation and shipping in port areas, the number of processors with which HERMESv3_BU can be executed is limited to the number of infrastructures considered in the domain of study (i.e. 47 airports and 48 ports in this example). Traffic is by far the most time-demanding pollutant sector due to a combination of several facts:

• i.

  a large amount of sources – the sector considers a road network with more than 100 000 road links and includes 491 vehicle categories;

• ii.

  multiple emission processes – the sector estimates emissions derived from hot- and cold-exhaust gases, tyre, brake and road wear, resuspension, and gasoline evaporation;

• iii.

  meteorologically dependent functions – three of these processes are estimated using meteorologically dependent parameterizations (i.e. cold-exhaust and gasoline evaporation take into account outdoor temperature, while resuspension considers precipitation information).

Figure 9 Computational times (s) obtained from the scalability test performed with HERMESv3_BU over Spain. Each line represents the amount of time needed to complete the emission estimation process of each pollutant sector as a function of the number of processors used (x axis, represented with a base 2 logarithmic scale).

The execution of each sector is performed in parallel (not sequentially). Therefore, the total execution time of HERMESv3_BU for the emission calculation process is equal to the sector's largest execution time. Following with the proposed example, if the user wants to perform the emission calculation process in less than 8 min, a total of 298 processors should be used with the following distribution: 256 for traffic, 16 for point sources, 12 for livestock, 4 for residential and commercial combustion, 4 for crop fertilizers and 1 for each one of the other sectors. The computational time needed to gather and write all the estimated emissions into a NetCDF output file is very low compared to the emission estimation time (i.e. between 13 and 25 s depending on the number of writing processors selected, not shown).

5 Conclusions

This paper presents HERMESv3_BU, a stand-alone and open-source emission model that estimates high-resolution spatial and hourly bottom–up anthropogenic emissions for air quality modelling. The tool combines state-of-the-art estimation approaches with local activity and emission factors along with meteorological data.

The main characteristics of HERMESv3_BU are as follows:

• multiple map projections and model grids – the model supports emission calculations in regular and rotated lat–long, Lambert conformal conic and Mercator map projections. Users can freely define their working domain introducing information on the starting x–y coordinates, number of grid cells in each direction and physical size of the grid cells.

• multiple pollutant sources – the model includes bottom–up emission estimation methodologies for several anthropogenic emission sources, including point sources, road transport (hot- and cold-exhaust, wear and resuspension, gasoline evaporation), agriculture (i.e. fertilizer application, livestock, crop operations), residential and commercial combustion, and other mobile sources (i.e. shipping in port areas, recreational boats, agricultural machinery, landing and take-off cycles at airports). Users can estimate emissions for each individual source or a combination thereof.

• subsetting of pollutant categories and regions – for each pollutant sector, users can select individual pollutant categories to be considered during the emission calculation process (e.g. 491 vehicle categories for road transport, 28 crop types for fertilizer application). Similarly, users can define a region of interest within the working domain (i.e. administrative areas or self-defined polygons) for which the emissions should be calculated. These functionalities can become useful when studying the contribution of certain categories or regions to total emissions (e.g. diesel vehicles from a selected province) or when performing source attribution modelling studies.

• emission estimation approaches – HERMESv3_BU integrates estimation methodologies and emission factors that are mostly based on (but not limited to) the approaches reported by the European EMEP/EEA (2016) emission inventory guidebook. The model also includes meteorologically dependent functions to take into account the dynamical component of the emission processes, both in terms of spatial and temporal allocation (e.g. plume rise calculations for point sources, effect of temperature and wind speed on the volatilization of agricultural NH₃).

• applicability – HERMESv3_BU is designed so that it can be applicable to any European country or region where the required input data are available. Global and regional state-of-the-art datasets are considered in the model to decrease the amount of information that needs to be provided by the user, including land use data (CORINE Land Cover inventory), livestock distribution (Gridded Livestock of the World), population maps (Global Human Settlement Layer) and soil properties (World Soil Information database). The application of the tool to non-European regions could be also possible, but it would involve reviewing in much more detail the emission factors databases and estimation methodologies that are proposed by default in the model.

• spatial operations – HERMESv3_BU contains a set of GIS functionalities that helps users manipulate and generate georeferenced data files related to emissions modelling (i.e. Esri shapefiles and Geotiff rasters). The operations included in the model allow developing automatically individual spatial surrogates or remapping spatial data from one gridded domain to another, among other things.

• emission outputs compatible with multiple air quality models – the output NetCDF emission files follow the convention of multiple chemical transport models, including CMAQ, WRF-CHEM and MONARCH. In the particular case of road transport, the resulting link-level emissions can be also used as input for the R-LINE Gaussian dispersion model.

• parallel implementation – the parallelization of the emission calculation and writing processes allows scaling up the execution time of the model, which may be relevant when using the model for air quality forecasting applications.

Several emission outputs obtained with HERMESv3_BU are provided in this paper and compared to other existing emission datasets to illustrate its potential. Future work will focus on performing a full evaluation of the emission results obtained for Spain through comparisons with other inventories (e.g. EDGAR, EMEP) and using air quality modelling.

HERMESv3_BU currently experiences some limitations. The first one is the non-inclusion of a submodule to estimate NMVOC emissions related to the use of solvents. For this source category, activity data are very uncertain and typically not recorded in statistics (e.g. domestic and industrial use of solvent products). For now, and as an alternative, emissions reported by other existing inventories will be used through the application of HERMESv3_GR. The second limitation is related to the suitability of the proposed emission estimation approaches for its application over different regions. For instance, in the case of residential combustion emissions, HERMESv3_BU assumes that wood combustion only occurs in rural areas, which is the case for Spain. Nevertheless, this assumption may be wrong for other countries such as Norway (López-Aparicio et al., 2017). Similarly, for the agriculture sector it is assumed that the temporal distribution of fertilizer application is not affected by any policy restriction (i.e. closed periods). While this is again true for Spain, there are some European countries (e.g. Germany) in which the application of fertilizer is prohibited during certain periods (Backes et al., 2016). In this sense, it is expected that the emission approaches will be refined in future versions of HERMESv3_BU. Regarding greenhouse gases, estimated emissions are currently only related to combustion processes. Emissions from non-combustion sources (e.g. CH₄ emissions from enteric fermentation processes in the livestock sector) are governed by complex processes that need to be described using specific estimation approaches, such as the ones reported in the Intergovernmental Panel on Climate Change (IPCC) inventory guidelines (IPCC, 2006). The implementation of such processes in HERMESv3_BU is a task that we plan to address in future versions of the system. One last aspect to consider is the availability and accessibility of all the required input data needed to run the model, which may limit the usability of the tool in certain regions. It is important to note that some of the required data can be obtained from European homogenized databases such as Eurostat (https://ec.europa.eu/eurostat/, last access: February 2020). Moreover, the continuous growth of the open data movement can also help to overcome this limitation. One example for this is the Open Transport Map (http://opentransportmap.info/, last access: February 2020), an application for accessing daily traffic volumes for the main European road network.

HERMESv3_BU represents an effort to integrate and combine in a flexible and transparent way state-of-the-art methods for estimating high-resolution bottom–up emissions from multiple anthropogenic sources. The purpose of the model is to serve as an emission tool for multiple applications, including air quality research and environmental management.

Appendix A

Table A1 Relationship between crop types and CORINE Land Cover (CLC) land use categories.

Appendix B

Table B1 Classification of HERMESv3_BU input data files per pollutant source.

¹ Fonts are used to specify the file format of each dataset: shapefile, CSV, NetCDF, raster. ² These meteorological parameters are not provided by ERA5 and therefore need
to be derived from other models. ³ ERA5 is proposed since it is open data. Users can alternatively use the outputs from other meteorological models