Constraining the budget of NO_x and volatile organic compounds at a remote tropical island using multi-platform observations and WRF-Chem model simulations

Volatile organic compounds (VOCs) act as precursors to ozone and secondary organic aerosols, which have significant health and environmental impacts. They can also reduce the atmospheric oxidative capacity. However, their budget remains poorly quantified, especially over remote areas such as the tropical oceans. Here, we present high-resolution simulations of atmospheric composition over Réunion Island, located in the Indian Ocean, using the Weather Research and Forecasting model coupled with Chemistry (WRF-Chem). The coexistence and spatial heterogeneity of anthropogenic and biogenic emission sources in this region present a valuable but challenging test of the model performance. The WRF-Chem model is evaluated against several observational datasets, including proton transfer reaction mass spectrometry (PTR-MS) measurements of VOCs and oxygenated VOCs (OVOCs) at the Maïdo Observatory, Réunion Island (2160 m above sea level), in January and July 2019, representing austral summer and winter, respectively, and capturing the seasonal extremes for the region. While the primary goal of our study is to gain a better understanding of the (O)VOC budget at remote tropical latitudes, important model refinements have been made to improve the model performance, including the implementation of high-resolution anthropogenic and biogenic isoprene emissions, updates to the chemical mechanism, and adjustments to the boundary conditions. These refinements are supported by comparisons with PTR-MS data as well as with meteorological measurements at Maïdo; in situ NO_x and O₃ measurements from the air quality Atmo-Réunion network; Fourier transform infrared spectroscopy (FTIR) measurements of O₃, CO, ethane, and several OVOCs, also at Maïdo; and satellite retrievals from the TROPOspheric Monitoring Instrument (TROPOMI).

TROPOMI NO₂ data suggest that anthropogenic emissions, particularly from power plants near Le Port, dominate NO_x levels over the island. Both TROPOMI and in situ surface NO₂ comparisons are used to adjust the power plant emissions at Le Port. Surface ozone concentrations are overestimated by ∼6 ppbv on average, likely due to the neglect of halogen chemistry in the model, though other factors may also contribute. While modelled NO₂ over oceans is too low in summer when the lightning source is excluded, including this source results in model overestimations, as corroborated by comparisons with upper tropospheric NO₂ mixing ratios derived from TROPOMI using the cloud-slicing technique (Marais et al., 2021). The model generally succeeds in reproducing the PTR-MS isoprene and its oxidation products (Iox), except for a moderate underestimation (∼30 %) of noontime isoprene concentrations, and modelled concentration peaks near dawn and dusk, which are not seen in the observations. The ratio of Iox to isoprene (0.8 at noon in January) is fairly well reproduced by the model. The methanol and monoterpenes observations both suggest overestimations of their biogenic emissions, by factors of about 2 and 5, respectively. Acetaldehyde anthropogenic emissions are likely strongly overestimated, due to the lumping of higher aldehydes into this compound. Without this lumping, the modelled acetaldehyde would be underestimated by almost one order of magnitude, suggesting the existence of a large missing source, likely photochemical. The comparisons suggest the existence of a biogenic source of methyl ethyl ketone (MEK), equivalent to about 3 % of isoprene emissions, likely associated with the dry deposition and conversion of key isoprene oxidation products to MEK. A strong model underestimation of the PTR-MS signal at mass 61 is also found, by a factor of 3–5 during daytime, consistent with previously reported missing sources of acetic and peracetic acid.

1 Introduction

Volatile organic compounds (VOCs) play a key role in important atmospheric chemical processes. They are precursors of ozone and secondary organic aerosols, both of which are responsible for negative effects on human health (Jerrett et al., 2009; Pye et al., 2021) and the environment (e.g. Sicard et al., 2017), while also having climatic consequences (Mickley et al., 2004; Kanakidou et al., 2005; Shindell et al., 2006; Shrivastava et al., 2017). VOCs are mostly emitted by biogenic, anthropogenic, and pyrogenic sources. Biogenic VOC (BVOC) emissions due to terrestrial vegetation are affected by meteorological conditions (Guenther et al., 2006). VOCs undergo oxidation mainly by reaction with the hydroxyl radical (OH). OH is responsible for the removal of numerous pollutants in the atmosphere, initiating oxidation processes that usually transform airborne species into more oxygenated and, therefore, more water soluble compounds (Comes, 1994). The reaction between VOCs and OH can effectively deplete OH and diminish the oxidative capacity of the atmosphere, thereby increasing the lifetime of pollutants and greenhouse gases such as methane (Zhao et al., 2019). Following the initial reaction of a VOC with OH, the further degradation of oxidation products can lead to additional OH loss, especially in remote regions with low levels of NO_x (NO_x = NO + NO₂) (Di Carlo et al., 2004; Read et al., 2012; Travis et al., 2020), where reaction with NO is not the dominant sink for peroxy radicals (Logan, 1985; Atkinson, 2000).

Field measurements of OH reactivity in the remote troposphere reveal the existence of a “missing” OH sink, primarily attributable to unknown organic compounds, in particular over boreal (Sinha et al., 2010) and tropical forests (Nölscher et al., 2016; Pfannerstill et al., 2021), as well as over the remote marine boundary layer (MBL) (Thames et al., 2020). In addition, models underestimate the concentrations of several known oxygenated VOCs (OVOCs) that contribute to OH reactivity over the remote ocean, the most important being acetaldehyde (Travis et al., 2020). The high observed abundances of acetaldehyde were recently supported by unexpectedly high concentrations of peroxyacetic acid (PAA) measured over the remote ocean (Wang et al., 2019), since PAA is photochemically produced almost exclusively by acetaldehyde oxidation under low NO_x conditions.  The underestimation of acetaldehyde in models (e.g. Millet et al., 2010; Read et al., 2012; Travis et al., 2020) indicates a missing source, likely due to air–sea exchange and secondary photochemical formation from oceanic precursors (Singh et al., 2003; Millet et al., 2010; Wang et al., 2019). Acetone is a substantial source of HO_x radicals (HO_x = OH + HO₂) in the upper troposphere and lower stratosphere (UT/LS) (Müller and Brasseur, 1999; Wang et al., 2020). Like acetaldehyde, it is a compound strongly regulated by the ocean, but which appears to be better simulated by models (Wang et al., 2020), although observations suggest uncertainties in its sea–air exchanges, continental emissions, and photochemical production and photodissociation rates (Fischer et al., 2012; Wang et al., 2020). Formic and acetic acid make up more than half of rain acidity in the remote atmosphere (Keene and Galloway, 1984), but their global budget remains poorly understood, and the significant discrepancies between modelled and measured distributions indicate large missing sources of both species in both polluted and remote regions (e.g. Paulot et al., 2011; Stavrakou et al., 2012; Millet et al., 2015; Khan et al., 2018). Although methanol is the most abundant tropospheric non-methane compound, there are still uncertainties in its source apportionment (e.g. Jacob et al., 2005). Biogenic emissions are the largest source of methanol over continental areas (Stavrakou et al., 2011), and secondary photochemical sources appear to be the main contributor to methanol abundances over remote oceanic areas (Bates et al., 2021). Finally, although methyl ethyl ketone (MEK) is much less abundant than acetone, it is also more reactive (Brewer et al., 2020). Besides diverse sources, including photochemical production and anthropogenic emissions, there is evidence of significant biogenic and oceanic sources of MEK that can affect its overall abundance, especially in more remote locations (Yáñez-Serrano et al., 2016; Brewer et al., 2020). Recently, the deposition of isoprene oxidation products on vegetation, and the subsequent conversion to MEK and other products, has been proposed to be the largest contributor to the MEK budget at the global scale (Canaval et al., 2020), although this estimate relies on relatively few measurements.

In this study, we compare multiple chemical observational datasets from the remote Réunion Island, in the southern Indian Ocean, against regional atmospheric composition simulations using the high-resolution Weather Research and Forecasting model coupled with chemistry (WRF-Chem). More specifically, we make use of proton transfer reaction – mass spectrometry (PTR-MS) measurements of VOC and OVOC concentrations (Verreyken et al., 2021), performed at the high-altitude site of Maïdo Observatory (21.1° S, 55.4° E, 2160 m a.s.l.). Despite its small size, Réunion Island has diverse emission sources and is influenced by both oceanic and continental emissions (Baray et al., 2013). For these reasons, Réunion Island is an area of continuous interest, where long-term observations as well as dedicated measurement campaigns are conducted and used to validate large-scale models (e.g. Vigouroux et al., 2012; Callewaert et al., 2022) and analyse factors influencing local atmospheric composition (e.g. Rocco et al., 2020, 2022; Verreyken et al., 2020, 2021; Dominutti et al., 2022; Duflot et al., 2022). Due to these multiple influences and the pronounced orography of the island, atmospheric chemistry modelling over Réunion is especially challenging. In particular, the significant role and strong spatial heterogeneity of anthropogenic and biogenic emissions (Verreyken et al., 2022) need to be considered. Equally, the steep topography of the island makes high-resolution data indispensable in forecasting local circulation patterns (El Gdachi et al., 2024). In this work, the WRF-Chem model is enhanced with high-resolution (1 km²) datasets for anthropogenic and biogenic isoprene emissions. The model is evaluated and further refined based on comparisons with meteorological and air quality in situ data, with Fourier transform infrared spectroscopy (FTIR) column measurements at Maïdo and spaceborne TROPOspheric Monitoring Instrument (TROPOMI) observations. Given the importance of meteorology for simulating transport, biogenic emissions, and photochemistry, the model is evaluated against meteorological observations at Maïdo. FTIR and PTR-MS observations of long-lived compounds are particularly useful to evaluate the background atmospheric composition and constrain the lateral boundary conditions of the regional model. TROPOMI column observations and air quality measurements are essential to test and constrain the emissions, particularly NO_x. The PTR-MS dataset at Maïdo is expected not only to validate the model and better constrain local emissions of important (O)VOCs, but it is also expected to help identify potential shortcomings that should be considered in future studies of atmospheric composition in environments similar to Réunion Island.

The set-up and configuration of the WRF-Chem model, the chemical mechanism, the initial and boundary conditions (ICBCs), and the emissions considered in the simulations presented in this study are described in Sect. 2.1–2.3. Section 2.4 presents the observational datasets used to evaluate the model, including meteorological observations, surface chemical concentration data (Sect. 2.4.1), the PTR-MS dataset of VOC and OVOC concentrations at Maïdo (Sect. 2.4.2), the FTIR column measurements at Maïdo (Sect. 2.4.3), and, finally, the spaceborne columns from TROPOMI (Sect. 2.4.4). Section 3 evaluates the model performance relative to the various measurement datasets. The results are recapitulated before the concluding remarks in Sect. 4. Complementary figures and statistics can be found in the Supplement.

2 Methodology

2.1 Simulation area

Réunion Island is a French overseas island located in the Indian Ocean, ∼700 km off the coast of Madagascar (Fig. 1). The island covers an area of 2512 km² (63 km long and 45 km wide), roughly spanning from −21.39° S to −20.87° S and 55.22° E to 55.84° E. Réunion is a volcanic island with mountainous orography, with a maximum altitude of 3070 m above sea level (Piton de la Fournaise). The region is largely covered by native vegetation (100 000 ha) and is mostly free of strong anthropogenic emission sources (Duflot et al., 2019). The main anthropogenic emission source is fossil-fuel combustion in the industrial and transport sectors (responsible for 87.1 % of energy generation, Praene et al., 2012). The three principal industrial emission hotspots are due to biomass power plants (Le Gol and Bois-Rouge) and a diesel power plant (Le Port) (Fig. 2).

Figure 1 Map of Africa and surrounding Indian Ocean indicating two model domains in blue and red (with 12.5 and 2.5 km horizontal resolution, respectively).

Réunion Island has a large concentration of endemic species (Myers et al., 2000). The island is largely dominated by vegetation, and the plant species distribution changes with altitude (Fig. 2 in Strasberg et al., 2005; Foucart et al., 2018; Duflot et al., 2019). Over 80 % of the human population is concentrated in the coastal regions.

2.2 WRF-Chem

WRF-Chem is used to calculate meteorological and chemical atmospheric processes (WRF-Chem; Grell et al., 2005). Model version 4.1.2 of WRF-Chem was used in conjunction with its preprocessing system, WPS (WRF Preprocessing System) version 4.1.

2.2.1 Model configuration

To minimize computational demand, the simulations were conducted at 12.5 and 2.5 km horizontal resolution in the parent and nested domains, denoted d01 and d02, respectively (Fig. 1). Note that a 2 km horizontal resolution was found appropriate for simulating FTIR and in situ observations in a model study of greenhouse gases over Réunion, also using WRF-Chem but with the chemistry turned off (Callewaert et al., 2022). The projection is set to Mercator, as is the recommended set-up for low-latitude simulations close to the Equator. For January and July 2019 30 d simulations were conducted, each starting on the first day of the month at 00:00 UT. The local time at Réunion is UT +4 h. January and July correspond to summer and winter in the southern hemisphere, although the difference in meteorology between the two seasons is relatively small due to the tropical climate of the island. January and July 2019 were relatively free of events interfering with data collection (weather, volcanic activity, maintenance, etc.; see Table 1 in Verreyken et al., 2021). Lightning NO_x emissions are ignored in the reference model simulation, but they are included in a sensitivity run, using the updated Price and Rind parameterization scheme based on cloud-top height (Price and Rind, 1992; Wong et al., 2013). Since Barten et al. (2020) note that the standard WRF-Chem settings lead to a large overestimation of lightning emissions, we also downscale the number of flashes (adopted flash rate factor of 0.1 and 0.02 for d01 and d02, respectively, to account for the difference in resolution), and the production of NO per lightning strike is set to 250 moles, instead of 500 in the standard setting (Barten et al., 2020).

Simulations were conducted using Silicon Graphics (SGI) high-performance computing (HPC), equipped with an Intel processor, using 72 cores (accounting for 5 h wall time for a 2 d simulation in a two-domain set-up). The simulations were conducted sequentially in 2 d intervals, whereby the initial chemical conditions of each run were obtained from the previous run, except in the case of the first day of the month. The meteorology was re-initialized at the start of each 2 d simulation. The physical parameterizations incorporated in the simulations are listed in Table 1.

Table 1 List of WRF-Chem physical parameterizations adopted for this study.

Figure 2 Map of Réunion Island with population density at 500 m × 500 m resolution. Key locations are denoted on the plot; cities (with 50,000+ population, based on the 2020 Population Census, https://www.insee.fr/fr/statistiques/fichier/4265439/dep974.pdf, last access: 22 October 2024) and the principal power plants. The blue triangle, representing the Maïdo Observatory, is the location of the PTR-MS and FTIR instruments. The locations of air quality monitoring sites are represented by red circles.

The land-use dataset utilized in WRF is provided by MODIS 21 class (Friedl et al., 2002; Hulley et al., 2016) and has a resolution of 30 s (roughly 1 km). The distribution of the MODIS land-use index at the resolution of the nested domain is displayed in Fig. 3.

The dominant surface type at Maïdo is woody savannas (8), which is characterized by a tree cover between 30 % and 60 % and a tree canopy higher than 2 m. There is spatial variability in the pixels adjacent to Maïdo, ranging across evergreen broadleaf forests (2), savannas (9), and grasslands (10), each with different specifications regarding major vegetation type and tree cover percentage. Note that the Model of Emissions of Gases and Aerosols from Nature (MEGAN) v2.04 module utilized to calculate biogenic emissions (Guenther et al., 2006) in WRF-Chem recognizes only four vegetation types: broadleaf trees (BT), needle leaf trees (NT), shrubs and brush (SB), and herbs, crops, and grasses (HB). Their spatial distribution is provided by an independent dataset developed specifically for Réunion Island (see further below).

Figure 3 Land use index used for defining surface features in WRF, following the classification of the MODIS Land Cover Type Product. The resolution shown is 2.5×2.5 km² (d02).

2.2.2 Chemical mechanism

The reference gas-phase mechanism used in the model simulations is the Model for Ozone and Related chemical Tracers, version 4 (MOZART-4) mechanism (Emmons et al., 2010), with the Kinetic PreProcessor (KPP) (Damian et al., 2002). Several updates to the mechanism were tested and implemented, as described further below. Aerosol chemistry is simulated using the Global Ozone Chemistry Aerosol Radiation and Transport model (GOCART) (Chin et al., 2002). The MOZART-4 mechanism was chosen in lieu of its more updated counterpart, MOZART-T1 (Emmons et al., 2020), for computational reasons, despite the improvements in the new version, namely in isoprene, monoterpene, and aromatic chemistry. Instead, the isoprene chemistry of the MOZART-4 mechanism used in the model has been updated to reflect recent mechanistic updates in OH recycling and formation of methyl vinyl ketone (MVK) and methacrolein (MACR) (Sect. 2.2.3).

2.2.3 Mechanistic updates

Isoprene

The MOZART-4 isoprene mechanism lacks the OH recycling mechanisms that were shown to occur through the reactions of isoprene hydroxy hydroperoxides (ISOPOOH) (Paulot et al., 2009) and the unimolecular reactions of isoprene peroxy radicals (ISOPO2) (Peeters et al., 2009, 2014; Wennberg et al., 2018; Müller et al., 2019). The reactions of ISOPO2 and ISOPOOH in MOZART-4 are amended based on the previously cited studies and adjusted to match the results of a more up-to-date mechanism (MAGRITTEv1.1, Müller et al., 2019) through box model comparisons at both high NO_x (1 ppb NO_x) and low NO_x (0.1 ppb NO_x). More specifically, as seen in Table 2, the MVK and MACR yields in the ISOPO2 reaction with NO are increased at the expense of the hydroxy aldehydes (HYDRALD), and OH is formed in reactions that play a dominant role at low NO_x (ISOPO2 + HO₂ and ISOPOOH + OH), in order to mimic OH recycling processes involving compounds that are missing from the mechanism. Note that the mechanistic changes were chosen to avoid the introduction of additional species to maintain a similar overall computational cost.

MEK oxidation

The mechanism of methyl ethyl ketone (MEK) is revised to avoid the overestimation of the photochemical production of acetaldehyde through the oxidation of MEK and higher alkanes, which are MEK precursors. The oxidation of MEK by OH forms an intermediate radical product (MEKO2 in MOZART-4), which leads to the formation of acetaldehyde in the presence of NO, with a unit yield in MOZART-4 (reaction MO1 in Table 2). MEKO2 actually represents a lumping of three isomers, denoted MEKAO2, MEKBO2, and MEKCO2, in the comprehensive Master Chemical Mechanism (MCM) v3.3.1 (Saunders et al., 2003; https://mcm.york.ac.uk/MCM/, last access: 1 March 2025). Only one of the three radicals leads to acetaldehyde in the presence of NO according to MCM, and the overall chemistry can be summed up as follows:

$$\begin{array}{}\text{(R1)}& {\displaystyle }\begin{array}{rl}\mathrm{MEK}& +\mathrm{OH}\to \mathrm{0.459}\mathrm{MEKAO}\mathrm{2}+\mathrm{0.462}\mathrm{MEKBO}\mathrm{2}\\ & +\mathrm{0.079}\mathrm{MEKCO}\mathrm{2}\end{array}\text{(R2)}& {\displaystyle }\begin{array}{rl}\mathrm{MEKAO}\mathrm{2}& +\mathrm{NO}\to \mathrm{0.24}{\mathrm{RO}}_{\mathrm{2}}+\mathrm{0.24}{\mathrm{HO}}_{\mathrm{2}}\\ & +\mathrm{0.76}\mathrm{HCHO}+\mathrm{0.76}{\mathrm{EO}}_{\mathrm{2}}+{\mathrm{NO}}_{\mathrm{2}}\end{array}\text{(R3)}& {\displaystyle }\mathrm{MEKBO}\mathrm{2}+\mathrm{NO}\to {\mathrm{CH}}_{\mathrm{3}}{\mathrm{CO}}_{\mathrm{3}}+{\mathrm{CH}}_{\mathrm{3}}\mathrm{CHO}+{\mathrm{NO}}_{\mathrm{2}}\text{(R4)}& {\displaystyle }\mathrm{MEKCO}\mathrm{2}+\mathrm{NO}\to \mathrm{HCHO}+{\mathrm{C}}_{\mathrm{2}}{\mathrm{H}}_{\mathrm{5}}{\mathrm{O}}_{\mathrm{2}}+{\mathrm{NO}}_{\mathrm{2}},\end{array}$$

where RO₂ denotes the acetonyl peroxy radical, CH₃C(O)CH₂O₂, and EO₂ denotes HOCH₂CH₂O₂. Note that intermediary MCM reactions leading to these final products were skipped, as they involved compounds that are not defined in the MOZART-4 mechanism. These compounds were assumed to react rapidly according to their major sink reaction in the presence of NO. Combining the previous equations leads to the updated MO1 reaction shown in Table 2. The oxidation of MEK at low NO is oversimplified in the MOZART-4 mechanism and leads to the same final products as at high NO. However, the reactions of the peroxy radicals, MEKAO2, MEKBO2, and MEKCO2, with HO₂ form ketohydroperoxides that are expected to photolyse rapidly, leading in part to very different products (enols) than in the high NO case (Liu et al., 2018).

Table 2 List of mechanistic changes adopted in the MOZART-4 mechanism.

2.2.4 Initial and boundary conditions

The model is initialized at the start of each run using input data from the Community Atmosphere Model with Chemistry (CAM-Chem; Lamarque et al., 2012; Emmons et al., 2020), which is a global model utilizing the MOZART-4 mechanism and providing output at 0.9° × 1.25° resolution. This allows for a direct representation of the initial and boundary conditions (ICBCs) in the run. For species where global Copernicus Atmosphere Monitoring Service (CAMS) reanalysis data (Inness et al., 2019) are available (NO_x, CO, O₃, HNO₃, H₂O₂, HCHO, peroxyacetyl nitrate or PAN, C₂H₆ and C₃H₈), we replace the CAM-Chem initial conditions with these higher resolution data (0.75° × 0.75°). Both inputs were taken at 6 h intervals. Adjustments to these ICBCs were made based on comparisons with ground-based and FTIR measurements, which are detailed in Sect. 3.4. More specifically, methanol and PAN concentrations were multiplied by 0.6. MEK and C_>3 alkanes (represented as the lumped compound BIGALK) were multiplied by 0.4 in both seasons, while the same factor was applied to acetone only in July. Ethane was increased by a factor of 1.6 in July and left unchanged in January.

2.3 Emissions

The source apportionment of the main gaseous pollutants and VOC precursors is summarized in Fig. 4, which provides the total emissions for every species over the island (average of January and July). In addition, Fig. 5 displays the emission distribution of the four most emitted VOCs over the island, namely isoprene, methanol, monoterpenes, and acetaldehyde, for the month of January. The assumptions and datasets used to derive these emissions are given in the following subsections.

2.3.1 Anthropogenic emissions

The horizontal resolution of global anthropogenic emissions inventories, such as the Emissions Database for Global Atmospheric Research (EDGAR) (0.1° × 0.1°) (Monforti Ferrario et al., 2022), is too coarse for air quality simulations over Réunion Island. The low-resolution dataset cannot accurately represent the transition from highly polluted areas (e.g. cities) to remote zones such as Maïdo. We utilized a 1×1 km² emissions inventory estimated for NO_x, SO₂, CO, and non-methane VOCs (NMVOCs) based on data from the local air quality agency (Atmo-Réunion; https://atmo-reunion.net/le-dispositif-de-surveillance, last access: 1 December 2023). The NMVOC species used in this study are listed in Table 3 below. The NMVOC speciation follows the Regional Atmospheric Chemistry Mechanism, Version 2 (RACM2) (Goliff et al., 2013). Traffic emissions were provided by Atmo-Réunion for the main roads of the island. Industrial emissions, mostly from cane sugar refining, rum distillation, and diesel-electric power production, were estimated as point sources. Agricultural emissions are distributed in cultivation areas.

The Atmo-Réunion emissions are complemented by the remaining sectors (namely shipping, residential burning, non-road transport, solvent processing, and waste management) from the EDGAR inventory, as well as by species not included in the high-resolution inventory (ammonia, black and organic carbon, and particulate matter). EDGARv6.1 (Monforti Ferrario et al., 2022) was used for the trace gases and particulate matter, and EDGAR Hemispheric Transport of Air Pollution (EDGAR-HTAP) v4.3.2 (Crippa et al., 2018; Huang et al., 2017) was used for the NMVOCs. These inventories were also used for all species and sectors in the rest of the model domain. Emissions from EDGARv6.1 reflect values from 2018, while emissions from EDGAR-HTAP v4.3.2 are based on values from 2012.

Temporal variations in the emissions are applied in accordance with Poraicu et al. (2023), based on EDGAR-specific temporal profiles (Crippa et al., 2020). As Réunion is a French territory, the temporal profile follows French specifications, which are defined based on French mainland regions. The temporal variations for some activities might not be representative of Réunion Island. For example, the seasonal temperature variations in mainland France are higher than in the tropics, with cold winters and hot summers, which determine different residential heating behaviours. The seasonal dependence of residential sector emissions was therefore omitted.

Table 3 MOZART-4 VOC precursors and correspondence with the species classes of the Atmo-Réunion and EDGAR-HTAP V4.3.2 inventories.

Figure 4 Percentage contributions of different emission sources (anthropogenic and biogenic) to the total emissions of different chemical compounds considered in this work. The total, shown on the right-hand y-axis, is calculated from the average of the two months. Results shown for the optimal run (R0), after adjustment of the anthropogenic and biogenic sources described in text (Sects. 2.2, 2.3 and 3).

Figure 5 Total emissions of isoprene, methanol, monoterpenes and acetaldehyde over the inner model domain at 2.5 km resolution (d02), averaged over the month of January 2019. The white circle indicates the location of the Maïdo Observatory.

Anthropogenic emission adjustments

Given the coarse resolution (0.1°) of the EDGAR inventory used for residential emissions, these emissions were redistributed spatially on the model grid using a distribution of population density specific to Réunion (https://public.opendatasoft.com/explore/dataset/population-francaise-par-departement-2018/table/?disjunctive.departement, last access: 8 December 2023). The population distribution is on a latitude–longitude grid at a resolution of 30 m and was regridded to the model resolution in the nested domain. In addition, since traffic activity occurs not only on highways but also in cities and generally where people live, the traffic sector was redistributed by assuming that 70 % of the emissions follow the population map, and 30 % of the emissions are distributed over the principal highways defined in the Atmos-Réunion inventory. The resulting high-resolution distribution is compared with EDGAR on Fig. 6.

The NO_x industrial emissions (primarily energy production) were redistributed based on the EDGARv6.1 emissions for this species and sector. Both the EDGARv6.1 dataset and TROPOMI NO₂ column data (see Sect. 3.5) suggest maximum emissions in the northwestern region of the island, around the Port-Est power plant. Therefore, the total industrial NO_x emissions were distributed amongst the 4 point source regions using percentage contributions based on the NO_x industrial EDGARv6.1 emissions dataset. EDGAR-HTAP v4.3.2 reports total VOC emissions that are three times higher than those in Atmo-Réunion, implying a significantly lower $\mathrm{VOC}/{\mathrm{NO}}_{\mathit{\text{x}}}$ emission ratio in the latter inventory. Although both inventories have similar spatial distributions, the Atmo-Réunion emissions are higher around the largest city (Saint-Denis) and lower in the other industrialized areas.

Figure 6 Anthropogenic emissions of (a–b) NO_x and (c–d) total VOC over the inner model domain at 2.5 km resolution (d02), averaged over the month of January 2019. (a) and (c) display low-resolution EDGAR emissions, while (b) and (d) show the high-resolution Atmo-Réunion dataset, with adjustments described in Sect. 2.3.1.

2.3.2 Biogenic emissions

Emissions from biogenic sources, mainly from vegetation, are calculated using MEGAN version 2.04. The emissions are calculated online (at the same time step as the model), based on the simulated meteorological fields and vegetation types defined by the land-use map. This model utilizes four general vegetation types (broadleaf trees, needleleaf trees, shrubs, and herbaceous land) together with standard emission factors and other parameters needed to estimate the emissions at each time step. MEGAN calculates emissions for 20 compounds/compound classes, based on which the emissions of 138 individual species can be estimated. The emissions of WRF-Chem species classes (e.g. monoterpenes) are obtained by lumping those individual compounds. The biogenic emission for each class (E) is calculated using

$$\begin{array}{}\text{(1)}& E=\mathit{\epsilon }\mathit{\gamma }\mathit{\rho };\end{array}$$

i.e., E is calculated as a product of the emission factor at standard conditions (ε), the emission activity factor (γ), and a factor accounting for production and loss within the plant canopy (ρ). The latter factor is not considered here, i.e. ρ=1. The standard emission factor (in µg m⁻² h⁻¹) is obtained from independent field and laboratory studies.

The dimensionless emission activity factor γ accounts for the dependence on environmental conditions such as the leaf area index (LAI); photosynthetic photon flux density, i.e. the amount of visible light at leaf level (γ_P); leaf temperature (γ_T); leaf age (γ_A); soil moisture (γ_SM); and CO₂ inhibition (γ_C). For more detail on the MEGAN algorithm, see Guenther et al. (2006, 2012).

MEGAN relies on a 1 km resolution global map of plant functional types (PFTs) and LAI. Since the standard database provided by the National Center for Atmospheric Research (NCAR) (https://www.acom.ucar.edu/wrf-chem/download.shtml, last access: 1 December 2023) has zero values for all these variables over Réunion, we replaced it with a detailed cartography of plant functional types and isoprene standard emission factors (100 m resolution) constructed from an extensive survey of natural habitats in Réunion Island (DEAL Réunion, 2025) and phytosociology studies (Strasberg et al., 2005). This is completed by LAI measurements for representative plant species (Duflot et al., 2019). The resulting distributions of plant functional types, the standard isoprene emission factor, and LAI are shown in Fig. 7 at the model resolution (inner domain). As seen in this figure, broadleaf trees are by far the dominant PFT over the island. On average, over the island (except bare soil), the isoprene emission factor is about 2900 µg m⁻² h⁻¹ (42.7 mol km⁻² h⁻¹). Around Maïdo Observatory, the emission factor ranges between 3000 and 6000 µg m⁻² h⁻¹. These values are much lower than the default values in the WRF-Chem model for the dominant PFTs (13 000 and 11 000 µg m⁻² h⁻¹ for broadleaf trees and for shrubs, respectively). This illustrates the importance of accounting for the local plant species distribution, when available.

Figure 7 (a–d) Plant functional type coverage (%) used to calculate biogenic emissions with MEGAN in WRF-Chem, (e) isoprene standard emission factor (mol km⁻² h⁻¹) and (f) leaf area index at 2.5 km horizontal resolution. The white circle represents the location of the Maïdo Observatory.

Except for isoprene, the emission parameters for each compound (or compound class) are assumed to be constant within each of the 4 vegetation types, despite possibly important variations. In this work, the MEGAN parameters for several compounds (methanol, monoterpenes, methyl ethyl ketone, and acetaldehyde/ethanol) were updated based on previous work and comparisons with the PTR-MS measurements at Maïdo, which will be further discussed in Sect. 3.3. More specifically, the methanol emission factor was decreased from 800 to 400 µg m⁻² h⁻¹ for broadleaf trees and herbaceous vegetation, in agreement with Stavrakou et al. (2011). The emission factors for monoterpenes were decreased by a factor of 5, and the light-dependent fraction (LDF) was increased to 0.9 (from initial values ranging between 0.05 and 0.8 for different monoterpenes). This high LDF value implies an emission algorithm similar to isoprene, with little emission at night. It might reflect the dominance of broadleaf forests and the marginal extent of coniferous trees in the island (Fig. 7), suggesting a minor role for temperature-controlled VOC emissions from storage pools (see Derstroff et al. 2017; Yáñez-Serrano et al., 2015). The MEGAN2.1 emissions of MEK were enhanced, since recent studies show an important biogenic contribution to its budget (Yáñez-Serrano et al., 2016). Since this biogenic source is thought to result from the partial conversion of isoprene oxidation products (MVK and ISOPOOH) deposited on plant leaves (Canaval et al., 2020), the biogenic source of MEK is assumed proportional to that of isoprene. Although this assumption might underestimate MEK emissions at night, since non-zero concentrations of isoprene oxidation products might be sustained during night-time (e.g. Langford et al., 2010), the error is likely small, since the dry deposition of those compounds is usually very slow at night (Nguyen et al., 2015). The scaling factor used for calculating the emissions has been adjusted as discussed in Sect. 3.3.

2.3.3 Biomass burning emissions

The fire emissions from the Fire INventory from NCAR (FINN) emissions inventory (Wiedinmyer et al., 2011) were zero over both domains in January and July 2019. FINN is a global, 1 km resolution dataset that provides daily estimates of biomass burning emissions. However, as with the biogenic emissions input files, the biomass burning emissions could be artificially missing in this inventory due to the small size of the island. In any case, there was no significant biomass burning activity detected at Maïdo in January and July 2019 (Verreyken et al., 2020). The fire emissions were therefore omitted in the simulations.

2.4 Measurements

2.4.1 In situ surface chemical observations

Common air pollutants are measured hourly at 18 air quality stations around Réunion Island. These stations are operated by Atmo-Réunion. The stations measuring NO₂, O₃, or CO are shown on Fig. 2. In the comparison of modelled and observed NO₂ concentrations, a correction factor is applied to account for known interferences in the NO₂ measurement (Lamsal et al., 2008). These interferences are due to NO_y reservoir compounds (such as HNO₃ and PAN) contributing to the measured signal, leading to higher NO₂ values than are actually present. The correction factor is calculated using NO₂, peroxyacetyl nitrate (PAN), and HNO₃ concentrations obtained from the model, and it is applied to the simulated NO₂, as detailed in Poraicu et al. (2023).

2.4.2 PTR-MS measurements

The PTR-MS instrument is used to measure VOC species (Hewitt et al., 2002). The dataset (Verreyken et al., 2021) includes 13 molecules, or, more precisely, 13 mass-to-charge ratios ($m/z$). The instrument (hs-PTR-MS, Ionicon Analytik GmbH, Austria) is located at the Maïdo Observatory high-altitude site (2160 a.s.l.) (21.079° S, 55.384° E). Further technical information on the PTR-MS measurements can be found in Verreyken et al. (2021). The major isoprene oxidation products at high NO_x, methacrolein (MACR) and methylvinyl ketone (MVK), are measured conjointly due to their identical molar mass (70 g mol⁻¹) and molecular formula (C₄H₆O). Isoprene hydroxy hydroperoxides (ISOPOOH), which are major isoprene oxidation products at low NO_x, are known to react heterogeneously in the PTR-MS instrument and decompose to either MVK or MACR (+HCHO) (Rivera-Rios et al., 2014). Therefore, the PTR-MS signal corresponding to MVK + MACR includes a contribution of ISOPOOH, which is significant at low NO_x conditions. Although the exact conversion efficiency of ISOPOOH into MVK + MACR is uncertain, a complete conversion is assumed here, and the corresponding signal will be referred to as Iox (isoprene oxidation product, MVK + MACR + ISOPOOH).

The PTR-MS signal for acetic acid (CH₃COOH) has interferences from other chemical species, namely glycolaldehyde, peroxyacetic acid, propanols, and ethyl acetate (Baasandorj et al., 2015). The contribution of propanol is very low (< 1 %) and is not considered. Protonated glycolaldehyde has the same mass-to-charge ratio as protonated acetic acid (61 $m/z$), while peroxyacetic acid and ethyl acetate were shown to fragment upon protonation with a resultant fragment ion (61 $m/z$) (Španel et al., 2003; Fortner et al., 2009). We therefore opt to compare the observations with the sum of the modelled concentrations of acetic acid, glycolaldehyde, and peroxyacetic acid, assuming a similar sensitivity for all three compounds. Ethyl acetate is not considered in the MOZART-4 mechanism, although this compound has both direct emissions and formation from the chemical oxidation of ethers (Orlando and Tyndall, 2010, and references within). Formic acid (46 g mol⁻¹), also measured by the PTR-MS, is not defined in the MOZART-4 mechanism, and it is thus not considered further in this study. The signal at $m/z$ 73 has a major contribution from MEK (Verreyken et al., 2021), but other compounds, namely butanal isomers and methylglyoxal (MGLY), also contribute (Yáñez-Serrano et al., 2016). The contribution of butanal is neglected here, as it is not calculated by the model. As the instrument's estimated sensitivity to MGLY is only a fraction (∼0.7) of the sensitivity to MEK (Koss et al., 2018), the PTR-MS concentrations will be compared to the sum MEK + 0.7× MGLY calculated using the model.

Benzene, toluene, and xylenes are all measured at the site. However, in the MOZART-4 mechanism, they are lumped and treated as toluene. As benzene, toluene, and xylenes have different chemical lifetimes (Atkinson, 2000) and different emission distributions, the PTR-MS observations of aromatics cannot be evaluated with the model. In situ meteorological data (2 m temperature, relative humidity, wind speed and direction, and solar radiation) are measured at and provided for the PTR-MS site (Maïdo Observatory) at the same temporal resolution as the chemical concentrations (2.7 min).

2.4.3 FTIR measurements

Ground-based FTIR measurements have been performed at Réunion Island (Maïdo and Saint-Denis) since 2002, initially on a campaign basis (Senten et al., 2008). Long-term measurements of VOCs and other compounds started in 2009 at Saint-Denis (Vigouroux et al., 2012). Since 2011, FTIR measurements in the prescribed Network for the Detection of Atmospheric Change (NDACC) spectral range (600–4500 cm⁻¹) are performed at Maïdo Observatory (Baray et al. 2013) using a Bruker 125 HR (since 2013).

As FTIR is a remote sensing technique that measures the absorption of solar light by atmospheric species along the line-of-sight (instrument–sun), the primary product of the FTIR retrievals is the total column of the absorbing gases. In addition, low resolution vertical profile information can also be derived based on the pressure and temperature dependence of the line shape. The choice of spectral windows and spectroscopic parameters is optimized for each target species, preferably within the whole NDACC FTIR community to ensure consistency within the network. Table S1 in the Supplement summarizes the main retrieval settings for the species used in this paper.

The HCHO retrieval has been optimized recently at more than 20 FTIR stations (Vigouroux et al., 2018), and it is now an official NDACC target species. Methanol is not an official NDACC species, and it is measured at only a few sites (see e.g. Wells et al., 2024) but not in a harmonized way. Ethane is an official NDACC species and has been harmonized within the network as described in Franco et al. (2015). At Maïdo, the humidity being important, the third window suggested for harmonization in Franco et al. (2015) is not used to avoid strong interference with water vapour lines. Carbon monoxide is also an NDACC target species, and as such, the windows and spectroscopy are harmonized and found in the NDACC InfraRed Working Group (IRWG) documentation (https://www2.acom.ucar.edu/irwg, last access: 1 November 2024). The existing IRWG retrieval strategies for official species are currently being reprocessed, with the aim of using improved settings. Ozone is the FTIR species with the most degrees of freedom for signal (DOFS = 4 to 5), allowing the retrieval of several independent partial columns: one in the troposphere and three in the stratosphere (Vigouroux et al., 2008).

Peroxyacetyl nitrate (PAN) is not an official target species. It has been measured at only a few stations (Mahieu et al., 2021), because its weak spectral signature makes it difficult to detect and retrieve. This is the first time that PAN time-series data from the Maïdo station are presented. From the two stronger PAN absorptions used in Mahieu et al. (2021), only the band at 1163 cm⁻¹ is useful at Maïdo, as the other one has a too low signal-to-noise ratio. As shown in Fig. S1, the spectral signature of PAN is much weaker than that of the other absorbing gases. Therefore, the interfering species must first be carefully taken into account. In addition to the gases shown in Fig. S1, HCFC-22  also interferes. We pre-retrieve HCFC-22 in a dedicated window (828.62–829.35 cm⁻¹), and then, for each individual spectrum, the retrieved profile of HCFC-22 is used as fixed values in the PAN retrieval. For each spectrum, profile retrievals are also made for H₂O, HDO, O₃, N₂O, and CFC-12, and they are used as a priori values in the PAN retrievals, in which the PAN and H₂O profiles are retrieved, while the other species are scaled from their a priori pre-retrieved profiles (with the exception of HCHC-22, which is not scaled but fixed).

2.4.4 TROPOMI observations

TROPOMI is a spaceborne instrument aboard the European Space Agency (ESA) Sentinel-5P (S5P) satellite (Veefkind et al., 2012), which monitors the global distribution of multiple air pollutants, including nitrogen dioxide, formaldehyde, ozone, and carbon monoxide. The S5P satellite is sun-synchronous, and the retrieved data has a daily global coverage with (typically) one TROPOMI measurement at ∼ 13:30 local time at a spatial resolution of 7×3.5 km² (updated to 5.5×3.5 km² in August 2019). The main species of interest for this study are NO₂ and HCHO. The measurement of multiple species is possible due to the large spectral range of the TROPOMI spectrometer, which covers the ultraviolet (UV), visible (VIS), near-infrared (NIR), and shortwave infrared (SWIR) ranges (ranging from 267 to 2389 nm). The instrument is a push-broom spectrometer that scans the Earth while the satellite moves and measures the composition of the atmosphere using differential optical absorption spectroscopy (DOAS). The light travelling from the Sun through the atmosphere is reflected back to space, where it is measured by the spectrometer. Molecules absorb photons at well-defined windows, depending on their molecular structure (primarily 405–465 nm for NO₂ and 328–359 nm for HCHO). Comparing the observed spectrum to a reference spectrum enables the retrieval of the slant columns through a fitting procedure involving multiple compounds that are active in the absorption bands of the species of interest. For both compounds, TROPOMI retrieves a slant column density (SCD) from the Level-1b radiance and irradiance spectra that represents the total amount of compound present along the effective solar light path (van Geffen et al., 2020).

In the case of NO₂, the total SCD combines both tropospheric and stratospheric distributions. The tropospheric slant column is obtained by subtracting the stratospheric contribution from the total SCD. This contribution is obtained from a global chemical transport model (TM5-MP, Williams et al., 2017). The tropospheric vertical column density (VCD) is obtained by dividing the slant column by an air mass factor (AMF) (Palmer et al., 2001) that is dependent on the vertical profile of the considered compound. AMFs are obtained from radiative transfer calculations and can introduce a large source of uncertainty in the VCD calculation (30 %–40 %, Lorente et al., 2017), especially in the presence of clouds. The vertical distributions of NO₂ and HCHO are taken from the global chemistry transport model TM5-MP (van Geffen et al., 2022a; Williams et al., 2017; De Smedt et al., 2018) at a resolution of 1° × 1°. In the HCHO VCD algorithm, a background correction is applied on a daily basis. Since the reference spectrum is obtained from Earth radiances in the equatorial Pacific, the slant column retrieved from the fit corresponds to an excess over the remote background, where the principal HCHO source is methane oxidation. The TM5-MP columns over the same region are therefore added to the vertical columns to account for this background (De Smedt et al., 2021).

Due to its role in AMF calculations, the choice of vertical profile influences VCD retrieval, and it is necessary to consider the vertical sensitivity of the TROPOMI instrument when comparing with other datasets. NO₂ VCD biases against independent measurement datasets decreased when the TM5-MP a priori profiles were replaced with model profiles at a higher resolution (e.g Tack et al., 2021; Judd et al., 2020; Douros et al., 2023) or with measured profiles (e.g. Dimitropoulou et al., 2020). The averaging kernels of the satellite data are applied to the model profiles to calculate a “smoothed” vertical column for comparison with the satellite retrieved columns (e.g. Boersma et al., 2016). The averaging kernel represents the sensitivity of the measurement to the tracer concentrations at different altitudes, weighted by the assumed vertical profile of the tracer (Eskes and Boersma, 2003). It is provided alongside the measurement in TROPOMI products.

We present comparisons between the reprocessed version (RPRO) of the NO₂ and HCHO retrievals from TROPOMI (v3.2) and simulated tropospheric columns. For both compounds, the averaging kernels are applied to the model profiles vertically interpolated to the TM5 vertical pressure grids. Both WRF-Chem and TM5 vertical pressure levels are calculated using hybrid sigma-pressure coordinates and the surface altitude, using values from WRF-Chem and TM5, respectively. Quality filtering (QF) follows Algorithm Theoretical Basis Document (ATBD) recommendations (QF > 0.75 for NO₂ and QF > 0.5 for HCHO) (van Geffen et al., 2022b; De Smedt et al., 2022).

In addition, we applied the oversampling technique to the NO₂ and HCHO data to gain insight into the fine-resolution distribution of those compounds (e.g. de Foy et al., 2009). This technique involves the long-term averaging of TROPOMI measurements on a very fine grid, taken here to be 0.01° × 0.01° (∼ 1×1 km²). The measurement from a given TROPOMI pixel is taken to apply to a circle defined by the centre of the pixel and a radius of 3.5 km. In this way, each 0.01° × 0.01° pixel accumulates > 500 measurements over the considered time period (May 2018–July 2022). This technique takes advantage of the variable offset of TROPOMI observations from day to day and achieves a high signal-to-noise ratio at high resolution, but it is not intended for direct comparison with the model, given the long-term averaging. It aims to present the average pollutant distribution at finer scales to inform about emission hotspots that are potentially lost when considering data over short periods of time.

3 Results and discussion

The optimal model set-up for this region was obtained through multiple sensitivity runs testing various changes in the emissions, lateral boundary conditions, and chemical mechanism used in the model. In this section, the reference run (R0) represents the simulation set-up incorporating all model updates. The following sections highlight the impact of the updates through model comparisons with the observations. The list of sensitivity runs portrayed in this study is provided in Table 4.

Table 4 Sensitivity runs conducted in this study, with the shorthand notation.

3.1 Evaluation against meteorological measurements

The model is evaluated against meteorological measurements at the Maïdo site in Fig. 8. Surface temperature and solar radiation are represented by the model for both January and July. The observed diurnal profile of temperature is well reproduced, except for a low night-time bias on many instances. On most days in both months, the midday modelled value matches the observation well. Exceptions occur, e.g. during the first days of January, when the model underestimates temperature and overestimates relative humidity and cloudiness. Errors in the WRF-Chem cloud parameterization were shown to impact the model performance for many physical and chemical processes (Zhao et al., 2012; Berg et al., 2015; Ryu et al., 2018). Furthermore, cloudy periods (when observed solar radiation is lower, e.g. on 8–9 January) are typically not well represented. The higher solar radiation fluxes in January (maximum of about 1000 W m⁻² in January compared to 750 W m⁻² in July) might enhance evaporation of the ocean and other water bodies, leading to a slightly higher humidity in the warmer month, as indicated by both model and measurement data in Fig. 8. Nevertheless, seasonal variations in meteorological parameters are well reproduced by the model. January is warmer (by 5.5 K and 6.4 K according to the measurements and the model, respectively), more humid (by 14 % and 11 %) and has higher radiation fluxes (by 51 W m⁻² and 65 W m⁻²) than the month of July. On average, the wind speed is slightly overestimated in both seasons. It is generally close to the observations, except during several predicted high wind episodes (e.g. on 7–9 January and 24 July) that are not present in the observations. The predominant wind direction is correctly predicted (mostly westerly winds during January, and an alternance of westerly, northerly, and easterly winds during July), but the model often fails to reproduce short-term variability in wind direction. For example, the sporadic occurrence of easterly winds (270°) in January is missed by the model, whereas sporadic southerly winds are simulated in January and July but are seldom observed. The observed circulation likely results from the competition between overflowing trade winds and meso-scale dynamics induced by anabatic thermal flow coupled with upslope transport to the Maïdo Observatory, a direct consequence of orography (Duflot et al., 2019; El Gdachi et al., 2024). This will directly affect transport patterns and model comparisons with chemical measurements over the site. Seasonal averages are given in Table S2.

Figure 8 Evolution of observed and modelled surface meteorology at Maïdo Observatory for (a–e) January and (f–j) July. The observations are shown in black, while the cyan line represents the WRF-Chem results. The wind direction is represented in degrees, relative to north.

3.2 Evaluation against air quality station measurements

Figure 9 presents the monthly averaged distributions of key chemical compounds calculated by the model, over Réunion Island, for January and July 2019. These distributions will be useful for the discussion of model comparisons with in situ and satellite observations (see next subsections). The strong heterogeneity of anthropogenic emissions and their high intensity over urban and industrial centres lead to strong contrasts between the chemical regimes over different regions of the island. For example, the NO_x mean surface concentration ranges between 0.12 ppb at the most remote part of the island to typically 5–15 ppb over most urban/industrial areas and over 30 ppb near the Le Port thermal plant. As seen in Fig. 9, other important compounds such as OH and O₃ are strongly impacted by anthropogenic emissions. Over a large portion of the island (southwestern part as well as a small portion of the northern coast), but also over a large oceanic region surrounding almost the entire island (especially to the northwest due to the influence of Le Port emissions), OH levels are enhanced, primarily because NO_x promotes the conversion of HO₂ to OH (Spivakovsky et al., 2000) through the reaction

$$\begin{array}{}\text{(R5)}& {\mathrm{HO}}_{\mathrm{2}}+\mathrm{NO}\to \mathrm{OH}+{\mathrm{NO}}_{\mathrm{2}}\end{array}$$

As a result, OH increases as NO_x increases when NO ranges between ca. 10 ppt and ca. 500 ppt (Logan et al., 1981). At very high NO_x levels ($>\sim \mathrm{5}$ ppb NO_x), OH decreases as NO_x increases, mainly due to the sink of HO_x (= OH + HO₂) due to the radical termination reaction

$$\begin{array}{}\text{(R6)}& \mathrm{OH}+{\mathrm{NO}}_{\mathrm{2}}\to {\mathrm{HNO}}_{\mathrm{3}}\end{array}$$

which explains the depletion of OH levels over the main NO_x hotspots. NO_x has a relatively weak effect on the patterns of ozone distribution, except for the clear titration effect near Le Port and the main cities, such as Saint-Denis and Saint-Pierre. Formaldehyde is clearly strongly affected by anthropogenic activities, as seen from the high surface concentration levels around the NMVOC emission areas (Fig. 6). The distribution also partly reflects the complex impacts of NO_x on OH, since the NMVOC oxidation processes leading to HCHO formation are mainly driven by their reaction with OH. Note that the NO_x concentrations rapidly decrease with altitude. Therefore, the OH depletion effect due to NO_x seen in Fig. 9 becomes rapidly negligible at higher altitudes. Therefore, the main effect of anthropogenic NO_x emissions on the oxidative conditions above the island is a strong enhancement of OH concentrations, except over localized hotspots, very close to the surface. NO_x plays a major role in the OH budget, and inaccurate predictions of its abundances can affect the model comparisons for VOCs and their oxidation products, in particular at Maïdo. Besides the insights provided by TROPOMI NO₂ column data (see Sect. 3.5), network measurements of in situ concentrations of NO₂, NO_x (= NO + NO₂), and O₃ are available mostly in polluted areas—in cities and in the vicinity of industrial sources (Fig. 2). Due to representativeness issues, caution is required when comparing model results with measurements often obtained very close to strong pollution sources. For this reason, model underestimation of NO_x levels is to be expected at many stations. A useful indicator is the [NO_x] $/$ [NO₂] ratio. At photochemical steady state (PSS), it is given by

$$\begin{array}{}\text{(2)}& {\left({\displaystyle \frac{\left[{\mathrm{NO}}_x\right]}{\left[{\mathrm{NO}}_{\mathrm{2}}\right]}}\right)}_{\mathrm{PSS}}=\mathrm{1}+{\displaystyle \frac{J_{{\mathrm{NO}}_{\mathrm{2}}}}{k_{\mathrm{1}}\cdot \left[{\mathrm{O}}_{\mathrm{3}}\right]+\mathrm{\dots }}},\end{array}$$

where $J_{{\mathrm{NO}}_{\mathrm{2}}}$ is the photolysis rate of NO₂, and k₁ (= $\mathrm{1.95}\times {\mathrm{10}}^{-\mathrm{14}}$ molec.⁻¹ cm³ s⁻¹) is the rate of the reaction

$$\begin{array}{}\text{(R7)}& {\mathrm{O}}_{\mathrm{3}}+\mathrm{NO}\to {\mathrm{NO}}_{\mathrm{2}}+{\mathrm{O}}_{\mathrm{2}}\end{array}$$

During the night, the expected value of the ratio is unity. However, close to a pollution source, photochemical steady state is not achieved, since directly emitted NO can travel over a few hundred metres within its chemical lifetime, of the order of several minutes for typical night-time ozone levels (a few ppb; see Fig. 10). Furthermore, ozone is titrated to even lower levels in the direct vicinity of strong sources. Given the model resolution and the distribution of NO_x sources (Fig. 6), strong and systematic ozone titration occurs only in the region of Le Port (Figs. 2 and 9), with its strongly emitting power plants. Elsewhere, the night-time $\left[{\mathrm{NO}}_x\right]/\left[{\mathrm{NO}}_{\mathrm{2}}\right]$ ratio is below 1.1 in the model. In the measurements, however, the observed ratio is usually ∼1.3, based on monthly averaged night-time concentrations, and it even exceeds a factor of 2 at several stations (Table S3), namely the stations 10–12 on the southern coast of the island (see Fig. S2 for station locations). These high values suggest the presence of very close NO_x sources that are unresolved by the model.

Figure 9 Averaged modelled surface concentrations of NO_x, OH, O₃, and HCHO (ppb) in January (a–d) and July (e–h), obtained from simulation R0.

Figure 10 presents a comparison of observed and modelled surface concentrations of NO_x and O₃ for a 30-d simulation in January 2019. The stations were divided into two classes—the two stations near the Le Port power plants (LP) and the other stations, excluding the least representative stations, i.e. those for which the observed night-time $\left[{\mathrm{NO}}_x\right]/\left[{\mathrm{NO}}_{\mathrm{2}}\right]$ ratio exceeds a factor of 2 (Table S3). Figures 10 and 11 display the averaged performance for the two classes of stations. Table S4 in the Supplement provides a statistical analysis for each station using the available observations, including the correlation coefficient, root mean square error (RMSE), and index of agreement.

Figure 10 Evaluation of the average modelled surface NO₂, NO_x, and O₃ in January, from simulations S1 (cyan) and R0 (red), against in situ network observations (dashed line) (a–b) in the region of Le Port (LP) and (c–e) at the other stations, for which the observed night-time ratio $\left[{\mathrm{NO}}_x\right]/\left[{\mathrm{NO}}_{\mathrm{2}}\right]$ < 2 (stations 3–9 of Table S3; see Fig. S2). The model NO₂ concentrations have been corrected to account for interferences in the measurements using the corresponding modelled concentrations of PAN and HNO₃ (Poraicu et al., 2023).

Figure 11 Average diurnal profile of NO₂, NO_x, and O₃ in January, from simulations S1 (cyan) and R0 (red), shown for the Le Port (LP) stations and at the other stations for which the observed night-time ratio $\left[{\mathrm{NO}}_x\right]/\left[{\mathrm{NO}}_{\mathrm{2}}\right]$ < 2. Note that ozone measurements were not available at LP stations. The black dashed line represents the observations at the same stations.

The LP stations (Fig. 2, Terrain de Sel and Centre Pénitentiaire) are strongly impacted by the NO_x emissions due to the power plants located near Le Port, on the northeastern coast of the island (Fig. 9). However, the observed night-time NO_x concentrations at those stations are largely overestimated by the model. During daytime, a large overestimation is also found for simulation S1, which adopted the original estimation of the NO_x emissions due to the Le Port power plants. On average, over the 30 d period, the model overestimation for NO_x (NO₂) against the LP measurements, by a factor of 3 (4) during daytime in the S1 simulation, is strongly reduced (to a factor of 1–1.4) when the power plant emission is reduced by a factor of 5 (run R0) (Fig. 11). This finding suggests that this reduction in the power plant emission is justified, consistent with the comparison of NO₂ columns with TROPOMI data (Sect. 3.5). The large model overestimation of night-time NO₂ and NO_x at the LP stations stands in contrast with the fairly good agreement at the other stations. A possible reason might be that all emissions are injected at surface level in the model, whereas the actual injection heights of power plant emissions might be above the (usually stable) night-time boundary layer. The impact of the 5-fold reduction in NO_x emissions on the night-time NO₂ levels at the two LP stations is small (−25 %), despite their close proximity to the power plants. This is explained by the fact that NO (not NO₂) is emitted by the power plants, and high emissions lead to O₃ titration through the reactions of NO_x with O₃. The much stronger ozone titration in the S1 run compared to R0 explains why NO is much less rapidly converted to NO₂ in S1, offsetting the higher emission. The titrating effect of the strong NO_x emissions on the ozone distribution is clearly seen in Fig. 9, especially at and near Le Port and to a lesser extent at other industrialized centres.

By comparison, the effect of this emission reduction is small at the other stations, except during the morning traffic peak (04:00 UT i.e. 08:00 LT, Fig. 11). The model correlates very well with the observed NO₂ and NO_x time-series data, although the model performance is lower on the first and last days of the month, coinciding with unusually large errors in meteorological parameters at Maïdo (e.g. solar radiation and wind direction; see Fig. 8). The average diurnal cycle of concentrations shows a good model agreement for NO₂ and NO_x during night and day, except for a moderate overestimation of NO₂ (Figs. 10 and 11). This discrepancy between the model agreement for NO_x and NO₂ is likely due to the model overestimation of surface O₃ by more than 10 µg m⁻³ on average, leading to an overestimated [NO₂] $/$ [NO] ratio (Reaction R7). Although the S1 run achieves a better agreement with O₃ data at the morning peak hour, the performance of the two simulations is similar at other times of the day, and the shape of the diurnal profile is better represented by the R0 run. The reasons for the ozone overestimation (by ca. 12 µg m⁻³ on average) are unclear, but it could be due to the neglect of halogen chemistry in the WRF-Chem version used in this study. For example, recent studies indicate that the inclusion of halogen (Cl, Br, I) chemistry in a regional or global model might deplete near-surface ozone levels by a much as 7 ppb (∼14 µg m⁻³) in the marine tropical troposphere (Badia et al., 2019; Caram et al., 2023). This reduction of ozone is partly due to a shortening of its lifetime – primarily attributed to iodine chemistry – and partly due to a reduction in ozone production, which is a consequence of the depleting effect of halogens on NO_x levels (Caram et al., 2023; Iglesias-Suarez et al., 2020; Sherwen et al., 2016). Since this influence of halogens is ignored in our model, the good agreement of the model with in situ NO_x measurements (Figs. 10 and 11) might mean that our NO_x emissions are actually too low, although a definitive assessment is not possible at this stage. Although the absence of halogen chemistry likely contributes to ozone overestimation, other factors such as uncertainties in emissions, model resolution, and vertical mixing may also influence the bias.

3.3 Comparisons with PTR-MS measurements

Figures 12 and 13 show modelled VOC concentrations (from simulations R0, S2, S3, S4, and S5) against observations from Maïdo Observatory for the months of January and July 2019. The averaged diurnal cycles of observed and simulated concentrations are shown in Fig. 14. For most species, the general evolution is well represented by the model. Both the observed and modelled diurnal cycles display a pronounced daytime maximum for almost all species. This general behaviour reflects the dominance of daytime sources for most compounds but also the weaker influence of surface emissions during the night, when the Maïdo station is located in or near the free troposphere (Verreyken et al., 2021).

3.3.1 Formaldehyde

Both the model and the measurements display a pronounced diurnal cycle of formaldehyde. The first and last days of each month show the best agreement with measurements, matching both day- and night-time values closely. On the other days, there is a consistent underestimation, especially at night, which can reach a factor of 2. This underestimation is unexpected, given the model's overestimation of FTIR columns (Sect. 3.4) and the fact that night-time observations at Maïdo reflect mostly free tropospheric air. Previous model calculations suggest very little contrast between daytime and night-time formaldehyde (HCHO) columns at Maïdo, due to the long photochemical lifetime of this compound during the night (Stavrakou et al., 2015). Formaldehyde remains mostly stable between the different sensitivity runs, with similar statistics (Table S5). Being mainly produced from the oxidation of other VOCs, including methane, both the simulated and observed formaldehyde concentrations are highest in January (austral summer) due to higher solar radiation fluxes and biogenic VOC emissions at that time.

3.3.2 Methanol

The simulated methanol from the standard run (R0) generally matches the observations very well, except for a slight average positive bias, which is very similar in both January and July (0.2 ppb; see Table S5). This good agreement, and the much higher bias (∼ 0.6 ppb) seen in the S3 run that adopted unadjusted biogenic emission factors (800 µg m⁻² h⁻¹ for all PFTs), appear to validate the halving of the emission factor for the dominant PFTs in run R0 (see Sect. 2.3.2). The lower biogenic emissions are in line with the MEGAN recommendation (Stavrakou et al., 2011). The boundary condition adjustment is also validated, since the simulation that adopted the unadjusted values (S2) also displayed a higher bias (0.3 ppb) and root mean square error (RMSE) than the R0 simulation. The budget of methanol over the island is largely dominated by biogenic sources, since anthropogenic emissions represent only ca. 1 % of the surface methanol emissions in the model (Fig. 4), and the photochemical production of methanol from the reaction of methylperoxy radicals (CH₃O₂) with themselves and with other peroxy radicals (Jacob et al., 2005) is only a minor source (Müller et al., 2016; Bates et al., 2021). Note that the formation of methanol from the reaction of CH₃O₂ with OH (Archibald et al., 2009), which was only recently shown to be a significant source (Bossolasco et al., 2014; Müller et al., 2016; Bates et al., 2021; Caravan et al., 2018), is not included in the model.

Although the R0 simulation performs very well on most days, large positive biases are found during specific periods (e.g. 29–30 January and 2–9 July). These discrepancies are likely not related to the emissions, since other compounds (e.g. Iox and acetaldehyde) are similarly overestimated during the same period, which suggests an issue with the representation of meteorology and possibly tracer transport. Cloudiness (in January) and relative humidity (in both months) are overestimated at Maïdo during those periods of low model performance (Fig. 8).

3.3.3 Isoprene and monoterpenes

Isoprene is generally well reproduced in WRF-Chem during both months, giving credence to the standard emission factor distribution used in the model. The seasonal variation is correctly predicted, with a factor of 3 higher concentrations in January compared to July in both the model and the observations. This difference is primarily due to the higher temperatures and radiation fluxes during January (Table S2). Model underpredictions of isoprene are found on days with very low simulated solar radiation due to excessive cloudiness, in particular, on 1–3 and 28–30 January and 27, 29, and 30 July (Figs. 12 and 13). Conversely, on days with high observed cloudiness and overestimated modelled radiation fluxes, e.g. on 7–8 July, the simulated isoprene is too high. These patterns reflect the strong influence of solar radiation on isoprene emissions in MEGAN. Duflot et al. (2019) and Verreyken et al. (2021) report that westerlies are generally associated with higher isoprene abundances at Maïdo, compared to the easterlies, presumably because of the closer proximity of vegetation west of Maïdo. However, in this study, the easterlies were frequently associated with low solar radiation fluxes due to cloudiness, especially in January (1–3 and 27–30 January), and to a lesser extent on 9 and 12 July, such that the influence of wind direction on isoprene cannot be established for the period considered here. There is essentially no variability in isoprene concentration between the different sensitivity runs, except for simulation S4, which predicts slightly higher concentrations around midday (+5 % at noon), due to the much weaker OH recycling in isoprene oxidation in the original MOZART mechanism compared to the modified mechanism used in R0 (see Sect. 2.2.3). While the average midday isoprene concentration is about 30 % too low in the model, possibly due to a small underestimation of the isoprene emission factor, the concentrations in early morning (05:00–06:00 UT i.e. 09:00–10:00 LT) and late afternoon (17:00–18:00 LT) are too high (Fig. 14). This feature, also found for other biogenic compounds such as methanol, Iox (MVK + MACR + ISOPOOH), and monoterpenes, might be due to insufficient boundary layer mixing.

The model performance for the isoprene oxidation products Iox (MVK + MACR + ISOPOOH) is similar as for isoprene. On many days on which isoprene is too low or too high, so are its oxidation products. There are exceptions to that pattern, such as the first and last days of January, when isoprene is strongly underestimated due to excessive cloud cover, whereas Iox is not. This difference might be related to wind transport errors in the model, as suggested by the poor model performance for wind direction on those days, or it might be due to the longer lifetime of Iox compared to isoprene (Verreyken et al., 2021), implying that Iox is influenced by more distant sources.

Examination of the average diurnal cycle of Iox concentrations (Fig. 14) shows that, while the PTR-MS data show an almost complete disappearance of Iox during the night, WRF-Chem calculates low, but non-negligible concentrations, especially in winter. This discrepancy is even larger in the sensitivity simulation S3. While isoprene remains relatively unchanged between the sensitivity runs, Iox shows important changes in the average diurnal shape between the R0, S3 (unadjusted biogenic emissions), and S4 runs (unadjusted chemistry). As seen in Table 2, the high NO_x yield of the sum MVK + MACR is much higher (+ 50 %) in the updated chemical mechanism of R0 than in the MOZART-4 mechanism used in S4, which explains the higher Iox concentrations in R0 (by about 10 %). Conversely, the simulation S3 predicts significantly higher Iox levels than R0—by ∼ 30 % near midday in July and up to a factor of 5 during the night. The improved diurnal cycle of R0 is due to the strong decrease in the monoterpenes emission factor and their increased light-dependence factor in the R0 simulation, leading to reduced emissions, especially during the night. Indeed, the MOZART-4 mechanism includes a source of MVK and MACR originating from the reaction of lumped monoterpenes with ozone, OH, and NO₃:

$$\begin{array}{}\text{(R8)}& {\displaystyle }{\mathrm{C}}_{\mathrm{10}}{\mathrm{H}}_{\mathrm{16}}+{\mathrm{O}}_{\mathrm{3}}\to \mathrm{0.7}\mathrm{OH}+\mathrm{MVK}+\mathrm{MACR}+{\mathrm{HO}}_{\mathrm{2}}\text{(R9)}& {\displaystyle }{\mathrm{C}}_{\mathrm{10}}{\mathrm{H}}_{\mathrm{16}}+\mathrm{OH}\to \mathrm{TERPO}\mathrm{2}\text{(R10)}& {\displaystyle }{\mathrm{C}}_{\mathrm{10}}{\mathrm{H}}_{\mathrm{16}}+{\mathrm{NO}}_{\mathrm{3}}\to \mathrm{TERPO}\mathrm{2}\text{(R11)}& {\displaystyle }\begin{array}{rl}& \mathrm{TERPO}\mathrm{2}+\mathrm{NO}\to \mathrm{0.1}{\mathrm{CH}}_{\mathrm{3}}{\mathrm{COCH}}_{\mathrm{3}}+{\mathrm{HO}}_{\mathrm{2}}+\mathrm{MVK}\\ & +\mathrm{MACR}+{\mathrm{NO}}_{\mathrm{2}}.\end{array}\end{array}$$

Since monoterpenes are not a source of MVK or MACR (e.g. https://mcm.york.ac.uk/MCM/, last access: 1 March 2025), this artificial source of Iox causes a model overestimation of Iox in both simulations, but especially in simulation S4, and even more during the night and in July, when lower radiation fluxes cause a stronger decrease in isoprene emissions compared to monoterpenes. Both the RMSE and mean bias of R0 are strongly reduced relative to run S4, and an even better agreement would likely be achieved without the artificial source of MVK and MACR from monoterpenes in the model.

The observed ratio of Iox to isoprene concentrations (0.74, based on concentrations averaged between 12:00 and 16:00 LT) is reproduced fairly well by the R0 run in January (0.82), whereas a higher ratio (1.05) is derived from run S3, due to its higher monoterpenes emissions and therefore higher Iox production from their oxidation. In July, the Iox to isoprene ratios are more dissimilar (0.62 based on observations, 1.0 and 1.54 for the R0 and S3 runs), due to the strong reduction of isoprene emissions and the higher share of monoterpenes to the Iox budget in winter compared to summer.

The much improved agreement for monoterpenes after the biogenic emissions adjustments (Table S5 and Fig.  14) appears to validate the lower emissions, especially during the night. A precise assessment is difficult given the large difference between the PTR-MS concentrations derived from the two signals ($m/z$ 137 and 81, corresponding with the protonated monoterpenes and their main fragments, respectively). These differences might result from e.g. temporal variations in monoterpene distribution or the contributions of other compounds to the $m/z$ signal. The LDF correction improved the temporal correlation between the model and the observations through the strong reduction of night-time concentrations, although, as for isoprene, the model simulates early morning and late afternoon concentration peaks that are not seen in the observations. On average, the modelled values in R0 fall between the two measured signals; i.e. they overestimate (by 8–17 pptv) the monoterpenes detected at $m/z$ 137 and underestimate (by 2–11 pptv) those determined based on the $m/z$ 81 signal.

3.3.4 Methyl ethyl ketone and methylglyoxal ($m/z$ 73)

Although both methyl ethyl ketone (MEK) and methylglyoxal (MGLY) contribute to the 73 $m/z$ signal, MEK is dominant according to the model results. As seen in Fig. S3, the modelled MGLY concentrations are highest near local noon and account for at most 25 % (15 %) of the total signal in January (July) in simulation R0, taking into account that the PTR-MS is less sensitive to MGLY (by a factor of 0.7) compared to MEK (Sect. 2.4.2). The observed 73 $m/z$ signal displays a pronounced diurnal cycle (Fig. S34), with daytime concentrations about 3–5 times higher than during night-time. Although the model simulation without any direct MEK emissions (run S5) reproduces this large night/day difference in January, this run strongly underestimates the observations, and it underestimates the amplitude of the diurnal cycle in July (Fig. 14). The photochemical production of MEK in the model originates exclusively from the oxidation of the surrogate compound higher alkanes (BIGALK) by OH. BIGALK emissions include the sources of all C₄ alkanes, but its oxidation mechanism is that of n-C₄H₁₀, a well-known major precursor of MEK (Yáñez-Serrano et al., 2016; Jenkin et al., 1997; Sommariva et al., 2011). Among the higher alkanes, only n-butane and 3-methyl propane are significant MEK precursors (Sommariva et al., 2011; https://mcm.york.ac.uk/MCM/, last access: 1 March 2025). Therefore, since n-butane and 3-methyl propane emissions make up only a fraction of total BIGALK emissions, of the order of 34 % (Stavrakou et al., 2015), the MEK photochemical production in the model is likely overestimated. MGLY has no direct source in the model, but it is photochemically produced, mainly through the oxidation of isoprene and other BVOCs (Mitsuishi et al., 2018). Given its short lifetime (∼ 1 h), primarily due to photolysis (Fu et al., 2008), MGLY displays a pronounced diurnal cycle with a noon maximum. The molar yield of MGLY through isoprene oxidation is of the order of 0.25 (Fu et al., 2008), and it is unlikely to be much underestimated. Therefore, the high MEK/MGLY mixing ratios observed during the day, ca. 0.07 and 0.05 ppb in January and July, respectively, are best explained by the presence of a substantial biogenic source of MEK. This source has been taken to be equal to 3 % of the biogenic isoprene emissions in run R0 (see Sect. 2.3.2). The diurnal shape of the modelled MEK (run R0) in this run is similar to the observations, except that the model underestimates the observed diurnal amplitude in July (Fig. S3). Furthermore, in January, the modelled concentrations rise too early in the morning and decrease too rapidly in the afternoon, suggesting that the biogenic emissions of MEK are delayed compared to those of isoprene. This delay is of the order of 2 h, and it is qualitatively in line with the proposal that MEK is released through the uptake of isoprene oxidation products (MVK and the 1,2-ISOPOOH isomer) by vegetation and their subsequent conversion and re-emission as MEK and other compounds (Cappellin et al., 2019; Tani et al., 2020; Canaval et al., 2020). The PTR-MS signal for Iox (the sum MVK + MACR + ISOPOOH) is slightly delayed (by 1–2 h) compared to isoprene, as shown in Fig. S4  and the delay for MVK is expected to be longer than for Iox, due to the lower rate of reaction of MVK with OH, compared to MACR and ISOPOOH (e.g. Müller et al., 2019). The adopted ratio between MEK and isoprene biogenic emissions (3 % on a mass basis) is larger than 1.5 % value derived by Canaval et al. (2020) based on enclosure measurements on grey poplar trees and field eddy-covariance measurements at two forested sites. This discrepancy could be partly due to uncertainties and natural variability in MVK/ISOPOOH deposition velocities and conversion rates to MEK in plant leaves. In addition, the good model agreement of the R0 run could possibly be achieved with a lower MEK/isoprene emission ratio, e.g. if the contribution of MGLY to the PTR-MS signal is higher. In addition, the transport of chemical compounds from source regions (e.g. cities) to Maïdo is imperfectly represented in the model due to its relatively coarse resolution (2.5 km), which may cause errors in the diurnal cycle of advected species, particularly for MEK, which has a significant anthropogenic component (Fig. 4) (Bon et al., 2011; Brito et al., 2015). Finally, note that the model ignores the potentially significant oceanic source of MEK, which is however very uncertain (Brewer et al., 2020). The S4 run leads to higher MEK/MGLY mixing ratios than R0 around midday (Fig. 14), due to differences in the isoprene chemical mechanism (Sect. 2.2.3); the original MOZART-4 mechanism (used in S4) has a higher yield of C₅-hydroxyaldehydes (HYDRALD), which generate slightly more MGLY than other isoprene oxidation products such as MVK and MACR (Emmons et al., 2010).

In July, the model overestimates both night-time and daytime observations in run S2, which has unadjusted lateral boundary conditions. This overestimation motivated the adjustment (multiplication by 0.4) of ICBCs for MEK (Sect. 2.2.4). This adjustment has a much lesser effect in January, due to its shorter photochemical lifetime in summer compared to winter.

3.3.5 Acetaldehyde

Regarding acetaldehyde (CH₃CHO), WRF-Chem largely underestimates PTR-MS concentrations in January, whereas the model is essentially unbiased on average in July (−0.009 ppb) but correlates very poorly with the observations (r=0.39, Table S5). The primary sources of acetaldehyde over the island, including anthropogenic and biogenic direct emissions and photochemical production, are detailed in Table S6 in the Supplement. The temporal correlation improves markedly (to ca. 0.75 in January and July) when the direct emissions of acetaldehyde and other VOCs are excluded (run S5), even though the mean bias increases. This finding suggests that direct emissions of acetaldehyde are overestimated, while the estimate of photochemical production might be too low. Acetaldehyde emissions in the model are largely due to the residential sector (Fig. 4). Aldehydes are the main NMVOC class released by the residential sector, and these emissions are mainly due to biomass burning (Huang et al., 2017). The model considers the combined emissions of all the higher (C_≥2) aldehydes as anthropogenic emissions of CH₃CHO. However, acetaldehyde accounts for only a small fraction of the total aldehyde emissions from boilers (7 %–17 %) (Macor and Pavanello, 2009) and the heating of cooking oils (at most ∼1 %) (Katragadda et al., 2010; Takhar et al., 2023). The anthropogenic emissions of CH₃CHO used in the model are therefore strongly overestimated, possibly by one order of magnitude or more. Without this lumping, we expect that the model would underestimate the PTR-MS data by almost one order of magnitude, suggesting that other sources of acetaldehyde are likely strongly underestimated. This is in line with previous assessments based on in situ data at remote marine locations (e.g., Wang et al., 2019; Travis et al., 2020). For example, at Cape Verde in the Atlantic Ocean, the CAM-Chem model was found to underestimate in situ concentration measurements by a factor 30 (Read et al., 2012), and the discrepancy was attributed to both oceanic acetaldehyde emissions and photochemical production from unknown, relatively long-lived precursors. Direct oceanic emissions might contribute significantly to the signal at Maïdo (Millet et al., 2010; Read et al., 2012), given the relative proximity of the ocean (∼20 km) and the relatively long lifetime of CH₃CHO (2–9 h, with [OH] = (2–10) ×10⁶ molec. cm⁻²; see Fig. 9). Photochemical production of acetaldehyde in the model results mainly from the oxidation of alkanes, alkenes, and ethanol (Millet et al., 2010). The total production of CH₃CHO calculated using the model (R0 simulation) over Réunion Island averaged over January and July amounts to 19, 0.92, 14, 13, 4.2, 3.86, and 18 t per month for ethane, propane, propene, higher alkenes (BIGENE), BIGALK, MEK, and ethanol, respectively. Unfortunately, only MEK and ethane were measured at Maïdo, and large model underestimations cannot be excluded for any of the other precursors. For example, the average mixing ratios of propane and BIGALK calculated by the model at Maïdo are only ca. 20 and 30 ppt, respectively, whereas long-term in situ measurements at other southern tropical marine sites (Seychelles and Ascension) are considerably higher, of the order of 100 ppt in January and 50 ppt in July for propane, and about 150–250 ppt in January and 100–200 ppt in July for the total concentration of measured butanes and pentanes (n-C₄H₁₀, i-C₄H₁₀, n-C₅H₁₂, i-C₅H₁₂) (Helmig et al., 2021). Biogenic emissions of acetaldehyde and ethanol are included in the model but are very uncertain (Millet et al., 2010) and might be underestimated. As suggested by Travis et al. (2020), oceanic emissions of various acetaldehyde precursors (primarily ethanol, alkanes, and alkenes) might partially explain the negative model biases against aircraft in situ data across the globe and particularly at remote locations and in the free troposphere. These emissions are currently ignored in the model simulations.

3.3.6 Acetone

Due to its relatively long lifetime (several weeks) (Jacob et al., 2002; Fisher et al., 2012), acetone is sensitive to the lateral boundary conditions, especially in July, when solar radiation and OH levels (Fig. 9) are at their minimum. The large positive model bias of the S2 run in July is reduced by a factor of 10 after adjusting boundary conditions (R0 run). Acetone is relatively insensitive to this change in January, likely due to its shorter lifetime and higher biogenic emissions in summer. Acetone is sensitive to a removal of its direct emissions (S5), which have a small anthropogenic component (∼ 3 t per month) and a dominant biogenic contribution (55 t per month, Fig. 4). The secondary source of acetone from photochemical production is comparatively weak over the island (17 t per month on average for the two months). The principal precursors are BIGENE, BIGALK, propane, and monoterpenes, with productions amounting to 8.5, 3.4, 3.5 and 1.7 t, respectively. As in the case of acetaldehyde, the contribution of propane and BIGALK oxidation is likely strongly underestimated. In agreement with previous studies (e.g. Wang et al., 2020), the direct biogenic source of acetone is dominant over the secondary production due to monoterpenes and anthropogenic VOCs. This source is poorly constrained and therefore very uncertain, as the acetone emission factor is constant for all PFTs in MEGAN, except for grasses and crops (Guenther et al., 2012). The ocean–atmosphere exchange of acetone plays a substantial role over the remote troposphere (Wang et al., 2020), but it is neglected in our study. The diurnal cycle of acetone mixing ratios is underestimated by the model during both months (Fig. 14). Since photochemical production plays only a minor role for this compound, this underestimation might be due to either a misrepresentation of transport patterns (Fig. 8) or an underestimation of daytime biogenic emissions. The LDF for biogenic acetone is only 0.25 in the model, implying limited diurnal variability of the emissions. Comparison of the model performance for acetone (Fig. 12) with the meteorological evaluation (Fig. 8) does not suggest a clear link between the modelled diurnal cycle and meteorological parameters, in particular wind direction.

3.3.7 Acetic acid, peracetic acid, and glycolaldehyde ($m/z$ 61)

The $m/z$`61 signal exhibits a pronounced diurnal shape in January and July, with a clear daytime maximum, and night-time values approximately 5 times lower than the midday peak (Fig. 14). Although this signal is mainly attributed to acetic acid, there are interferences, mainly from glycolaldehyde (GLYALD) and peracetic acid (PAA). Here, this measurement is compared against the sum of acetic acid, GLYALD, and PAA from WRF-Chem. The contributions of the three species to the combined modelled concentration are displayed in Fig. 15. The model strongly underestimates the magnitude of the signal, by a factor of 3–4, but it replicates the observed diurnal and seasonal variations. All three species have their highest concentrations during daytime in January, when photochemical activity and biogenic emissions are at their highest. Acetic acid has only a weak biogenic source (5 tons/month on average over the island, Fig. 4), in agreement with previous studies (e.g. Paulot et al., 2011), while all three compounds have significant photochemical production in the model. GLYALD is produced at a high yield (∼ 0.3 at high NO_x) by the oxidation of isoprene by OH and has a short chemical lifetime (a few hours during the day), which explains why its night-time concentrations are very low (Fig. 15). Given the large isoprene emissions, GLYALD is the main contribution to $m/z$ 61 during the day according to the model. The oxidation of monoterpenes also leads to GLYALD formation in the MOZART-4 mechanism, which explains the higher GLYALD concentrations of the S3 run, which uses unadjusted (e.g. 5-times higher) monoterpenes emissions.

The contributions of the two acids show less diurnal variation than GLYALD, particularly in July. Photochemical production of these two compounds over the island amounts to 19.5 and 62.4 t per month for acetic acid and PAA, respectively, on average for the two months. This is mainly due to the reaction of the acetylperoxy radical (CH₃CO₃) with HO₂, represented in the model as

$$\begin{array}{}\text{(R12)}& \begin{array}{rl}& {\mathrm{CH}}_{\mathrm{3}}{\mathrm{CO}}_{\mathrm{3}}+{\mathrm{HO}}_{\mathrm{2}}\to \mathrm{0.75}\mathrm{PAA}+\mathrm{0.25}{\mathrm{CH}}_{\mathrm{3}}\mathrm{COOH}\\ & +\mathrm{0.25}{\mathrm{O}}_{\mathrm{3}}\end{array}\end{array}$$

This reaction alone accounts for 79 % and 94 % of the total production of acetic and peracetic acid over the island, respectively. The model ignores the radical-forming channel,

$$\begin{array}{}\text{(R13)}& {\mathrm{CH}}_{\mathrm{3}}{\mathrm{CO}}_{\mathrm{3}}+{\mathrm{HO}}_{\mathrm{2}}\to {\mathrm{CH}}_{\mathrm{3}}{\mathrm{CO}}_{\mathrm{2}}+\mathrm{OH}+{\mathrm{O}}_{\mathrm{2}}\end{array}$$

(followed by CH₃CO₂ (+O₂) → CH₃O₂ + CO₂), which was shown to account for about half of the total reaction (Jenkin et al., 2007). The total rate constant of the reaction, $\mathrm{1.4}\times {\mathrm{10}}^{-\mathrm{11}}$ molec.⁻¹ cm³ s⁻¹ at 298 K in MOZART-4, is also likely underestimated by a factor of ca. 1.5 compared to more recent estimates (Gross et al., 2014). These approximations might contribute to a moderate overestimation of acetic and peracetic acid production from the reaction, although the probable underestimation of alkanes and other VOCs noted previously likely implies an underestimation of the production rate of CH₃CO₃ radicals. In addition, many reactions of CH₃CO₃ radicals with other organic peroxy radicals are neglected in the mechanism, although they might contribute significantly to acetic acid production (Khan et al., 2018). Furthermore, the rate constant of the reaction of peracetic acid with OH was recently shown in a combined experimental and theoretical study (Berasategui et al., 2020) to be much lower (factor of ∼ 30) than previously determined, implying a much longer lifetime in the troposphere.

The model underestimation of the $m/z$ 61 signal is consistent with (1) the suggestion that there are underestimated or unknown sources of acetic acid, based on model comparisons with in situ measurements from aircraft and surface stations (Paulot et al., 2011; Khan et al., 2018), and (2) observations of high PAA concentrations over remote oceans during the Atmospheric Tomography Mission (ATom) (Wang et al., 2019). Given that acetaldehyde is a known major precursor of PAA, box model calculations by Wang et al. (2019) show that these high levels of PAA are fully consistent with the (simultaneously observed) acetaldehyde measurements. As noted previously, there is abundant evidence of a large missing source of acetaldehyde over the oceans.

Figure 12 Time series of observed and modelled concentrations of PTR-MS measured species in January. All concentrations are shown in ppb. The model results are shown for simulations R0, S2, S3, S4 and S5 (see Table 4). Both monoterpene measurements are depicted, differentiated by solid ($m/z$ 81) and dashed black lines ($m/z$ 137).

Figure 13 Same as Fig. 12, for July 2019.

Figure 14 Average diurnal cycle of organic compounds measured at Maïdo Observatory and corresponding modelled concentrations from the 30 d runs defined in Table 4 (except S1). The concentrations are shown in ppb. Both monoterpene measurements are depicted in the same figure, differentiated by solid ($m/z$ 81) and dashed black lines ($m/z$ 137).

Figure 15 Averaged diurnal cycle of the $m/z$ 61 signal in January (a) and July (b) and comparison with the modelled concentrations of acetic acid, GLYALD, and PAA, which are all assumed to contribute to the observed signal. The PTR-MS sensitivity is assumed identical for the three compounds.

3.4 Evaluation against FTIR columns

The model is evaluated against FTIR data in Fig. 16. The model values are sampled at the same times as the measurements and smoothed using the averaging kernels and a priori profile of the FTIR retrieval. The results of two simulations are shown: R0, the standard run with all model updates, and S2, identical to R0 but without the adjustments to the ICBCs. The results of the other simulations differ only marginally from those of the R0 run and are therefore not shown. ICBC adjustments were made only for CH₃OH, PAN, and C₂H₆; for the other compounds, the discrepancies are generally small.

For HCHO, a model overestimation is noted during both months, amounting to 8 and 2×10¹⁴ molec. cm⁻² in January and July, respectively. In January, this bias exceeds the reported systematic uncertainty of the FTIR column ($\sim \mathrm{2.5}\times {\mathrm{10}}^{\mathrm{14}}$ molec. cm⁻²) (Vigouroux et al., 2018) and is not likely due to ICBC errors for HCHO, given the short lifetime of this compound (a few hours). The monthly averaged retrieved columns (3.2 and 1.61×10¹⁵ molec. cm⁻² in January and July, respectively) are similar to the reported FTIR columns at Maïdo for the same months in previous years (Vigouroux et al., 2018, 2020). The weak sensitivity of the modelled columns to changes in local VOC emissions (e.g. from run S5) suggests that the model overestimation is either due to HCHO precursors that are not well represented in the model or background HCHO, i.e. primarily the contribution of methane oxidation, above ∼ 2 km altitude (the altitude of Maïdo). The addition of lightning NO emissions (run S6) aggravates the bias, as it increases the HCHO column (+ 24 % in January; see Fig. S5). The HCHO enhancement due to lightning NO is mostly located in the mid-troposphere, between 4 and 11 km (Fig. S6), and it is primarily caused by the increase in OH resulting from the added NO_x emissions in this altitude range, as the OH radical promotes methane oxidation and its associated HCHO production. At lower altitudes, OH levels are less impacted, whereas at higher altitudes, HCHO is more strongly influenced by the deep convection of low-altitude pollutants to higher altitudes (Bozem et al., 2017; Fried et al., 2008).

Figure 16 Time series of FTIR-measured and modelled HCHO, CH₃OH, C₂H₆, CO, PAN, and tropospheric O₃ columns in January (a–f) and July (g–l) 2019. The symbols denote the individual FTIR measurements and corresponding WRF-Chem results, while the solid lines (red for R0, cyan for S2) show the 24 h averages.

Given the longer atmospheric lifetimes of the other compounds shown in Fig. 16, the model performance for these species is very dependent on their ICBCs, obtained from either CAM-Chem (for CH₃OH) or CAMS (for the other compounds). Without any adjustment, the model agreement with the observed CO and tropospheric O₃ columns is reasonable, with biases generally well below 10 % for CO and 20 % for O₃. The simple ICBC adjustments, presented in Sect. 2.2.4, generally succeed in reducing the gap between the model and observations for the other compounds. Since the ICBC adjustments consist of a single scaling factor for each species, the short-term variability in the observations remains poorly reproduced by the model after adjustment, e.g. for CH₃OH. Nevertheless, the model correlates well with the daily-averaged data for O₃ (Pearson's correlation coefficient r=0.6 and 0.7 in January and July) and to a lesser extent for C₂H₆ (r=0.7 and 0.3).

The monthly averaged FTIR PAN columns, ca. 1×10¹⁵ molec. cm⁻², are lower, but of the same order of magnitude as those typically observed above the Jungfraujoch station in the Swiss Alps during summer (up to 4×10¹⁵ molec. cm⁻²) (Mahieu et al., 2021). They are also significantly lower than those observed by the Infrared Atmospheric Sounder Interferometer (IASI), typically > 2×10¹⁵ molec. cm⁻² above oceanic areas around Réunion Island (Franco et al., 2018), although this difference is likely partly due to the high altitude of the station. The spatial distribution displayed by IASI suggests that the long-range transport of African emissions has a substantial effect on PAN levels in the region. Although PAN is short-lived in the lower troposphere due to fast thermal decomposition, it is much longer-lived in the cold upper troposphere, where it can be transported over long distances. Due to this long lifetime, a large fraction of the column lies at those high altitudes, which explains why the modelled PAN column at Maïdo responds strongly to the ICBC change between simulations S1 and R0. Note that this does not imply that ICBCs have a strong influence on PAN and NO_x close to the surface, where the PAN lifetime is short. The modelled columns of the S1 run overestimate the FTIR data by about a factor of 2, providing the justification for the adjusted ICBCs for PAN in simulation R0. However, both FTIR data and the R0 run underestimate the IASI columns by a factor of 2 or more. The acknowledged uncertainties in the PAN retrievals by FTIR and IASI suggest that the ICBC adjustment for PAN should be considered very uncertain.

3.5 TROPOMI

The distribution of the NO₂ and HCHO columns obtained from the oversampling of TROPOMI retrievals over a period of more than 4 years is displayed in Fig. 17. Despite the much lower NO₂ columns retrieved over and around the island (< 1.5×10¹⁵ molec. cm⁻²), in comparison with e.g. Europe (typically (2–10) ×10¹⁵ molec. cm⁻²) (Poraicu et al., 2023; Lange et al., 2023), sharp gradients are observed, reflecting the strong heterogeneity of NO_x emissions. A very well-localized maximum (1.5×10¹⁵ molec. cm⁻²) is found at the precise location of the Le Port power plants (Fig. 17), whereas the other thermal power plants, located near Sainte-Suzanne (Bois Rouge) on the northern coast, and near Saint-Pierre (Le Gol) on the southern coast, are not clearly detected. The emissions from the Le Port power plants (and to a lesser extent, from the cities of Le Port, Saint-Denis and Saint-Paul) generate a broad hotspot of NO₂ columns exceeding 1×10¹⁵ molec. cm⁻², which extends over land as well as over the ocean over distances of the order of 5–10 km. This distribution confirms the dominance of Le Port power plant emissions over those from the other point sources and justifies the redistribution of industrial emissions based on the EDGAR 6.1 dataset described in Sect. 2.3.1. The NO₂ column values around Maïdo are significantly lower, of the order of 0.5×10¹⁵ molec. cm⁻². Even lower columns are observed in the southeastern part of the island, around the volcano (Piton de la Fournaise).

Figure 17 Oversampled NO₂ and HCHO columns from TROPOMI, covering 2018–2022. The resolution is 0.01° × 0.01°.

The oversampled distribution of TROPOMI HCHO is noisier compared to NO₂, but it also shows well-defined features and gradients. The largest columns (∼ 5×10¹⁵ molec. cm⁻²) are found in a large area along the northwestern coast and, in particular, near the Le Port power plants. Since the power plants are not a large VOC source, and since the HCHO distribution shows little correlation with the anthropogenic VOC emission distribution (Fig. 6), these high columns clearly do not reflect the distribution of VOC emissions. Instead, they are likely due to the influence of anthropogenic NO_x, leading to higher OH levels and therefore to enhanced HCHO production rates from the oxidation of long-lived hydrocarbons, especially methane. Methane oxidation is by far the largest source of HCHO at the global scale (Stavrakou et al., 2009), and to a large extent also above Réunion Island, even though non-methane hydrocarbons are the dominant source over most continental areas. This dominance of the background is the main reason for the many common features shared by the HCHO and NO₂ distributions in Fig. 17, such as the high values to the west and northwest of the island, and the lower columns to the east and southeast and over several specific areas in the island interior. Besides this likely dominant contribution of the methane background to HCHO, the TROPOMI distribution also shows evidence of elevated columns due to isoprene emissions, particularly along the eastern coast (see Fig. 5).

The model is evaluated against NO₂ columns from TROPOMI in January and July 2019 in Fig. 18. The domain shown on Fig. 18 encompasses both Réunion and Mauritius Island, located to the Northeast of Réunion, at a distance of ca. 200 km. The results of two simulations are shown: the reference run without lightning (R0) and the simulation with lightning emissions (S6). The modelled columns of run S1 (with higher NO_x emissions from the Le Port power plants) are given in the Supplement (Fig. S7). The TROPOMI averaging kernels were applied for both compounds, and the model averages were calculated from values sampled at the same times as the TROPOMI monthly averages. Both modelled and retrieved columns were regridded to 0.1° resolution.

Both the model and TROPOMI data display higher NO₂ columns on Mauritius Island compared to Réunion. Based on EDGARv6.1 emissions data, the main hotspot on Mauritius is due to the energy sector, with several fossil fuel power plants located in Port-Louis, the capital city. The model comparison with TROPOMI strongly suggests that the EDGAR NO_x emissions are overestimated. These emissions generate a large plume extending towards the southwest (in January) or the west (in July), but their impact on Réunion appears to be limited, based on the distributions shown in Fig. 18.

Over Réunion, both the R0 and S6 model simulations reproduce the main spatial patterns observed in the measurements, namely, the NO₂ maximum at Le Port and over adjacent regions, primarily the eastern and (to a lesser extent) the southeastern coast in January, and the northern coast, in July. A secondary NO₂ maximum is also seen along the southern coast near Saint-Pierre and especially near the thermal power plant of Le Gol. The model also reproduces the minimum seen in the eastern part of the island in January and near the southeastern extremity of the island. Overall, both R0 and S6 overestimate the observations in winter, but R0 underestimates during summer (−11 % and +64 % for the island-averaged column of R0 and S6 in January, respectively, and + 25 % for both in July). The discrepancy is, however, locally much higher, e.g. near Maïdo in July (factor of ∼ 1.8 at the two pixels nearest to Maïdo), likely due to the export of pollution from the Le Port area. Caution is warranted when comparing model and satellite data over remote locations, as the satellite columns are very noisy, despite their spatial (0.1°) and monthly averaging, especially over the mountainous island interior. The TROPOMI NO₂ uncertainties are of the order of 5×10¹⁴ molec. cm⁻² (see Fig. S8); i.e. they are of the same order as the retrieved columns over most of the island except the strongest anthropogenic hotspots. The model correlates better with the data in July (r=0.72 for R0 over the inner domain; see Fig. 18) than in January (r=0.43), likely due to the higher values in July. The model reproduces the observed seasonal variations in the NO₂ columns, with higher values in July (winter) compared to January (summer), in response to changes in OH radical concentrations (Fig. 9). The R0 run provides a significantly better match with TROPOMI than the S1 run, which strongly overestimates the columns along the western and southwestern coast (Fig. S7) by a factor of 2 as well as over the ocean, west, and southwest of the island. These overestimations essentially disappear in the R0 run, which assumes a 5-fold reduction in NO_x emissions due to the Le Port power plants.

The causes for the slight NO₂ model overestimation in the R0 run in July are unclear. The NO_x emissions could be overestimated, the NO_x lifetime could be too long, or the TROPOMI NO₂ columns might be too high. Unfortunately, most TROPOMI NO₂ validation studies were conducted at mid-latitudes, in regions with strong pollution sources (Verhoelst et al., 2021; Poraicu et al., 2023; Lange et al., 2023; van Geffen et al., 2022a). The previously cited studies have generally reported the slope (s) and intercept (i) of linear regressions of the type C = i + sC^′, where C denotes the TROPOMI column and C^′ the co-located independent measurement. Most validation studies reported positive values for the intercept i and slope values s lower than unity, suggesting that TROPOMI underestimates high values and overestimates very low values. More studies are needed to better characterize the potential biases of TROPOMI NO₂ in remote areas.

Figure 18 Monthly-averaged NO₂ columns from TROPOMI and WRF-Chem (simulations R0 and S6) for January (a–c) and July (d–f), regridded to 0.1° resolution. The WRF-Chem averages account for TROPOMI sampling times and for the averaging kernel corresponding to each overpass. The white line denotes the boundary of the inner model domain (d02).

Implementation of lightning emissions (run S6) leads to a strong column enhancement in January, while its impact is negligible in July. The S6 run worsens the model overestimation over land and leads to overestimated columns over the ocean. For example, the average column from run S6 over an oceanic area north of Réunion is 6.3×10¹⁵ molec. cm⁻², a factor of 1.75 above the TROPOMI average of 3.61×10¹⁵ molec. cm⁻², whereas the R0 run is a factor of 2 too low (1.83×10¹⁵ molec. cm⁻²). The lightning emissions flip the seasonality of NO₂ columns over the ocean, with higher values predicted in January, while both TROPOMI and the R0 run suggest a wintertime maximum. In January, lightning NO_x is responsible for increasing the average mixing ratio of NO₂ in the upper troposphere by a factor of 4 compared to the R0 simulation (Fig. S6). Note that the default WRF-Chem settings would lead to much larger lightning emissions since, following Barten et al. (2020), we decreased the number of flashes by a factor of 10, and we also reduced the NO production per flash to 250 moles. Nevertheless, the S6 run substantially overestimates the upper tropospheric (UT) NO₂ mixing ratios obtained from TROPOMI using the cloud slicing technique (Marais et al., 2021; Horner et al., 2024). Indeed, the TROPOMI-based UT mixing ratios (between 180 and 450 hPa) are typically 40–50 pptv above the Indian Ocean around Réunion Island during December 2019–January 2020 (Marais et al., 2021), whereas the model averages over this region are 28 and 120 pptv in the R0 and S6 runs, respectively. Both the model comparison with the tropospheric columns (Fig. 18) and with the UT mixing ratios based on TROPOMI indicate that the lightning source in the S6 run is too high, despite the reduction of flash count and number of moles per strike. In July, lightning emissions are largely insignificant, with a less than 2.5 % increase in NO_x in the upper troposphere.

Much like for TROPOMI NO₂, the model evaluation against TROPOMI HCHO data shows a relatively good agreement with regards to spatial representation (Fig. 19), even though the spatial correlation coefficients are very low (e.g. r=0.2 for the R0 run in January), largely due to the high noise in the data. In January, relatively high columns (of the order of 8×10¹⁵ molec. cm⁻² in the R0 run, slightly less in the TROPOMI columns) are found to the southwest of the two islands, as well as along the western coast of Réunion. These patterns mirror the NO₂ distribution, which also shows enhancements in those areas (Fig. 18). Nevertheless, the HCHO columns from both R0 and S6 show overestimations when compared to the observations. For example, TROPOMI shows values close to 6×10¹⁵ molec. cm⁻² along the northwestern coast of the island in January, slightly lower than the modelled values close to 8×10¹⁵ molec. cm⁻². The discrepancy is close to the uncertainty associated with the retrieval product, 2×10¹⁵ molec. cm⁻² for the systematic component (trueness) and ca. 1×10¹⁵ molec. mcm⁻² for the random part (precision), when accounting for the number of measurements used in the shown averages (see Fig. S9). This model overestimation is consistent with the model overestimation against FTIR HCHO columns at Maïdo (see Sect. 3.4). As pointed out previously, the main source of HCHO in this region is methane oxidation. Based on the model output, this production is estimated at about 3.0 and 1.2 Gg per month over the land area of Réunion Island for January and July, respectively. The overestimation of HCHO levels in the model, in comparison with TROPOMI and FTIR data, can probably be explained by an overestimation of OH levels in the model simulations. Note that the addition of lightning in S6 increases OH levels and therefore worsens the overestimation of HCHO columns in January. There might be several causes for the HCHO overestimation. NO_x might be too high, as suggested by the moderate overestimation of NO₂ columns against TROPOMI data, and this could impact OH. In addition, the omission of halogen chemistry in the model likely also leads to OH overestimation (Sherwen et al., 2016). Finally, it was recently found out that the absorption of UV radiation by water vapour is more significant than previously assumed and slows down the rate of ozone photolysis and therefore the production and concentrations of OH radicals in the lower troposphere (Prather and Zhu, 2024). This process is particularly efficient at tropical latitudes where water vapour is most abundant. Note, however, that this effect might be partially counterbalanced by the reduction of HCHO photolysis rates, following the decrease of UV radiation levels.

The HCHO production due to methane oxidation might be compared with the contribution of isoprene, a major source of HCHO at the global scale. Taking an average isoprene flux of 0.25 Gg month⁻¹ over the island (Fig. 4), and assuming 2.5 HCHO molecules produced for one molecule of isoprene (Stavrakou et al., 2009), the HCHO production from this source is estimated to be 0.28 Gg month⁻¹, i.e. almost an order of magnitude lower than the contribution of methane oxidation (see previous paragraph). The secondary source of HCHO from anthropogenic VOC oxidation might be significant, but as the lifetimes of these compounds are much longer than that of isoprene, the HCHO production over the island is diluted due to mixing and is partially exported away from the island.

While the model overestimation of HCHO columns with respect to TROPOMI might be partially due to model errors, the TROPOMI columns are very low and noisy, and therefore close to the detection limit of the instrument. Taking the detection limit as three times the precision, most frequently in the range (0.5–1.5) × 10¹⁵ molec. cm⁻² for 0.1° pixels over the island (Fig. S9), and taking into account that the TROPOMI precision was found to be underestimated by a factor of about 2 (1.6–2.3) in comparisons against FTIR data from a wide network of stations (Vigouroux et al., 2020), almost all TROPOMI monthly averages over the island (except around the Le Port hotspot) fall below the detection limit, which explains the high noise seen on the TROPOMI HCHO maps (Fig. 19).

Figure 19 Monthly-averaged formaldehyde columns from TROPOMI and WRF-Chem (simulations R0 and S6) for January (a–c) and July (d–f), regridded to 0.1° resolution. The WRF-Chem averages account for the TROPOMI sampling times and for the averaging kernel corresponding to each overpass. The white dotted line denotes the borders of the inner model domain (d02).

4 Conclusions

The WRF-Chem model has been used to compute the atmospheric composition in a region surrounding Réunion Island in January and July 2019. The oversampled TROPOMI NO₂ distribution demonstrates the predominant influence of anthropogenic emissions on NO_x abundances over the island. In particular, the power plants near Le Port cause by far the strongest NO_x hotspot over the island, implying that their emissions largely exceed the emissions of any other point source or major city over Réunion. The model (R0 simulation) is only moderately biased against TROPOMI NO₂ columns (−11 % in January, + 25 % in July) and reproduces the main observed patterns. To reach this agreement, the high-resolution (1 km²) anthropogenic emission inventory used as the input to the model has been adjusted with regard to the repartition of energy sector emissions among the different power plants (in accordance with EDGAR v6.1), and furthermore, the NO_x emissions due to the Le Port power plants were downscaled by a factor of 5. This reduces the industry/energy sector emission estimates from Atmo-Réunion and EDGAR v6.1 by factors of 3.8 and 1.9, respectively. The emissions from the other power plants could not be constrained based on TROPOMI data. The 5-fold reduction here is a crude adjustment, but it is confirmed by the model comparisons with both TROPOMI NO₂ and in situ NO_x measurements in the direct vicinity of the Le Port power plants. The fairly good model agreement with in situ NO_x measurements at other air quality stations adds credibility to the anthropogenic NO_x emission estimates for the other sectors (mostly traffic). In situ O₃ concentrations are overestimated by the model at the air quality stations, by ca. 6 ppb on average. This is likely due to the neglect of halogen (Cl, Br, I) emissions and chemistry in the model. These compounds were shown to deplete O₃ substantially in the tropical marine troposphere, by about 7 ppb, in line with previous work.

Over the ocean, the model comparisons with TROPOMI (for January) show the importance of lightning as a source of NO. Accounting for this source (S6 run), or neglecting it (R0) leads to either overestimated or underestimated modelled NO₂ columns in comparison to TROPOMI over the ocean; the comparison with upper tropospheric NO₂ mixing ratios based on cloud-sliced TROPOMI NO₂ data follows the same trend and confirms the significant impact of lightning in this region in January. However, the lightning emissions are strongly overestimated by WRF-Chem with its current parameterization, even when the flash count is reduced by an order of magnitude in the model.

The model performance against PTR-MS VOC measurements is species-dependent, and the model comparisons prompted several adjustments to the MEGAN model-calculated emissions. Most importantly, biogenic emissions of methanol and monoterpenes were downscaled by factors of about 2 and 5, respectively, and the light-dependent fraction for monoterpenes was increased to 0.9 (see Sect. 2.3.2) to provide a fair match with the measurements. In addition, the chemical mechanism of the model was updated to better account for OH recycling in isoprene oxidation and for MEK degradation mechanisms. The model performs quite well for isoprene and its oxidation products (Iox), but their diurnal shape displays unrealistic peaks in the early morning and late afternoon (a feature also seen for other compounds), likely due to an underestimation of vertical mixing in the model. The ratio of Iox to isoprene (ca. 0.8 around noon in January) is fairly well reproduced. This provides some reassurance regarding the model performance, since the ratio is strongly dependent on oxidative conditions, i.e. on OH levels.

Despite a good reproduction of the seasonal variation and diurnal shape of the formaldehyde concentrations observed by PTR-MS, the model underestimates the measurements by up to 25 % during the day and a factor of 2 during the night. This stands in contrast with the model overestimation of HCHO column measurements by FTIR (also at Maïdo) and TROPOMI for reasons unclear. The discrepancy when compared with the column measurements could be due to an overestimation of the background HCHO production, primarily due to methane oxidation, resulting from overestimated OH levels in the free troposphere. The underestimation against PTR-MS could possibly indicate direct HCHO emissions; however, the night-time measurements at Maïdo are primarily influenced by the free troposphere. More work will be needed to elucidate the factors influencing formaldehyde concentrations around Maïdo.

The model fails to reproduce both the seasonal cycle of acetaldehyde (higher values observed in summer compared to winter) and its diurnal variation (strong midday peak). The anthropogenic source of CH₃CHO in the model is almost certainly too high, due to the lumping of higher aldehydes into this compound. Lower anthropogenic emissions and a much higher photochemical production would be needed to match the observations. The modelled concentration of C_>3 alkanes is very likely underestimated, as shown by the evaluation of the model against measurements in similar environments. Therefore, as proposed in previous studies, a large part of the missing production of acetaldehyde may be from the oxidation of alkanes, alkenes, and alcohols, presumably released by the oceans.

MEK being the largest contribution to the 73 $m/z$ signal from the PTR-MS, the large observed daytime mixing ratio suggests the presence of a significant biogenic source of MEK, in line with the proposal that the deposition of isoprene oxidation products (MVK and 1,2-ISOPOOH) on vegetation is followed by their partial conversion and re-emission as MEK. The ratio between MEK and isoprene emission that achieves the best agreement with the data is 3 % (mass basis) higher than the value of this ratio derived by Canaval et al. (2020) based on more direct emission measurements for different environments. The difference might be due to natural variability and/or model uncertainties related e.g. to other contributions to the observed signal.

Clearly, more work will be needed to understand and quantify the budget of key compounds (O₃, OH, NO_x, VOCs, and OVOCs) around Réunion Island and similar environments. From the modelling perspective, efforts should be made to improve the representation of NMVOC speciation, with more explicit chemical mechanisms (e.g. for monoterpenes and aromatics) and more detailed emission inventories (e.g. for higher aldehydes). The advection of chemical compounds should be improved, e.g. with finer resolution modelling to better represent transport in a mountainous environment, and the vertical mixing of the WRF-Chem model (especially during night-time) should be addressed in future model studies. The ocean–atmosphere exchange of important VOCs should be implemented in the model, as it is a well-recognized source of several OVOCs and their precursors (alkanes, etc.). The ocean is also a large source of halogens that have far-reaching impacts on ozone, OH, NO_x, and VOC levels, and their emissions and chemistry should be implemented in WRF-Chem. Efforts should also be made to improve the representation of the boundary conditions of the regional atmospheric models. FTIR column data were used in this work to adjust the boundary conditions for several compounds; other species should be adjusted based on current knowledge relying on campaign data and network in situ measurements.