Dynamics of volatile organic compounds in a western Mediterranean oak forest

Abstract Volatile organic compounds (VOCs) are emitted from many sources and have important implications for plant fitness, ecological interactions, and atmospheric processes, including photochemistry and ozone formation. Forest ecosystems are strong sources of biogenic VOCs. We aimed to characterize forest below-canopy VOC mixing ratios, monitored by Proton Transfer Reaction Mass Spectrometry (PTR-MS), at Montseny Natural Park, a Mediterranean forest 50 km from the Barcelona urban area. Measurements were taken every 2 min during six months around the maximum emission period of summer. All VOCs had diel cycles with higher mixing ratios during the day, but different patterns over time. Monitored VOCs were grouped as biogenic, oxygenated, or aromatic compounds. Additionally, a positive matrix factorization analysis identified four emission profiles that were attributed to photochemical VOC production, biogenic emissions, mixed VOC emission sources, and traffic emissions. Even though the biogenic source was the strongest source profile at the site, we found a strong influence of anthropogenic air masses infiltrating the forest canopy and altering the biogenic air masses at the site.


H I G H L I G H T S G R A P H I C A L A B S T R A C T
• We analysed the below canopy volatile organic compounds of a Mediterranean forest of the Iberian Peninsula for six months. • Biogenic VOC emissions dominated at the forest, but there was strong influence from anthropogenic sources from elsewhere. • We identified photochemical VOC production, biogenic, mixed VOC, and traffic emissions sources. • We show how the atmospheres of forested ecosystems could be substantially affected by anthropogenic VOC sources.

Introduction
Biogenic sources emit larger quantities of volatile organic compounds (VOCs) than anthropogenic sources worldwide (Kelly et al., 2018). Atmospheric interactions amongst these two groups, however, are key to understanding atmospheric chemistry at a regional level. Biogenic VOCs (BVOCs) emitted by plants have diverse roles at multiple scales, from cellular protection and defence at the foliar level, to chemical signalling at the regional level, to influencing rainfall at the ecosystem scale (Laothawornkitkul et al., 2009;Yáñez-Serrano et al., 2020).
Biogenic and anthropogenic VOCs profoundly affect biosphereatmosphere interactions by altering atmospheric reactivity, processes of aerosol growth, and cloud formation and thus the radiative balance (Kulmala et al., 2013;Pöschl et al., 2010). VOCs are reactive compounds that can be rapidly oxidised by available oxidants, hydroxyl radicals during day, ozone during day and night and to NO 3 radicals at night. VOCs play key roles in atmospheric chemistry, one of which is the formation of tropospheric ozone (under sufficient NOx abundance), which is phytotoxic and contributes to global warming (Jia et al., 2019). VOCs are also precursors of secondary organic aerosols. The gas-phase oxidation of VOCs can form new particles by nucleation (Kirkby et al., 2016), or the oxidised VOCs can condense onto pre-existing particles, undergo heterogeneous reactions on particle surfaces, be processed in clouds, or experience further atmospheric degradation and deposition (Hallquist et al., 2009). Aerosols have a large impact on the solar-radiation balance locally, regionally, and globally, either directly by scattering Kulmala et al., 2013) and absorbing solar radiation (Boucher et al., 2013) or indirectly by influencing cloud and rain formation (Li et al., 2017). Aerosols can also increase global primary production by up to 25% for some forest ecosystems via diffuse radiative fertilisation (Cirino et al., 2014;Koren et al., 2012;Rap et al., 2018). The presence of VOCs therefore influences the associated albedo by modifying the number and size of cloud condensation nuclei (CCN) and cloud droplets, affecting how the clouds reflect and absorb light (Boucher et al., 2013) and thus the radiation balance (Andreae and Crutzen, 1997;Artaxo et al., 2009;Ehn et al., 2014;Sena et al., 2013) and hydrological cycle (Sheil, 2018).
Mediterranean forest ecosystems have high BVOC emissions with an important ecological role, and together with high solar radiation and the influence of anthropogenic sources enhance photochemistry favoring the production of ozone and aerosols. In fact, ozone production in the western Mediterranean Basin increases with altitude (Díaz- de-Quijano et al., 2009;Ribas and Peñuelas, 2006), has already damaged plants (Díaz-de-Quijano et al., 2009;Gimeno et al., 1995;Sanz et al., 2000), and often exceeds the standard threshold for ozone (Filella and Peñuelas, 2006a). Characterizing and understanding the volatilome of these Mediterranean forests is important, not only because of the strong implications to atmospheric chemistry, but also due to the strong ecological and biological implications, including positive impacts on human health . This characterization is thus particularly important below canopies, because sub-canopy plants and animals, including humans, are exposed to VOCs.
The objective of this study was to characterize the composition and dynamics of VOCs inside a forest, below the canopy of a holm oak forest representative of Mediterranean montane forests (Terradas, 1999). This location was chosen specifically because it has footprints from important local biogenic emissions and distant anthropogenic emissions (Seco et al., 2011b;Seco et al., 2013).

Sampling site
VOCs were measured at Masia Mariona (41 • 43 ′ 43.29"N, 2 • 26 ′ 24.35"E; 422 m a.s.l.) in the Biosphere Reserve of Montseny Natural Park (northeastern Iberian Peninsula, Spain). This site is 43 km northwest of the city of Barcelona and 20 km from the Mediterranean coast. The site is surrounded by a holm oak (Quercus ilex) forest and is a characteristic Mediterranean rural site with a typical Mediterranean climate, including warm summers, temperate winters, and irregular rates of precipitation . This forested area ranges in altitude from 415 to 550 m. Air circulation is strongly determined by the topography of the valley Wind speeds are highest from the south when the sea breeze develops during daytime, entering the valley from the south. Westerly winds correspond to intense advection from the north and northwest channelled into the valley and abate during the night ( Fig. 1) (Seco et al., 2011b).

Atmospheric variables
Meteorological variables were obtained from different areas. Daily air temperature at 08:00, maximum and minimum temperature, and precipitation were obtained from a manual meteorological station at Masia Mariona. We also obtained meteorological data from two other stations of the Servei Meteorològic The Puig Sesolles station is 1668 m a.s.l., so we used the data for air temperature from Tagamanent (1030 m a.s.l.), because this data set had magnitudes similar to the Masia Mariona daily values. Radiation data were obtained from the Puig Sesolles station, because data were not available from the Tagamanent station. Some discrepancies between the actual and reported meteorological conditions at the site are thus expected. These approximations to the local conditions at the site are nevertheless reliable. The location of the different sites can be observed in Fig. 1.

PTR-MS sampling methodology
VOC mixing ratios were measured from 26 June to 15 November 2019 using a Proton Transfer Reaction Mass Spectrometer (PTR-MS, Ionicon Analytic GmbH, Innsbruck, Austria). The inlet was attached to a tree trunk 1.5 m above the ground in a holm oak forest near a trekking path. A 10-m isolated Teflon tube (OD 1/4 in.) attached to a pump (JUN-AIR, Benton Harbor, USA) were used to draw air to a house containing the PTR-MS system. The PTR-MS was operated at standard drift-tube conditions (2.2 mbar, 40 • C, 600 V, and 127 Td) (Lindinger and Jordan, 1998). A catalytic converter (custom made with platinum pellets heated to 380 • C) was used to convert ambient VOCs to CO 2 . The monitored compounds in the study are presented in Table 1. The frequency of the measurements was 2 min, and data was averaged in 30 min. The background signal for each compound was measured each hour throughout the measurement period. Background signals were interpolated during the measurements. Limits of detection (LODs) were calculated as 2σ of the background averages (Table 1). Humiditydependent calibrations (using bubbled zero air to dilute the standard, regulated as closely as possible to the ambient humidity) were performed at several dilution steps using a gravimetrically prepared multicomponent standard (Riemer Environmental Inc., Miami, USA) that included all monitored compounds. The PTR-MS technique is not compound-specific, so compounds with the same mass as those reported above may interfere with the acquired signal. We have nevertheless chosen a compound assignment ( Table 1) that is well established for most of the compounds (Yáñez-Serrano et al., 2021). The measured signals at m/z 75, 107, and 135 were attributed to the compounds in the calibration cylinder, i.e. methyl acetate (MA), benzaldehyde, and p-cymene, but other compounds have likely contributed to these masses during ambient measurements.

Description and application of positive matrix factorization (PMF)
Receptor models can use measurements of ambient VOC to identify and quantify the VOC sources by analysing the changes in species correlations with time and identifying the optimal solution that explains all observed constituents (Sun et al., 2020). PMF is a receptor model that uses multivariate factor analysis based on weighted least-square fits. PMF uses realistic estimates of errors to weight the data by enforcing non-negative constraints during factor computation and is widely used to identify and quantify the main sources of atmospheric pollutants (e.g. Bari et al. (2015); Belis et al. (2014); Crippa et al. (2014)). PMF has been widely applied to VOC source apportionment and profile contribution, spanning urban/rural environments and ship campaigns (Baudic et al., 2016;Bourtsoukidis et al., 2020;Brown et al., 2007;Leuchner et al., 2015;Morino et al., 2011). The main advantage of PMF is that it can provide source profiles and contributions with no prior knowledge of VOC emission profiles. We applied PMF to the 30-min data samples of the VOC measurements collected at Masia Mariona in summer and autumn for identifying and quantifying the main sources of VOCs. PMF5.0 software (https://www.epa.gov/air-research) was used in this study.
The mathematical background of PMF analysis is comprehensively described elsewhere (Paatero and Tapper, 1994). Briefly, this statistical method uses a mass-balance equation, which in the receptor model is expressed as: where X ij is the concentration of species j measured in sample i, G ik is the species contribution of source k to sample i, F kj (also reported as "source profiles") is the fraction of species j from source k, E ij is a residual associated with the concentration of species j measured in sample i, and p is the total number of sources. The goal of the model is to reproduce the x ij matrix by finding values for the G ik and F kj matrices for a given number of sources. The values of the G ik and F kj matrices are adjusted until a minimum Q, "the loss function", for a given number of sources is found (Paatero et al., 2002). PMF solves the receptor modelling equation by minimising the loss function, Q, based on the uncertainty of each observation: Fig. 1. Google Earth satellite photograph showing the locations of the meteorological stations whose data were used in this study.
where σ ij is an estimate of the uncertainty for species j in sample i, n is the number of samples, and m is the number of species. More information about PMF analysis can be found in Supplementary materials.

Bivariate polar plots
The PMF analysis was complemented by bivariate polar plots with conditional probability functions (CPFs) to gain insights into the regional and local sources in the area. CPF analysis is very useful for identifying the dominant wind direction and wind-speed interval of high mixing ratios for a given probability of doing so. CPF was applied to find the direction of source contributions from the PMF source factors. Bivariate CPF polar plots represent the number of events with a concentration larger than the 90th percentile as a function of both wind speed and direction: where m θ,j is the number of samples in wind-direction sector θ and windspeed sector j with a concentration larger than the 90th percentile, and n θ,j is the total number of samples with the same wind direction and speed interval. The resultant CPF polar plots represent the probability that high concentrations of a pollutant were from a particular wind direction and speed, providing information about the sources in the area under consideration. A full explanation of the development and use of the bivariate case of a CPF is described elsewhere (Uria-Tellaetxe and Carslaw, 2014). Polar bivariate plots with CPF analyses were all performed in R (Carslaw and Ropkins, 2012), and the wind data were derived from Puig Sesolles. The study area is surrounded by complex terrain that includes mountains and the Mediterranean Sea, so local winds play an important role in the prevailing atmospheric conditions and the regional/local sources (Matthaios et al, 2017(Matthaios et al, , 2018. We investigated the impact of such complex sea/mountain flow interactions in the source-profile contributions by dividing the CPF bivariate polar plots into daytime and night-time in summer and autumn.

Statistics
Statistical analysis used R software (R Core Team, 2018) and Sig-maPlot software (Systat, San Jose, USA). A Pearson product-moment correlation was used to group the various monitored VOCs though a heatmap with dendrogram. Correlations were considered significant at p < 0.05. Note that the averaging for wind direction was vectorial.

Results
The monitored VOCs had higher mixing ratios in July and August than September and October (Fig. 2). This pattern represents the seasonality of the site, with air temperature, solar radiation, evapotranspiration, and ozone levels decreasing from July to October (Fig. 2).
We grouped the seasonal dynamics of the monitored VOCs based on the results of a dendrogram and heatmap (Fig. 3). We used this visualization technique to show hierarchical clustering based on the VOCs correlation values. Given the clustering obtained, we thereby clearly identified biogenic, oxygenated, and aromatic groups of VOCs. The biogenic group was then separated into two subgroups. The first subgroup contained isoprene, MVK+MACR, and acetone, likely associated with the oxidation of isoprene and monoterpene, because MVK+MACR and acetone are oxidation products of isoprene and monoterpenes, respectively. The other subgroup contained monoterpenes and p-cymene (supporting the compound assignment of signal m/z 135 to p-cymene, which is an aromatic monoterpene, thus expected to be related to monoterpene dynamics. The oxygenated group was separated into three subgroups, one for methanol, acetonitrile, and acetaldehyde, one for methyl acetate and MEK, and one for only DMS. All correlations were significant (p < 0.05). The clustering of the oxygenated group does not imply that these compounds do not originate from biogenic sources. In contrast, it implies that their atmospheric dynamics are less related to compounds such as monoterpenes and isoprene.
We split all available data into three periods to explore the seasonal changes of the VOC mixing ratios: summer (28 June to 9 August 2019),

Fig. 2.
Complete time series of the monitored VOCs for the measurement period. All VOC units are parts per billion. Temperature, radiation, rain, ozone concentration, wind direction and wind speed are also reported. The different colours represent the periods chosen for averaging. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) late summer (30 August to 1 October 2019), and autumn (1 October to 3 November 2019) (see Table 2). Diel cycles for solar radiation, air temperature, ozone mixing ratios, and wind direction clearly varied seasonally, with decreasing magnitudes from summer to autumn (Fig. 4). Solar radiation averaged 800 W m − 2 in summer, decreasing to 450 W m − 2 in autumn. Air temperature differed between periods by approximately 5 • C. Summer temperatures were highest at noon, averaging 27 • C, whereas autumn temperatures averaged only 15 • C. Ozone levels had distinct diel profiles in summer. Levels at night in summer were approximately 55 ppb, with a minimum at 09:00. This level was considerably higher than the 37 ppb of ozone at night in summer measured at La Castanya in 2009 (Seco et al., 2011b). The ozone mixing ratios then peaked at 14:00, at an average of 65.4 ppb. This peak shifted to 15:30 in late summer and autumn. Wind patterns at the La Castanya research station indicated distinct diel profiles for each period. Summer diel cycles shifted from northwest to southwest during daylight, coinciding with the maximum wind speed, with a strong shift at 18:00 back from a more northernly direction. The diel cycles for wind direction shifted more strongly in late summer than in summer and autumn. These meteorological parameters, however, were notably not measured directly at the study site but in nearby areas.

Biogenic compounds: isoprene, total monoterpenes, p-cymene, MVK+MACR, and acetone
The biogenic compounds identified in this study were isoprene, total monoterpenes, p-cymene, MVK+MACR, and acetone. This grouping was based on the dendrogram, the correlation heatmap, and the similar patterns of temperature and radiation. Isoprene and total monoterpenes had a correlation coefficient of 0.78. p-Cymene, an aromatic monoterpene, had a correlation coefficient of 0.61 with isoprene and a lower coefficient with monoterpenes of 0.52. MVK+MACR was strongly correlated with isoprene (r = 0.9) and slightly less strongly with monoterpenes (r = 0.69). This biogenic group had the strongest seasonal changes, with higher mixing ratios in summer and the lowest mixing ratios in autumn, following the same pattern as air temperature and solar radiation (Fig. 4). These changes were nearly negligible for acetone and p-cymene, which had much higher mixing ratios in late summer (0.5 ppb in late summer and 0.2 ppb in summer). Mean isoprene mixing ratios reached 0.8 ppb at midday in summer, whereas the midday ratio in autumn was 4-fold lower, at 0.2 ppb. Isoprene diel cycles increased at sunrise, peaked at 12:00, and sharply decreased after sunset. Monoterpenes had a mean mixing ratio at midday of 0.25 ppb in summer, decreasing to 0.07 ppb in autumn. Monoterpene mixing ratios peaked in the morning in summer, but the peaks were less robust in late summer and autumn, occurring later in the day because sunrise was later. These ratios peaked again at 19:30. The diel cycles of MVK+MACR were very  similar to that of isoprene, with higher mixing ratios in the summer. The MVK+MACR summer diel cycle had two peaks, one at noon and one at 18:00. Interestingly, the acetone mixing ratios did not vary with season, peaking at 3.5 ppb at midday for all periods (Fig. 5).

Oxygenated compounds: methanol, acetaldehyde, methyl ethyl ketone, and methyl acetate
The grouping of the oxygenated compounds was selected based on the strong correlation between methanol and acetaldehyde (r = 0.72). The correlation coefficients for acetaldehyde with MEK and MA were 0.64 and 0.65, respectively. Interestingly, MEK and MA were very strongly correlated (r = 0.92).
Methanol had the highest mixing ratios of all VOCs measured, similar to other nearby studies in the same forest (Seco et al., 2011b). Methanol had the highest mixing ratios in summer, peaking at 8 ppb, and very similar mixing ratios in the late summer and autumn, up to 3 ppb. The summer diel cycle peaked in the morning at 10:00, decreased at midday, peaked again at 18:00, and finally decreased to night-time levels, when methanol remained constant at 2.3 ppb. The diel cycles of the methanol mixing ratios were slightly higher in late summer than autumn (3.2 and 2.7 ppb, respectively; Fig. 6) but were generally lower than in summer.
The acetaldehyde and methanol diel cycles were quite similar, except for the peak at 10:00, which was lower for acetaldehyde. Both methanol and acetaldehyde mixing ratios decreased at midday in summer, although this decrease was not visible in late summer and autumn. MEK and methyl acetate had very similar diel cycles and magnitudes. The mixing ratios were highest in summer, peaking at 0.5 ppb for both compounds. The diel cycles of both compounds decreased at midday. Differently to the other OVOC species, these two compounds had higher mixing ratios in autumn than late summer, whereas the diel cycles in late summer peaked at midday, which shifted to 14:00 in autumn (Fig. 6).
We monitored m/z 63, which we calibrated to DMS, although the acetaldehyde water cluster may have strongly contributed to this mass. The correlation between DMS (or m/z 63) and acetaldehyde had a coefficient of 0.65, suggesting that the signal for DMS was not fully associated with acetaldehyde. This was further evident in the diel cycles. DMS had higher mixing ratios during the day and distinct diel cycles in summer, when DMS peaked at 10:00 and decreased at 18:00, similar to Fig. 4. Half-hourly average diel cycles of solar radiation at the Puig Sessolles station (a); air temperature at the Tagamanent station (b); and ozone mixing ratios (c), wind direction (d), and wind speed (e) at the Vilar de la Castanya station in summer (red circles), late summer (amber asterisks), and autumn (brown triangles). Error bars represent standard errors. Note that the diel cycles of wind direction and speed are expressed as hourly averages rather than half-hourly averages for better visualization. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.) the oxygenated compounds, where diel cycles peaked at 18:00 in late summer and winter. The LOD for DMS, however, was 0.3 ppb, so only values above this level could be assessed. The mixing ratios for DMS were higher when sea breezes developed. The DMS and acetaldehyde summer diel cycles were notably quite similar, so interference may have been be higher; the diel cycles in late summer and autumn are different to those of acetaldehyde, with a peak of DMS mixing ratios at 18:00 (Fig. 6f). Lastly, acetonitrile was assigned to the group of oxygenated VOCs. Acetonitrile correlated surprisingly well with acetaldehyde (r = 0.88) and methanol (r = 0.68). Average mixing ratios (0.15 ppb) were highest in summer, and the diel cycles peaked at 10:00 and decreased at 18:00.

Aromatics: benzene, toluene, and benzaldehyde
Aromatic compounds had the highest correlations amongst themselves, even though the correlations were generally weak. The similarity was highest between toluene and benzaldehyde (r = 0.69). Benzene correlated to about the same degree with toluene and benzaldehyde (r = 0.47). Benzene also had a correlation coefficient of 0.44 with p-cymene, an aromatic monoterpenoid. Interestingly, toluene was correlated strongly with the oxygenated compounds MEK (r = 0.66) and MA (r = 0.68).
The diel cycles of benzene had similar mixing ratios amongst the three periods, averaging 0.06 ppb. The mixing ratios started increasing with solar radiation (at 05:00) and peaked at 10:00, when the ratios plateaued until 18:00, when they decreased. Mixing ratios for toluene were highest in autumn, peaking at 0.33 ppb. The diel cycles differed but had similar magnitudes in summer and late summer. The cycles decreased at midday in summer but not late summer. Benzaldehyde mixing ratios were similar in magnitude but had different diel cycles in autumn and late summer, peaking in the morning in late summer and in the afternoon in autumn. The mixing ratios for the three periods were lowest in summer (Fig. 7).

PMF and bivariate polar plots
Additionally to the heat map with dendogram, a PMF analysis was performed to assess the VOC signatures of the different sources. This type of receptor model uses ambient measurements to apportion the observed VOC mixing ratios to signature sources (i.e. factors) by assessing changes in the correlations amongst the VOC species over time and finding the optimal solution to explain the VOC variability. This PMF analysis of the Montseny air masses identified four emission sources (Fig. 8), whose diel cycles are shown in Fig. 9. Therefore, depending on the atmospheric conditions a single compound can have different emission sources and the contribution of the compound to the identified factors can be depicted in the form of percentage and mixing ratios for each of the assigned factors (Fig. 8).
Factor 1 illustrates an emission source containing high levels of acetonitrile, a tracer of burning biomass, and other oxygenated compounds, such as methanol, acetaldehyde, and acetone. In fact, the diel cycles of this factor resembled those of methanol, acetaldehyde, and acetonitrile. This factor was assigned to photochemical VOC production. The similarity between weekdays and weekends indicated the regional nature of this photochemistry, occurring on any day of the week. Factor 2 illustrates an emission source with higher percentages of isoprene, products of isoprene oxidation, and terpenes and resembled the diel cycles of isoprene and monoterpenes. This factor dominated the mixing ratios and was termed as emissions of biogenic VOCs. Factor 3 had high mixing ratios of acetone, benzaldehyde, and p-cymene and low concentrations of methanol, acetaldehyde, isoprene, MVK, and methyl acetate. The diel cycles of this factor were very similar to those of acetone and p-cymene in late summer and autumn. The high abundance of benzaldehyde, an anthropogenic tracer, and the slight decrease on weekdays suggested that this source was likely associated with a mixed profile of VOC emissions, with VOC atmospheric oxidation, aged anthropogenic VOCs, and local BVOCs. Lastly, factor 4 was attributed to sources of traffic emissions due to its high content of aromatics and oxygenated compounds such as MEK and methyl acetate and its low content of biogenic compounds such as monoterpenes and isoprene  Figs. 8 and 9), together with the differences between weekdays and weekends. The differences between weekdays and weekends were clearer in late summer and autumn, tending to have higher mixing ratios on weekdays. The use of cars increases at weekends in summer, when many people take their vacations, may explain the less clear difference in the summer diel cycle.
Bivariate polar plots with wind speed and direction were used to assess the origin of the air masses in summer and to further aid the identification of the various emission profiles. Factor 1 could then be associated with photochemical oxidation in situ, due to its high levels in oxygenated compounds and the local origin of the source profile. The biogenic nature of factor 2 was further corroborated by its extremely local origin, with no effect at night. The bivariate polar plots aided the identification of factor 3, which had two possible origins, one more local and another slightly more distant Lastly, the bivariate plot for factor 4, containing aromatics, had the highest concentrations, with the highest wind speeds from the south, corresponding to the pollution from Barcelona. Interestingly, a small local source was from the northeast at night, possibly due to traffic emissions from the AP7 highway (Fig. 10).

VOC groups in the Montseny Mediterranean forest
This study presents the longest data set for ambient VOCs in the Iberian Peninsula, a step towards understanding VOC emission sources and gas-phase atmospheric chemistry in the western Mediterranean Basin.

Biogenic VOCs
Biogenic VOCs varied the most seasonally due to changes in the physiological activity of the vegetation and consequent VOC emission (Seco et al., 2011b). Monoterpene mixing ratios in our study were consistent with monoterpene seasonality reported by Seco et al. (2011a) and with the monoterpene emissions previously reported for holm oak, with lower emissions in winter (Llusia et al., 2012). The magnitude of monoterpene emissions in our study was nevertheless 5-fold lower than in previous studies, where monoterpenes had higher mixing ratios than isoprene (Seco et al., 2011b;Seco et al., 2013). This decrease may have been due to the slightly different forest compositions, with more cropland and mixed forest areas at our study site  compared to the La Castanya research station where Seco et al. (2011a) obtained their measurements. The monoterpene mixing ratios peaked in the morning at the start of local emission when sunlight hit the forest. This peak and the peak at 19:30 may also have been due to decreased levels of OH radicals as compared to midday OH levels, decreasing the rate at which monoterpenes are being oxidised (Harrison et al., 2001).
Previous nearby measurements of isoprene suggested low levels of biogenic emissions (Seco et al., 2011b), because the forest is dominated by holm oak, which emits isoprene at a rate less than 5% of its rate of monoterpene emission . In contrast, we found that the isoprene diel cycle was identical to the radiation cycle. Emissions from Erica arborea, a shrub species and a strong isoprene emitter Kesselmeier and Staudt, 1999), at our site was a possible origin of such isoprene dynamics. Additional local sources of isoprene, including vehicle exhaust (Borbon et al., 2001) and the industrial activity at Sant Celoni and its surroundings, 7 km from the site, may nonetheless have contributed to the mixing ratios during the development of the sea breezes. Such additional sources must have been local (maximum distance of 11 km) given the lifetime of isoprene, about 1.5 h (Atkinson, 2000). Myrcene, which fragments into the isoprene signal in the PTR-MS, may have also contributed to this signal (Tani, 2013;Yáñez-Serrano et al., 2021), because Q. Ilex emits this monoterpene (Kesselmeier et al., 1996).
MVK+MACR and acetone were classified as biogenic by the correlation heatmap and the dendrogram and had diel cycles typical of oxygenated compounds (Fig. 6). We expected that these two major intermediate products from the photochemical oxidation of isoprene and acetone would be more related to photochemistry, but some other processes could also potentially account for their biogenic classification. Biogenically, isoprene oxidation is one of the dominant contributors to MVK+MACR (Yáñez-Serrano et al., 2015), and we identified a very strong correlation between isoprene and MVK+MACR, suggesting their possible source as a product of isoprene oxidation. However, MVK+MACR have also been reported to be directly emitted by vegetation (Canaval et al., 2020;Jardine et al., 2011), including Q. ilex (Holzinger et al., 2000), and also to be taken up by plants (Canaval et al., 2020). The biogenic signal of MVK+MACR may thus represent a contribution of in situ and ex situ isoprene oxidation and direct plant emission and uptake. An anthropogenic contribution is also likely to be present, especially during periods when the advection of air masses from Barcelona is strong. Acetone had a similar combination of sources (Seco et al., 2007). Even though acetone was selected as a biogenic compound, the lack of strong seasonality may represent a more regional source, such as the photochemical oxidation of both anthropogenic (Jacob, 2002) and biogenic (Sommariva et al., 2011;Wisthaler et al., 2001) VOCs. The acetone summer diel cycle was very similar to that of methanol (Figs. 5 and 6), suggesting a possible combination of biogenic and regional-transport origins, which is further supported by the emissions of acetone reported from Q. ilex (Holzinger et al., 2000). The diel cycles in late summer and autumn also suggested transport from other regions, because the diel profile had the same profile as the source of photochemical VOC production. Additionally, propanal may have contributed to this signal (Seco et al., 2013;Yáñez-Serrano et al., 2021).

Oxygenated VOCs
Possible sources from atmospheric transport and oxidation (Seco et al., 2007) are dominant, but the strong seasonality of the methanol diel cycles, together with the summer morning peak, suggested also a biogenic component. In fact, methanol emissions by Q. ilex have not only been reported but have been found to trigger isoprenoid emissions (Roger Seco et al., 2011). Photochemical oxidation may have played an important role in the dynamics of acetaldehyde (Seco et al., 2007), similar to methanol, although a biogenic source was evident from the sharp decrease in acetaldehyde mixing ratios after sunset, particularly in summer. The seasonality of MEK and methyl acetate was less closely associated with the change in magnitude than with the change in the pattern of the diel cycles. Both compounds have been attributed to biogenic and anthropogenic sources (Andreae, 2019;Canaval et al., 2020;Fasbender et al., 2018;Holzinger et al., 2000;Yáñez-Serrano et al., 2016). Lastly, acetonitrile originates mostly from the burning of biomass , anthropogenic sources such as vehicle exhaust (Holzinger et al., 2001), and biogenic sources (Bange and Williams, 2000;Nyalala et al., 2011). The pattern of the acetonitrile diel cycle in our study coincided with the mesoscale circulation of mountain and sea breezes that transport air masses from Barcelona and its industrial area to Montseny Natural Park, driven by solar radiation (Seco et al., 2013). Montseny Natural Park had no major fires during our study, and chimneys and heaters were not expected to be used at this time, so we concluded that acetonitrile did not originate from the burning of biomass. The high percentage of acetonitrile attributed to regional photochemical VOC production, and the difference between weekday and weekend levels suggested that the acetonitrile originated from vehicle exhaust and other regional sources of pollution. DMS had the lowest atmospheric mixing ratios, which could suggest a marine regional origin, because the sea breezes transport DMS inland (Charlson et al., 1987;Jardine et al., 2014). Such air masses first pass by the Barcelona urban and industrial areas, thus depleting the DMS in these air masses.

Aromatic VOCs
Recent evidence suggests that plants can be strong sources of aromatics (Fasbender et al., 2018;Misztal et al., 2015), but their diel cycles resemble background atmospheric sources, with no changes in magnitude amongst the seasons. The diel patterns and magnitudes were similar to those at a nearby site reported by Seco et al. (2011a) and much lower than those measured at the Barcelona urban area (Filella and Peñuelas, 2006b;Seco et al., 2013). The ratio of benzene to toluene in an air parcel tends to increase over time, because toluene is scavenged more quickly than benzene by its reaction with OH radicals and thus has a shorter photochemical lifetime than benzene (Seco et al., 2013). We found this pattern in the summer and late summer (ratio of 0.35) in the morning and a ratio of 0.9 in the evening, but the ratio in autumn remained constant at 0.25. The lower ratios in the morning may have been associated with photochemical processing, and the higher ratios in the afternoon and evening may have been associated with a stronger effect of anthropogenic air masses containing vehicle exhaust (because the exhaust has a ratio of benzene to toluene of 0.5), in combination with the decrease in the inputs of benzene and toluene to the site from the sea breezes.

Source profiles
We monitored VOC concentrations below the canopy where we expected a strong influence of biogenic sources, because the inlet was below the canopy of a holm oak forest and therefore closer to biogenic sources. Even though we identified a dominant source profile of biogenic VOC emissions, anthropogenic air masses infiltrating the forest canopy strongly influenced the biogenic air masses at the site. This influence was clearly demonstrated by the PMF results identifying a biogenic VOC source profile, but also a photochemical VOC oxidation, mixed sources, and residual traffic emissions source profiles.
We identified a mixed VOC emission source profile (i.e. factor 3) combining distant and local origins, as shown by the bivariate polar plots (Fig. 10). This source could be anthropogenic, given the difference between weekday and weekend emissions, the strong presence of benzaldehyde, the footprints from the Barcelona urban and industrial areas, and the footprint of the AP7 highway at night. The source may also have had a biogenic origin from the strong presence of p-cymene and acetone, which can be emitted from biogenic sources. Lastly, secondary formation from atmospheric oxidation could also have contributed to this source profile, such the formation of acetone from the atmospheric oxidation of monoterpenes (Sommariva et al., 2011). These processes suggest that this source was likely associated with a mixed emission profile, with atmospheric VOC oxidation, local and distant anthropogenic VOCs, and local BVOCs.

Air mixing
The bivariate CPF polar plots indicated that at least factors 3 (mixed VOC emissions) and 4 (traffic emissions) had distant sources, because the effects of wind direction and speed were strong at the site, bringing air masses from elsewhere to infiltrate the canopy. This infiltration can be aided by the structure of the forest canopy at the site that allows good ventilation and the mixing of air between the canopy and the atmospheric boundary layer, as reported for other oak forests (Kalogridis et al., 2014). The atmospheric circulation and wind direction at the site were mainly dominated by the topography and height of the boundary layer. Winds from the south and southwest have been enriched with Barcelona urban and industrial pollution, but winds from the south, southeast, and east have a marine origin. The measurement site, however, was between the Besòs and Tordera basins, so the topography affected the sea breezes, favouring southern and southwestern directions. This influence was evident in the anthropogenic VOC diel cycles (benzene, toluene, and benzaldehyde) following sea breeze. Additionally, the diel cycle mixing ratios of the weekdays for factor 4, assigned to traffic emissions, may be related to the traffic rush hour peaks from Barcelona city and surrounding, as the diel cycles for factor 4 had the same peaks but with approximately 2 h lag time, potentially attributable to the time for air transport from Barcelona city to the site. The wind direction after midday was from the southwest, which abated at night. Air masses from the Barcelona urban and industrial areas are therefore transported inland due to the sea breezes and topography, increasing VOC loading in the region of Montseny Natural Park in the middle of the day, as previously reported for a nearby site (Seco et al., 2013). More local emissions from the city of Sant Celoni and the AP7 highway may have also contributed to the local polluted air masses at the study site.

Photochemistry
These wind patterns could have profound implications for photochemistry, because bringing pollutants to the site alters atmospheric oxidation (Liu et al, 2016(Liu et al, , 2018. Products of VOC oxidation can be formed at the site or transported to the site from their emission sources (Barcelona urban and industrial areas). This photochemistry process was observed from the four source profiles that demonstrated a residual traffic source, with high levels of oxygenated and aromatic compounds, originating in the Barcelona urban and industrial areas. The mixing ratios of some oxygenated compounds, namely acetaldehyde, MVK, MA, and acetonitrile, but particularly MVK+MACR, decreased after midday. This decrease has been previously reported and was attributed to photochemical destruction (Atkinson, 2000;Warneke et al., 2003) and to the effects of dilution caused by an increase in the mixing depth during the day (Filella and Peñuelas, 2006b). These patterns could also be attributed to the diel pattern of traffic emission, with one peak coinciding with the morning rush hour and another corresponding to the evening rush hour, because vehicle exhaust contains oxygenated VOCs (Inomata et al., 2016).
We chose MVK+MACR and isoprene as examples of local photooxidation. We calculated the ratio of MVK+MACR to isoprene to assess the photochemical destruction of isoprene, even though other sources may have contributed to the MVK+MACR signal. The average ratio was 0.7, within the range (0.3-0.75) for other forested ecosystems (Kalogridis et al., 2014;Yáñez-Serrano et al., 2015). Isoprene is quickly oxidised during daylight, primarily by OH to MVK+MACR and formaldehyde (Liu et al., 2013). This oxidation is fastest in the summer from 09:00 to 11:00, when ratios are higher (0.75); and then decreased to 0.5 during the afternoon. This decrease could indicate that isoprene oxidation was the dominant source of MVK+MACR in the morning but that the dominant source shifted to previously oxidised air masses transported in the afternoon. Oxidation is dominated by NO 3 radicals at night when biogenic isoprene emissions abate (even though reaction rates are low and the nocturnal boundary layer tends to trap species in the shallow layer near the ground (Biesenthal et al., 1997)), so the ratio increased. The relatively high ratio of MVK+MACR to isoprene could indicate a higher yield of isoprene oxidation products in rural forested regions, such as those in the western Mediterranean Basin, where NOx levels are higher than at other remote forest sites (Kalogridis et al., 2014;Fig. 10. Polar plots of bivariate conditional probability functions (CPFs) in summer for the four factors identified in the PMF analysis. Liu et al, 2016Liu et al, , 2018Yáñez-Serrano et al., 2015). This high ratio likely affected the formation of ozone in the area, because higher VOC emissions with high NO X levels may have been transported from Barcelona, particularly in winter (Seco et al., 2013). In fact, high ozone concentrations are already damaging vegetation in the western Mediterranean Basin (Díaz-de-Quijano et al., 2009;Gimeno et al., 1995;Sanz et al., 2000).

Atmospheric implications of observed VOC patterns
Biogenic VOCs are more reactive than anthropogenic VOCs, and isoprene and monoterpenes particularly produce high SOA yields (Xavier et al., 2019). The emissions and dynamics of forest VOCs thus have important implications for the formation and accumulation of aerosols. The anthropogenic influences in these processes also disrupt aerosol yields and products. For example, the formation of new particles from the same level of precursors is more pronounced under mixed anthropogenic-biogenic conditions compared to pure anthropogenic or biogenic conditions (Debevec et al., 2018). Furthermore, extreme photochemical nucleating bursts at midday in spring under clear skies have been attributed to the increase in biogenic VOC emissions at this time of year (Cerro et al., 2020).
The atmosphere in Montseny Natural Park is changing over time. First, a biomic shift in the Montseny Mountains has already been observed and was attributed to global change (Penuelas and Boada, 2003). Second, the long-range transport of air pollutants in the western Mediterranean Basin has decreased and is being replaced with a local to regional increase in the number of sources of pollutants (Cerro et al., 2020). Climate change and anthropogenic actions could thus further alter the atmospheric mixing ratios of trace gases at Montseny Natural Park. This influence has important implications for the formation of ozone and consequent plant fitness by affecting the atmospheric oxidation of tropospherically relevant compounds such as acetone, and has implications for the formation of SOA, because changes in BVOC emission and oxidation will lead to changes in SOA dynamics.

Conclusions
We report a six month time series of VOC mixing ratios for Montseny Natural Park, a Mediterranean forest in the Iberian Peninsula. We measured VOC concentrations below the forest canopy at a height of 1.5 m, where the effect of biogenic sources was dominant, but found a strong influence of anthropogenic VOCs transported to the forest by sea breezes. A PMF analysis identified four emission sources, which we attributed to photochemical VOC production, biogenic VOC emissions, mixed VOC emissions, and traffic emissions. Biogenic VOC emissions dominated early in the morning, but polluted air masses arrived at the site and interacted with the forest air when the direction of the wind changed due to the development of sea breezes altering photochemistry. Thus atmospheres of forested ecosystems can be substantially affected by anthropogenic VOC sources.

Data availability
Data is available upon request to corresponding author.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.