A new implementation of FLEXPART with Enviro-HIRLAM meteorological input, and a case study during a heavy air pollution event

ABSTRACT

We integrated Enviro-HIRLAM (Environment-High Resolution Limited Area Model) meteorological output into FLEXPART (FLEXible PARTicle dispersion model). A FLEXPART simulation requires meteorological input from a numerical weather prediction (NWP) model. The publicly available version of FLEXPART can utilize either ECMWF (European Centre for Medium-range Weather Forecasts) Integrated Forecast System (IFS) forecast or reanalysis NWP data, or NCEP (U.S. National Center for Environmental Prediction) Global Forecast System (GFS) forecast or reanalysis NWP data. The primary benefits of using Enviro-HIRLAM are that it runs at a higher resolution and accounts for aerosol effects in meteorological fields. We compared backward trajectories generated with FLEXPART using Enviro-HIRLAM (both with and without aerosol effects) to trajectories generated using NCEP GFS and ECMWF IFS meteorological inputs, for a case study of a heavy haze event which occurred in Beijing, China in November 2018. We found that results from FLEXPART were considerably different when using different meteorological inputs. When aerosol effects were included in the NWP, there was a small but noticeable difference in calculated trajectories. Moreover, when looking at potential emission sensitivity instead of simply expressing trajectories as lines, additional information, which may have been missed when looking only at trajectories as lines, can be inferred.

KEYWORDS:

• Atmospheric and chemical transport modelling
• trajectory and particle dispersion modelling
• severe air pollution episode
• FLEXPART
• Enviro-HIRLAM

1. Introduction

FLEXPART (FLEXible PARTicle dispersion model) is a Lagrangian model that is used for simulating atmospheric trajectories, transport, and dispersion of particles and air parcels (Pisso et al., 2019). FLEXPART can be run either in forward mode, which generates atmospheric trajectories and tracks atmospheric transport and dispersion after a release, or in backward mode, which traces an air parcel or set of particles arriving at a point of interest back to its source. In addition to drawing a trajectory as a single line, FLEXPART can calculate potential emission sensitivity, also called source-receptor relationship (Pisso et al., 2019; Seibert & Frank, 2004). In backward runs, the potential emission sensitivity provides a footprint of emission source areas.

FLEXPART has a multitude of applications (Pisso et al., 2019). Examples include air quality analysis (e.g. Madala et al., 2016; Wang et al., 2023; Zhu et al., 2020), tracking wildfire smoke (e.g. Paris et al., 2009), atmospheric chemistry studies (e.g. Xavier et al., 2022), risk assessment for the dispersion of a release of hazardous material such as chemicals or radionuclides (e.g. ElShafeey et al., 2023), and tracing a measured compound back to its probable source (e.g. Sun & Wang, 2014). These applications are invaluable and important to assessments of impact on human health and environment, safety, policy- and decision-making, and scientific research.

FLEXPART requires meteorological input from a numerical weather prediction (NWP) model. The publicly available version of FLEXPART (available at http://flexpart.eu) supports NWP input from the European Centre for Medium-Range Weather Forecasts (ECMWF) Integrated Forecast System (IFS), or from the U.S. National Center for Environmental Prediction (NCEP) Global Forecast System (GFS). Either forecast or reanalysis data from the GFS or IFS models can be used, specifically the Final Operational Model Global Tropospheric Analyses (FNL) from the GFS model, and the ECMWF Re-Analysis Interim (ERA-Interim) or the newer Fifth Generation ECMWF Re-Analysis (ERA5, also abbreviated EA5) from the IFS model.

The GFS and IFS are both global weather forecast models and have many similarities, but being developed by different organizations, they also have differences. Numerous studies have been published that compare the two models. In the context of air quality, Zuo et al. (2023) provides a comprehensive comparison of ERA5 versus FNL reanalysis datasets for satellite-based retrieval of particulate matter pollution in China. The Zuo et al. (2023) study evaluated ERA5 and FNL, along with two other reanalysis datasets (the National Aeronautics and Space Administration Modern-Era Retrospective Analysis for Research and Applications, version 2 (MERRA2), the GEOS-FP Atmospheric Data Assimilation System) with focus on selecting the best data for specific regional applications. The goal was to provide a standard for data selection when choosing between different reanalysis datasets. Results from the study showed that in the context of air quality in China, ERA5 was the most accurate for many of the primary meteorological parameters, and of the four datasets studied, FNL had the greatest uncertainty. Thus, while complementing each other, it is useful to look at results from both input models rather than drawing too many conclusions based on a FLEXPART trajectory created using only one model.

A table listing the historical versions of FLEXPART and a description of changes can be found on the Roadmap section of the FLEXPART website: https://www.flexpart.eu/wiki/FpRoadmap.

In addition to the standard model version that can read ECMWF and NCEP meteorological data, external developers have created modified versions of FLEXPART, which can utilize meteorological input from other NWP models, such as the Weather Research and Forecast model (WRF; Brioude et al., 2013), the Consortium for Small-scale Modeling (COSMO; Katharopoulos et al., 2022), the Applications of Research to Operations at Mesoscale model (AROME; Verreyken et al., 2019), the Fifth-Generation Penn State/National Center for Atmospheric Research Mesoscale Model (MM5; Grell et al., 1994), and the Norwegian Earth System Model (NorESM; Cassiani et al., 2016). Links to these external variants of FLEXPART can be found at: https://www.flexpart.eu/wiki/FpOthermetinput.

In this project, we have developed a method and added new functionality to FLEXPART so that it can read and process meteorology from an additional NWP model, the Environment-HIgh Resolution Limited Area Model (Enviro-HIRLAM). Enviro-HIRLAM is a seamless, online integrated NWP and atmospheric chemical transport modelling system, which simultaneously simulates meteorology and atmospheric composition. It is based on the hydrostatic NWP model known as HIRLAM (Undén et al., 2002), which has been used for research and operational forecasts since 1990. Beginning in the 2000s, the chemistry and aerosol dynamics, including feedback of aerosols and chemical composition on meteorology, have been integrated into the HIRLAM NWP code, and the result became the Enviro-HIRLAM modelling system. Enviro-HIRLAM was originally developed by the Danish Meteorological Institute in collaboration with the HIRLAM consortium and universities/research organizations (Baklanov, 2008; Baklanov et al., 2017). Since spring 2017, the Enviro-HIRLAM research, development, application, and science education activities have been led by the University of Helsinki. Figure 1 shows a workflow of the new functionality in this project.

Figure 1. Workflow of enviro-HIRLAM+FLEXPART, IFS-ERA5+FLEXPART, and GFS-FNL+FLEXPART simulations we generated and compared in this study, including three models (enviro-HIRLAM, IFS and GFS), types of model runs, resolutions, output, steps of preprocessing, trajectory calculations and analysis, as well as types of meteorological input to FLEXPART. Inputs to enviro-HIRLAM, including initial and boundary conditions, and the data assimilation of observations are described in detail in section 2.2.3. Remark: ICs = initial conditions, BCs = boundary conditions, CAMS = Copernicus Atmosphere Monitoring Service, DA = data assimilation, EIs = emission inventories; EH0-EH3 = enviro-HIRLAM model meteorology. Colored arrows (green, cyan, blue, red) indicate enviro-HIRLAM meteorology as input for FLEXPART model corresponding nests (mother domain, nest 1, nest 2, nest 3).

There are two novelties of using Enviro-HIRLAM meteorological fields over IFS or GFS. The first is that Enviro-HIRLAM accounts for aerosol direct and indirect effects and atmospheric composition when simulating meteorology. This is particularly useful for air quality studies. The second novelty is that it uses a downscaling chain that can be run at much higher resolution than GFS and IFS, and users can set up and run their own custom-defined high-resolution domains over a region of interest.

Air quality in the megacity of Beijing, China is an area of research focus, especially with the development of a new comprehensive measurement station (Liu et al., 2020). There are several recent studies that use FLEXPART for air quality studies in the Beijing region, for instance, Wang et al. (2023), Hakala et al. (2022), Guo et al. (2020), Bei et al. (2020), and He et al. (2020). The latter two of these studies used WRF meteorology as input for FLEXPART, with the WRF-FLEXPART model described in Brioude et al. (2013), the most recent version being version 3.3.2, released in 2017. The introduction of Enviro-HIRLAM meteorology for the newest version of FLEXPART can complement these air quality studies with its detailed and high-resolution meteorology as input for trajectory calculations.

Our goal is to be able to use the new combined Enviro-HIRLAM + FLEXPART modelling system in air quality studies, and we will trace observed pollutants back to their source regions and simulate chemical processes which occur in airmasses along the trajectories, eventually arriving at and being observed by the measurement station. Thus, it is important for us to accurately calculate the atmospheric backward trajectories.

Effects of aerosols are especially important during severe, widespread haze events, where heavy aerosols directly and indirectly influence regional and urban meteorology (Baklanov et al., 2016). For example, the North China Plain (NCP), including the Beijing region, is particularly affected by patterns of severe haze episodes in autumn and winter (Zhao et al., 2019; Zheng et al., 2016; Zheng et al., 2015). Aerosols from regional haze can create a feedback effect by cooling air underneath the aerosols and lowering the boundary layer height, further enhancing the severity of the air pollution episode (Zhao et al., 2019). Therefore, when choosing an NWP model as input to FLEXPART for an application involving air quality studies, it may be useful to consider accounting for aerosol effects. In particular, when using trajectory analysis to trace emissions, accounting for aerosols in the meteorology used to calculate the trajectory can make a difference in tracing and identifying sources of air parcels, and thus affecting inferences one might make about air pollution sources and precursors.

Additionally, Enviro-HIRLAM can be run using a downscaling chain, using multiple nested domains with higher resolutions than what GFS or IFS models offer. Downscaling, also called nesting in the context of FLEXPART input, is defined as using a high-resolution model run to dynamically interpolate across a coarser resolution run of the model. For Enviro-HIRLAM, the high-resolution run is done on a smaller geographical extent than the coarse-resolution run, and the output generated from the coarse-resolution run is used as boundary conditions for the high-resolution run. The purpose of nesting is to be able to simulate or resolve small-scale features that depend on large-scale effects, without the computational expense of running the model at high resolution over a large geographical area. For example, local-scale meteorology modelled at high-resolution needs boundary conditions based on mesoscale and synoptic-scale meteorology. Nesting is an effective compromise between high-resolution modelling and efficient use of computing resources (Baklanov & Nuterman, 2009).

In the case of analysis with FLEXPART, ERA5 input, for example, has a horizontal resolution of up to 0.25° × 0.25° (Hersbach et al., 2020), and FNL has a resolution of 1° × 1° (National Centers for Environmental Prediction/National Weather Service/NOAA/U.S, 2000). For this project, we started with generating Enviro-HIRLAM simulations over a large geographical domain with 0.25° × 0.25° horizontal resolution, with subsequent downscaling chains of 0.15° × 0.15°, 0.05° × 0.05°, and finally down to resolution of 0.025° × 0.025° (roughly 2.5 × 2.5 km) over the Beijing metropolitan area and the surrounding region (see Figure 2 for a map of domains used in this project). For this project, using nested domains offers us the opportunity to resolve regional, subregional, and urban scale effects on meteorology, as well as effects of terrain which may become smoothed out in a coarse-resolution model run. In the high-resolution domain, we can also resolve impact of direct and indirect aerosol effects in more detail than we could in coarse-resolution runs. Moreover, high-resolution, when compared with low-resolution, can be especially valuable when studying the contribution of emissions along a trajectory simulated by FLEXPART, for application in air quality studies, when using high-resolution emission information. Higher horizontal spatial resolution is also very important in an urban location like Beijing, which is surrounded by complex terrain.

Figure 2. Map of the Enviro-HIRLAM nested model domains used in this study, with horizontal resolutions as follows, from largest to smallest: 0.25° (mother domain; green), 0.15° (cyan), 0.05° (blue), and 0.025° (red). The innermost nested domain is centered over the greater Beijing metropolitan area. Beijing is marked with the red cross.

The motivation for this project is to apply the FLEXPART model for atmospheric chemistry, aerosol processes and interactions and air quality modelling studies, especially during heavy haze events in megacities. In this study, we present a newly developed method for the combined Enviro-HIRLAM + FLEXPART modelling system and describe the method of pre-processing and integrating Enviro-HIRLAM meteorological output into FLEXPART. We compare the calculated airmass trajectories from the combined modelling system with those calculated from meteorological inputs from FNL and ERA5. In the future, we and other users can apply Enviro-HIRLAM + FLEXPART for studying air pollution episodes in China or any worldwide site.

The primary goal of this paper is to introduce our new implementation of FLEXPART which reads and processes Enviro-HIRLAM meteorology as input. Because the motivation for this work is to study air quality and haze formation, we also provide an example of how FLEXPART can be applied to air quality analysis, including a case study of a haze episode where we show some example atmospheric trajectories created with FLEXPART. The goal is not to try to prove that Enviro-HIRLAM meteorology is a better input to FLEXPART compared to meteorology from IFS and NFS models.

The FLEXPART trajectories used in this case study will be used in an upcoming modelling study on air quality to understand the formation, atmospheric chemical and physical processes, and lifecycle of the haze episode that occurred in Beijing and the surrounding region.

Our comparison of atmospheric trajectories will demonstrate the effects of different meteorological inputs when calculating atmospheric trajectories, and it underlines the necessity of understanding the strengths and weaknesses of trajectory analysis.

2. Methodology

2.1. FLEXPART

The predecessor of FLEXPART is the FLEXible TRAjectory model (FLEXTRA), which was developed in the early 1990s and was used to calculate mean wind trajectories (Stohl et al., 1995). The first version of FLEXPART, developed in the mid-1990s based on FLEXTRA, used ECMWF forecast model data as input and interpolated it onto its own Cartesian height levels. In this model, a user could simulate the release of a predetermined number of particles (for example 50,000 or 100,000 individual particles). This model version was evaluated with the European Tracer Experiment (ETEX) and was published by Stohl et al. (1998). FLEXPART has since been continually developed. For instance, Stohl and Thomson (1999) developed density correction and boundary layer physics in the model. The ability to use nested input was introduced in version 3 (mentioned in Stohl, Seibert, et al. (2005)). Seibert (2001) introduced the first convection scheme into version 4. Stohl et al. (2002) developed updated convection and turbulence schemes, and it introduced cluster analysis to generate an average trajectory path – what we now visualize as an atmospheric trajectory in FLEXPART. Seibert and Frank (2004) developed the source-receptor relationship that is still in use in the model today. Major updates were made in version 6 (Stohl, Forster, et al., 2005; Stohl, Seibert, et al., 2005), which included domain filing and the ability to define output nests in the model. Also in 2005, the ability to use NCEP’s GFS forecast data as input was introduced in version 6.4 (mentioned in Pisso et al. (2019) and on the FLEXPART website). Support for the newer GRIB2 file format for meteorology was added in version 8.2 in 2010 (mentioned on the FLEXPART website).

The most recent version of FLEXPART publicly available on the website, www.flexpart.eu, is version 10.4.3, which was released in November 2019 (Pisso et al., 2019). We chose to use this version because it is currently the latest publicly available version of FLEXPART. We downloaded the code for this project on 4 August 2021.

Technical details of FLEXPART, including a description of the trajectory calculations, atmospheric transport equations, vertical motion, convective and turbulent schemes, and handling of orography, see Pisso et al. (2019), Stohl, Forster, et al. (2005), Stohl et al. (2002), Stohl and Thomson (1999), and Stohl et al. (1998).

We ran FLEXPART on CSC’s Puhti HPC; technical details of the computing environment and modules we used are in Appendix A.

2.1.1. FLEXPART runs using various meteorological model input

In this study, we performed FLEXPART runs using as input meteorological data from the ECMWF Re-Analysis, version 5 (ERA5), NCEP GFS Final Operational Model Global Tropospheric Analyses (FNL), and Enviro-HIRLAM (Figure 1).

We obtained the ERA5 datasets using the Flex_extract v.7.1.2 software package (Tipka et al., 2020), with 1-hourly temporal resolution, and 137 vertical levels. The ERA5 datasets used in this study were created with the IFS model run at spectral truncation T639 (which corresponds to 0.28° resolution at the equator), and Flex_extract allows the datasets to be extracted at any resolution of choice. For this study, we obtained global datasets at 0.25° × 0.25° horizontal resolutions, which is the highest extraction resolution available.

Next, we downloaded FNL datasets directly from the National Research Data Archive (RDA) website: https://rda.ucar.edu/datasets/ds083.2. The FNL datasets available from this website have 6-hourly temporal resolution, 31 vertical levels, and 1° × 1° spatial resolution. Then we ran FLEXPART using these ERA5 (IFS+FP) and FNL (GFS+FP) datasets as input (Figure 1).

Next, we ran Enviro-HIRLAM in several modes with aerosol modes turned on and off, and we created a downscaling chain with four nested domains. Details of our Enviro-HIRLAM runs are described in Section 2.2. We then used the resulting meteorological fields from Enviro-HIRLAM as input to our newly developed version of FLEXPART. The results of FLEXPART calculations with these different meteorological inputs were compared to each other for evaluation.

2.1.2. Modifications to FLEXPART

To integrate and process the Enviro-HIRLAM meteorological output as input for FLEXPART, we have made modifications to the FLEXPART source code. In our implementation, we wanted to keep support for NCEP and ECMWF inputs, so that one compiled executable could be used for any of the available possible inputs, not requiring separate compilations. The high-level FLEXPART routines will detect which model the input meteorology is from and will automatically branch into either the old or new subroutines, based on the input provided (IFS, GFS, or Enviro-HIRLAM), and the respective subroutines will read and process the meteorology accordingly. As much as possible, we added new code in new subroutines separate from existing subroutines and kept new code independent from existing code to have minimal impact on existing program behavior. Thus, Enviro-HIRLAM functions like a new branch in the main code, for example, with calls to subroutines related to reading and processing Enviro-HIRLAM datasets placed inside if-else blocks that previously called into equivalent subroutines for reading and processing either NCEP or ECMWF produced datasets.

We tested the new code using ERA5 and FNL and compared it to the current version available online, and the results were identical, thus proving that our new additions are completely back-compatible with the current publicly available version of FLEXPART.

A detailed description of the modifications made to FLEXPART version 10.4 code can be found in Appendix B.

2.1.3. FLEXPART ensemble over a box

In addition to generating FLEXPART trajectories at the point of interest, which is the measurement site – Beijing University of Chemical Technology, Aerosol and Haze Laboratory, we defined a 20 km × 20 km horizontal box with the measurement site located in the center, and calculated trajectories at the four corners of the box as well as in the center of the box. The rationale for this is to test the sensitivity of the exact location, which shows how much impact a small horizontal difference of arrival location has over the course of the backward-trajectories.

2.2. Enviro-HIRLAM

In this project, we ran Enviro-HIRLAM version 7.2 (currently the latest available version) on a high-performance supercomputer (HPC) maintained by the Finnish Centre for Scientific Computing (CSC). Details of the computing environment are in Appendix A.

No modifications to the Enviro-HIRLAM source code were made in this project. All necessary code modifications were done in the FLEXPART source code. Postprocessing the Enviro-HIRLAM output is done with a Linux shell script, which calls the Climate Data Operators (CDO; Schulzweida, 2022). The script used for this project is included in Appendix C.

2.2.1. Model domains

We started with setting up Enviro-HIRLAM over a continental, large-scale geographical domain, which is large enough to enclose all of the 4-day backward trajectories arriving in Beijing during our case study. The outer domain (see Figure 2), with 0.25° resolution, covers most of Eurasia, extending from west-central Europe to eastern China and the Yellow Sea of the Pacific Ocean, with north-south extents from the Arctic Ocean toward the Middle East and South Asia. From Beijing, this domain extends more westward than eastward because during our case study, all trajectories generally arrived from the western directions, which is consistent with the prevailing winds in midlatitudes. In addition to being large enough for all of our backward trajectories, this domain is large enough to capture synoptic-scale meteorological processes.

The first downscaled nested domain, with 0.15° resolution, extends over China, Mongolia, and the Korean Peninsula. This domain also covers the synoptic scale in terms of meteorology, and certain features, such as fronts, can be more precisely resolved in this higher resolution. Moreover, this domain covers approximately the entirety of China.

The next downscaled domain, with 0.05° resolution, extends over the North China Plain (NCP) and southeastern Mongolia through the Yellow Sea. For purposes of air quality studies, NCP is one of the most significant geographical areas of interest because pollution from sources in the region could be transported into the Beijing metropolitan area (Zhao et al., 2013; Zheng et al., 2015). At this resolution, there are more insights into simulated regional/subregional and mesoscale meteorological effects that may not be apparent in larger and coarser runs.

Finally, the smallest domain with resolution of 0.025°Covers the Beijing-Tianjin-Hebei urban region. This domain was chosen to resolve small-scale meteorological features and effects on urban and local scales.

Enviro-HIRLAM utilizes a rotated system of coordinates, and hence, the domains shown above all use South Pole latitude of −80° and South Pole longitude of 0°. Our implementation of FLEXPART requires that all nested domains use the same South Pole latitude and longitude as the parent domain.

2.2.2. Aerosol effects and reference run

We ran Enviro-HIRLAM in four different modes: (i) Reference (REF), which has no aerosol effects included (in other words, a control run, with the aerosol module turned off); (ii) direct aerosol effects (DAE); (iii) indirect aerosol effects (IAE); and (iv) combined direct and indirect aerosol effects (CAE). We separately ran FLEXPART simulations for our evaluation, using output from each of these four modes and four resolutions as nested meteorological input to FLEXPART trajectory calculations.

The impact of aerosols on the atmosphere is widely known (e.g. Cheng et al., 2015; Wang et al., 2019, 2020; Zhang et al., 2017; Zhao et al., 2019). Aerosols are important in numerical weather prediction and climate related modelling/studies as there is growing concern on their impact on population’s health. The presence of aerosols in the atmosphere influence aerosol concentrations near the surface through changing meteorology. In particular, aerosols change optical properties of clouds, and hence, impact formation and development of clouds and precipitation at various scales. They impact through direct scattering and absorption of incoming solar radiation as well as through trapping of outgoing long-wave radiation.

The modes are defined in the Enviro-HIRLAM model experiment setup control file, where direct aerosol effects and indirect aerosol effects can be enabled or disabled. For more details, see Baklanov et al. (2017).

We hypothesize that during the haze episode in our case study, aerosol effects – both direct and indirect – will affect the meteorology and subsequently affect the atmospheric trajectories calculated by FLEXPART. On the other hand, during the clean period after the pollution episode, we expect that the trajectories simulated by FLEXPART with meteorology from the four modes of Enviro-HIRLAM will not be as significantly different from each other because there was lower aerosol loading to affect the meteorology over the greater Beijing metropolitan area.

2.2.3. Initial and boundary conditions, and observational data assimilation

For running the Enviro-HIRLAM modelling system, we extracted meteorological fields from ERA5 datasets (based on IFS model simulations) and atmospheric composition (aerosols and gaseous components) through the Copernicus Atmosphere Monitoring Service (CAMS) service. We extracted initial and boundary conditions at 3-hour time intervals utilizing the ECMWF’s Meteorological Archival and Retrieval System (MARS), downloaded through ECMWF’s ecgate. Moreover, using MARS, we also extracted as well as assimilated meteorological observations. These observation datasets are in Binary Universal Form for Representation of meteorological data (BUFR) format and are assimilated by Enviro-HIRLAM model. Note, the BUFR datasets contain observations every 3 hours, starting from 00 UTC. Furthermore, inventories of emissions from global biomass burning are from the Integrated System for vegetation fires (IS4FIRES; Sofiev et al., 2009). Emissions of gases, volatile organic compounds (VOCs), and selected aerosols are from Evaluating the Climate and Air Quality Impacts of Short-Lived Pollutants (ECLIPSE, 5th version; Stohl et al., 2015). Shipping emissions also include dimethyl sulfide (DMS; Nightingale et al., 2000).

Before any meteorological or NWP model tries to predict the future state of the atmosphere, it needs to know the current state of the atmosphere. Data assimilation consists of reading observational datasets that contain information about the current state of the atmosphere, and then interpolating this information into the grids used by the model. Such observational data can include, for example, measurements from weather stations, soundings (weather balloons), satellites, aircraft measurements, radar observations, and ship measurements. For an informative general overview of the concept of data assimilation in geosciences, see for example Carrassi et al. (2018). The method of data assimilation in Enviro-HIRLAM follows that of HIRLAM, as described in Randriamampianina and Storto (2008).

2.2.4. Chemistry and aerosol dynamics schemes

Enviro-HIRLAM uses several complex gas- and particle-phase chemistry schemes. Version 7.2, used in this project, uses the M7 (seven modes) aerosol microphysics module (Vignati et al., 2004) combined with the ECHAM5-HAM (Stier et al., 2005, 2006) for removal processes. Full details of the chemistry and aerosol dynamics schemes and modules in Enviro-HIRLAM are described in Baklanov et al. (2017) and references therein.

2.2.5. Model output

Output data from Enviro-HIRLAM are stored in GRIB (General Regularly-distributed Information in Binary format) files. GRIB is a World Meteorological Organization (WMO) standard for storing output from NWP forecasts and re-analysis data (World Meteorological Organization WMO, 2020). Both ECMWF and NCEP use GRIB format for their NWP models, which means FLEXPART already incorporates the necessary tools for reading GRIB files. Therefore, support for reading Enviro-HIRLAM output as input to FLEXPART does not require any modifications for reading the datasets (i.e. no new libraries are needed to be included in FLEXPART compilation).

The FLEXPART working directory contains a text file called AVAILABLE, which provides FLEXPART with a list of input meteorology files to use for its calculations. This file contains a list of Enviro-HIRLAM output files. For details of the file structure and example of an AVAILABLE file for Enviro-HIRLAM output, see Appendix C.

2.2.6. Merging of the model output files

There are several types of output files from each Enviro-HIRLAM run, with output files saved at every 3-hour interval. At each output time, there are three files generated which contain the meteorological variables needed for FLEXPART: forecast (*.fc), diagnostics (*.md), and verification (*.ve), which must be pre-processed and merged into a single file before reading it into FLEXPART. A detailed description on how the files are merged, along with an example script, is provided in Appendix C.

2.3. Case study

For our case study, we ran FLEXPART in backward mode, with the arrival of trajectories at a measurement site in Beijing, China during a severe regional air pollution episode from 14 to 15 November 2018. Our case study includes a severe air pollution episode followed by a synoptic weather system which rapidly brought in a new airmass with clean air. For evaluation, we compared 4-day (96-hour) backward atmospheric trajectories arriving at a measurement site in Beijing with 5 arrival times, two times during the haze episode (T1 and T2), one time during the weather shift (T3), and two times after the weather system has passed (T4 and T5). Note that Beijing local time (LT) is 8 hours later than UTC.

The evaluation times are as follows:

• T1: 2018-11-14, 12 UTC (20:00 LT) – during the haze episode.

• T2: 2018-11-14, 18 UTC (02:00 LT) – near the end of the episode, when aerosols and haze were the heaviest.

• T3: 2018-11-15, 00 UTC (08:00 LT) – during the synoptic weather shift.

• T4: 2018-11-15, 06 UTC (14:00 LT) – after the synoptic weather shift, when the air was clean and some high clouds lingered.

• T5: 2018-11-15, 12 UTC (20:00 LT) – well after the synoptic weather shift, when air is clear and cloudless, and aerosol loading at the measurement site was comparatively low.

We expect that FLEXPART will show airmasses sourced from a different region after the weather shift, and this also allows us to evaluate how different models predict the change in weather conditions. We also hypothesize that because the heavy aerosol loading during the haze episode (T1–T2) will have both direct and indirect aerosol effects feeding back on the local meteorology, there will be a difference in trajectories calculated using Enviro-HIRLAM with aerosol effects included (DAE, IAE and CAE) version the Enviro-HIRLAM reference run.

In our evaluation, we used 500 m arrival heights because at night, the boundary layer height is very low, and choosing an arrival height above the nighttime boundary layer avoids influence of local emission sources at nighttime. Moreover, this height avoids the influence of tall buildings surrounding the site, which could interfere with the observations.

2.3.1. Measurement station in Beijing, China

The site of interest is the Aerosol and Haze Laboratory (AHL), an atmospheric and air chemistry research station at Beijing University of Chemical Technology (BUCT). This station (hereafter, BUCT-AHL) is in the west-central part of Beijing, near the 3^rd Ring Road, ca. 9 km west-northwest of the center of the city. The exact location is 39°56’30“ N, 116°17’50” E. We chose this test case in part because our research groups at the University of Helsinki and BUCT have been working on several studies involving the analysis of severe air pollution that occurred in Beijing and Northern China in November 2018, with focus on this measurement station. More details about this station are described in Liu et al. (2020).

The air pollution episode that occurred between 10 and 15 November is of particular interest to study because of its intensity, with PM_2.5 values reaching as high as 250 µg/m³ during the peak of the episode. This episode, with its heavy aerosol loading, is a good example for testing combined Enviro-HIRLAM + FLEXPART modelling system because it can demonstrate the usefulness of accounting for aerosol effects on different scales and model resolutions when generating weather and air quality forecasts.

The output grid from FLEXPART must be wholly inside the Enviro-HIRLAM largest domain (i.e. with a resolution of 0.25° × 0.25°. In FLEXPART’s terminology, this is the mother domain. Note that it is not necessary to restrict the output grid to the size of any of the high-resolution nests. If the position of a particle is inside one of the nests, FLEXPART will calculate parameters based on meteorological input from the corresponding nested input. If the position of a particle is outside the nested domains but inside the mother domain, it will calculate parameters based on meteorological input from the mother domain. For a full description of how FLEXPART determines whether to use meteorology from nested input or from the mother domain input, see Stohl, Forster, et al. (2005) and Section 2.1. Therefore, we chose a longitude range of 2.25º–132º E and a latitude range of 29º–62º N, and we selected an output grid of 0.1° × 0.1° horizontal resolution. This is defined in FLEXPART’s OUTGRID file, and for consistency, we used the same output parameters throughout all of our FLEXPART calculations.

2.3.2. Timeline of synoptic meteorological events

The air pollution/haze episode was set up by persistent high pressure over the NCP region, which set in around 12 UTC on 9 November. This synoptic pattern continued through 14 November, with pollution gradually increasing over time. This is consistent with typical pollution episodes in the region, described in Zheng et al. (2015). On 13 November, there was heavy aerosol loading across the NCP, which is also visible in satellite imagery (Figure 3). On 14 November, a cloud band approached from the west, and the daytime satellite imagery on 14 November showed that the clouds and aerosols became indistinguishable. In the early hours of 15 November (starting around 22 UTC 14 November), a synoptic weather feature moved through and relatively rapidly changed the conditions. By 16 November, the clear air throughout the region was evident in satellite imagery.

Figure 3. Visible satellite imagery of the north China plain, taken from the visible infrared imaging radiometer suite (VIIRS). Beijing is marked with the red cross. The images were captured shortly after noon local time on each day as the satellite passed overhead (i.e. 13 Nov. 04:38 UTC/12:38 LT; 14 Nov. 04:20 UTC/12:20 LT; 15 Nov. 04:03 UTC/12:03 LT; 16 Nov. 05:22 UTC/13:22 LT).

Figure 4 shows a time-series of meteorological parameters (air temperature, relative humidity, atmospheric pressure, precipitation, visibility, wind speed and direction) observations at BUCT-AHL, along with PM2.5 concentration, which serves as a general measure of air pollution during the haze episode.

Figure 4. Time-series of meteorology (air temperature, relative humidity, atmospheric pressure, precipitation, visibility, wind speed and direction) and PM_2.5 concentration from 12 UTC 14 November (T1) through 12 UTC 15 November 2018 (T5), observed at the BUCT-AHL station in Beijing, China.

Meteorological instruments at the BUCT-AHL station recorded a very distinct change in meteorological and air quality conditions starting around 22 UTC 14 November (6:00 am LT on 15 November) which lasted for a short period (as seen in Figure 4). Air temperature measured at the station increased by about 2°C (to approximately 9°C), atmospheric pressure rose by over 3 hPa to 1018 hPa, relative humidity dropped from 80% to 40%, visibility improved from less than 1 kilometer to over 20 km, and a wind direction shift from west to north was observed. Additionally, a rapid drop in concentrations of PM2.5 from nearly 200 µg/m3 to less than 10 µg/m3 was observed during the case study period.

These observations show some characteristics of a frontal passage over the measurement station, even though it is not a well-defined cold front described in Ahrens and Henson (2021). Surface weather maps obtained from the Japan Meteorological Agency showed a significant change in the atmospheric pressure and weather pattern between 14 and 15 November, even though they do not definitively mark a front (these charts can be found in Appendix D).

According to the surface charts, at 00 UTC 14 November, there was a slack pressure gradient in the Beijing area, with an area of low pressure to the north. By 18 UTC on the 14th, the low-pressure system had moved eastward and the trailing trough was just south of Beijing. The pressure gradient in the Beijing area increased as the interior high pushed a bit eastward, and eventually a pressure gradient moderately changed over Beijing. Such development was enough to clear out the cold pool and bring cleaner air from the northwest, including downslope flows from the mountainous regions, where air is almost always clean (Wang et al., 2019). Therefore, this synoptic shift and subsequent change in weather conditions is consistent with cold air pool mixout and the rapid cessation of the air pollution episode.

The haze episode serves as a good test case for evaluating the combined modelling system because the heavy aerosol loading will affect the local meteorology through the aerosol effects on radiation (as shown in Wang et al. (2020) and Zhao et al. (2019)). This would subsequently affect the calculation of backward trajectories. Therefore, we would expect that during the heavy haze, there would be a difference between FLEXPART trajectories based on meteorological fields calculated from the Enviro-HIRLAM simulations including aerosols (DAE, IAE, and CAE) vs. without aerosols (REF). Based on the satellite imagery in Figure 3, the haze is seen without clouds around T1. Thus, we would expect DAE to contribute more to the total aerosol effects than the IAE. At the end of the pollution episode, closest to T3, the clouds are the heaviest, and we would expect that IAE would play a bigger role in the trajectory calculations than DAE.

Additionally, the relative humidity during this haze episode and the subsequent drop of humidity during the mixout period is consistent with previous studies on such haze episodes (Cheng et al., 2015; Zhao et al., 2019; Zheng et al., 2016; Zheng et al., 2015). Because of the high relative humidity during the event, we predict that the indirect aerosol effect is likely to play an important role in the overall aerosol effects, in addition to direct aerosol effects.

During the clean period after the mixout period, when there is low aerosol load, we would expect a less significant difference between trajectories that use output from Enviro-HIRLAM runs with aerosol effects versus the reference run.

Additionally, the mixout period (T3) is also an interesting test case because of the expected shift in wind direction and change in airmass. When calculating the backward trajectories before and after the passage, we would expect the trajectories to originate from noticeably different regions.

3. Evaluation

3.1. Comparison of FLEXPART trajectories with different meteorological inputs

Figure 5 shows atmospheric backward trajectories using the inputs described in Section 2.1 plotted for the case study mentioned in Section 2.3.

During the stagnant conditions of the haze episode (T1; Figure 5(a)), the trajectories calculated with Enviro-HIRLAM’s meteorology show atmospheric transport from farther away regions compared to trajectories calculated using FNL and ERA5 datasets. However, based on this figure alone, it is difficult to discern a difference between the different trajectories for this arrival time.

Figure 5. Spatiotemporal positions of atmospheric backward trajectories calculated by FLEXPART using different meteorological inputs (from model runs: IFS-ERA5 (extracted at 0.25° horizontal resolution), GFS-FNL, and enviro-HIRLAM (EH), with full downscaling chain, without aerosol effects (REF) and with aerosol effects (CAE)) arriving at the BUCT-AHL measurement station at height 500 m above ground level with arrival times of a.) T1; b.) T2; c.) T3; d.) T4; and e.) T5. In addition to the center point over the station, an ensemble representing a 20 km × 20 km box (described in Section 2.1.3) is plotted in lighter colors, with the four corners (northeast, northwest, southeast, and southwest). Vertical profiles of trajectories, as a function of altitude (km) vs. distance (km) to the station along trajectories, are also plotted above each subplot map, together with the underlying terrain. The numbers on the trajectories indicate the time (in days) before arrival at the station.

At T2 (Figure 5(b)), the trajectories calculated with Enviro-HIRLAM meteorology again reach farther away compared to the trajectory calculated using ERA5 meteorology. The trajectory generated with FNL has already begun to change, indicating that FNL is already forecasting the beginning meteorological system passage over the station at this time. This difference could be partly because of FNL’s coarser temporal resolution (6 hours versus 3 hours for Enviro-HIRLAM or hourly for ERA5). Being at the latter end of the air pollution episode, the observation that FNL differs most from the others agrees with the finding of Zuo et al. (2023).

Figure 5. (Continued).

The trajectories with different meteorological inputs differ most during the passage of the synoptic meteorological system (T3 and T4; Figure 5(a,b)). Most likely, this is because the different meteorological models forecast the passage of the weather system at slightly different times, depending on their scales, processes accounted, physical parameterizations used, spatial and temporal resolutions, etc.

After the meteorological system passed (T5; Figure 5(e)), there is once again more agreement between the trajectories calculated with different inputs. During this time with stronger westerly winds, there is good agreement between the trajectories for the first day backwards in time, but after that, the trajectories calculated with Enviro-HIRLAM deviate from the trajectories calculated with FNL and ERA5.

We will look further into the differences and spread between trajectories in Figure 6.

In Figure 6, we see that there is more agreement between the trajectories (calculated based on different input datasets – ERA5, FNL, and Enviro-HIRLAM CAE) arriving during the haze episode on 14 November (with daily averages plotted as solid black lines), compared to the trajectories arriving the next day (Nov 15^th), especially during the 36–48 hours closest to arrival. However, the spread increases considerably when tracing back in time for more than two days. Along with more spread, we see more constantly increasing spread in trajectories arriving on 15 November. Since the distance between two trajectories is related to the length of the trajectory, it can be useful to consider the relative spread, which we define as the Euclidean distance between two trajectories at any given hour, divided by the distance to the station along the reference trajectory at that time (Figure 6(b)). The relative differences show more spread on 15 November compared to the 14^th during the 24–48 hours closest to arrival, but less or equal spread further back in time. This is because on the 15^th the airmasses were arriving with a large-scale air flow, and the trajectories were generally straighter and longer, as opposed to the shorter trajectories arriving on the 14^th with the more stagnant air flow. The relative differences often show extreme values very near the station, where the distances travelled are small, and therefore even small deviations between trajectories produce large relative values.

Figure 6. Differences in atmospheric trajectories. Each colored line represents the difference between two trajectories, arriving at the same time but calculated using different meteorology. The six panels compare all three FLEXPART realizations (input from ERA5, FNL, or enviro-HIRLAM with combined aerosol effects (CAE)) with each other. The label on the right shows the reference trajectory; the axis on the left shows the distance from comparison trajectory to reference trajectory (panel a.) or relative difference, defined as distance to station (along trajectory) versus distance to reference trajectory (panel b.) solid black lines show the average for 14 November (during the haze episode) and dashed black lines show the average for 15 November (after the haze has cleared out).

Figure 6(a) shows that the further back in time the model goes, there are greater differences between trajectories using different meteorology input. In the same way that running a forecast model will lead to more uncertainty with longer forecast length. What is interesting, though, is that in the cases shown in Figure 6(a), there is more agreement within about 36 hours of arrival, but beyond that, the trajectories calculated with different model meteorology diverge significantly from each other. This can be seen as a minimum point in the relative differences in Figure 6(b).

Figure 6. (Continued).

Figure 6 also shows that trajectories calculated with meteorology from FNL and ERA5 datasets were more similar to each other than to trajectories calculated with meteorology from Enviro-HIRLAM. Of the two global models, results using Enviro-HIRLAM were closer to results with FNL. Perhaps the most likely reason for the larger difference between trajectories from Enviro-HIRLAM and the other two is due to the higher resolution of the Enviro-HIRLAM model simulations, especially inside the highest resolution nested domain. Moreover, Enviro-HIRLAM is a regional model run over a selected geographical area, whereas the GFS and IFS models cover the whole globe. This would result in different scale meteorological effects modelled in Enviro-HIRLAM which may be smoothed, parameterized, and/or not resolved in the GFS and IFS models due to the lower resolution.

Another possible cause for the differences between trajectories is the vertical structure or the number of vertical layers used in the models. As noticed in Figure 5, the calculated trajectories based on Enviro-HIRLAM meteorology have a slightly different vertical profile and considering different wind speeds and directions at different heights, this could result in different spatial positions of trajectories.

Furthermore, the difference could also be partly attributed to the aerosol effects when trajectories are calculated with Enviro-HIRLAM CAE meteorology. In addition to comparing trajectories using input from three different meteorological models, we would also like to compare trajectories using input from Enviro-HIRLAM runs with different aerosol modes (DAE, CAE, IAE, and REF). These are shown in Figure 7.

Figure 7. FLEXPART trajectories calculated based on enviro-HIRLAM meteorology with combined aerosol effects (CAE), direct aerosol effects (DAE), indirect aerosol effects (IAE), and control/reference (REF) without aerosol effects included. Arrival times are as follows: a.) T1; b.) T2; c.) T3; d.) T4; and e.) T5. The numbers on the trajectories indicate time (in days) before arrival at the station.

Figure 7. (Continued).

Figure 7. (Continued).

We found small but noticeable differences between the trajectories when the different aerosol modes were used. During the haze episode at T1, when there is heavy aerosol loading but only a few visible clouds, trajectories calculated using DAE and CAE meteorology were very similar to each other, which differed from trajectories with the IAE and REF. When aerosols are actively interacting (absorbing, reflecting) with sunlight, the role played by the direct aerosol effect makes sense, especially during the daytime.

At T2, all of the trajectories appear almost identical. Since this is overnight the day before arrival, this can be somewhat expected, at least for the first half day, because there is no sunlight to be reflected by aerosols.

At T3 during the shift in meteorological conditions, the trajectories using IAE and CAE were most similar, differing from the DAE and REF. This is when the cloud banks are heaviest in the region, as seen in the satellite imagery in Figure 3, and thus it makes sense that there is a strong influence of the IAE on the trajectories. In this case, the sunlight is not interacting directly with the aerosols, but rather it interacts with the clouds that are formed with involvement of the aerosol particles. Namely, the clouds, which are white, reflect sunlight back to space, resulting in less sunlight reaching the surface and thus different temperatures and meteorological conditions. These results are consistent with the prediction we made in Section 2.3.2.

As the meteorological conditions continued to change, at T4, we found that the IAE and DAE seem to offset each other. Finally, at T5, late evening after the weather system has passed over the measurement site, there is a noticeable difference in the calculated trajectories using inputs with different aerosol effects included, but it is unclear whether direct or indirect aerosols played a more significant role. Notably, at T5, the trajectory based on CAE is in between the trajectories based on DAE and IAE, which varies from the REF trajectory.

3.2. Potential emission sensitivity

In addition to comparing the trajectories, we compared the potential emission sensitivity (hereafter, emission sensitivity) of the FLEXPART simulations. This is also called source-receptor relationship (Pisso et al., 2019; Seibert & Frank, 2004). One way to look at it is like a “plume” (particles-dispersion) backward in time instead of a single line (trajectory-pathway). For many purposes, visualizing a trajectory as a line could be sufficient, for example for general insight of the origin of an airmass. Additionally, plotting a trajectory as a line is more readable than creating a 3-dimensional visualization over time and space, which is why a single continuous connecting line is commonly used to represent a trajectory (Stohl et al., 2002). However, plotting a trajectory as a single line has also limitations. In particular, when we see clear differences between FLEXPART simulations using different meteorological input, we should be careful before drawing solid conclusions about the potential source region of an airmass origin based on a single mean trajectory plotted as a line. This is even more important when, for instance, using trajectories as a basis of determining emission sources. For example, in our air quality applications, emission sensitivity over an area may be more useful than a single line trajectory, especially when there may be individual point sources that contribute to air pollution, which could be missed by the single-line approach.

Figure 8 shows emission sensitivity, seen in this figure as a series of shaded plumes, for the backward trajectories used for this study. The plumes show the source receptor relationship (SRR) values described by Pisso et al. (2019), calculated in the FLEXPART simulations. This figure shows the sum of SRR values in all output levels up to 1500 m above the ground. The SRR can also be interpreted as sensitivity to emissions with regard to the geographical location from where the airmasses which arrive at the measurement station are originally sourced.

Figure 8. Spatiotemporal distribution of dispersion of particles (modelled by FLEXPART using GFS-FNL, IFS-ERA5, and enviro-HIRLAM CAE meteorological inputs), calculated as emission sensitivity (in units of seconds), arriving at the BUCT-AHL measurement station with arrival times at: a.) T1; b.) T2; c.) T3; d.) T4; and e.) T5. The emission sensitivity is related to the residence time of dispersed particles in the atmosphere, and it gives an idea of the source area of potential emissions of pollutants that could eventually arrive at the station.

In Figure 8, we get a different visualization for potential source regions of airmasses, this time spread over the geographical area rather than simply along single trajectories. Interestingly, there is a much larger geographic/spatial spread in the residence times when the Enviro-HIRLAM meteorology is used compared to meteorology from FNL or ERA5 datasets. This, along with the different placement of the particles, affects the mean trajectories that are calculated. This demonstrates the importance of looking at a distribution of released particles, rather than drawing conclusions only from a single line.

The difference is especially noticeable at 00 UTC on 15 November (T3). When looking at the plumes using FNL and ERA5 meteorology (Figure 8(c)), and likewise when looking at trajectories as lines (Figure 5(c)), it appears that the airmass source is solely from the northwest. However, when looking at the plumes using Enviro-HIRLAM meteorology (Figure 8(c)), we get a slightly different outcome: although the majority of the source footprint is also from the northwest, a good portion of the emission area is also from the southwest. In this case, the mean trajectory, which is in between the two separate plumes, does not give a good representation of either plume, and might even misrepresent the true story.

Figure 5. (Continued).

Additionally, the plumes based on input from all three meteorological models have the most spread in dispersion during the weather system passage (T3–T4), implying that there will be larger uncertainty in the calculated trajectory. This is likely due to the fast-changing conditions, and thus, a small difference between two different particles’ paths results in greater change as they are projected further back in time.

To offer better quantification of the emission sensitivity, we attributed obtained results into quadrants around BUCT-AHL, with the percentage of source from each quadrant shown in Table 1.

Table 1. Percent of emission sensitivity values calculated by FLEXPART (based on meteorology from 3 different models, including Enviro-HIRLAM model output for 4 aerosol modes), divided into quadrants (NW, NE, SW, and SE) around the BUCT-AHL station. Darker shading represents higher values.

Analyzing the source region based on quadrant is especially significant in air quality analysis in Beijing because air from the northwest is from low-populated regions, where the source airmasses tend to be cleaner; on the other hand, the region south and southwest of Beijing contains heavy industry, and air originating from this area tends to be more polluted (Hakala et al., 2022; Wang et al., 2019; Zhang et al., 2017). Using the data from the trajectory compiled into quadrants, as shown in Table 1, provides useful insight on regional emission sources, as opposed to single point emissions along a line.

For example, at T2, we can see from Figure 8 and Table 1 that regardless of which NWP model meteorology we use as input in FLEXPART, most of emissions is from the southwest, and considering that this quadrant is with multiple emission sources, we would expect the airmass arriving at the station to be influenced by the emission sources from that region. If we looked only at a single trajectory calculated by FELXAPRT (based on FNL meteorology) at T2, we would have guessed the airmass was from the west. In this case we might have completely underestimated the influence of emission sources from the southwest quadrant, and hence, we might have incorrectly assumed that clean airmass arrived at the station. This demonstrates importance in estimation of emission sensitivity vs. different meteorology input in the FLEXPART model trajectory calculations.

Therefore, it is important to look at all possibilities before analyzing and interpreting trajectories to draw solid conclusions. Using a different choice of meteorology or approach to visualization of different types of data might tell a different story.

3.3. Comparison of enviro-HIRLAM with nesting vs. without nesting

Finally, we decided to do an analysis of running FLEXPART using meteorology from only the 0.25° resolution Enviro-HIRLAM run and comparing it against FLEXPART using Enviro-HIRLAM with the full nesting/downscaling chain (i.e. with different resolutions ranging from 0.25 to 0.025 degrees). Note, meteorology is based on the CAE mode of Enviro-HIRLAM run. We used the same FLEXPART output grid from both Enviro-HIRLAM runs. Results from this comparison are shown at different times in Figure 9.

As expected, there are somewhat more details in the spatial distribution of atmospheric trajectories calculated based on Enviro-HIRLAM meteorology from a full downscaling chain.

During the low-wind period during the haze episode (T1–T2), as shown in Figure 9(a,b), trajectories using the full downscaling change for input traverse a slightly longer path compared to the trajectories using only 0.25° horizontal resolution input. Several studies (e.g. Giannaros et al., 2017; Kalverla et al., 2019; Schlager et al., 2019; Solbakken & Birkelund, 2018) have concluded that high-resolution limited area models, such as HARMONIE, HIRLAM, WRF, and any other NWP models, have a tendency to slightly overestimate wind speeds during low wind conditions, and this is especially true in complex terrain. This could also influence our case study when using the high-resolution input.

Figure 9. Atmospheric trajectories calculated by FLEXPART without nesting (based on enviro-HIRLAM CAE meteorology at 0.25° horizontal resolution) and fully nested (i.e. mother domain with all three higher-resolution nests & based on enviro-HIRLAM-CAE meteorology at 0.25°, 0.15°, 0.05° and 0.025° horizontal resolutions, respectively) at arrival times at: a.) T2; b.) T2; c.) T3; d.) T4; and e.) T5. The numbers on the trajectories indicate time (in days) before arrival at the station.

Figure 9. (Continued).

Figure 9. (Continued).

However, the differences in using the full downscaled chain versus only the 0.25° horizontal resolution domain are not significantly different, especially at the end of the case study period (T5) when the clear air is present, and the trajectories are almost identical and cross over each other. In this case, if we are interested in general meteorological patterns as opposed to looking at very specific features, e.g. meteorology around specific terrain, then perhaps for many use-cases, only the 0.25° resolution input would be sufficient for the FLEXPART runs, without the need for the extra computation expense of using nests.

This result is somewhat contrary to our initial hypothesis that there would be a significant difference between trajectories calculated by FLEXPART using only the coarse-resolution meteorology (i.e. without nesting) versus trajectories calculated using the full nesting. The differences turned out to be smaller than expected. At the ends of the trajectories (4 days backward in time) the differences between the trajectories calculated with fully nested input versus no nests, were between 50 and 200 km. On one hand, this similarity in trajectories gives confidence in the obtained results, in particular that the different resolution meteorology datasets have agreement with each other in terms of meteorology, and it also gives confidences in FLEXPART’s algorithm for calculating trajectories within nests. Despite the comparable results, even a difference of 50 km in trajectory path might be significant, depending on the use case. For example, when looking at source emissions for air quality modelling, a 50 km difference in trajectory calculation could lead to different results. Therefore, future users should consider also these differences when deciding on whether or not to use nesting in their model setup. For some use cases, such as preliminary evaluation and analysis, calculation based on low-resolution meteorology (for example, 0.25°) might be suitable, but the full nesting (including several resolutions) would be preferable for in-depth air quality studies.

4. Conclusions and future work

In this study, we applied the Enviro-HIRLAM seamless, online integrated meteorology – atmospheric composition modelling system’s meteorology as an input option for the FLEXPART particle dispersion and trajectory model. For our selected case study of a severe pollution episode on 14 November 2018 in Beijing, China, we found that accounting for aerosol direct and indirect effects in simulation meteorological fields showed a small – yet not insignificant – effect on the FLEXPART calculated trajectories. This is important, for instance, when using atmospheric trajectories to determine sources of airmasses and emissions in air quality studies.

Moreover, the results show there is a significant difference in trajectories simulated by FLEXPART when different input meteorology from different NWP models is used. This is evident during the stagnant, low-wind conditions prior to the cold front passage, and it is especially noticeable during the sudden change in synoptic weather conditions. Most likely, this is because meteorological models simulate the timing of the weather system passage over the measurement site slightly differently.

Elaborating on the aerosol effects, we found that had effects on the resulting trajectory on not only the polluted time (as expected) but even on the cloudless, non-polluted time. Perhaps the aerosols have more effects on clear days than originally expected because clouds exist regardless of aerosols, but without clouds, there was more impact by the aerosols. We speculate that on the cloudy time during the pollution episode, the cloud cover overshadowed the aerosol effects, and the aerosol direct and indirect effects were negligible in the cloudy case.

We also found when analyzing the four points in a box around the BUCT-AHL station, rather than a single arrival point, that a small difference in target place affected the trajectory and emission source footprint. This highlighted the usefulness of looking at trajectories for more than just a single arrival point but also nearby arrival points, for a more comprehensive analysis.

Finally, we demonstrated that using nesting (downscaling from regional to subregional and urban scales) from lower to higher resolutions resulted in slightly different meteorological input to FLEXPART which subsequently calculated different trajectories. Although the differences may be small, they could be significant for in-depth air quality analysis, especially when looking at emission sources.

This study also highlights the importance of avoiding over-reliance on a single trajectory when performing meteorological analysis involving airmass sources. Every model and type of model run with various parameters/settings has its own strengths and weaknesses, and each has places where they are beneficial, depending on the desired case study. Our conclusion is not to claim that one input or model setup is more accurate or reliable than another, but rather, we point out the differences as our new work offers a new option to use Enviro-HIRLAM meteorology as input to FLEXPART. In some cases, it may be useful to look at results from multiple different model runs (e.g. using different inputs or different model setups), and compare the results of the different runs.

We believe that FLEXPART trajectory simulations generated using Enviro-HIRLAM meteorology as input can complement FLEXPART simulations using input datasets from ECMWF or NCEP, and comparing the results can offer a slightly different picture than looking at results based only on one model as input.

Now that we have developed the tools to integrate and process the Enviro-HIRLAM meteorological output as input into the FLEXPART model, we plan to use the resulting FLEXPART trajectories in SOSAA (model to Simulate the concentration of Organic vapors, Sulfuric Acid, and Aerosols), which is a column model developed at the University of Helsinki to study atmospheric aerosols and chemistry Chen et al. (2021). One of our applications of SOSAA is to analyze and understand severe haze episodes. We will run SOSAA along trajectories modelled by FLEXPART, using emissions along the trajectories to include in the simulations of chemistry and aerosol processes. Because of the heavy aerosol loads during the air pollution episodes, our hypothesis is that using meteorological output from Enviro-HIRLAM, rather than GFS or IFS models, will improve the quality of the SOSAA simulations because of the Enviro-HIRLAM includes chemistry and aerosol feedback effects when performing meteorology – atmospheric composition forecasts. We also believe that the high-resolution nested domain focused on the Beijing urban area will improve the simulation of chemical processes in the area during the time of low and variable wind speeds, because lower-resolution models are not capable of resolving smaller-scale features of meteorological processes.

We are also planning to use this newly developed approach in other cities around the globe. In addition to focusing on air quality, the combined Enviro-HIRLAM + FLEXPART modelling system can be useful for other studies that involve trajectory analysis and emission sources. Moreover, such modelling system will lead to better understanding of atmospheric chemistry and aerosol effects on meteorology, which can be upscaled to larger-scale climate and Earth system models.

Thus, this work presented the developed, tested and applied the Enviro-HIRLAM-FLEXPART modelling system, which uses meteorology generated at different resolutions and domains as input for atmospheric trajectory calculations. It has been our first step for enlarging, improving, and enhancing studies on air quality and impact of aerosols on meteorology on different scales.

Potential applications of the combined model system include but are not limited to: air quality analysis, tracking smoke and pollution, emissions of air pollution, evaluation of atmospheric hazards, source apportionment of observed atmospheric chemical pollutants, and many other applications. These applications are invaluable to human health, safety, and scientific research.

Acknowledgements

We wish to acknowledge CSC – IT Center for Science, Finland (www.csc.fi/csc), for computational resources, technical support and advice used in this project.

We also wish to acknowledge the FLEXPART development team (FlexTeam, Univ. of. Vienna; flexteam.univie.ac.at), especially Andreas Plach and Lucie Bakels, for technical support, advice, and assistance on setting up and running FLEXPART.

CSC (Atos BullSequana) HPC & ECMWF (CRAY-XC) HPCs were used. The ECMWF boundary conditions, meteorological and air quality observations/datasets are utilized in this project.

The work is supported by Enviro-PEEX(Plus) on ECMWF (www.atm.helsinki.fi/peex/index.php/enviro-peex_plus), Enviro-HIRLAM on CSC-HPC projects, AoF ACCC, H2020 CRiceS.

Additionally, we acknowledge Prof. Jim Steenburgh, with whom we had a brief conversation about the meteorological conditions in Beijing, and he provided insight into the synoptic meteorology timeline during the case study period.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

1. New FLEXPART code

The FLEXPART code used for this publication has been uploaded to Zenodo, available at https://doi.org/10.5281/zenodo.8300429. The code will also be available in a repository on the FLEXPART website at https://doi.org/10.5281/zenodo.8300429, and new versions will be published here. We suggest potential users who wish to use the Enviro-HIRLAM-FLEXPART modelling system use the version that will be published on the FLEXPART website. We also plan to work with the FLEXPART team to integrate this into a future version of the public FLEXPART version. We suggest that users refer to the official FLEXPART website for the latest version.

2. Enviro-HIRLAM post-processing script

An example script to select and merge the Enviro-HIRLAM output for use in FLEXPART is in Appendix C of this manuscript.

3. Enviro-HIRLAM model code

The Enviro-HIRLAM modelling system is a community model. The source code is available for non-commercial use (i.e. research, development and science education) upon agreement through contact with Alexander Mahura ([email protected]), Bent Sass ([email protected]) and Roman Nuterman ([email protected]). Documentation, educational materials and practical exercises are available from http://hirlam.org and hirlam.org/index.php/documentation/chemistry-branch, and Young Scientist Schools (netfam.fmi.fi/YSSS08, www.ysss.osenu.org.ua; aveirosummerschool2014.web.ua.pt and https://megapolis2021.ru).

Supplementary Material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/20964471.2024.2316320

Additional information

Funding

Funding for this project came in part from the Jenny and Antti Wihuri Foundation project, with the grant for “Air pollution cocktail in Gigacity. ”Funding was also received from the Research Council of Finland (formerly the Academy of Finland, AoF) project 311932 and applied towards this project. Partially, funding included contribution from EU Horizon 2020 CRiceS project “Climate relevant interactions and feedbacks: the key role of sea ice and snow in the polar and global climate system” under grant agreement No 101003826; and AoF project ACCC “The Atmosphere and Climate Competence Center” under grant agreement No 337549. Additionally, internal funding at the University of Helsinki’s Lahti University Campus was used for work on this project.

Notes on contributors

Benjamin Foreback

Benjamin Foreback is a PhD student in atmospheric science at the Lahti campus of the University of Helsinki. He has bachelor’s degrees in physics and atmospheric science from the University of Utah and a master’s degree in atmospheric science from the University of Helsinki. His research includes modelling and observations of air quality, with a focus on severe air pollution episodes in megacities. His anticipated graduation is early spring 2024.

Alexander Mahura

Alexander Mahura is University Researcher at Institute for Atmospheric and Earth System Research, University of Helsinki (UH), Finland since 2017; received BSc in 1991, MSc in 1998, and PhD Phys&Math Atm.Sci. in 2002; worked in Russia, USA, Austria, France, Denmark, and Finland. His research interests include online integrated multi-scales and -processes modelling of meteorology and atmospheric composition, atmospheric chemistry, atmospheric boundary layer processes, numerical weather prediction, fine-scale road weather and birch pollen forecasting, statistics for data analysis and post-processing models output, environmental impact and risk assessment.

Petri Clusius

Petri Clusius has worked with atmospheric modelling since his bachelor studies in 2017, and has focused heavily on new particle formation in boreal forests. His master thesis work was the development of the new Atmospherically Relevant Chemistry and Aerosol Box Model – ARCA box. In his PhD studies he is now focusing on the contributing role of biogenic hydrocarbons emitted from forests to aerosol concentrations and their effect to cloud condensation nuclei number concentration. In this work he is using the ARCA-model together with its 1D version SOSAA and developing it to a Lagrangian trajectory model.

Carlton Xavier

Carlton Xavier is a post-doctoral researcher at Lund university and SMHI who’s work focuses primarily on the formation and growth if secondary aerosols in polar regions. To accomplish this, he uses a detailed chemical and aerosol column model ADCHEM which is run along Lagrangian trajectories generated using Flexpart. The overarching aim of his research is to understand the radiative effects of these marine aerosols on polar climate.

Metin Baykara

Metin Baykara is a researcher at the Istanbul Technical University, Eurasia Institute of Earth Sciences, Department of Climate and Marine Sciences. He is also part of the University of Helsinki, Institute for Atmospheric and Earth System Research, Multiscale Modeling group.

Michael Boy

Michael Boy is a Professor of Atmospheric Science and Applied Mathematics. The position is shared between the University of Helsinki, INAR and the Lappeenranta University of Technology, Department of Computational Engineering. His main scientific interests are in atmospheric chemistry and aerosol dynamics and how to apply machine learning and artificial intelligence techniques in these areas.