The Influence of Saturation on the Surface Structure of Mixed Fatty Acid-on-Water Aerosol: A Molecular Dynamics Study

The reactivity of C=C groups in aerosols towards incoming species is highly influential to their chemical evolution and thus plays an important role in determining the properties of these aerosols and their impact on a variety of atmospheric processes. Reactions between aerosol components and gas-phase radical species often occur at the gas-liquid interface of the aerosol and thus the availability of different groups at this interface is central to determining their reactivities towards these species and the rates of these reactions. Here we look at model aerosol systems, C18 fatty acids on water, and carry out molecular dynamics simulations to determine how the presence of ‘inert’, fully saturated stearic acid molecules affects the accessibility of alkene groups within oleic acid molecules to ozone radicals. A range of stearic acid:oleic acid ratios have been studied, and a methodology has been developed in order to grow the organic layer in a random and stepwise manner that as closely as possible mimics growth in the atmosphere. The surface presence of HC=CH was found to undergo a near-linear decrease as the stearic component in the slab increased, however, the coverage of other groups was found to vary in a less linear fashion and there was an increase in the overall ordering of the organic components as the stearic acid concentration increased. It was concluded that the presence of fully saturated fatty acids is unlikely to significantly increase the rate of oxidation of unsaturated species at the surface of atmosphereic aerosols, however, it cannot be ruled out that differences in the overall structure of the aerosols with varying compositions could affect the rate of ozonolysis of these within the bulk of the organic layer.

The Influence of Saturation on the Surface Structure of Mixed Fatty Acid-on-Water Aerosol: A Molecular Dynamics Study: Impact Statement Oleic acid on water has often been used as an idealised aerosol system for experimental studies. The HC=CH group within oleic acid is known to be highly reactive towards ozone, with cleavage of this bond an important process in aerosol aging. Knowledge about the accessibility of the C=C bond to incoming ozone species is thus vital to developing accurate models about the evolution of aerosols in the atmosphere. This molecular dynamics study investigates how the presence of a second, ozoneinert, organic species (stearic acid) affects the structure of oleic acid-on-water aerosol and the accessibility of C=C groups to incoming ozone molecules.

Introduction
Aerosols play an important role within a variety of atmospheric processes, however, despite the large body of work that has been carried out to understand their origins and roles within the atmosphere, their contribution to the radiative forcing of the atmosphere remains uncertain. 1- 8 Modelling the structure and reactions of aerosols remains one of the greatest challenges to atmospheric chemists, due to the complexity and wide range of their compositions. [9][10][11] The structures of aerosols vary widely as a result of the environment in which they were produced and the aging processes that they undergo. Marine aerosols, for example, are often based around a saline core and will contain minerals found within seawater, as well as organic compounds originating from marine biological processes, and their derivatives. 9,12 Aerosols that those that come into being close to industrial areas often contain large amounts of sulfurous and nitrogen-based components, 13,14 whilst those that grow over areas of forest often contain significant amounts of isoprene or monoterpenes. 5,[15][16][17] In addition to this there are a vast number of processes that can occur to alter the structure and composition of aerosols. These include the uptake of atmospheric vapours into the aerosol, the coalescence of smaller particles to form larger ones, and chemical reactions occurring on interaction of the aerosol species with atmospheric radicals, such as O 3 , OH and Cl. [18][19][20][21][22][23] Chemical aging of aerosols, in particular, can affect their hydrophilicity and therefore their ability to act as cloud condensation nuclei and thus contribute to radiative forcing. 20,[24][25][26] Common proxies for aqueous marine aerosol are systems of oleic acid on water, as oleic acid is found within atmospheric aerosols and C=C double bonds are known to be highly reactive towards atmospheric radicals. [20][21][22]27,28 Oleic acid is released into the atmosphere both from biological sources 29 and via human emissions, particularly from cooking processes. [30][31][32] The oxidative cleavage of oleic acid aerosol at the HC=CH site by ozone, a common atmospheric radical, has been particularly well documented in laboratory experiments. 25.28, 31,33-43 These studies give an insight into the mechanisms and kinetics of the oleic acid-ozone reaction in idealised aerosol mimics.
Oleic acid-based aerosol mimics, however, usually either contain only oleic acid or oleic acid and water. Whilst such studies give an important insight into the rates and mechanisms of individual reactions that may occur at the aerosol-air interface, it is important to note that in real-world situations aerosols are not present with a single organic species. The organic component of aerosols is complex and varies from particle to particle. The presence of other species at the surface of an aerosol particle may affect the accessibility of the oleic acid C=C group to incoming radicals, and thus alter the rate of its reaction towards species such as ozone. The accessibility of alkene double bonds within oleic acid and other fatty acid species could potentially be enhanced or decreased by the presence of other species. 35,44 Indeed, there has previously been significant discrepancies noted between the rates of reaction of oleic acid molecules with ozone within atmospheric aerosols compared to aerosol mimics in laboratory studies. It has been suggested that the presence of other species within aerosols that may interfere with the ability of ozone to interact with the oleic acid unsaturated groups, contributing to these differences. 27,37 As such expanding studies of fatty acids-on-water aerosols to include multiple organic species represents an important step in gaining a greater understanding of the workings of aerosols within real-world systems. 35,[42][43][44] The most widely accepted model for fatty acid-on-water aerosol particles is known as the reverse micelle model, 45 whereby the aerosol has a water droplet core, surrounded by a 'shell' of fatty acid molecules, with their carboxylic acid groups pointed inwards to allow hydrogen bonding to the water core, and primarily non-polar groups exposed to the atmosphere. The exact proportions of polar and nonpolar groups at the surface, however, depends on the packing of the molecules within the organic layer. In the case of oleic acid and another, 'inert' species, it is feasible that another species could disrupt the ordered packing of the unsaturated fatty acid molecules, creating gaps between them, allowing ozone molecules to more easily reach alkene groups that would otherwise be buried. It is also possible that 'inert' molecules could compete with oleic acid to bond to the water core, thus leaving more 'free' oleic acid that may give groups different surface preferences than for the bound oleic acid. Alternatively, it is also possible that oleic acid molecules could become buried by the other species, thus making their reactive groups less accessible to gas-phase oxidants.
Molecular dynamics (MD) simulations are a popular method for modelling the structure of aerosol particles. The majority of MD studies of aqueous aerosol mimics have focussed on a single organic species, either pure or in conjunction with water, 19, 46-54 as opposed to looking at mixtures of species and how these interact with each other and the water core. In addition to this the focus is often on the overall structure of organic-on-water particles, as opposed to specifically what is present at the air-aerosol interface. It is thought, however, that reactive gaseous species such as ozone are often more likely to react at the surface of aerosol samples, rather than penetrating deeply into the bulk. 33,35,43,55 The nature of the functional groups that are present at the interface is therefore very important in determining the rates of reactivity and of aerosol aging.
In this work we establish an MD method for simulating aqueous aerosols, using a slab of water that is extended infinitely in two dimensions, and allowing organic species to attach to these one by one, in a way that mimics aerosol growth in the real world as closely as possible. We use this methodology to probe the random growth of monolayers containing different ratios of oleic acid and stearic acid ( Fig.1) and then investigate the differing interfacial structures of the aerosols at these different compositions. Stearic acid has been chosen as a model 'inert' species, as it is found within the atmosphere 30,56 and has also been investigated before in experimental studies of mixed organic aerosols, with oleic acid. 35,43 We place particular emphasis on probing which groups are present at the surface, with a view to determining whether the presence of stearic acid in the monolayers could affect the rates of surface-based reactions of the oleic acid component with incoming species. This represents one of a small number of studies of aerosols involving multiple organic species and also one of few computational studies that try to mimic the growth of aerosols in a stepwise, random fashion, as opposed to predefining the positions and orientations of the organic species on the water core. It thus gives a more unbiased look at the growth of multispecies aerosols that would be a step closer to the processes in the atmosphere.

Methods
Simulations in this work are carried out using a procedure designed to mimic, as closely as possible, the growth of aerosols within the atmosphere, with organic molecules being added one by one and allowed a short time to equilibrate before further molecules were added. All simulations were carried out using the GROMACS software package (version 2020.2). 57,58 The fatty acids were modelled using the OPLS-AA force field, 59 with parameters being obtained from the LigParGen server. [60][61][62] Water was simulated using the extended simple point charge (SPC/E) model 63 and constrained using the SETTLE algorithm. 64 A cubic box of edge dimensions 6.5 nm was filled with water molecules and equilibrated at 298 K, before being expanded in one axis (designated the z direction) to 15 nm, creating areas of vacuum above and below the water molecules. The total number of water molecules in each simulation was 8925. Molecules of acid were added individually, with an energy minimisation and an NVT (constant Number of molecules, Volume and Temperature) equilibration step carried out on the whole system after each addition to generate an organic-on-water slab (see Fig. 2). Energy minimisation employed the steepest-descent algorithm and was For each of the compositions, molecules were inserted into a simulation box in groups of ten according to the desired ratios of the two molecules in the final mixture. For example for the 7:3 oleic:stearic sample, 7 oleic acid molecules were added, followed by 3 of stearic acid, then a further 7 of oleic acid and so on until a total of 260 molecules, with energy minimisation and equilibration after the addition of each individual molecule. 260 molecules were chosen as this equated to a monolayer coverage for both pure oleic acid and pure stearic acid-on-water simulations. Having reached 260 molecules in the organic layer the sample was then subjected to a 10 ns equilibration followed by a 10 ns production run, both under NVT conditions, using a velocity Verlet alogorithm 66,67 and a time step of 0.5 fs. Analysis was carried out on the final 6 ns of each production run, in order to ensure that only the fully equilibrated system was sampled. The simulations all employed and particle-mesh Ewald long-range electrostatics and 3D periodic boundary conditions. All simulations were run at 298.15 K. The ESI (S2 and S3) contains a detailed discussion of the equilibration process.

Results and Discussion
The stepwise addition of the fatty acid molecules allows for the slabs to grow as 'naturally' as possible, with there not being any initial biases about the positioning of the acid molecules either with relation to the water core or with respect to each other. For each of the compositions studied four slabs were built independently of each other. The random nature of the insertions of molecules meant that each of the initial growth processes was different, prior to equilibration in order to find the energetically preferred locations of different species and functional groups. Because of the random nature of where the inserted molecules were added within the unoccupied areas of the periodic box, the distributions of the two species over the upper and lower surfaces of the slabs would not have necessarily been equal, however, the composition as a total across the two surfaces was the same for a given slab of a stated oleic:stearic ratio. This allowed for the possibility that regions of higher surface concentration of stearic or oleic acid could develop by chance, as may well be the case within a real-world system, whilst still allowing for multiple slabs of a given total composition over which statistical analysis could be carried out. The use of an infinite slab to represent the aerosol core, as opposed to the alternative method used in many studies, 19,46,47,49,68 allowed the samples to mimic larger aerosols particles where curvature issues are much less significant, but which would be too computationally expensive to simulate without the use of periodic boundary conditions. For a fuller discussion of curvature effects see the ESI section S4. Fig. 2 shows side a side-on view of one of the samples, containing equal number of oleic and stearic acid molecules. The image has been produced using the Visual Molecular Dynamics (VMD) software 69 and the atoms are color-coded based on the functional group to which they belong. The periodic box has also been marked on in black. The relatively ordered structure of the slab can be seen clearly, with the central water core (cyan) surrounded by a monolayer coverage of fatty acid molecules. These are primarily directed with their COOH groups (dark blue) directed towards the water molecules, in order to maximise hydrogen bonding with these. The methyl ends (red) of the molecules are directed towards the areas of vacuum above and below the slab. This arrangement is typical for all the compositions studied.
In order to learn more about the composition of the surface that is accessible to an incoming gaseous species, it is informative to look at the top-down VMD images of the slabs. Fig. 3 shows such images for all compositions studied, with atoms color-coded by the species that they belong to. All atoms are displayed as spheres with the van der Waals radius of the element of interest. Water molecules are shown in cyan, oleic acid in red and stearic acid in black. It can be seen that for all of the samples shown there is very little water visible, supporting the use of this number of fatty acid molecules to obtain a monolayer coverage (see ESI section S5 for details on how the number of organic molecules required for a monolayer coverage was estimated). The oleic and stearic acid molecules tend to be distributed fairly evenly across the surface of the slab, instead of forming clumps of exclusively one species. This means that even at relatively low oleic acid concentrations, unsaturated molecules are able to effectively disrupt the ordering of the stearic acid molecules, which for the higher stearic acid concentrations are able to pack in a very ordered fashion due to their straight but fairly flexible, fully saturated chains. The coverage and distributions of different species at the organicvacuum interface can give us insights into the surface structure of aerosols, however, in order to ascertain how these will interact with other species present in the atmosphere, it is vital to probe which functional groups are present at the surface. While it is possible for gaseous species to dissolve into the bulk and undergo reactions over longer time scales, there is a higher likelihood for gaseous species to interact with groups that are closer to the surface as they will encounter these more readily. As such as enhanced surface presence of certain groups above their stoichiometric ratios in the sample as a whole can significantly impact on the rates of particle aging. In addition, many of the gaseous species that will come into contact with aerosols are highly reactive. This means that they are likely to react within their first few collisions and therefore primarily with surface available species. Fig. 4 shows top-down VMD images of the same simulations as for Fig. 3, however this time the atoms are color-coded based on the functional group to which they belong, regardless of the species. Methyl groups are shown in red, carboxylic acid groups in blue, CH 2 in green, HC=CH in black and water molecules in cyan. On inspection, it can be seen that for all concentrations the surface is overwhelmingly covered by methyl and CH 2 groups, with very little HC=CH and virtually no COOH or H 2 O present. This is as would be expected for a monolayer system in which the vast majority of acid molecules are hydrogen bonded to the water core using their COOH groups, with the methyl groups directed outwards into the vacuum. The sheer number of CH 2 groups and the fact that these are also present close to the methyl end of the chain means that they are also visible in significant amounts at the interface.
A more quantitative method for determining the surface presence of different groups is the use of the solvent accessible surface area (SASA) method. 70,71 This calculates the area of a sample that is accessible to species approaching it from outside. These probe species are represented as a hard sphere of a given radius, which is allowed to move around the surface of the sample. The probe records which atoms it comes into contact with and the total surface area contact that it has with these. Surface area contact is defined as instances at which the surface of the probe is in contact with the surface of a second atom without any part of the volume of the probe overlapping the volume of a third atom. SASA analysis was carried on all of the samples that were created, and across both interfaces of these. The probe was chosen to have a radius of 0.22 nm, making it similar in size to an ozone molecule. After the SASA analysis was carried out in GROMACS the resulting surface areas for individual atoms could be grouped to give total surface area coverages of particular functional groups, as shown in Fig. 5a. While, as mentioned above, the ozone will not exclusively react at the interface, SASA gives an insight into the groups which O 3 will first encounter on interaction with the aerosol surface and therefore be at a greater chance of reacting with. SASA results are known to be similar for  Please do not adjust margins Please do not adjust margins smaller probe sizes and therefore relevant to the size of other reactive radicals, such as Cl or OH. 72 The results shown in this figure are, for each composition, averages over four independently generated slabs, each with interfaces at the top and bottom. For all SASA measurements included in this paper a thresholding parameter has been employed. This dictates a minimum surface area contact that is needed between the probe and a given atom before it can be counted as being truly part of the surface, and not in a pore in the bulk of the sample. The threshold value used was 0.05 nm 2 , with the reasoning behind this choice being discussed in our previous work. 72 Analysis was carried out at times >4 ns of the final production run, in order to ensure that the samples were fully equilibrated prior to this.
From the SASA analysis presented in Fig. 5a it can also be clearly seen that there is a large preference for methyl groups to be present at the surface, with these representing over 60% of the surface area coverage, despite there only being one methyl group per C18 acid molecule. This seems like an even greater methyl coverage than would be suggested by looking at Fig 4 alone and is indicative of the methyl groups protruding from the surface, allowing for areas on the side of these groups to interact with the probe species. As the VMD images shown are purely a 2D projection from the top down this would not be visible from these pictures alone.
Additionally, in Fig. 5a COOH is barely visible on the scale shown, supporting the idea that these parts of the molecules are primarily as close as possible to the water core (see also S7 of the ESI). There is a small amount of HC=CH present at the interface. This is highest for the slab containing only oleic acid in its organic layer and decreases as these oleic acid molecules are replaced by the fully saturated stearic acid molecules. Figures 5b and 5c show the plots of the mole fraction of oleic acid within the slab against the surface area coverage for selected individual functional groups. The error bars represent one standard error of the mean across the four production runs for a given composition. COOH has not been included as its surface area coverage is so small for all compositions that any trends are not thought to be statistically significant. It can be seen from the plot in Fig. 5b that the majority of compositions follow a linear trend for HC=CH coverage, which would suggest that for the systems studied the accessibility of HC=CH groups to incoming ozone molecules is simply a function of the proportion of species containing this group within the sample, and that the presence of the stearic acid is not otherwise changing the surface structure of the slabs in such a way as to hinder or increase the availability of these groups at the organic-vacuum interface. This supports experimental work by King et al., 43 which suggests that stearic acid is largely a bystander in the reaction of monolayers of oleic/stearic acid mixtures with ozone.
Potential differences in the structures of slabs of different compositions can be probed by examining changes in the distance between carbons C1 and C18 at the opposite ends of the chains with differing oleic:stearic ratios (Fig 6). These distances have been measured using gmx distance (included within the GROMACS software package). Distances shown are averages over both interfaces of four slabs per composition and over times within the trajectory >4 ns. Error bars represent one standard deviation of the mean and the points plotted represent distances that have been ordered into 0.05 nm bins. It can be seen that for all compositions the average C1-C18 distance of stearic acid is slightly greater than for oleic acid. As discussed above, this is likely due to stearic acid being able to extend more in order to pack more efficiently, due to its lack of planar HC=CH group. For lower concentrations (0-40%) of stearic acid the distributions of the C1-C18 distances for both stearic and oleic acid seem to be fairly constant in terms of both peak distance and spread, however, on increasing to higher concentrations of stearic acid there is a notable increase in the average C1-C18 distance and decrease in the spread of bond distances for stearic acid, supporting the idea that at these concentrations stearic acid is able to pack more efficiently and there is a corresponding preference for chains to be more fully extended. Interestingly there is also a corresponding trend in oleic acid C1-C18 distances at these coverages, indicating that the increased packing of the organic layer as a whole also forces a greater extension of the oleic acid molecules towards their greatest possible lengths.
Another indicator of the differing surface structures of organic aerosols with different compositions is the change in tilt angles as the ratios of different species within the organic layer are varied (Fig.  7). The tilt angle is defined here as the angle that the C1-C18 vector of the organic molecules make with the z axis (surface normal for a completely flat surface). We measure here the tilt angle for each of the species separately and display them on the same axes for comparison. Angles have been calculated using the gmx gangle program that is included within the GROMACS software. In addition to averaging over four slabs of the same composition and over times >4 ns of the 10 ns production runs, the results have also been averaged over the positive and negative z directions, thus taking into account molecules that are attached to the upper and lower faces of the slab. Whilst this means that we cannot distinguish between a molecule that has its COOH group pointing upwards or downwards on a given face of the slab the vast majority of molecules will be oriented with their COOH groups pointing closer to the water core than their methyl groups (as discussed above). It is therefore only the magnitude of the angle relative to the z axis that is important here and not its direction (clockwise versus anticlockwise). As before error bars represent one standard error or the mean. Angles calculated have been grouped into 1° bins.
Looking at Fig. 7 there are some striking differences between the overall angular distributions for the organic molecules at different compositions. For slabs of higher oleic acid content (>50% oleic acid) the angular distributions are very similar for oleic acid and stearic acid molecules within slabs of a given composition. The angular spreads for these slabs are relatively broad and peak at angles ~20-30°. On increasing the stearic component further the distribution shows significant changes, becoming narrower and shifting to higher angles, with the peak at 1:9 oleic acid:stearic acid being around 50°. The narrowing of the distributions indicates that there is a greater ordering of the chains due to the ability of the less bent saturated molecules to pack together in a closer and more organised fashion. In general, the angular distributions of the stearic acid molecules for slabs containing greater than 50% stearic acid in their organic layer are narrower than the oleic acid angular distributions, indicating that within these slabs the regions of stearic acid have a greater degree of ordering than the rest of the slab as a whole. There is, however, a corresponding narrowing of the oleic acid distributions, suggesting that the tighter packing of the stearic acid molecules also enforces a CH 3 Fig. 7 Distributions of the angles that the C1-C18 vectors make with the surface normal (z axis) for different ratios of oleic:stearic acid. Plots are averages over the of four slabs, and over the positive and negative z directions. Only data at times >4 ns into a 10 ns production run have been included. Error bars give one standard error of the mean. Distances have been binned into 1° bins. Fig. 6 Distributions of C1-C18 distances for different oleic:stearic ratios. Plots are averages over the top and bottom of four slabs, and at times >4 ns into a 10 ns production run. Error bars give one standard error of the mean and their small size indicates the consistency of these results. Distances have been binned into 0.05 nm bins greater degree of order on those oleic acid molecules that are present within the same organic layer. The SASA analysis above shows that the HC=CH presence at the interface scales approximately linearly with the oleic acid presence in the slab as a whole, suggesting that stearic acid is largely a bystander in whether or not oxidation reactions of oleic acid take place at the aerosol-atmosphere interface. However, the increased fraction of stearic acid may cause the organic layer of the slab to be more viscous. This may prevent the penetration of O 3 further into the organic layer and slow down the rate of HC=CH oxidation within the bulk. 73,74

Conclusions
Modelling the surface structures of aerosols and their interactions with other atmospheric species is a hugely important task for atmospheric chemists and yet it is made incredibly complex by the vast numbers of possible compositions of these aerosols. Here idealised aerosols mimics consisting of a water core covered by a monolayer of different ratios of saturated and unsaturated fatty acids have been simulated. A methodology has been developed in order to mimic the natural and largely random growth of organic layers on aerosols in the atmosphere and analysis has been carried out in order to learn more about the structures and surface compositions of these aerosol mimics.
Changing the oleic:stearic acid ratio within the organic layer of the aerosol mimics significantly changed the overall ordering of molecules within these. An increase in the stearic acid component at values of >50% stearic acid led to the distributions of C1-C18 distances for both the stearic acid and the oleic acid molecules becoming narrow and peaking at higher values, showing that increased concentrations of the less bent saturated acids lead to a greater extension, and probable closer packing of the chains in the sample as a whole. There was also a greater tendency for chains of both species to be tilted away from the surface normal at higher concentrations of stearic acid, and angular spreads were much smaller, again supporting the idea that there is a greater degree of ordering amongst slabs that contain a greater proportion of the saturated fatty acids.
Fatty acid molecules are found to bond to the water core using their carboxylic acid groups, with their methyl groups extending out into the vacuum. This leads to an overwhelming preference for methyl groups at the vacuum-organic interface, as well as a significant surface coverage of CH 2 , due to the number of these present per molecule and their proximity to the methyl ends of the chains. These surface preferences were quantified using SASA analysis and it was found that there was a near-linear relationship between the number of unsaturated groups present within the sample and the number of these present at the surface. On the other hand, methyl and CH 2 coverage at the interface varied in a non-linear way with oleic:stearic ratio within the organic layer. On increasing the stearic content of the sample there was little change in the surface coverage of these groups at coverages of <50% stearic acid, however at coverages of >50% stearic acid there was a near-linear decrease in surface CH 2 presence and a corresponding increase in surface CH 3 presence. This aligns with the tighter packing of stearic acid, which leaves proportionally fewer CH 2 groups at the surface. It is therefore thought that the presence of ozone-inert species, such as stearic acid, does not affect the rate of surface reactions that ozone can undergo with HC=CH, however, the changes to the overall structure and packing of the organic layer may lead to a decrease in the rate of bulk oxidation within this layer.
Gaining more knowledge about the structure of organic aerosols in general and the organic-atmosphere interface of these in particular is of great value when trying to generate accurate models of processes involving aerosols within the atmosphere. This represents one of only a small number of molecular dynamics studies of organicon-water aerosols containing more than one species within the organic layer and to our knowledge the first to look at mixtures of saturated and unsaturated fatty acids, despite the fact that these have been used as model systems for the study of the reactions of aerosols with ozone. It therefore provides insight into how the presence of different species of similar size but differing functionality and geometry can affect the overall structure of aerosols, and represents a step towards more accurate modelling of real-world aerosols systems.

Conflicts of interest
There are no conflicts to declare.