Estimating rate coefficients for the reactions of ethers + OH from atmospheric to combustion temperatures: an extension of the electrotopological state method

Abstract

The reactions of ethers with OH present a challenge to estimation techniques. These hydrogen abstraction reactions exhibit both positive and negative temperature dependences according to temperature regime and degree of halogenation. Positive temperature dependence relates to reactants overcoming energy barriers in the primary process of hydrogen abstraction, and is readily described using a recently developed structure–activity relationship based upon the electrotopological state (E-state). Conversely, we attribute the negative temperature dependence to the formation of van der Waals (vdW) pre-reactive complexes that facilitate reactions through quantum tunnelling, for which we have developed a simple parameterization. This is implemented in the current version of the E-state code in a modular fashion, leaving the original parameterization for alkanes and haloalkanes unmodified. We find that this treatment combined with the a priori information provided by the E-state affords, to some extent, the decoupling of these opposing behaviours leading to a robust prediction of temperature-dependent rate coefficients with OH, k(T), compared with existing methodologies. This approach, named E-stateXvdW, provides site-specific, temperature-dependent estimates implemented in a freely available open-source Python programme that can estimate large batches of compounds automatically and at negligible computational cost over a broad temperature range (∼190–1500 K). This method applies to all acyclic non-halogenated and halogenated ethers of composition C_nO_≤n−1H_mX_2n+2−m, and provisionally, all cyclic ethers of composition C_nO_≤n−1H_mX_2n−m, where X represents any combination of F, Cl, Br and I, n = [1, 2, …] and m = [1, 2, …, 2n + 2]. This approach may be useful to modellers and environmental scientists who wish to assess the fate of industrially important etheric molecules, as well as those chemists with an interest in rationalizing the reactivity of molecules based on their structural characteristics.

---

Environmental significance

Ethers have many industrial and commercial uses and can enter both the atmospheric and combustion environments where they will be further oxidized. Ether oxidation is initiated primarily by the hydrogen-abstraction reactions of the OH radical. The rates of these reactions and their temperature-dependences are strongly affected by substitution patterns within molecules. Accordingly, we develop a site-specific temperature-dependent structure–activity relationship, based on the electrotopological state and a purpose-built parameterization of van der Waals complex formation, that is able to capture the behaviour of hydrogen-abstraction reactions of ethers between ∼190 and 1500 K. We implement this algorithm into an automated open source Python programme that allows atmospheric and combustion scientists to estimate the reactivity of these molecules with convenience.

Introduction

Ethers have a variety of industrial and commercial applications that will ultimately result in their release into the environment. Indeed, as societies transition from fossil fuel feedstocks towards more carbon neutral biomass alternatives, this will provide convenient access to a range of ethers and other oxygenated molecules.^1,2 This, combined with their favourable material properties has led to interest in the use of ethereal compounds as fuels, lubricants, solvents and other specialty chemicals.^2–4 Halogenated ethers are an industrially important subset of the ethers that have several specialist applications such as vapour degreasing, foam blowing, refrigeration and dry etching,⁵ in which roles they may replace ozone-depleting substances such as the chlorofluorocarbons and the hydrochlorofluorocarbons.

Because of their volatility, ethers readily enter the atmosphere, where their main loss process will be reaction with the OH radical. The site-specific kinetics of these reactions will determine their longevity and the secondary products that are formed, defining their environmental impact. Although ethers are a class of compound whose oxidation kinetics are relatively well-characterized,⁶ their site-specific reactivity is less understood. Furthermore, reactions with OH play a major role in combustion systems, and therefore, it is also necessary to determine the kinetics and mechanisms of these reactions at high temperature such that their suitability as fuels may be assessed.^7–9

Under both atmospheric and combustion conditions, ethers are expected to undergo hydrogen abstraction by reaction with OH radicals, leading to a water molecule and a carbon-centred radical that will rapidly form a peroxy radical upon reaction with molecular oxygen. Ethers with several reaction sites will lead to a variety of products, each of which will play roles as secondary pollutants in the atmosphere, or as reactive intermediates and products in combustion environments. An example is given in Fig. 1, which summarizes some of the potential differences between these environments, although it is noted that theory suggests that some of the hydroperoxide-forming hydrogen shift reactions of importance to combustion systems may also be significant under atmospheric conditions.¹⁰

Fig. 1 A simplified scheme of the different oxidation pathways that are available in the reaction of an asymmetrical ether, methyl ethyl ether, initiated by OH radicals in the presence of NO_x. The number of possible oxidation products grows rapidly with the number of unique hydrogen abstraction sites, affecting secondary pollution and combustion characteristics alike. For clarity, stable end products have been outlined.

Where experimental data is lacking, it is advantageous to estimate rate coefficients such that their environmental impact can be properly assessed. Several structure–activity relationship (SAR) methods are available for this purpose, each based on the group-additivity method of Atkinson.¹¹ Among these is the work of Kwok and Atkinson (1995),¹² which is implemented in EPA software package, EPI Suite;¹³ Jenkin et al. (2018),¹⁴ which provides estimates for the automated atmospheric chemical mechanism generation software, GECKO-A;¹⁵ and Carter (2021),¹⁶ which forms the basis for estimates in the SAPRC model.¹⁷ Kwok and Atkinson (1995)¹² can be applied to all ethers of CHOX composition (where X = F, Cl, Br and I) over a broad temperature range (∼230–1000 K); Jenkin et al. (2018)¹⁴ is restricted to CHO composition over a narrower temperature range (∼240–400 K); and Carter (2021)¹⁶ can be used to estimate rate coefficients of CHO ethers over a small range in temperature (270–330 K). At the time of writing, EPI Suite represents the only publicly available automated method to provide estimates of the halogenated ethers based on the approach of Kwok and Atkinson (1995).^12,13

Clearly, the extent to which an SAR covers the parameter space of OH + ether reaction kinetics varies between these approaches. A given SAR might reproduce its training set more precisely by narrowing its domain of applicability in both temperature range and composition. However, this may come at the cost of generality and predictive qualities. The merits and demerits of such approaches shall be discussed below.

The method described in the current work is a modular extension of the E-states approach that was recently applied to alkanes and haloalkanes.¹⁸ This extended method, named E-stateXvdW, now applies to all saturated etheric compounds of CHOX composition over a large temperature range of ∼190 to 1500 K.

Methodology

The E-state represents a fundamental property that describes substituent effects upon hydrogen environments in molecules, thus accounting for the electronic character of individual atoms within these molecules and their respective distances from one another. In doing so, the E-state therefore provides algorithmically generated a priori information on the electronic properties of every substitution within a given compound. This method has been demonstrated to provide accurate estimates using an efficient, automated tool, relying on fewer fitting factors than existing methods. The details to calculate the unmodified rate coefficient, k_base, using the original E-state approach for saturated alkanes and haloalkanes can be found in the previous work of McGillen et al. (2024).¹⁸

To test the performance of this method, a dataset of ethers and haloethers was compiled from the kinetic database of McGillen et al. (2020),⁶ together with some additional references that were either too recent or found to be missing from the original database.^4,19–21 These data are summarized in Table 1 and provided in an Excel file available with this paper (see SI). These include total rate coefficients and (where available) temperature dependence parameters for a variety of cyclic, acyclic substituted (i.e. of CHOX composition) and unsubstituted ethers (i.e. of CHO composition). Note that the Br-containing compound of the dataset represents a measurement of a single species (1-bromo-2-ethoxyethane) and that there are no data for I-containing ethers. However, the inclusion of Br- and I-containing ethers is justified on the grounds that the E-state method of McGillen et al. (2024)¹⁸ was previously trained on these atom types in haloalkanes, and is extended to ether functionalities in this work without further modification. This illustrates one advantage of this approach, which, by calculating a priori electronic effects, can provide estimates beyond the immediate experimental training set.

Table 1 Number of experimental rate coefficients available for each type of ether (where X = F, Cl, Br), from the most recent AtmVOC database V3.1.⁶ (see SI-1) unless otherwise specified.^4,19–21 The training set is composed exclusively of acyclic compounds, while estimates are provided for both acyclic and cyclic ethers

Number of molecules in the dataset
Substitution types	Acyclic	Cyclic	Total
CHO	38	19	57
CHOX (X = F, Cl, Br)	64	2	66
Total	102	21	123

---

The present treatment assumes that the overall temperature-dependent rate coefficient, k_total(T) in cm³ molecule⁻¹ s⁻¹, for hydrogen abstraction reactions of OH + etheric molecules can be described using the following modified Arrhenius expression:

where i refers to each abstraction site in the molecule, κ represents the effective tunnelling coefficient, defined in eqn (2), A is the pre-exponential factor, B is the ratio of the activation energy (in J mol⁻¹) to the gas constant (in J mol⁻¹ K), and n is the curvature term associated with the i^th respective hydrogen abstraction site. Except for κ, the remaining terms of this equation can be estimated using a previous parameterization based on the electrotopological state, a full description of which can be found in McGillen et al. (2024), resulting in the base rate coefficient, k_base,i.¹⁸ This parameterization was trained on a diverse set of alkanes and haloalkanes, and therefore describes a variety of substitution effects, spanning the electron-donating to the strongly electron-withdrawing. We hypothesize that this provides a reasonable estimate of the reactivity of each hydrogen abstraction site in the absence of a significant quantum tunnelling effect. The parametrization in eqn (1) was trained and tested against all experimental data available for etheric molecules as presented in Tables 1 and SI-1.

Fig. 2 shows the estimates made using the unmodified room-temperature parameterization (k_base,298) of McGillen et al. (2024), from which it is evident that it is systematically underestimated for the majority of etheric molecules for which data are available (notable exceptions include subsets of the cyclic ethers and highly fluorinated ethers). This general behaviour is hardly surprising, since oxygenated molecules are considered to form van der Waals pre-reactive complexes that serve to increase the efficiency of hydrogen abstraction reactions through quantum tunnelling,^22–26 suggesting that it is also necessary to parameterize κ_i in order to correctly parameterize k_total(T) in eqn (1).

Fig. 2 Estimated room-temperature rate coefficients, k₂₉₈, for the reaction of OH with etheric molecules using the E-state method from McGillen et al. (2024) compared with experimental data. Data are taken from the database of McGillen et al. (2020),⁶ unless otherwise specified. Circles and triangles represent acyclic and cyclic ethers respectively. Data are colour coded according to their composition. This treatment tends to systematically underestimate the majority of rate coefficients for this set of molecules, with the exception of cyclic and some highly fluorinated molecules.

To achieve this, we first considered the possible van der Waals complexes that can form in ether + OH reactions, with examples shown in Fig. 3. Similar to Smith and Ravishankara (2002),²⁷ we considered 5- and 6-membered cyclic structures, the rationale being that smaller ring sizes would be too strained and larger ring sizes would have a higher entropic cost, both of which are expected to inhibit the formation/stability of a pre-reactive complex. As shown in Fig. 3, the number of possible vdW complexes increases with the size and complexity of a given ether.

Fig. 3 The various van der Waals pre-reactive (vdW) complexes that are considered in the parameterization of the effective tunnelling parameter, κ. Monocyclic structures with one ether bridge, e.g. (1)A, (1)B, (1)C, (2)A, (2)C, and (2)E, and polycyclic structures with two ether bridges, e.g. (2)B, (2)D, and (2)F are included in the current version. Relatively simple structures such as methyl ethyl ether (1) yield fewer possible vdW complexes compared with larger, more oxygenated molecules such as diisopropoxymethane (2). The structures include 5-membered, e.g. (1)A, (1)B and 6-membered, e.g. (1)C cyclic vdW complexes.

This is visualized further in Fig. 4, showing how this parameterization of vdW complexation is affected by carbon number and oxygen content. Note, however, that only the complexes where the OH is associated with the H being abstracted are considered to enhance the abstraction rate coefficient. For example, in the simple case of (1), each of the three sites of OH attack is enhanced by one (and only one) of the vdW complexes shown.

Fig. 4 The dependence of the number of possible unique van der Waals complexes considered upon carbon and oxygen number for a fully enumerated list of acyclic alkyl ethers, C_nO_≤n−1H_2n+2, where 2 ≤ n ≤ 8, generated using the MOLGEN 5.04 software package.²⁸ As carbon number increases, the maximum number of van der Waals complexes increases rapidly. To reduce the complexity of our algorithm, we prioritize certain structures as described in the text. The same figure is presented for the complex types shown in Fig. 3 (cyclic and bicyclic 5-membered and 6-membered van der Waals complexes), in SI-2.

In the present parameterization, κ_i is calculated for the i^th reactive site (positioned ≤2 carbon atoms away from the ether bridge) and the j^th oxygen bridge in the molecule, using the following eqn (2):

where n represents the number of complexes involved, possessing a value of ≤4, α_j and β_j are fitting factors associated with the j^th vdW complex configuration (one of the six possible 5-membered or 6-membered complexes forming on the primary, secondary or tertiary abstraction sites), and S(O_j) is the E-state of the j^th oxygen bridge that participates in the vdW complex, as calculated in McGillen et al. (2024), using the additional intrinsic value of I(–O–) = 3.500 taken from Hall et al. (1991).²⁹ A total of four corrections are considered, with a maximum of two per 5- or 6-membered complex configuration. Priority is given to: (1) bicyclic formations over monocyclic ones; and (2) among configurations exceeding two possibilities, to those involving oxygen atoms with the highest E-state value. If a reaction site i is involved in more than two 5-membered complexes and two 6-membered complexes, the resulting κ_i value would be: κ_i = κ(5-membered, highest E-state) + κ(5-membered, second highest E-state) + κ(6-membered, highest E-state) + κ(6-membered, second highest E-state).

It must be noted that in the case of alkanes and haloalkanes, the original approach of McGillen et al. (2024) is applied,¹⁸ which corresponds to κ being equal to unity in all cases, this is also the case for atoms at a distance of >2 carbon atoms away from the ether bridge(s). We refer to κ as an “effective” tunnelling parameter in the sense that it can include non-tunnelling effects that may serve to enhance (e.g. via vdW complex formation) or decrease the rate coefficient (e.g. entropic effects that limit ring formation or persistence, or possible steric repulsion effects). As observed in Fig. 2, some rate coefficients for highly fluorinated species appear to be overestimated by the original E-state treatment (k_base). This would indicate that in such cases, additional effects on the reactivity must be accounted for, resulting in optimized values of κ dropping below unity. Ideally, this inhibition effect would be described separately to the vdW enhancing effect. However, this is beyond the scope of the present work, and furthermore, may not be possible with the current available data.

It is known from transition state theory that the tunnelling parameter is determined primarily by the imaginary frequency of the transition state, which defines barrier curvature and width, and therefore the tunnelling probability. Viegas and Jensen³⁰ took the approach of describing this phenomenon through intramolecular electronic interactions (analysed via interacting quantum atoms), where they calculated transition state bond distances, the imaginary frequency, and the resulting reactivity of a selection of partially fluorinated ethers. These intramolecular interactions were found to be a dominant control on the reactivity in these systems. Our approach represents an alternative approximation, which describes these parameters in a simplified way based only upon a description of the vdW pre-reactive complexes. We can't expect perfect agreement between these methods, given that Viegas and Jensen effectively decouple several components of this process compared with our less specific effective tunnelling parameter. Furthermore, there are some differences between k_base and their equivalent term. Nonetheless, a qualitative comparison between the techniques is interesting. One could argue, for example, that to some extent we capture the conformational complexity inherent in these systems by considering various possible van der Waals complexes on a given site. Furthermore, we could hypothesize that by accounting for the E-state of the oxygen bridge, we are in fact describing a feature that is expected to be proportional to the transition state bond distances/energies. Ultimately, in our current parameterization, we don't have the luxury to decouple any of these terms, and so we are left with the “effective” tunnelling parameter.

The formulas used to calculate the Arrhenius parameters (A, B, and n) for predicting site-specific temperature-dependent rate coefficients given by eqn (1), based on the original E-state method, are provided in Tables SI-2 and SI-3 (McGillen et al., 2024)¹⁸ and in the SI-3 of this work in a copy/paste format for use in Excel or other programs. Optimized values of α_j and β_j are presented in Table 2. A flow diagram describing this algorithm is provided in Fig. 5. Note that several alternative algorithms were tested, including allowing for only monocyclic and also the possibility of including tricyclic complexes, but we selected the present approach based on general good performance and comparative simplicity. Detailed calculations for the example molecules in Fig. 3 are provided in SI-4. In all cases, this algorithm represents a modular upgrade to the original E-state code, leaving the previous parameterization unchanged. It should be noted that future upgrades of the original algorithm may be needed for further improvements. Our rationale for keeping k_base unchanged is that it allows us to observe the effects of structure and substitution in etheric molecules compared with their alkane and haloalkane analogues.

Table 2 Optimized values for the ether correction factors used in eqn (2). Values were optimized against the experimental measurements of the training set (acyclic compounds)

αj	αprim	αsec	αtert
5-Membered vdW complex	2.55762	0.23939	1.40790
6-Membered vdW complex	3.43778	1.14463	0.15977

---

βj	βprim	βsec	βtert
5-Membered vdW complex	694.76597	1225.65567	4.58926
6-Membered vdW complex	402.32975	739.16523	930.22576

---

Fig. 5 Flow diagram of the method presented in this work. ^aThis calculation also applies to cyclic structures, except epoxides, which are treated as non-etheric compounds.

Cyclic structures may exhibit different steric behaviour due to their lack of free rotation, which may necessitate specific treatment depending on ring size. For this reason, they are excluded from the training set, though the predictions are included in the plots and analyses. In the parametrization, they are treated similarly to their acyclic counterparts – with the exception of epoxides, whose rate coefficients tend to be overpredicted by our method. These appear not to behave like regular ethers, which may result from their highly strained ring structure, although only a small amount of data is available for testing at the time of writing. Therefore, they are treated as non-etheric compounds in the calculations. However, we note that two epoxides for which temperature-dependent kinetics are available (1,2-butylene oxide and cyclohexene oxide) both demonstrate a subtle negative temperature dependence at lower temperatures, which may indicate a minor contribution of vdW complexes.^23,24 Estimates of other cyclic compounds are in reasonable agreement with experimental data for most of the compounds using these assumptions, although there is a slight tendency for overestimation with the current treatment.

In all cases except for epoxides, the E-state value of the ether bridge was found to contain useful information for describing vdW complex formation. According to eqn (2) and Fig. 6, lower values of S(O_j) (associated with the highly halogenated ethers) suppress vdW complex formation, resulting from smaller values of κ_i (κ_min = 0.07), and a more halocarbon-like behaviour. Conversely, higher values of S(O_j) (associated with alkyl ethers) promote vdW complex formation, increase the size of κ_i, and therefore result in more vdW-influenced behaviour with a stronger negative temperature-dependent component (κ_max = 187). Within this training set, E-state values of S(O_j) range from 1.685 to 6.060 compared with an unperturbed value of 3.5, highlighting the extent to which the ether bridge is modified by substitution in this suite of etheric molecules. From Fig. 6, κ is shown to vary depending on temperature, electronic environment and the position of the reactive site of the molecule relative to the ether bridge(s) (see eqn (2) and Table 2).

Fig. 6 Frequency distribution of the effective tunnelling parameter, κ, for all reactive sites in acyclic compounds. The values of κ range from 0.07 to 187. The reactive sites of each molecule are chosen as examples here since they represent the lowest or highest extremes of κ throughout the dataset.

Results

The performance of the E-stateXvdW approach was tested against the available experimental database, and compared with other methods (where available) at room temperature and between 190 and 1500 K. While E-stateXvdW (at all temperatures) and Kwok and Atkinson (1995) (at room temperature) can be conveniently calculated through publicly available programmes using SMILES input, this is not currently the case for Jenkin et al. (2018)¹⁴ and Carter (2021).¹⁶ The Jenkin et al. (2018) estimates were obtained manually using the SAR for k₂₉₈ and k(T), and validated against k₂₉₈ when available. The Carter (2021)¹⁶ estimates were obtained using MechGen at room temperature.^31,32 Likewise, it was necessary to post-process the output of EPI Suite, in order to obtain estimates of k(T) using the information contained in Kwok and Atkinson (1995).¹² Estimates of each of these methods at room temperature and between 190 and 1500 K, where applicable, are provided in SI-1 and SI-5, summarized in file SI_k298_kT_values.xlsx in the SI.

Room-temperature estimates

Statistics for the performance of the various estimation methods are calculated using common parameters and eqn (3), and results are summarized in Table 3. These predictive methods (k_pred) were evaluated by calculating the following statistical values: (1) the percentage of predictions within a factor of 2 of experimental data (k_exp), (2) the coefficient of determination, R², (3) the standard deviation, SD, and (4) the mean unsigned error, defined as:where n is the total number of values in the dataset. The percentage of unique predictions, defined as non-identical rate coefficient estimates, is also shown and gives a measure of the degeneracy in each technique. Linear regression was applied on a logarithmic basis, as the datasets cover several orders of magnitude (∼10⁻¹⁸ to 10⁻¹⁰ cm³ molecule⁻¹ s⁻¹).

Table 3 Statistics of the performance of each method to predict room-temperature rate coefficients as defined in eqn (3). Note that these statistics are calculated based on the scope of each method (e.g. Jenkin et al.¹⁴ and Carter¹⁶ apply only to the CHO compounds)

		% Ratio ≤ 2	R^2	SD	MUE	% of unique predictions
k298	E-stateXvdW	65%	0.99	0.30	1.96	100%
	EPI Suite/Kwok and Atkinson (1995)^12	63%	0.99	0.40	2.76	89%
	Jenkin et al. (2018)^14	94%	0.99	0.17	1.37	98%
	Carter (2021)^16	87%	0.99	0.27	1.76	98%

---

We found that the E-stateXvdW approach to etheric substitution demonstrates a clear improvement in estimation accuracy compared with the original E-state method (cf. Fig. 2 and 7). A generally good performance is observed for the acyclic molecules (upon which the method is trained), and reasonable agreement for the cyclic ethers (other than epoxides, where the rate coefficient tended to be overpredicted), which were excluded from the training set (65% of predictions within a factor 2 of experimental data, MUE of 1.96 and R² = 0.99 with a standard deviation of 0.30, see Table 3 and eqn (3)). This improvement is achieved using the 12 fitting parameters shown in Table 2, but note that this doesn't include the original fitting parameters that were optimized for estimating k_base(T) in the original work of McGillen et al. (2024) for the C_nH_mX_2n+2−m set, and remained unmodified in the present work.¹⁸ These are given in Tables SI-2 and SI-3 in the SI. We consider that this treatment is justified based on the absence of all etheric molecules in the original optimization.

Fig. 7 Estimated room-temperature rate coefficients, k₂₉₈, for the reaction of OH with etheric molecules using the E-stateXvdW method compared with experimental data. These data are taken from the database of McGillen et al. (2020),⁶ unless otherwise specified (see Tables SI-1 and SI-5). Circles and triangles represent acyclic and cyclic ethers respectively. Data are colour coded according to their composition. The predictions are in good agreement with the experimental data for the majority of this set of molecules, with a clear improvement compared to the unmodified E-states method, k_base, (cf. Fig. 2). The major outlier in this dataset is an epoxide, which is expected based on the discussion in the methodology.

Comparisons with other SAR approaches from the literature are shown in Fig. 8 and Table 3. As mentioned above, the domains of applicability are different for each method, and only the approach of Kwok and Atkinson (1995)¹² and E-stateXvdW can be applied to all CHOX molecules in this dataset. In general, all methods provide reasonable estimates for the majority of compounds in the dataset where the methods are applicable. However, it is noted that Kwok and Atkinson (1995)¹² find molecules of intermediate halogenation particularly challenging, where major overestimates up to a factor of 30 are encountered. For a similar reason, even though the MechGen system^31,32 can make predictions for halogenated compounds, these estimates are judged to be out of scope of the Carter (2021) method. In this case, the predictions exhibit a higher degree of degeneracy (see percentage of unique predictions of 90% in Table 3). This suggests that this method is unable to describe the effect of substitution on the ether bridge itself, which is integral to the E-stateXvdW approach (see eqn (2)). The estimates of Jenkin et al. (2018)¹⁴ and Carter (2021)¹⁶ are found to be generally excellent (Table 3), although it must be stressed that these techniques only apply to the CHO compounds of this dataset, and therefore a limited subset (∼46%) of the more reactive ethers considered here.

Fig. 8 Estimated room-temperature rate coefficients, k₂₉₈, for the reaction of OH with etheric molecules, as predicted by the E-stateXvdW method, the Kwok and Atkinson (1995)¹² approach (implemented in EPI Suite¹³), the Carter (2021)¹⁶ method, and the Jenkin et al. (2018)¹⁴ method, compared against experimental data. These data are taken from the database of McGillen et al. (2020),⁶ unless otherwise specified. Circles and triangles represent acyclic and cyclic ethers respectively. Improvement in accuracy and applicability is observed using the E-state method corrected for ethers, compared to other available SARs.

Temperature-dependent estimates

Temperature-dependent rate coefficients were estimated with E-stateXvdW and EPI Suite, which represents an update and extension of the Kwok and Atkinson (1995) method. The output of EPI Suite was applied over a wide temperature range (∼190–1500 K) for all CHOX compounds in the dataset (Fig. 9) using the temperature-dependent treatment found in Kwok and Atkinson (1995). The behaviour of Jenkin et al. (2018) outside its stated temperature range is examined in Fig. SI-6.

Fig. 9 Temperature-dependent predictions of k(T) based on the E-stateXvdW method (panel a, b and c), the Kwok and Atkinson (1995)¹² approach as implemented and updated in EPI Suite¹³ (panel c, d and e), and the Jenkin et al. (2018)¹⁴ method (panel b, e and f), plotted against all available temperature-dependent experimental data (panel a, d and f) or against each other. The extrapolation of Jenkin's method over the temperature range not covered by the original model is shown in additional plots provided in SI-6.

As with the room-temperature estimates, each of the methods provides reasonable agreement across its specified temperature range, with the E-stateXvdW and EPI Suite tools showing broadly similar estimation accuracy compared with the temperature-dependent experimental kinetic dataset (see Fig. 9a, respective MUE of 1.88 for E-stateXvdW and 1.96 for EPI Suite, see eqn (3)). The comparisons are estimated based on a set of 23 temperature bins evenly spaced in inverse temperature-space between 195 and 1443 K, where data and estimates were available.

Similar to room temperature, the method by Jenkin et al. (2018)¹⁴ provides accurate estimates, but for a limited number of CHO compounds (57 compounds), and within a narrower temperature range than other techniques (∼250–400 K) (see Fig. 9b, e and f). When extrapolated beyond its working temperature range, this method becomes increasingly inaccurate towards higher and lower temperatures (see Fig. SI-6).

Fig. 10 shows a comparison of the ratio between estimates and experimental values at room temperature (panel a) and over the full temperature range studied (panel b). The E-stateXvdW method demonstrates better performance, achieving complete dataset coverage (100%) with reasonably low errors (r values < 11), showing close agreement with experimental data. In contrast, the method of Jenkin et al. (2018) successfully estimates all covered values with r < 5, suggesting good accuracy; however, it only covers a maximum of 46% of the dataset (i.e. only the CHO molecules are included) (see Fig. 10). The E-stateXvdW and Kwok & Atkinson (1995) methods show similar levels of accuracy, though E-state appears to perform slightly better, yielding fewer high-ratio outliers while still achieving full dataset coverage. By contrast r values for Kwok and Atkinson (1995) reach up to a factor of 30.

Fig. 10 Comparison between cumulative percentage of the ratio between predicted and measured rate coefficients, defined as r = (k_pred/k_exp)^a, with a = 1 if k_pred > k_exp and −1 if k_pred < k_exp, for the 4 estimation methods (E-stateXvdW, EPI Suite,¹³ Jenkin et al. (2018)¹⁴ and Carter (2021)¹⁶) (a) at room temperature, and (b) for the temperature range ∼190–1500 K. Note that the predictions of Jenkin et al. (2018)¹⁴ are strictly within the applicability range of ∼240–400 K – extrapolated value were not included. Note that Jenkin makes predictions for only 46% of the dataset, which is why the cumulated percentage is no greater than that amount.

In SI-7, each method is compared to experimental data and to the other approaches, together with the ratio between these values. These plots allow a more detailed comparison between each estimate.

Discussion

The current version of E-stateXvdW has been shown to perform well compared with existing techniques for estimating the temperature-dependent OH rate coefficients of etheric molecules in terms of estimation accuracy, domain of applicability (composition and temperature) and ease of calculation. This new parameterization is based around a concept that separates the contribution of direct hydrogen abstraction that has been acquired through training on the alkane and haloalkane dataset versus complex-mediated abstraction, and allows for the estimation of an effective tunnelling coefficient, κ, as well as the site-specific reactivity of a given reactive site eqn (4).

The percentage contribution of site-specificity of the j^th unique reactive site in the molecule, %k_j(T), is quantified as follows:

where i is the index of the i^th respective hydrogen abstraction site, N is the number of identical reactive sites, n is the total number of carbon atoms bearing an H atom in the molecule. The rate coefficient on the unique reactive site j, k_j, is the sum of the rate coefficients on all N identical reactive sites i, k_i. For example, for diisopropoxymethane (example 2 in Fig. 3), the site-specificity of site j = 1 is the ratio of the sum of site-specific rate coefficients for the 4 identical methyl groups on atoms i = 0, 2, 7 and 8 (see SI-4 for more details) to the total rate coefficient, k_total. Fig. 11 and 12 show how site-specificity of k(T) changes with temperature between ∼190 and 1500 K.

Fig. 11 Temperature dependences of the percentage contribution of site-specific rate coefficients, k_j, and effective tunnelling coefficients, κ_i, for a simple system, methyl ethyl ether. Site-specific reactivity is provided as a percentage, see eqn (4), on the left-hand axis, the tunnelling coefficient, κ_i(T) is shown on the right-hand axis. Values at 298 K are shown adjacent to the red dashed line, the percentage contribution of k_j(T) and κ_i(T) values are coloured according to reactive site. The structures of the vdW complexes are labelled as in Fig. 3. The overall contribution of k₂ is dominant at 190 K, but k_j becomes more evenly distributed at 1500 K, becoming more proportional with the number hydrogens contained within a given reactive site.

Fig. 12 Temperature dependences of the percentage contribution of site-specific rate coefficients, k_j, and effective tunnelling coefficients, κ_i, for a simple system, diisopropoxymethane. Site-specific reactivity is provided as a percentage (see eqn (4)) on the left-hand axis, the tunnelling coefficient, κ_i(T) is shown on the right-hand axis. Values at 298 K are shown adjacent to the red dashed line, the percentage contribution of k_j(T) and κ_i(T) values are coloured according to reactive site. The structures of the vdW complexes are labelled as in Fig. 3. At 190 K, the dominant contributions are those of k₁ (∼40%) and k₃ (∼50%), but k₁ becomes increasingly important as the temperature increases (∼80% at 1500 K). Due to priority rules in the parametrization, the rate coefficients on sites i are only accounting for specific complexes (grey ones are not accounted for).

For the simple asymmetrical ether, methyl ethyl ether (see Fig. 11), the secondary alpha site, k₂ decreases from ∼80% to ∼10% across this temperature range. This reflects the comparatively low activation energy for this secondary site, and even though k₁ and k₃ are significantly enhanced by the strongly negatively temperature dependent κ_i values, this is insufficient to compensate for the reduction in reactivity at lower temperatures imparted by the higher activation energies of these primary sites. Furthermore, even though the bond-dissociation energies (BDEs) of the C–H bonds are predicted to be smaller in the k₁ methyl compared with the k₃ methyl (95.2 vs. 101.3 kcal mol⁻¹),³³ our algorithm estimates the reactivity of both these sites is almost equal over the whole temperature range. This contrasts with approaches in the combustion community,^7,8 where BDE is taken to be the primary determinant in product branching ratios, although it is noted that such approaches may become more applicable in higher temperature regimes.

The diether, diisopropoxymethane, presents a different picture (see Fig. 12). Firstly, the predicted vdW complexes are of a polycyclic nature, compared with the monocyclic complexes of methyl ethyl ether. Secondly, the primary groups (k₁) are the dominant reactive site at all temperatures, which is partly the result of the comparatively large number of hydrogen atoms in these primary groups, which corresponds to 75% of the C–H bonds within the molecule. As with methyl ethyl ether, the secondary site, k₃, becomes more important at lower temperatures, aided by a combination of a lower activation energy and higher values for κ₃. In both cases, our algorithm predicts that κ tends towards relatively low values (∼1) as the temperature approaches 1500 K, in-line with the expectation that at these higher temperatures, the kinetic energy will easily exceed the binding energy of a hydrogen-bonded van der Waals complex, rendering it short-lived and ineffective at facilitating a reaction. Furthermore, both cases demonstrate some of the variation in the reactive sites that can be expected for etheric and similar molecules in their interactions with OH.

Concerning the range in composition and temperature, there are various strategies that SARs may utilize with respect to this parameter space. In this instance, we take a more inclusive approach, whereby all ethereal compounds of CHOX composition are estimated over a temperature range of ∼190 to 1500 K, which is similar in its scope to the approach taken by Atkinson and co-workers.¹² By contrast, the more recent work of Jenkin et al. (2018),¹⁴ in its approach to ethers, is restricted to the CHO compounds over a narrow temperature range of ∼240 to 400 K. It is logical that Jenkin et al. (2018)¹⁴ and Carter (2021)¹⁶ focus upon CHO compounds, which are dominant in atmospheric organic reactive fluxes; likewise, their emphasis on lower temperatures, which better reflect the conditions of Earth's atmosphere. Nevertheless, we contend that inclusion of the halogenated molecules in SAR development has several potential advantages. Firstly, the capacity of halogen substitutions to modify the properties of neighbouring atoms and bonds presents a diverse information-rich training set for testing, challenging and refining an SAR technique. Secondly, based on our E-state calculations, it appears that neighbouring groups containing halogenated substitutions have a potent impact upon the electronic character of ether bridges, thus affecting their ability to participate in vdW complex formation, and hence affect reactivity. Thirdly, despite their smaller contributions to the reactive flux, many of the CHOX compounds are of great industrial and commercial importance, and a detailed understanding of their photochemical degradation is necessary for understanding their environmental footprint.^34,35 Similarly, the application of an SAR to an extended temperature range provides a much more stringent test of its underlying parameterizations, and a generally good performance in this regard provides additional confidence when it is applied outside the immediate parameter space of its training set.

We find that the current implementation of E-stateXvdW is generally robust for a broad variety of chemical structures and temperatures, which is achieved using the original parameterization of McGillen et al. (2024),¹⁸ combined with a new parameterization of the hydrogen-bonded vdW structures that are thought to mediate the OH–ether reactions. These parameterizations each have a physical basis, and utilize the a priori information provided by the E-state. This allows us to interpret the role of vdW complexes in OH chemistry, indicating that a combination of lower temperatures and electron-donating substitutions on the ether bridge play important roles in enhancing the stability of vdW complexes. Furthermore, this effect appears to be minor in the case of tertiary hydrogen sites, but becomes significant for primary and secondary sites (see Table 2 and Fig. 12), which suggests that complex formation on tertiary sites, and the associated shutting off of free rotors carries a prohibitive entropic cost. It is possible that combining this current approach with more detailed calculations,³⁰ would enable us to decouple several of the factors contributing to ether hydrogen abstraction and therefore provide a more detailed physical description in future work.

By contrast, the group-additivity methods rely on purely empirical fitting factors, F(X) corresponding to various types of substitutions. This necessitates the grouping of substitutions, which leads to a commonly encountered drawback of the group-additivity methods, which have a tendency towards degeneracy in their estimates,^18,36 whereby a range of substitution patterns yield an identical estimated rate coefficient. Therefore, even when an SAR can be applied to a wide range of compositions such as Kwok and Atkinson (1995),¹² the extent to which these estimates map onto unique structural configurations remains unclear. Nevertheless, it is notable that our method and that of Kwok and Atkinson (1995)¹² provide broadly similar results for k(T) despite using very different approaches. We assume that at least some of this similarity originates from the training set for the base rate coefficients of Atkinson's method deriving from the simpler, non-oxygenated molecules, which, paired with a functional form for the temperature-dependent F(X) factor, E(X) that tends towards unity at higher temperatures, affords a method of decoupling the effect of complex-mediated vs. direct hydrogen abstraction. The main differences that are observed between our method and theirs relate to the treatment of partially fluorinated ethers, which can be significantly overestimated in the latter, and which appears to relate to the modification of the ether bridge itself by halogenated substituents. It appears possible that some of the general limitations of Atkinson's group-additivity method that have been noted by Bethel et al. (2001)³⁷ with respect to multifunctional alcohols could originate from this lack of description of the oxygen atoms in this method. There are various similarities between ethers and alcohols, and the increasing importance of vdW complexes at lower temperatures mirrors what has been observed from purely empirical work on the butanol isomers,^25,38 where the contribution of beta-sites is found to overtake the alpha-sites at low temperatures for the n and i isomers, and where a transition from positive to negative temperature dependence is noted for the t isomer as a consequence of vdW formation.²⁵ This indicates that the E-stateXvdW approach may be of utility for this functional group also.

Conclusions

We present a novel, procedurally generated approach to estimating the site-specific temperature-dependent rate coefficients of all saturated etheric molecules of CHOX composition (s.s. C_nO_≤n−1H_mX_2n+2−m and C_nO_≤n−1H_mX_2n−m) over a temperature range of ∼190–1500 K. This method utilizes the a priori information contained within the electrotopological state²⁹ – which can describe the reactivity of direct hydrogen abstraction reactions,¹⁸ combined with a newly conceived parameterization of the effect of van der Waals pre-reactive complexes upon the effective tunnelling coefficient, κ. Our approach is distinct from other SAR methods in this decoupling of direct and vdW complex-mediated hydrogen abstraction mechanisms in a physically interpretable way. We expect that this type of calculation can avoid some of the long-standing challenges associated with the non-additive effects on the reactivity of molecules with multiple functionalities,³⁹ and its physico-chemical basis is expected to provide a more reliable estimate of the reactivity of molecules with unusual or complicated combinations of substituents.

In comparison to previously published methods,^12,14,16 we find that E-stateXvdW provides a more accurate and consistent representation of the complex and highly variable temperature-dependences exhibited by the ether family in its reactions with OH, which points to a general robustness and validity in our parameterization of these systems. The physical basis for this SAR contrasts with existing approaches, which tend to use fitting parameters to describe the substitution effect, and which are optimized purely to reproduce experimental rate coefficients, without a fundamental chemical description of the substitution.

Unlike the other SAR methods with which we compare, the E-stateXvdW method provides a fully automated open-source programme for temperature-dependent rate coefficient estimation, based on the commonly used SMILES line notation (see Python file enclosed, SI-8). Consequently, it is ready at the point of publication to be applied to complex atmospheric and combustion models as necessary.

In summary, we find that this new approach offers several advantages over the traditional group-additivity techniques. Nevertheless, it should be noted that there remains scope for further improvement. For example, the original parameterizations of McGillen et al. (2024)¹⁸ form the basis for the direct abstraction contribution in E-stateXvdW, but these can be refined in the future as more data and/or improved fitting techniques become available, which, in turn, would allow us to refine the parameterization of the vdW component. Furthermore, now that the principle has been demonstrated in its application to OH + alkanes, haloalkanes, ethers and haloethers – a challenging suite of reactions whose rate coefficients span 6 orders of magnitude – it appears likely that E-stateXvdW can be extended effectively to other families of compounds, which we anticipate can improve the general accuracy of kinetic estimates in more explicit simulations of both atmospheric⁴⁰ and combustion environments.^41,42

Author contributions

Lisa Michelat: methodology, software, data curation, writing, reviewing and editing, visualization. Max R. McGillen: conceptualization, methodology, data curation, writing, reviewing and editing, supervision. John J. Orlando, William P. L. Carter: data curation, writing, reviewing and editing, funding.

Conflicts of interest

There are no conflicts to declare.

Data availability

The Python code for the E-stateXvdW programme can be found in <code.zip> together with necessary ReadMe and input files.

The data supporting this article have been included as part of the supplementary information. Supplementary information: Tables SI-1 and SI-5 (also as an Excel file <SI_k298_kT_values.xlsx>), which compose all of the experimental and predicted data associated with the parameterization. See DOI: https://doi.org/10.1039/d5ea00165j.

Acknowledgements

LM, MRM, JJO and WPLC thank the Coordinating Research Council (CRC) for their partial support of the project through contract A-108.

References

1. J. Clark, T. Farmer, A. Hunt and J. Sherwood, Int. J. Mol. Sci., 2015, 16, 17101–17159.
2. J. E. Rorrer, A. T. Bell and F. D. Toste, ChemSusChem, 2019, 12, 2835–2858.
3. T. A. Semelsberger, R. L. Borup and H. L. Greene, J. Power Sources, 2006, 156, 497–511.
4. J. D'Souza Metcalf, R. K. Winkless, C. Mapelli, C. R. McElroy, C. Roman, C. Arsene, R. I. Olariu, I. G. Bejan and T. J. Dillon, Atmos. Chem. Phys., 2025, 25, 9169–9181.
5. W. Tsai, J. Hazard. Mater., 2005, 119, 69–78.
6. M. R. McGillen, W. P. L. Carter, A. Mellouki, J. J. Orlando, B. Picquet-Varrault and T. J. Wallington, Earth Syst. Sci. Data, 2020, 12, 1203–1216.
7. Z. Serinyel, M. Lailliau, G. Dayma and P. Dagaut, Fuel, 2020, 263, 116554.
8. M. Belmekki, B. R. Giri, D. Liu and A. Farooq, Energy Fuels, 2021, 35, 16075–16085.
9. J.-T. Chen, A. A. E.-S. Mohamed, P. Wang, Y. Zhai, S. S. Nagaraja, C.-W. Zhou and H. J. Curran, Proc. Combust. Inst., 2024, 40, 105685.
10. S. Wang and L. Wang, Phys. Chem. Chem. Phys., 2016, 18, 7707–7714.
11. R. Atkinson, Chem. Rev., 1986, 86, 69–201.
12. E. Kwok and R. Atkinson, Atmos. Environ., 1995, 29, 1685–1695.
13. EPA, EPI Suite™-Estimation Program Interface, https://www.epa.gov/tsca-screening-tools/epi-suitetm-estimation-program-interface, accessed 17 January 2020.
14. M. E. Jenkin, R. Valorso, B. Aumont, A. R. Rickard and T. J. Wallington, Atmos. Chem. Phys., 2018, 18, 9297–9328.
15. B. Aumont, S. Szopa and S. Madronich, Atmos. Chem. Phys., 2005, 5, 2497–2517.
16. W. P. L. Carter, Atmosphere, 2021, 12, 1250.
17. W. P. L. Carter, Gateway to the SAPRC-16 Mechanism Generation System, https://mechgen.cert.ucr.edu/.
18. M. R. McGillen, L. Michelat, J. J. Orlando and W. P. L. Carter, Environ. Sci.: Atmos., 2024, 4, 18–34.
19. P. Dagaut, R. Liu, T. J. Wallington and M. J. Kurylo, Int. J. Chem. Kinet., 1989, 21, 1173–1180.
20. K. Tokuhashi, K. Takizawa and S. Kondo, J. Phys. Chem. A, 2024, 128, 1134–1141.
21. K. Stemmler, W. Mengon and J. A. Kerr, Environ. Sci. Technol., 1996, 30, 3385–3391.
22. D. E. Heard, Acc. Chem. Res., 2018, 51, 2620–2627.
23. H. El Othmani, Y. Ren, Y. Bedjanian, S. El Hajjaji, C. Tovar, P. Wiesen, A. Mellouki, M. R. McGillen and V. Daële, ACS Earth Space Chem., 2021, 5, 960–968.
24. H. El Othmani, Y. Ren, A. Mellouki, V. Daële and M. R. McGillen, Chem. Phys. Lett., 2021, 783, 139056.
25. M. R. McGillen, M. Baasandorj and J. B. Burkholder, J. Phys. Chem. A, 2013, 117, 4636–4656.
26. C. Bänsch, J. Kiecherer, M. Szöri and M. Olzmann, J. Phys. Chem. A, 2013, 117, 8343–8351.
27. I. W. M. Smith and A. R. Ravishankara, J. Phys. Chem. A, 2002, 106, 4798–4807.
28. R. Gugisch, A. Kerber, A. Kohnert, R. Laue, M. Meringer, C. Rücker and A. Wassermann, in Advances in Mathematical Chemistry and Applications, Bentham Science Publishers, 2015, vol. 1, pp. 113–138.
29. L. H. Hall, B. Mohney and L. B. Kier, J. Chem. Inf. Model., 1991, 31, 76–82.
30. L. P. Viegas and F. Jensen, J. Phys. Chem. A, 2020, 124, 3460–3470.
31. W. P. L. Carter, J. Jiang, Z. Wang and K. C. Barsanti, Geosci. Model Dev., 2025, 18, 8461–8483.
32. W. P. L. Carter, J. Jiang, J. J. Orlando and K. C. Barsanti, Atmos. Chem. Phys., 2025, 25, 199–242.
33. P. C. St. John, Y. Guan, Y. Kim, S. Kim and R. S. Paton, Nat. Commun., 2020, 11, 2328.
34. M. R. McGillen, Z. T. P. Fried, M. A. H. Khan, K. T. Kuwata, C. M. Martin, S. O'Doherty, F. Pecere, D. E. Shallcross, K. M. Stanley and K. Zhang, Proc. Natl. Acad. Sci. U. S. A., 2023, 120, e2312714120.
35. O. J. Nielsen, M. P. Sulbaek Andersen and T. J. Wallington, Atmos. Environ., 2025, 343, 120953.
36. L. Michelat, A. Mellouki, A. R. Ravishankara, H. El Othmani, V. C. Papadimitriou, V. Daële and M. R. McGillen, ACS Earth Space Chem., 2022, 6, 3101–3114.
37. H. L. Bethel, R. Atkinson and J. Arey, Int. J. Chem. Kinet., 2001, 33, 310–316.
38. M. R. McGillen, G. S. Tyndall, J. J. Orlando, A. S. Pimentel, D. J. Medeiros and J. B. Burkholder, J. Phys. Chem. A, 2016, 120, 9968–9981.
39. B. Ervens, A. Rickard, B. Aumont, W. P. L. Carter, M. McGillen, A. Mellouki, J. Orlando, B. Picquet-Varrault, P. Seakins, W. R. Stockwell, L. Vereecken and T. J. Wallington, Atmos. Chem. Phys., 2024, 24, 13317–13339.
40. Z. Peng, J. Lee-Taylor, H. Stark, J. J. Orlando, B. Aumont and J. L. Jimenez, Atmos. Chem. Phys., 2021, 21, 14649–14669.
41. B. Rotavera and C. A. Taatjes, Prog. Energy Combust. Sci., 2021, 86, 100925.
42. F. Battin-Leclerc, E. Blurock, R. Bounaceur, R. Fournet, P.-A. Glaude, O. Herbinet, B. Sirjean and V. Warth, Chem. Soc. Rev., 2011, 40, 4762–4782.

---