Rapid oxidation of phenolic compounds by O₃ and HO^●: effects of the air–water interface and mineral dust in tropospheric chemical processes

Environmental media affect the atmospheric oxidation processes of phenolic compounds (PhCs) released from biomass burning in the troposphere. To address the gaps in experimental research, phenol (Ph), 4-hydroxybenzaldehyde (4-HBA), and vanillin (VL) are chosen as model compounds to investigate their reaction mechanism and kinetics at the air–water (A–W) interface, on TiO₂ mineral aerosols, in the gas phase, and in bulk water using a combination of molecular dynamics simulation and quantum chemical calculations. Of the compounds, Ph was the most reactive one. The occurrence percentages of Ph, 4-HBA, and VL staying at the A–W interface are ∼ 72 %, ∼ 68 %, and ∼ 73 %, respectively. As the size of (TiO₂)_n clusters increases, the adsorption capacity decreases until n > 4, and beyond this, the capacity remains stable. A–W interface and TiO₂ clusters facilitate Ph and VL reactions initiated by the O₃ and HO^•, respectively. However, oxidation reactions of 4-HBA are little affected by environmental media because of its electron-withdrawing group. The O₃- and HO^•-initiated reaction rate constant (k) values follow the order of A–W_Ph > TiO_2 VL > A–W_VL > A–W_4-HBA > TiO_2 4-HBA > TiO_2 Ph and TiO_2 VL > A–W_Ph > A–W_VL > TiO_2 4-HBA > TiO_2 Ph > A–W_4-HBA, respectively. Some byproducts are more harmful than their parent compounds, so they should be given special attention. This work provides key evidence for the rapid oxidation observed in the O₃/HO^• + PhC experiments at the A–W interface. More importantly, differences in the oxidation of PhCs by different environmental media due to the impact of substituent groups were also identified.

1 Introduction

Biomass burning, stemming from natural wildfires and human activity, significantly contributes to atmospheric particulate matter (PM). Biomass burning is a primary source of approximately 90 % of the global primary organic aerosols (POAs) and releases a substantial quantity of organic pollutants (Ito and Penner, 2005; Chen et al., 2017, 2023). Biomass burning is to blame for about 62 % of the total annual emissions of about 8.0 Tg of black carbon and 93 % of the total annual emissions of about 33.9 Tg of organic carbon worldwide (Bond et al., 2004). Emissions from biomass combustion are one of the primary sources of atmospheric and particle pollutants that negatively affect human health, air quality, and climate (Reid et al., 2005; Yao et al., 2016). One of the three main types of biopolymers responsible for the formation of biomass is lignin (Sun et al., 2011), which is also the polymeric organic molecule most abundant in plants (Lou et al., 2010; Soongprasit et al., 2020). Pyrolysis of lignin releases phenolic compounds (PhCs) into the air, including phenols, phenolic aldehydes, and methoxyphenols. By mass, these PhCs make up between 21 % and 45 % of the aerosol composition (Hawthorne et al., 1989; Diehl et al., 2013; Liao et al., 2020; Soongprasit et al., 2020). Methoxyphenols are one of the potential tracers that can be found in atmospheric wood smoke pollution, with the emission rate ranging from 900 to 4200 mg kg⁻¹ fuel (Hawthorne et al., 1989; Rogge et al., 1998; Simoneit, 2002; Chen et al., 2017). Evidence shows that the oxidation processes of PhCs can result in the formation of secondary organic aerosol (SOA) (Yee et al., 2013; Jiang et al., 2023). Hence, it is imperative to explore the effects of PhCs when exposed to atmospheric oxidants.

After being released into the atmosphere, PhCs will be oxidized by ozone (O₃) and hydroxyl radicals (HO^•). Both are significant contributors to SOA (Arciva et al., 2022). The homogenous oxidation of PhCs has been the emphasis of previous studies (Henry et al., 2008; Yee et al., 2013; Liu et al., 2019; Arciva et al., 2022). In the past decade, researchers have investigated the kinetics and reaction mechanisms of gas-phase interactions of PhCs with O₃ and HO^• (Kroflič et al., 2018; Smith et al., 2016; J. Sun et al., 2021; Y. Sun et al., 2021; Liu et al., 2022). Furthermore, they have investigated the hydroxylation, ring opening, and oligomerization processes of PhCs in the atmospheric liquid phase, with a focus on the potential environmental toxicity and climatic effects of these events (Ma et al., 2021; Liu et al., 2022; Arciva et al., 2022; Carena et al., 2023).

However, there is a dearth of specific data as well as of explanations of the mechanisms involved in the atmospheric oxidation of PhCs at the air–water (A–W) interface. The atmosphere contains a high concentration of aqueous aerosols and water microdroplets (Zhong et al., 2019; Guzman et al., 2022). The oxidation of PhCs can rapidly occur at the A–W interface (Rana and Guzman, 2022c). The term “water surface catalysis” denotes the phenomenon where chemical reactions happen at a faster rate at the A–W interface compared to the bulk phase (Lee et al., 2015a, b; Yan et al., 2016; Banerjee et al., 2017). In chemical engineering, titanium dioxide (TiO₂) is an essential photoactive component found in atmospheric mineral dust (Sakata et al., 2021; Wang et al., 2023). The interaction between PhCs and TiO₂ is continuous (Grassian, 2009; Rubasinghege et al., 2010; Shang et al., 2021), despite the relatively low prevalence of TiO₂ mineral particles (comprising 0.1 % to 10 % by mass). Therefore, it is essential to investigate the disparity between the oxidation reaction mechanisms and kinetics of PhCs at the A–W interface and mineral dust particles.

Increasing the number of constituents on the aromatic ring would affect the reactivity and lead to complex compounds after reaction addition and/or open-ring pathways. Phenol (Ph), 4-hydroxybenzaldehyde (4-HBA), and vanillin (VL) are typical lignin pyrolysis products (Jiang et al., 2010; Kibet et al., 2012). Thus, we selected Ph, 4-HBA, and VL as model compounds to present comprehensive mechanistic information at the A–W interface, on TiO₂ clusters, in the gas phase, and in bulk water, using a combination of molecular dynamics simulation and quantum chemical calculations. Rate constants were calculated within a wide temperature range in various environmental media. Additionally, computational toxicology was employed to evaluate the ecotoxicological impact of PhCs and their transformation products.

2 Methods

2.1 Molecular dynamics simulation

All of the molecular dynamics simulations were carried out by utilizing the GROMACS 2019 package, which included the AMBER force field. Parametrization of the Ph, 4-HBA, and VL was accomplished using the GAFF force field in conjunction with RESP charge calculations performed at the M06-2X/6-311++G(3df,2p)//M06-2X/6-31+G(d,p) level. The TIP3P water model was utilized so that individual water molecules may be represented (Jämbeck and Lyubartsev, 2014).

2.1.1 Properties of Ph, 4-HBA, and VL at the A–W interface

Considering the significance of the interfacial behavior of Ph, 4-HBA, and VL at the A–W interface, the properties of these three substances were initially examined by focusing on the A–W interface. Figure S1a in the Supplement depicts a rectangular box that has dimensions of 4 × 4 × 9 nm³ and has a z axis that is perpendicular to the A–W contact. This box was used for all simulations. A water box that is too small may cause the central PhC molecules to be too close to the interface region, leading to inaccurate results. Conversely, opting for a water box that is too large can lead to an unnecessary waste of computational resources. To begin the process of constructing the initial configurations, a water slab measuring 4 × 4 × 4 nm³ was positioned at the coordinates (2, 2, 4.5 nm) of the center of mass (COM). Because the rest extension along the z axis of the box was sufficiently large (2.5 nm³), it was possible to steer clear of the intersection of two A–W interfaces. Prior to the formal simulation, six Ph molecules were randomly selected and placed in a vacuum above the water box for 150 ns of NVT molecular dynamics simulation. The purpose of simulating 150 ns is to capture the fundamental molecular dynamics that occur on this timescale, such as bond formation, conformational changes, and interaction events. The results show no significant π–π interactions or formation of hydrogen bonds between the Ph molecules. To simplify the model, this was followed by simulations of individual molecules. Ph, 4-HBA, and VL were each placed in their own compartment at the coordinates (2.0, 2.0, 7.75 nm) for each system in order to simulate the behavior of these molecules in the A–W interface region of cloud/fog drops and aerosol liquid water (ALW). To begin, the three different systems were optimized to use the least amount of energy possible. After that, NVT molecular dynamics simulations were carried out for a total of 150 ns.

2.1.2 Umbrella sampling simulations

In Fig. S1b, the molecule of Ph, 4-HBA, or VL was placed inside the box (their COM is (2.00, 2.00, 6.00 nm)), which is located directly 2.00 nm away from the COM of the water slab. The distance between the COM of Ph, 4-HBA, or VL and that of the water slab was used as the definition for the reaction coordinate (Fig. S1). The weighted histogram analysis approach, also known as WHAM, can be used to calculate the free-energy profiles of Ph, 4-HBA, or VL when they transition from the gas phase into bulk water (Kumar et al., 1992; Hub et al., 2010); details about WHAM are in Text S1 in the Supplement.

2.1.3 Radial distribution function

Estimating the strength of hydrogen bonds (HBs) between specific atoms can be done with the help of a tool known as the radial distribution function (RDF). Text S2 has an explanation of the characteristics of the RDF and the coordination number.

2.2 Density functional theory (DFT) calculations

In this work, all structural optimization and energy calculations were accomplished by utilizing the Gaussian 16 program (Frisch et al., 2016). Calculated at the CCSD(T)/cc-pVDZ, CBS-QB3, B3LYP/6-311+G(d,p), MP2/6-311+G(d,p), and M06-2X/6-311+G(d,p) levels, Cao et al. (2021) found that M06-2X/6-311++G(3df,2p)//M06-2X/6-31+G(d,p) is reliable for PhCs in the gas phase. After analyzing the stability of the wave function, the method we used is reliable. Therefore, all calculations for gas-phase reactions are performed at this level. Text S3 contains a description of the additional calculated details. The frequency correction factor (0.967) has been taken into account. Multiwfn (Lu and Chen, 2012) was used to construct the electron density map. This program integrates Visual Molecular Dynamics (version 1.9.3) (Humphrey et al., 1996) in order to conduct an analysis of the electrostatic potential (ESP) and the average local ionization energy (ALIE).

2.3 IRI analysis

Interaction region indicator (IRI) analysis (Lu and Chen, 2021) was used to determine the chemical bonds and weak interactions of Ph–4-HBA–VL adsorbed to TiO₂ clusters (the details are in Text S4).

2.4 Kinetic calculations

Text S5 contains an explanation of the kinetic calculation methods.

3 Results and discussion

3.1 Enrichment of Ph, 4-HBA, and VL at the A–W interface

3.1.1 The uptake of gaseous PhCs at the A–W interface

Figures S1 and S2 illustrate the relative distributions of water, O₃, and PhC molecules (Ph, 4-HBA, and VL) in the A–W interface system along the z axis. Hydroxyl radicals (HO^•) are primarily situated at the A–W interface contact, with the potential to diffuse through the water slab interior (Roeselová et al., 2004). Figure 1a displays the variation in water density along the z-coordinate distance from 0 to 9 nm, categorizing three zones: the A–W interface (2.25 to 2.79 and 6.21 to 6.75 nm), air (0 to 2.25 and 6.75 to 9 nm), and bulk water (2.79 to 6.21 nm). This method accurately determines the interfacial range (Zhang et al., 2019; Shi et al., 2020). According to location definitions, the O₃ percentage distribution was as follows: 26 % at the A–W interface, 72 % in the air, and 2 % in bulk water (Fig. 1b). Figure 1c depicts molecular dynamics (MD) trajectories of Ph diffusion through the water slab from the air region over a 150 ns period. Ph is distributed in the air (8 %) and bulk water (20 %), with the majority at the A–W interface (72 %) (Fig. 1d). The majority of 4-HBA and VL molecules are located at the A–W interface, constituting 68 % and 73 % of the total locations as presented in Fig. S2.

Figure 1 (a) Relative concentration distributions in the A–W system along the z axis; (b) probability of O₃ at the A–W interface, in gas phase, and in bulk water; (c) MD trajectories of Ph diffusion through the water slab over a 150 ns period; (d) probability of Ph at the A–W interface, in gas phase, and in bulk water.

In Fig. 2a, we observe the three key processes involving PhCs (Ph, 4-HBA, or VL) diffusing into the water slab from the air region: (i) the mutual attraction of gaseous Ph, 4-HBA, or VL; (ii) the uptake of PhCs (Ph, 4-HBA, or VL) at the A–W interface; and (iii) the hydration reaction of PhCs (Ph, 4-HBA, or VL) in the bulk water. Figure 2b displays the free-energy profile of the trajectories as Ph, 4-HBA, and VL transition from the air into the bulk water (see Text S6 for calculations details). The ΔG_{gas→interface} values are −0.22 kcal mol⁻¹ for the Ph–A–W (phenol–air–water) system, −0.45 kcal mol⁻¹ for the A–W_4-HBA (4-hydroxybenzaldehyde at the air–water interface) system, and −0.20 mol⁻¹ for the A–W_VL (vanillin–air–water) system. This finding is consistent with previous studies about per- and polyfluoroalkyl substances (PFASs) at the A–W interface (Yuan et al., 2023). These values suggest that it is thermodynamically favorable for PhCs to approach the interfacial water molecules. Figure S3 illustrates typical snapshots from the trajectories of PhCs (Ph, 4-HBA, or VL). Initially, one molecule of Ph, 4-HBA, or VL was placed in the center of the water box, with an equivalent COM distance of 2 nm between the PhCs and the air phase. Subsequently, the PhCs moved closer to the interface, leading to adsorption at the A–W interface. During the adsorption process, the H atom of the phenolic hydroxyl group binds to the oxygen atom of the H₂O molecules at the A–W interface, forming H bonds and preventing its return to the bulk water. This property allowed the phenolic hydroxyl groups on PhCs to effectively adhere to the A–W interface, consistent with the experimental observations using steady-state interfacial vibrational spectra (Kusaka et al., 2021) and Fourier transform infrared (FTIR) imaging micro-spectroscopy (Guzman et al., 2022). Based on these findings, compared to the number of PhC molecules distributed in the gas phase and in bulk water, the location where air and water meet exhibits an increase in the number of PhC molecules.

Figure 2 (a) Three key processes for the reaction of gaseous PhCs (Ph, 4-HBA, or VL) with the water drops; (b) free-energy change profile of gaseous PhCs (Ph, 4-HBA, or VL) approaching the bulk water.

3.1.2 Interface properties of PhCs

Introducing more hydrophilic functional groups increases the characteristic angle α and β of PhCs at the interface, allowing for more secure adsorption at the water–air interface. The interaction between H_PhCs and O${}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ is the primary factor influencing the stability of PhCs at the interface. The coordination numbers (N) of H_Ph-OH−O${}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$, H_4-HBA-OH−O${}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$, and H_VL-OH−O${}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ are 2.68, 2.51, and 2.09, respectively. The number of functional groups attached to the benzene ring affects the N value; a higher number of functional groups leads to a lower N value. The reason for this is that aldehyde and methoxy are strong electron-withdrawing groups, which will reduce the conjugation effect between the benzene ring and the hydroxyl group and make the hydrogen atom on the hydroxyl group partially positively charged, thus weakening the hydrogen bonding ability with water molecules. See Text S7 for the interface properties of PhCs.

3.2 Adsorption of Ph, 4-HBA, and VL by TiO₂ clusters

The placement of PhCs on TiO₂ clusters significantly impacts adsorption energies (Bai et al., 2020). The adsorption capacity of pollutants on cluster surfaces is a key factor influencing degradation efficiency (Qu and Kroes, 2006). The primary mechanism of C atom adsorption to (TiO₂)_n (n = 1–4) clusters occurs at a range of 2.57 to 2.61 Å and involves interaction between the H_−OH atom and the O${}_{{\mathrm{TiO}}_{\mathrm{2}}}$ atom, as seen in Fig. 3a. Hydrogen bonds can be formed between the H_−OH atom and the O${}_{{\mathrm{TiO}}_{\mathrm{2}}}$ atom (1.80–2.61 Å), improving the adsorption capacity. In contrast, Ph adsorption to (TiO₂)_n (n = 5–6) clusters, ranging from 2.08 to 2.09 Å, is primarily due to interaction between the Ti atom and O_−OH atom. For a detailed description, see Text S8.

Figure 3 Adsorption details of PhCs on (TiO₂)_n (n = 1–6) clusters: (a) structure of Ph adsorption on the (TiO₂)_n (n = 1–6) surface; adsorption energy of (b) Ph, (c) 4-HBA, and (d) VL on (TiO₂)_n (n = 1–6, unit: kcal mol⁻¹); and (e) interaction region indicator (IRI) analyses of Ph on the (TiO₂)_n (n = 1–6) surface.

Adsorption energy a metric of adsorption capacity, is illustrated in Fig. 3b–d for Ph, 4-HBA, and VL on (TiO₂)_n (n = 1–6). TiO₂ exhibits the highest adsorption capacity for Ph ($\mathrm{\Delta }{G}_{\text{ad}}=-$72.35 kcal mol⁻¹) (Fig. 3b). The adsorption energy values of TiO₂ and (TiO₂)₃ for 4-HBA and VL are −45.32 (Fig. 3c) and −102.46 kcal mol⁻¹ (Fig. 3d), respectively. Physisorption energy ranges from −1.20 to 9.56 kcal mol⁻¹ (Nollet et al., 2003); thus the adsorption process in this study is spontaneous chemical adsorption. However, the capacity of TiO₂ to adsorb VL is significantly higher than that to adsorb Ph and 4-HBA. Figure 3b–d show that the adsorption capacity falls as the size of TiO₂ clusters increases when n ≤ 4. In contrast, the adsorption capacity remains constant when n > 4. IRI measurements of Ph on the (TiO₂)_n surface (Fig. 3e) reveal Ph–TiO₂ hydrogen bonds (H_Ph–O${}_{{\mathrm{TiO}}_{\mathrm{2}}}$ bonds) and their electrostatic and dispersion effects. Benzene C atom of Ph exhibits sp² hybridization, meaning it forms one σ bond and one π bond. The sp² hybridization of benzene for Ph explains its limited interaction with TiO₂ clusters and accounts for the substantial adsorption energy. Similar interactions occur with 4-HBA and VL (Fig. S7). Hydrogen bonds form between the H_−CHO atom of 4-HBA or VL and the O${}_{{\mathrm{TiO}}_{\mathrm{2}}}$ atom, despite the presence of the H_Ph atom.

3.3 Continuous-oxidation mechanisms

3.3.1 O₃- and HO^•-initiated reactions

PhCs, once released into the atmosphere, undergo several processes, including adsorption on mineral aerosol surfaces, accumulation at the A–W interface, dispersion in bulk water within liquid droplets, and oxidation reactions initiated by atmospheric oxidants (Lin et al., 2017). Reactions inside the aqueous particle have also been a hot topic of interest in recent years (Tilgner et al., 2021; Go et al., 2023; Zhang et al., 2024; Rana et al., 2024), so we also focused on the process by which phenolic compounds enter the interior of the droplet. This section delves into the detailed mechanisms and characteristics of these reactions. At the M06-2X/6-311++G(3df,2p)//M06-2X/6-31+G(d,p) level, the structures with the minimum free energy for Ph, 4-HBA, and VL have been determined (Fig. S10). In the case of VL, a significant reduction in molecular energy is observed due to the formation of a powerful intramolecular hydrogen bond with a length of 2.09 Å between the H and O atoms near the methyl group. Moreover, the lone-pair electrons of oxygen atoms can additionally form p–π conjugations with the π electrons of the phenyl ring, further reducing the overall energy of VL in gas phase. The statistical charts of calculated Δ_rG and ΔG^‡ values for O₃- and HO^•-initiated reactions are displayed in Figs. 4 and S8, and further details are available in Tables S1–S4​​​​​​​.

Figure 4 Statistical charts of calculated (a–d) Δ_rG and (e–h) ΔG^‡ values for O₃-initiated reactions and (i–l) Δ_rG and (m–p) ΔG^‡ values for HO^•-initiated reactions.

O₃ is a major oxidant in the atmosphere, with high concentrations in the troposphere of 9.85 × 10¹¹ molec. cm⁻³ (Tomas et al., 2003; Pillar-Little et al., 2014). Investigating the fate of PhCs in the presence of O₃ is essential (Pillar-Little et al., 2014; Rana and Guzman, 2020). The ozonolysis of PhCs involves the synthesis of primary ozonide, the formation of active Criegee intermediates (CIs), and the disintegration of CIs (Rynjah et al., 2024). The O₃-initiated reactions of Ph, 4-HBA, and VL involve radical adduct formation (RAF) channels on their benzene ring ($R_{{\mathrm{O}}_{\mathrm{3}}-\mathrm{RAF}}$1–$R_{{\mathrm{O}}_{\mathrm{3}}-\mathrm{RAF}}$6), highlighted in red in Fig. S10. Figure 4a–d show that the ozonolysis pathways $R_{{\mathrm{O}}_{\mathrm{3}}-\mathrm{RAF}}$ are exergonic, indicating their spontaneity. The average ΔG^‡ values for the ozonolysis of Ph, 4-HBA, and VL are ranked as Ph > VL > 4-HBA. The following is a list of the average values for the ozonolysis of Ph, 4-HBA, and VL, as illustrated in Fig. 4e–h: Ph is superior to VL and 4-HBA, with the exception of TiO₂ clusters. Figure 4e illustrates that the average value of ΔG^‡ for O₃ + Ph reactions at the A–W interface is 15.38 kcal mol⁻¹, the lowest value out of the three PhCs. The average ΔG^‡ values for the ozonolysis of Ph, 4-HBA, and VL are as follows: VL (13.95 kcal mol⁻¹) < Ph (24.70 kcal mol⁻¹) < 4-HBA (25.16 kcal mol⁻¹) on TiO₂ clusters (Fig. 4f). The average ΔG^‡ values for O₃ + VL reactions in the gas phase are the highest among the four different environmental media (23.28 kcal mol⁻¹) shown in Fig. 4g. Comparing the phenolic oxidation in each of these four environmental media (bulk water, interface, TiO₂ clusters, and gas phase) reveals that the A–W interface is more conducive to the ozonolysis of Ph, whereas TiO₂ clusters are more conducive to the ozonolysis of VL. The effect of solvation on ΔG^‡ is predominantly caused by the hydration of the phenolic OH group, as this is the part of the molecule being solvated. However, the presence of water molecules in the region around the phenyl group has been shown to have a considerable influence on the ΔG^‡ values.

Figure 5 Δ_rG and ΔG^‡ values of (a–b) O₃-initiated reactions and (c–f) HO^•-initiated reactions at different reaction positions.

HO^•, known as the “atmospheric detergent”, is another significant atmospheric oxidant (Atkinson, 1986; Zhang et al., 2020). The worldwide mean tropospheric concentration of HO^• is roughly 11.3 × 10⁵ molec. cm⁻³ (Lelieveld et al., 2016). For this reason, elucidating the reaction mechanism underlying HO^• + PhCs reactions in the troposphere is of the utmost importance. HO^•-initiated reaction pathways of Ph, 4-HBA, and VL include RAF, hydrogen atom abstraction (HAA) channels from the benzene ring (R_HAAben1–R_HAAben6), and the substituent group (R_HAAsub7–R_HAAsub9). Previous research (Gao et al., 2019) has shown that the process of single electron transfer (SET) does not significantly contribute to the HO^•-initiated reactions examined. Once the hydroxyl adducts or H₂O are formed, significant heat (4.21–30.28 kcal mol⁻¹) is released (Figs. 4i–l, S8a–d and i–l; detailed data in Table S3), indicating high thermodynamic feasibility. The average ΔG^‡ values for HO^•-initiated reactions (Figs. 4m–p, S8e–h and m–p) are lower than those for O₃-initiated reactions. Routes R_HAAben make a minimal contribution to HO^•-initiated reactions. At the A–W interface, the interaction energies are in the following order: VL (3.52 kcal mol⁻¹) < Ph (4.52 kcal mol⁻¹) < 4-HBA (9.50 kcal mol⁻¹). Additionally, the ΔG^‡ value for Ph is the lowest at −0.97 kcal mol⁻¹ in the pathways $R_{\mathrm{RAF}-{\mathrm{HO}}^{\mathrm{•}}}$ (Fig. 4m). Among the three aromatic compounds, the $R_{\mathrm{RAF}-{\mathrm{HO}}^{\mathrm{•}}}$ routes of VL on TiO₂ clusters are most likely to take place (Fig. 4n). When compared to HO^•-initiated reactions of aromatic compounds in the gas phase (Fig. S9e) or bulk water (Fig. S9f), the process of Ph + HO^• reactions at the A–W interface is accelerated, whereas the process of VL + HO^• reactions is accelerated by TiO₂ clusters. These findings are in agreement with the ozonolysis findings. The same guidelines can be used for routes R_HAAsub (Figs. 4o, p and S8g, h) and R_HAAben (Fig. S8m–p). The following is a ranking of the average ΔG^‡ values for routes R_RAF−HO• in the gas phase or bulk water: Ph < 4-HBA < VL. As a result of having the lowest ΔG^‡ values among all HO^•-initiated reaction mechanisms, routes R_RAF are the most advantageous of all the possible reaction mechanisms. In light of this, each and every route of R_RAF−HO• and R_HAAsub will be dissected in detail.

Figure 5 shows the Δ_rG and ΔG^‡ values of O₃- and HO^•-initiated reactions at various reaction locations. These reactions are almost entirely exothermic, with a close correlation between Δ_rG values and ΔG^‡ values. The ΔG^‡ values for the Ph + O₃ reactions shown in Fig. 5a are the lowest among the three compounds, ranging from −0.97 to 7.86 kcal mol⁻¹. Exergonic and spontaneous addition reactions took place at the C1–C2 and C3–C4 locations of Ph and VL, respectively. Because of their low ΔG^‡ values, the C1–C2 and C2–C3 sites of O₃-initiated reactions for 4-HBA are advantageous. Their values are 21.76 and 22.03 kcal mol⁻¹, respectively. The C1–C2 location of 4-HBA is activated to a greater extent at the A–W interface in comparison to the gas phase and bulk water. However, the ΔG^‡ values of O₃ + Ph reactions on TiO₂ clusters are significantly greater than those of the A–W interface (12.86–18.10 kcal mol⁻¹) at 24.30–25.34 kcal mol⁻¹. The VL + O₃ reactions on TiO₂ clusters are favorable at the C2–C3 and C4–C5 locations (the ΔG^‡ values are 11.42 and 11.14 kcal mol⁻¹, respectively, Fig. 5b). This can be explained by the fact that the electron cloud has a greater propensity to congregate in C2–C3 and C4–C5, respectively. In addition, the p orbitals of the methoxy and hydroxy groups are conjugated to the benzene ring, which offers a powerful electron-donating conjugation effect (Aracri et al., 2013). Because of this, the oxidation of aromatic molecules is thermodynamically more favorable than the oxidation of the aldehyde group. This is consistent with previous studies that electron density influences the oxidative activity of PhCs (Rana and Guzman, 2022a). Clearly, the ΔG^‡ values of HO^•-initiated reactions (−0.97 to 13.46 kcal mol⁻¹) in Fig. 5c–f are lower than those of O₃-initiated processes (11.14 to 27.83 kcal mol⁻¹) at different points in the A–W interface and TiO₂ clusters. This can be seen by comparing the values to each other. At the A–W interface, the most advantageous position for the phenol hydroxyl group to be in for Ph/4-HBA/VL + HO^• reactions is the ortho-position (Fig. 5c). OESI-MS, which stands for online electrospray ionization mass spectrometry, was also able to identify the hydroxylation product known as 3,4-dihydroxybenzaldehyde (Rana and Guzman, 2020). In Fig. 5d, the ortho- and meta-sites of phenol hydroxyl are the most favorable positions for Ph + HO^• and 4-HBA + HO^• reactions, respectively, on the TiO₂ clusters. On the other hand, all of the VL sites on the TiO₂ clusters are advantageous. At the A–W interface and on the TiO₂ clusters, the abstraction of hydrogen atoms follows the order of H_−CHO atom > H${}_{-{\mathrm{OCH}}_{\mathrm{3}}}$ atom > H_−OH atom in Fig. 5e and f. This can also be explained by the ALIE values of these atoms listed in the same order of H_−CHO atom (11.67–11.74 eV) > H${}_{-{\mathrm{OCH}}_{\mathrm{3}}}$ atom (14.06 eV) > H_−OH atom (15.46 eV), as shown in Fig. S8.

Figure 6 Subsequent reaction mechanisms of important intermediates (IMs) (unit: kcal mol⁻¹) (a) in the gas phase (g)/bulk water (l) and (b) at the A–W interface.

3.3.2 Generation and degradation of key products

For the purposes of this discussion, the primary atmospheric fate of the selected aromatics was considered to be their reactions with O₂ (typically mediated by reactive intermediates or catalytic processes) and O₃. Figures 6 and S10 illustrate the subsequent reaction mechanisms of intermediates (IMs). IM_1–2 was produced using the pathway that offered the best conditions for the HO^•-initiated reaction of Ph. As can be seen in Fig. 6a, the addition of O₂ to the C3 sites of the C₆H₅O radicals results in the formation of C₆H₅O-OO radicals with no barriers in either the gas phase or the bulk water. This is a desirable outcome. For the transformation of the C₆H₅O₂-OO radicals that were created, the ring closure reaction to form C₆H₅O₂-OO-d is an attractive option. However, it must overcome an energy barrier of 18.83 kcal mol⁻¹ in the gas phase or 13.67 kcal mol⁻¹ in bulk water. The C₆H₅O₂-OO-d₁ radical, which was produced by the C₆H₅O₂-OO-d reaction, interacts once more with O₂. In the atmosphere, these Criegee intermediates also may undergo bimolecular reactions with NO_x (Sun et al., 2020). Malealdehyde (P1) is what should mostly result from the reaction of the C₆H₅O₂-OO-d₁ radical with NO. However, during this process, it still needs to overcome an energy barrier of 49.5 (in the gas phase) or 50.83 kcal mol⁻¹ (in the bulk water) to generate the C₆H₅O₂-OO-d₃ radical; as a result, the further transformation of the formed C₆H₅O₂-OO-d₂ should continue very slowly. Pyrocatechol (P2) is the primary product generated in the gas phase and bulk water when the H atom of the C₆H₅O₂-OO radical is displaced. P2 generates o-semiquinone radicals via pathways R_HAA by HO^• or O₃, which in turn generate oligomers (Guzman et al., 2022). This results in the formation of brown organic carbon in atmospheric aerosols. At the A–W interface, a sequence of hydroxylation products, including pyrocatechol (P2), benzene-1,2,3-triol (P3), and benzene-1,2,3,4,5-pentaol (P4), are generated through hydroxylation processes rather than by a single SET (ΔG^‡ = 111.79 kcal mol⁻¹). It is difficult for P4 to form benzene-1,2,3,4,5,6-hexaol because it is difficult for hydrogen transfer reactions to occur (ΔG^‡ = 34.32 kcal mol⁻¹). These hydroxylation products have been detected by experimental means (Pillar-Little et al., 2014; Pillar-Little and Guzman, 2017; Rana and Guzman, 2020). The HO^• abstracts a hydrogen atom from the hydroxyl group of catechol, forming a C₆H₅O₂ radical and a water molecule. Due to the widespread presence of NO₂ in the environment, it adds to the C₆H₅O₂ radical at the ortho-position of the extracted hydrogen atom through an addition reaction. Subsequently, a hydrogen transfer reaction occurs, resulting in the formation of 4-nitrobenzene-1,2-diol (P_2-a). This computational result validates the previous experimental hypothesis by Finewax et al. (2018). The P_2-a subsequently transforms into benzoquinone, maleic acid, fumaric acid, acetic anhydride, acetic acid, and formic acid or is directly mineralized into carbon dioxide and water (Chen et al., 2015). In order to gain a more comprehensive understanding of the reaction mechanism at the A–W interface, the major product (the C₇H₅O₂ radical) for pathways R_HAA of 4-HBA was also taken into consideration. According to Fig. S11a, the addition of HO^• to the C7 sites of the C₇H₅O₂ radical can occur without any obstructions. The overpowering of the 18 kcal mol⁻¹ barrier resulted in the formation of the hydroxylation products (4-hydroxybenzoic acid (P5), 3,4-dihydroxybenzoic acid (P6), 2,3,4-trihydroxybenzoic acid (P7), and 2,3,4,5,6-pentahydroxybenzoic acid (P8)). There was found to be one transition route for the continued ozonolysis of the hydroxylation products that were produced in P6. The C2–C3 site of P6 to create P6-5O₃ (ΔG^‡ = 16.59 kcal mol⁻¹) has the lowest activation energy of all the available paths for the relevant reactions (Fig. S11b). This corresponds to a value of 16.59 kcal mol⁻¹. When the ΔG^‡ values of the breakage of five-membered rings created by ozonolysis pathways are compared, one can reach the conclusion that the formation of IM_P6-5O₃-a is the most favored pathway. All of the hydrogen abstraction processes involving H₂O and IM_P6-5O₃-a have rather high energy barriers (32.93 kcal mol⁻¹). On the other hand, in Fig. S10a, the very low ΔG^‡ values (19.74–22.89 kcal mol⁻¹) of the –NO–O abstraction make it a desirable choice. Following a chain of ozonolysis reactions, the following products were obtained (Fig. S11c): ((2E,4Z)-2-formyl-4,5-dihydroxy-6-oxohexa-2,4-dienoic acid (P9), 2,3-dihydroxymalealdehyde (P10), and 2,3-dioxopropanoic acid (P11). Therefore, the product that was created, P10, may also be the product that was discovered through experimentation (mass-to-charge ratio ($m/z$) of 115) (Rana and Guzman, 2020). The VL subsequent reaction mechanism is demonstrated in Fig. S11d. The final oxidation products of VL are P12 ((2E,4E)-4-formyl-2-methoxy-6-oxohexa-2,4-dienoic acid), P13 (ethene-1,1,2-tricarbaldehyde), P14 (2-methoxy-2-oxoacetic acid), P15 (oxalaldehyde), and P16 ((E)-2-methoxy-4,5-dioxopent-2-enoic acid). The formation of these products could explain the biomass burning material for the formation of SOA (Rana and Guzman, 2022b).

Table 1 The available experimental and calculated reaction rate constants (k) values of O₃- and HO^•-initiated reactions at 298 K. Unit: cm³ molec.⁻¹ s⁻¹.

^a Calculated values of phenolic compounds at A–W interface. ^b Calculated values of phenolic compounds on TiO₂ clusters. ^c Calculated values of phenolic compounds in the gas phase. ^d Calculated values of phenolic compounds in the bulk water. ^e Experimental values of catechol in the gas phase. ^f Experimental values of guaiacol in the gas phase. ^g Experimental average $k_{{\mathrm{HO}}^{\mathrm{•}}}$ values of methoxyphenols in the gas phase.

3.4 Comparison with available experimental results

The rate constants (k) of the overall reaction under the temperature range of 278–318 K were computed based on acquired potential energy surfaces for the O₃-initiated and HO^•-initiated reactions of selected compounds. The results of these calculations are listed in Tables S5 and S6, respectively. The temperature dependences of the various k values for Ph, 4-HBA, and VL at the A–W interface and in bulk water are depicted in Fig. 7. At low values of k, there is a positive dependence on temperature. When the k values are raised to a certain degree, the temperature dependency seems to lose any significance it may have had before. The following is an order of the k values for O₃-initiated reactions: A–W_Ph > TiO_2 VL> A–W_VL > A–W_4-HBA > TiO_2 4-HBA > TiO_2 Ph (Fig. 7a). According to Fig. 7b, the k values of HO^•-initiated reactions are ordered as follows: TiO_2 VL > A–W_Ph > A–W_VL > TiO_2 4-HBA > TiO_2 Ph > A–W_4-HBA. In Fig. 7a and b, the k values of HO^•-initiated reactions are 100 times greater than those of O₃-initiated reactions. Table 1 is a listing of the experimental and estimated k values that are available for O₃-initiated and HO^•-initiated reactions at 298 K. According to the findings, the ozonolysis of Ph was promoted by the water–gas interface as well as by TiO₂ clusters, and the HO^•-initiated reactions of VL were promoted by TiO₂ clusters. However, the O₃/HO^• + 4-HBA reactions have the lowest k values among the three molecules when tested in a variety of environments. The estimated $k_{{\mathrm{O}}_{\mathrm{3}}+\mathrm{Ph}}$ values at the A–W interface are 11 orders of magnitude greater than those of catechol under dry conditions in the gas phase (Zein et al., 2015) when compared with the experimental data. Because it has a higher $k_{{\mathrm{O}}_{\mathrm{3}}}$ value, catechol, which is one of the main products of Ph's oxidation in the atmosphere, has a higher degree of reactivity than its parent compound (Table 1). The estimated value of VL is lower than the experimentally determined value of $k_{{\mathrm{O}}_{\mathrm{3}}}$ for guaiacol under dry conditions, which is (0.40 ± 0.31) × 10⁻¹⁸ cm³ molec.⁻¹ s⁻¹ in the gas phase (Zein et al., 2015). The difference between the predicted value of $k_{{\mathrm{HO}}^{\mathrm{•}}+\mathrm{VL}}$ is 1.14 × 10⁻¹⁰ cm³ molec.⁻¹ s⁻¹, and the average experimental value of $k_{{\mathrm{HO}}^{\mathrm{•}}}$ for methoxyphenols is just an order of magnitude. As a consequence, the findings of our calculations are reliable. Previous studies measured the second-order rate constants of guaiacylacetone + HO^• reaction to be (14–25) × 10⁹ M⁻¹ s⁻¹ at pH 5 and 6 for aqueous secondary organic aerosol, which is lower than our results (Arciva et al., 2022). This is because galactose reduces the steady-state concentration of HO^•. The reaction rate constants of PhCs increase with increasing pH, and we calculated the rate constants at pH 7 in bulk water (Ma et al., 2021). This study summarizes the O₃- and HO^•-initiated reaction sequences of three PhCs in different environmental media. The reaction sequences for O₃- and HO^•-initiated reactions of Ph and 4-HBA are identical in different environmental media, while VL shows slight variations. For O₃-initiated reactions, the reaction sequences are as follows: for Ph, A–W interface > bulk water > gas phase > TiO₂ clusters; for 4-HBA, bulk water > A–W interface > TiO₂ clusters > gas phase; and for VL, bulk water > TiO₂ clusters > A–W interface > gas phase. For HO^•-initiated reactions, the sequences are, for Ph, A–W interface ≈ bulk water > gas phase > TiO₂ clusters; for 4-HBA, bulk water > A–W interface > TiO₂ clusters > gas phase; and for VL, TiO₂ clusters > bulk water > A–W interface > gas phase. According to the atmospheric concentration of O₃, the atmospheric lifetime of Ph is the shortest (< 1 s) of the three PhCs at the gas–water interface or in bulk water, whereas 4-HBA and VL were oxidized more slowly than Ph (Smith et al., 2016). See Fig. S12 and Text S9 for ecotoxicity assessment.

Figure 7 Calculated rate constants for the initial reactions of Ph, 4-HBA, and VL with O₃ and HO^• at different temperatures (278–318 K) and 1 atm.

4 Conclusions

Combining molecular dynamic simulations (with the AMBER force field) and quantum chemical calculations (at the M06–2X/6–311++G(3df,2p)//M06–2X/6-31+G(d,p) level) methods has provided comprehensive insights into the surface properties of Ph, 4-HBA, and VL, as well as their reactions induced by O₃ and HO^•, in both homogenous and heterogenous environments. Here are some key findings from this research:

1. The free-energy well of Ph, 4-HBA, and VL favors the A–W interface as the compounds' preferred location, with occurrence percentages of approximately ∼ 72 %, ∼ 68 %, and ∼ 73 %, respectively. Ph and 4-HBA show a preference for the A–W interface over the air, with an energy difference of around 0.22 and 0.45 kcal mol⁻¹. The VL adsorbed on the TiO₂ clusters has a higher likelihood of remaining compared to VL adsorbed at the A–W interface.

2. The adsorption capacity of TiO₂ clusters decreases with increasing cluster size until n > 4. After that point, the adsorption capacity remains constant. Strong electrostatic attractive interactions and attractive dispersion effects occur between the benzene of the Ph and Ti atoms. Hydrogen bonds form between the atom of O${}_{{\mathrm{TiO}}_{\mathrm{2}}}$ and the H_−CHO group of 4-HBA or VL.

3. The O₃- and HO^•-initiated reactions for Ph and VL are facilitated by the A–W interface and TiO₂ clusters, respectively. For O₃-initiated reactions at the A–W interface, the C1–C2 position on the benzene ring is most favorable. Both in the A–W interface and on TiO₂ clusters, the total branching ratio for routes R_RAF and R_HAAsub is 72.68 %–100 %. For route R_HAAsub, the order is H_−CHO atom > H${}_{-{\mathrm{OCH}}_{\mathrm{3}}}$ atom > H_−OH atom.

4. The k values (in molec. cm⁻³ s⁻¹, at 298 K and 1 atm) of O₃-initiated reactions follow the order of A–W_Ph (5.98 × 10⁻⁷) > TiO_2 VL (3.30 × 10⁻¹⁵) > A–W_VL (1.27 × 10⁻¹⁷) > A–W_4-HBA (6.79 × 10⁻²³) > TiO_2 4-HBA (5.32 × 10⁻²⁴) > TiO_2 Ph (1.84 × 10⁻²⁴). The k values of HO^•-initiated reactions follow the order of TiO_2 VL (6.70 × 10⁻⁶) > A–W_Ph (2.69 × 10⁻⁶) > A–W_VL (1.73 × 10⁻⁷) > TiO_2 4-HBA (3.16 × 10⁻⁹) > TiO_2 Ph (3.17 × 10⁻¹⁰) > A–W_4-HBA (9.49 × 10⁻¹¹).

5. Toxicity risk assessment on aquatic species reveals that most of the reaction products are significantly less harmful than the parent compounds. However, products P1, P2, P3, P10, and P11 are more hazardous, and further investigation of their atmospheric fate is recommended.

Ph undergoes transformation to malealdehyde and catechol when exposed to O₃ or HO^• in the troposphere (Xu and Wang, 2013). When Ph $/$ VL is at the droplet aerosol interface, rapid oxidation to polyhydroxylated compounds occurs (Ma et al., 2021). VL eventually creates tiny molecule aldehydes and acids. This is consistent with experimental observations (Rana and Guzman, 2020). This oxidation process is accelerated when VL is encased in a mineral aerosol represented by TiO₂ clusters. Li et al. (2022)​​​​​​​ found that seasonal average concentrations of total nitrophenol compounds in particulate matter were comparable to those measured in the gas phase. However, the reactivity order of nitrophenols in the atmospheric compartments is water droplets > gas phase > particles (Vione et al., 2009). The formation of some low-molecular-weight acids and aldehydes (2,3-dihydroxymalealdehyde, 2,3-dioxpropanoic acid, etc.) confirms their association with the formation of SOA. It is recommended that enterprises producing lignin, such as those in the pulp and paper industry, or factories that employ lignin in the manufacturing of adhesives, rust inhibitors, color dispersants, diluents, or other similar products be constructed in regions with low relative humidity. It is recommended that treatment facilities that collect lignin pyrolysis products and recycle the byproducts be located in the surrounding area.