Since accurate simulation of temperature, wind speed and wind direction is critical for the accurate simulation of surface ozone and other trace gases, it is instructive to evaluate model performance of meteorology. Thus, in this section, we will first assess model skill in simulating observed meteorology before evaluating the model simulation of surface pollutants. Finally, we will examine the model scenarios.
3.1. Simulation of Observed Meteorology
January 2013 was an extreme month for Sydney in terms of both high temperatures and high rainfall with conditions so severe that the 2013 summer was also dubbed the ‘angry summer’ [
45]. The Australian Bureau of Meteorology (BoM) weather station at Sydney Observatory Hill showed an average recorded temperature of 27.6 °C for January, which was 1.8 °C warmer than the historical average for this location in this month. This was a result of prolonged and widespread extreme heat periods that affected most parts of Australia during January 2013 so that higher temperatures were recorded across the whole southeast of Australia.
Figure 3 shows maps of simulated sea-level pressure (SLP) and 2-m temperature from 7 January to 18 January 2013 at 3:00 p.m. local time. Generally, the simulated synoptic pattern shows a high pressure system over eastern Australia for the whole duration of the model run. This high pressure system, combined with hot northerly winds, caused the extreme temperatures that were observed on 8 January, 12 January and 18 January.
Figure 4a shows a time-series of mean observed and simulated temperature averaged over all 17 measurement sites. The observed time-series of mean temperature is captured well by the model with a correlation coefficient of 0.9, a mean bias (MB) of 0.2 °C and a normalised mean bias (NMB) of 0.8%. In general, as a mean across the 17 sites, the simulated temperatures are very similar (to within 1–3 °C) to the measured temperatures, with the exception of the night of the 12th of January when the model is much warmer (27 °C) than observed (20 °C).
In order to analyse model skill at simulating observed wind speeds, we used data from the BoM observation network because these are better located or sited (i.e., not sheltered) than the OEH sites (for which we found consistent underestimation for wind speeds compared to the BoM data). Across these eight sites, the observed mean wind speeds range from 1.5 m·s
−1 at night to a maximum of 7–10 m·s
−1 during the day, representing moderate wind speeds during the extreme heat periods. Not only does the model capture the magnitude of observed wind speed very well with correlation coefficient of 1.0, slope of 0.8, MB of 0.2 m·s
−1 and NMB of 3.9% (see
Figure 4b), but it also tracks the measured wind direction (
Figure 4c) with correlation coefficient of 0.8, slope of 0.8, MB of −5.7° and NMB of −3.8%. As is usually expected, the extreme heat periods were associated with winds from the warmer north or northeast interior. A time-series of measured (calculated from vertical meteorological soundings of the atmosphere using the method described in Seidel et al. [
46]) and simulated PBL heights (diagnosed by both the Yonsei University (YSU) and the Mellor-Yamada-Janjic (MYJ) PBL schemes) from Sydney airport (see
Figure 4d) shows that the PBL heights are deeper during 8 January and 18 January—extreme heat days. Both PBL schemes correctly predict deep PBL height for 8 January, but only the MYJ scheme predicts the deep PBL height observed on 18 January. Extreme PBL heights exceeding 4 km are reminiscent of those that were observed during the extreme heat events that led to the 2009 Black Saturday fires in Victoria [
47]. Across the region, a map of simulated PBL heights is shown in
Figure 5 for each of the days from 7 to 18 January 2013 at 3:00 p.m. local time. Clearly, the extremely warm days have a deeper mixed layer over land due to convection from the extreme heat, and a shallower PBL height over the Pacific Ocean to the east of Sydney. On the hot days, the region is characterised by offshore flow, which means that warm air from the land is travelling over the cooler ocean surface. The near-surface air is stable, which leads to a shallow PBL. For onshore flow, cool southerly air is travelling over warmer water (water further north is warmer) creating an unstable air mass (with air being heated from below) which leads to a deeper PBL. The model tracks the measured shortwave radiation (see
Figure 4e) very well. Although there is a reduction in radiation during the non-heatwave days as expected, there is no significant variation in simulated and measured radiation from one extreme heat event to another.
A companion paper [
48] provides more in-depth analysis of model comparison with meteorology.
3.2. Simulation of Observed Surface Trace Gases
Figure 6 presents a timeseries of simulated and observed ozone and NO
x mixing ratios averaged over 17 sites across Greater Sydney. Isoprene and formaldehyde were measured at only one site and thus are not shown here but are shown later in
Section 3.3.1. Across the 17 sites, the model generally captures the observed surface ozone mixing ratio but seems to underestimate the initial phase of the extreme heat period from 8 January to 11 January. For the extreme heat on 12th January, the model overestimates the observed mean peak mixing ratios of 50 ppbv by over 10 ppbv, but there is much closer agreement for the later extreme heat period on 18 January. Overall, the model simulates ozone with a slope of 0.8, a Pearson correlation coefficient (R) of 0.8 and MB of −3 ppbv. NO
x is more impacted by local sources. Therefore, it is not surprising that, at a resolution of 3 km × 3 km, the model struggles with simulation of NO
x (slope = 0.4, R = 0.6 MB = 0.3 ppbv and NMB = 3.8%).
The spatial distribution of simulated ozone mixing ratios from 7 January to 18 January 2013 at 3:00 p.m. local time is shown in
Figure 7. Also plotted on these maps are wind vector fields showing wind speed and wind direction. Simulated ozone mixing ratios range from 10 to 120 ppbv with background mixing ratios in the range of 10–40 ppbv over large parts of the region (over land and the Pacific ocean). During the extreme heat episodes of 8 January, 12 January and 18 January, plumes of high ozone mixing ratios are advected offshore (
Figure 7b,f,l) but over land the mixing ratios are lower (30–50 ppbv). The reduced mixing ratios over land are possibly due to a number of factors: higher wind speeds during the warm days combined with higher mixing layer height over land and the lack of ozone titration due to low NO
x and VOCs over the ocean. This is shown clearly during the cooler days when we see comparatively higher ozone further inland in areas (
Figure 7j,k) which have lower PBL height (see
Figure 5j,k). We also see that the simulated ozone plumes (
Figure 7b) coincide with simulated isoprene dips (
Figure 8). Various studies have shown strong gradients in ozone mixing ratios near forest canopies due to high fluxes of biogenic emissions (e.g., [
49,
50,
51]). This negative ozone sensitivity to isoprene can occur through various pathways. For example, it can be caused by the direct reaction between ozone and isoprene or indirectly, under NO
x-limited conditions, by reactions that remove NO
x to form organic nitrates [
52]. However, more likely, the isoprene dips are caused by direct titration through the ozonolysis reaction, the products of which include formaldehyde and acetyl peroxy radicals. The latter reacts with NOx to form PAN. The isoprene dips coincide with increase in both formaldehyde (not shown here) as well as PAN (shown in
Section 3.3).
The simulated surface isoprene maps (
Figure 8) indicate that the model may be overestimating isoprene mixing ratios. However, it is worth noting that there are very few measurements of isoprene in Australia. The highest ever observed isoprene mixing ratio is about 7 ppbv and these are normally measured in city locations and not close to the sources of isoprene. At the model resolution of 3 km, the grid cell in which the measurement site is located includes the local forest so that the model will see higher isoprene mixing ratios which in reality will be titrated out by ozone by the time the isoprene gets to the measurement site. Our simulations (not shown in the paper) of isoprene in the 9 km and 3 km domains show that the latter simulates about 5 ppb less isoprene than the former. A plot of hydroxyl radical (OH) concentrations (
Figure 9) shows OH concentrations in the range of 5 × 10
5 to 2.0 × 10
8 molecules/cm
3. As a comparison, OH measurements at Cape Grim are 3.5 × 10
6 molecules/cm
3 average maximum [
53].
3.3. Model Sensitivity Studies
Although this study has not looked at the direct influence of meteorology (such as PBL depth and stagnation events) on the elevated pollutant mixing ratios during the extreme heat, it is worth noting that these too can lead to elevated pollutant mixing ratios. However, in this study, all model sensitivity studies used the same meteorology and PBL scheme, so we only address the contribution of biogenic emissions and chemistry to air pollution episodes.
Extreme temperatures and extreme biogenic emissions were simulated during the extreme heat period in Greater Sydney.
Figure 10a,b show the simulated and monthly averaged temperature and isoprene emissions, respectively. As a mean across all 17 sites, we see that differences in temperatures between average and standard conditions were as high as 15 °C, which translates into isoprene emission differences of over 40 mole/km
2/h. Maps of the difference between simulated standard and average temperature and isoprene emissions for 3:00 p.m. (local time) on 18 January are shown in
Figure 11. The largest differences of up to 20 °C are found close to the shoreline (and therefore close to the city), with differences decreasing towards high altitude in the mountainous western and northwestern part of the domain (see topographical map in
Figure 1). The reasons for this spatial pattern are two-fold. Towards the mountains, the air is cooler and hence the difference from average temperature is smaller. Towards the cities, the urban heat island effect tends to exacerbate the intensity of heatwaves so that the difference from average conditions is more pronounced [
54,
55]. Even though the relatively cooler conditions in the Blue Mountains means lower isoprene emissions relative to the case under warmer conditions, the emissions of isoprene are still large since these areas are heavily vegetated. As shown in
Figure 11b, the difference between standard and averaged isoprene emissions simulated in January 2013 ranges from −180 to 300 mol/km
2/h (shown for 18 January only).
Figure 12 shows a time-series plot of the mean ozone mixing ratios from all 17 sites simulated with standard model conditions (scenario STD_ET) as well as simulated with the various sensitivity studies: with biogenic emissions averaged from January standard emissions (scenario AVG_E), with temperatures averaged from January standard temperatures (scenario AVG_T) and with both temperature and emissions averaged (scenario AVG_ET). As expected, the largest ozone differences are seen during the warm periods (8 January, 12 January and 18 January 2013) when the differences between standard and averaged conditions (for both temperature and emissions) are at their highest. Zooming in on the extreme heat days on 12 January and 18 January (see
Figure 13a,b), we see that the runs with averaged biogenic emissions (but standard temperatures) and averaged temperatures (but standard biogenic emissions) give similar ozone responses, with simulated mean ozone mixing ratios reduced by about 5 ppbv from the standard run. As expected, the run in which both the biogenic emissions and temperatures are averaged gives the largest ozone difference but, interestingly, this difference is approximately equal to the sum of differences from the separate runs in which only temperature or only emissions are averaged (i.e., STD_ET-AVG_ET = (STD_ET-AVG_E) + (STD_ET-AVG_T)).
Table 3 and
Table 4 show the hourly ozone differences between 6:00 a.m. and 5:00 p.m. on 12 January and 18 January 2013, respectively. On average (from 6:00 a.m. to 5:00 p.m.), we find reductions in ozone mixing ratios from the standard simulation of 3.8 and 3.2 ppbv for runs from AVG_E and AVG_T, respectively, on the 12 January, and 5.8 and 5.3 ppbv for AVG_E and AVG_T runs on 18 January. For AVG_ET runs, we find a reduction of 6.5 and 10 ppbv (relative to the standard run) on 12 January and 18 January, respectively. That the effect of temperature is as large as that of biogenic emissions is perhaps surprising, but it should be noted that on these extremely hot days the difference between simulated standard and average temperature was large (mean difference of 10 and 14 °C across all 17 sites on 12 January and 18 January, respectively). It should also be noted that in these model runs we have ignored implicitly the effect of temperature change on humidity, which will affect the OH concentrations and, therefore, photochemistry.
In the results presented above, we have assumed that there is no significant change in photolysis rates from aerosol-radiation feedback during the simulation period. As noted earlier in the introduction, the absorption and scattering of radiation by aerosols can affect radiation flux, which, in turn, affects photolysis rates. We have seen in
Section 3.1 that the simulated shortwave radiation during the extreme heat periods was largely similar. In order to assess the impact of aerosol-radiation feedback, we conducted further studies in which the model was run with aerosol-radiation feedback switched on.
Table 5 shows the simulated differences in average ozone (ppbv) between 6:00 a.m. and 5:00 p.m. for various scenarios for A (STD_ET-AVG_E), B (STD_ET-AVG_T) and C (STD_ET-AVG_ET) for runs with and without aerosol radiation feedback. Also shown is the mean simulated ozone for STD_ET run (D). The differences between runs with feedback and without feedback are on average 0.2 ppbv. This confirms that there is no significant change in photolysis rates resulting from aerosol-radiation feedback during the extreme heat period.
Figure 14 shows simulated differences in ozone between the standard run and each of the three scenarios (i.e., AVG_T, AVG_E and AVG_ET) on the three extreme heat days on 8 January, 12 January, and 18 January 2013. The maps show ozone plumes being advected from land to the Pacific ocean with the largest difference in ozone for the STD_ET-AVG_ET run. Although the spatial structures for the STD_ET-AVG_T and STD_ET-AVG_E runs are very similar, there are some salient differences. For example, on 8 January, running the model with average isoprene emissions results in more ozone (i.e., STD_ET-AVG_E O
3 is negative up to −2 ppbv) over parts of the land, whereas, for STD_ET-AVG_T, we see a positive difference of 0–2 ppbv over large parts of the land. This suggests that over land areas conditions are mostly NO
x-limited such that running the model with more isoprene removes NO
x, which results in less ozone production over land. We also see isolated positive ozone difference ‘plumes’ (i.e., less simulated ozone in average scenarios than for the standard scenario) over land, presumably from NO
x sources. The advected ozone difference (positive) plumes are more pronounced (i.e., more ozone in the standard case than in the sensitivity scenario) over the ocean where there is less isoprene (due to its short lifetime).
Due to the interconversion between O
3 and NO
2, it is instructive to consider a more conserved marker of photochemical production than O
3. Following the works of e.g., Levy et al. [
56] and Chou et al. [
57], we define total oxidant (O
x) as:
where NO
z = NO
y− NO
x i.e., all the total reactive odd nitrogen less NO
x. Differences in O
x between STD_ET and AVG_T and between STD_ET and AVG_E runs are shown in
Figure 15 and
Figure 16, respectively. The magnitudes in O
x differences are similar for the two sets of experiments, although larger for STD_ET-AVG_E, again signifying the importance of biogenic emissions to their contribution to total oxidant.
The reduction in ozone at the lower temperatures can be attributed to a number of factors, including the reduced photolysis rates at lower temperatures, the altered reaction rates for various photochemical reactions (which increase or decrease depending on the sign of the pre-exponential Arrhenius factor, whether endothermic or exothermic), and the reduced water vapour mixing ratio at lower temperatures. In the work of Sillman and Samson [
13] in which they studied the effect of temperature on ozone production, it was found that the change in reaction rates for the majority of the reactions did not impact the simulated ozone mixing ratio as much as the impact of reduced decomposition of PAN at lower temperatures. PAN is well known to represent a major sink for NO
x, especially in rural environments where NO
x is limited.
Figure 17 shows the difference in simulated PAN mixing ratios between the standard run and the average temperature run (STD_ET-AVG_T) for the three extreme heat days on 8 January, 12 January and 18 January 2013. As expected, the change in PAN is negative (up to −2.8 ppbv) so that there is more PAN at the lower (average) temperature than under normal (standard) temperature. The loss in PAN at the higher standard temperature (compared to average temperature) results in more ozone. An analogous plot is also shown in
Figure 18, this time showing the difference in PAN between the standard run and the averaged emissions run (STD_ET-AVG_E). In this case, there is higher PAN in the standard run than in the averaged emissions run. In both cases, the differences in simulated PAN concentrations tracks the differences in simulated ozone perfectly on the map, thus underlying the importance of PAN in regulating ozone mixing ratios downwind of precursors [
58].
As already mentioned, in this study, we have assumed no change to water vapour mixing ratios for the average temperature run relative to the standard run so, we can discount the water vapour effect in these simulations. This assumption is partially consistent with the results from Sillman and Samson [
13], who found the impact of water was negligible (simulated 54 and 53 ppbv ozone mixing ratios for runs with and without water corrections) for rural Michigan, although they found a very large water vapour effect for Detroit (75 and 107 ppbv with and without water correction, respectively). However, other studies have shown the opposite effect (negative for rural areas but mixed for polluted environments) of temperature on ozone via its effect on humidity [
59,
60]. Due to its proximity to the presence of biogenic sources in Greater Sydney, the atmospheric composition in Greater Sydney may be more similar to rural Michigan than Detroit. Future work will account for changes in water vapour mixing ratio to determine the impact on photolysis rates.
In our study, we have looked at the effect of temperature on biogenic emissions and not on anthropogenic emissions. Studies have shown that some anthropogenic emissions are also affected by temperature, such as NO
x emissions from soils (e.g., [
61]) (e.g., from nitrogen fertilisers). However, the temperature dependency of biogenic emissions is stronger.
3.3.1. Ozone Sensitivity to Biogenic Emissions
Figure 19 shows a comparison of measured and simulated O
3, NO
2, HCHO and isoprene mixing ratios at the MUMBA site (the only site where isoprene and HCHO were measured [
24]; see
Figure 1 for location). As we saw earlier, the model under-predicts ozone in the early phase of the study period from 8 January to 14 January, despite over-predicting isoprene. This could be due to NO
x-limited conditions as the model also generally under-predicts NO
2 as shown in
Figure 19d. Compared with the measured isoprene, the model over-predicts isoprene even for the case when isoprene emissions are halved. Unfortunately, during the extreme temperature episode of 18 January, the instruments failed due to the extreme heat and thus we do not have complete measurement data during this day. Running the model with halved isoprene emissions has the effect of decreasing simulated ozone during the extreme temperature days by as much as 30 ppbv (e.g., during the 18 January extreme heat period) below the standard run. The effect of halving biogenic emissions is only significant during the extreme heat days. For the non-extreme heat days, the simulated ozone for the various emission scenarios is not much different from the standard run. The ozone timeseries shown in
Figure 19a also shows that running the model with zero biogenic emissions has the effect of completely removing ozone episodes during the extreme heat days, which shows the impact of biogenics especially during periods of extreme temperatures. Timeseries of observed and simulated HCHO (a byproduct of isoprene oxidation,
Figure 19c) show increased mixing ratios during the extreme heat periods as expected from the increased photochemistry and with the increased isoprene emissions. It is interesting to note that, although the model over-predicts isoprene during the extreme heat period on 8 January, it still under-predicts both ozone and HCHO, which can be attributed to NO
x-limited conditions (see
Figure 19d) in the model. These NO
x-limited conditions are in agreement with a study by Linfoot et al. [
20] who analysed ozone exeedances between 1998 and 2003 using an Integrated Empirical Rate model. They found that the majority of the events were in the NO
x-limited regime especially for western Sydney.
Table 6 shows the regression statistics from the measured and simulated O
3, NO
2, HCHO and isoprene during the MUMBA campaign. For O
3, going from zero emission factor (EF) = 0.0) to standard (EF = 1.0) biogenic emissions increases the slope from 0.2 to 0.6 ppbv ppbv
−1 with a reduction in negative mean bias from −6 ppbv to −5 ppbv. The correlation coefficient is improved from 0.4 to 0.6. Comparing the half (EF = 0.5) and the standard (EF = 1.0) biogenic emissions runs, there is not much change in the sensitivity of ozone to the increase in isoprene, with the slope increased by 0.1 ppbv ppbv
−1 from 0.5 ppbv ppbv
−1 and mean bias increased by 0.6 ppbv. For HCHO, the agreement with observations improves as the biogenic emissions are increased, with the slope rising from 0.1 ppbv ppbv
−1 for zero biogenic emissions run to 0.9 ppbv ppbv
−1 for the standard run, and a negative MB reduction from −1.4 ppbv to −0.6 ppbv. Isoprene is over-predicted, with slopes of 1.6 ppbv ppbv
−1 and 4.0 ppbv ppbv
−1 and mean biases of 0.4 ppbv and 1.5 ppbv for the half and standard biogenic emissions runs, respectively. The model’s over-prediction of isoprene is in agreement with Emmerson et al. [
62], who found that the MEGAN model tends to over-predict isoprene emissions over southeast Australia by up to a factor of 6. This is thought to be due to the fact that the high emission factors in MEGAN v2.1 are estimated from young Eucalypts, which may have higher emission fluxes than the population of older trees generally found in Australian ecosystems.