Statistical modelling of air quality impacts from individual forest fires in New South Wales, Australia

.

measured by ozone levels, tends to occur with a high-pressure system to the east of Sydney with light north-westerly winds and a sea-breeze (Hart et al., 2006).
There are a variety of ways to improve our understanding of PM 2.5 from individual fires. Atmospheric dispersion models can predict the spread of particulates from fires based on modelled atmospheric dynamics and are routinely used in many countries 110 to guide burning operations and community warnings for HRBs. However, while evaluations of such systems are rare, existing evaluations indicate a poor to moderate agreement between predictions and observations (Yao et al., 2014;Saide et al., 2015), possibly because the local effects of HRBs are poorly captured by the models.
An alternative method is to relate air quality observations directly to real fires to calculate how far the smoke impact is likely to spread and under what conditions. Air quality measurements can be from ground-based stations or via satellite-based 115 measurements, e.g. aerosol optical thickness (Gupta et al., 2007). For ground-based measurements, studies have been done using monitors mostly stationed within ~10 km from HRBs (Pearce et al., 2012;Price and Forehead, 2021). Pearce et al. (2012) made 684 24-hour observations of PM2.5 by placing monitors around 55 forest HRBs. They found that PM 2.5 concentrations fell to near-background levels within 3 km of the fire perimeters. Price and Forehead (2021) made 5445 hourly observations of PM 2.5 with a combination of fixed and mobile monitors around 18 forest HRBs. They also found that PM 2.5 concentrations 120 had largely fallen to background levels by 3 km but this depended on weather conditions. One of the HRBs caused poor air quality at monitors more than 30 km away. These studies captured the local effects of the HRBs but did not explain why HRBs can impact air quality much further away.
Deploying air quality monitors to wildfires is difficult due to the large size of wildfires, unpredictable ignition and spread and the safety risks of working near an active wildfire. However, large permanent air quality monitoring systems can be used to 125 gather PM2.5 data for wildfires and HRBs, for example, the NSW Air Quality Monitoring Network. Here, we used historical fire and air quality data to identify the occasions when an individual fire was burning within 150 km of a monitor in the NSW Air Quality Monitoring Network from 2012 to 2021, and developed random forest models of PM 2.5 concentrations at individual monitors as a function of fire area, distance and weather conditions. Our aims were: 1) Improve understanding of the fire and weather conditions that influence smoke dispersal and PM 2.5 levels. 130 2) Develop predictive models of PM 2.5 output from individual forest fires, as a complement to physical models, to improve warnings. 3) Make inferences about potential changes in HRB protocols that could reduce PM 2.5 impacts.

Fire Data 135
Our study period was from February 2012 to September 2021 because this was when our main fire data set was available (see below). For the study period, we created a spatial dataset of forest fires that were actively burning within 150 km of air quality monitoring stations (AQS) maintained by the NSW Department of Planning and Environment (DPE) (Fig. 1). 150 km captures most of the eucalypt-dominated Blue Mountains that is subject to the majority of fire activity near Sydney. We assigned attributes of fire location, fire type (Hazard Reduction Burn (HRB) or Wildfire (WF)), date of fire activity and AQS name and 140 location. Each fire had at least one active date and most burnt on several days. As a fire could be within a 150 km buffer of multiple AQS, there was a separate row in our data for each fire and AQS combination. For our modelling, we used only cases where, for each AQS and day, only one fire was active within 150 km of the AQS. We did not analyse cases where multiple fires were burning on the same day near the same AQS as it was unclear which fire produced the smoke that reached the AQS.
We relied on two data sources to identify fire locations, type and active dates: NSW fire history GIS polygons (NPWS Fire 145 History -Wildfires and Prescribed Burns, 2022), maintained by NSW National Parks and Wildlife Service (NPWS), and VIIRS SNPP hotspots, downloaded from NASA's Fire Information for Resource Management System (Schroeder et al., 2014;Fire Information for Resource Management System (FIRMS)). VIIRS SNPP, which refers to the Visible Infrared Imaging Radiometer Suite -Suomi National Polar-orbiting Partnership, hotspots were available beginning 20 January 2012.
The fire history dataset is a spatial polygon dataset of the final burnt area of fires across NSW, which has attributes of fire 150 identity (name and number), fire type (HRB or WF) and start and end dates. We did not rely solely on the fire history to identify fire locations and dates because an initial inspection suggested some issues for our analysis. These included fires identifiable from VIIRS hotspots/images that were missing from the fire history; occasional errors in the start and end date recording; the final fire polygon being the combination of separate fires that eventually merged; and the data identifying only fire start and end date, not whether a fire was actively burning on each day between those dates (e.g. fires may have extinguished then 155 reignited on different days). Also, the data did not capture daily fire progression only the final boundary, meaning the location of fire activity on the first day (perhaps a few hectares) was not well represented by the final fire polygon (perhaps tens or hundreds of thousands of hectares), which was particularly an issue for WF.
We employed a process to map active fire dates and locations from clusters of VIIRS SNPP hotspots. We used VIIRS SNPP hotspots instead of MODIS as VIIRS are higher resolution (at nadir, 375 m vs. 1km for MODIS), thus can detect more hotspots 160 per fire than MODIS, which reduces the chance that an active fire is missed (Schroeder et al., 2014). The process to create hotspot clusters for each day for each AQS was to: 1. Extract all hotspots within 150 km of the AQS. 2. To focus on forest fires, remove hotspots that were in grassland or open woodland by removing hotspots with low foliage projective cover score (Gill, 2012;Gill et al., 2017). This measure of canopy density is equal to the 165 proportion of ground that the vertically projected area of the green foliage covers. We removed hotspots with a cover fraction of less than 0.25 so that our analysis only included dense woodlands, open forests and closed forest types (Specht and Specht, 1999). 3. Buffer each hotspot by 2.5 km and dissolve overlapping buffers into a single polygon, thus creating hotspot cluster polygons (Fig. 2). 170 4. Remove clusters that did not have at least five day or night hotspots. This was our minimum threshold for fire activity, as we wanted to exclude small fires such as burning heaps on farmland that can be detected by VIIRS. We also tested three as a minimum threshold, which produced similar but slightly less accurate models. 5. For each cluster, calculate the daily active fire area by intersecting the hotspot points with a 500m x 500 m grid (25 ha cells). The area assigned to each cluster was the number of unique intersecting cells x 25 ha. 175 6. Repeat the process for each combination of date and AQS.
Where a fire identified from the above process (a "VIIRS fire") intersected an NPWS fire history polygon between its start and end date, we assigned the fire name, number and type (HRB or WF) to the VIIRS fire. If multiple VIIRS fires intersected the same fire history polygon, we merged them into a single fire with the same attributes for analysis. If a VIIRS fire intersected multiple fire history polygons, we assigned the attributes from the fire history polygon with the largest overlapping area. NPWS 180 fire history polygons were excluded from analysis if either the start or end date was missing or a polygon had no intersecting VIIRS hotspots. If a VIIRS fire did not intersect a fire history polygon, we assigned the fire type based on the date: from October to February (inclusive) were WF and all other months were HRB. For each fire identified we added attributes of distance and direction from AQS to fire centroid (Fig. 2), i.e. the arithmetic mean of the hotspot coordinates, with a separate row for each fire and AQS combination (within 150 km). 185  overlapping buffers merged and hotspots given separate identity numbers based on which final buffer they fell within. Two separate fires are created here because of distinct fires where buffers don't overlap, i.e. greater than 2 buffer widths (5 km) apart. Fire 1 (grey) has > 50 hotspots, fire 2 (blue) has 5 hotspots. Note that fire area used in analysis was calculated via an intersection with a 500 m x 500 m grid, not buffer size (see methods).

PM2.5 Data 200
We modelled PM 2.5 (µgm -3 ) as a function of several environmental predictors. We downloaded all available PM 2.5 data (hourly averages) from the NSW DPE for the period 2012 -2021, which comprised 48 AQS. Data were available free online at https://data.airquality.nsw.gov.au/docs/index.html. We calculated mean PM 2.5 for each AQS for three six-hour periods: 1. Afternoon: 1400 to 1900 AEST inclusive. This period covered peak burning conditions in the afternoon and after sunset, although sunset and fire ignition times varied. 205 2. Night: 2100 to 0200 AEST inclusive. Covered the night period starting on the same day as the fire. 3. Morning: 0500 to 1000 AEST inclusive, next day after fire day. Captured early next morning conditions after the main periods of fire activity are likely to have ended, although some fires may have burnt through the night and smoke may still have lingered.
Note that there were some missing PM2.5 values in the data, which meant some summary afternoon/night/morning values had 210 < six records. However, > 98 % of records were summarised from >= four hourly PM 2.5 values.
We chose these times to represent different periods in the daily cycle that may have distinct smoke, weather and fire behaviour characteristics. All fires identified in the hotspot analysis were matched to AQS summary PM 2.5 for active days when the fire was within 150 km. Not all AQS had records for all years, as some were not operational until later in the study period ( Fig. 1).
Note that we modelled PM 2.5 observed at air quality stations, which would include primary and secondary PM 2.5. Secondary 215 PM 2.5 can be formed via atmospheric chemistry processes that transform emitted gases into particulates, with the processes influenced by factors including season, solar radiation, temperature and relative humidity (Cope et al., 2014;Fine et al., 2008).

Predictor variables
We sampled hourly weather variables at each AQS and each fire centroid from ERA5 weather grids, which is an atmospheric 220 reanalysis product with multiple weather variables and atmospheric levels available at 30 km spatial and hourly temporal resolution (Hersbach et al., 2018b, a) (Table 1). We calculated the mean weather values for both surface and upper atmospheric conditions (Table 1) for the afternoon, night and morning periods as described for PM2.5. We calculated additional variables describing the spatial relationship between the fire and each AQS. We used the AQS-to-fire direction and wind direction to calculate the percent of time-period where the surface wind was blowing directly to the AQS, with directly meaning ±22.5 225 degrees of the AQS-to-fire bearing. We also used the daily active fire area based on the intersection of hotspots and a 500 m by 500 m grid (area = N intersecting cells x 25 ha), as a predictor. We included a month variable (i.e. month of the active fire date) as a predictor variable to account for any seasonal variation in background PM2.5 levels . Month was represented as a cyclic variable, where the sine and cosine of the month (1-12) were both included in the modelling. We included the latitude and longitude of the AQS to account for spatial dependence, and fire type as a factor variable to account for differences not 230 captured by the weather/fire area variables. We also experimented with making separate models for each fire type (HRB model and WF model) for each time period but found that resulting accuracy statistics on the training and test sets were similar, so instead just used one model for each time period with fire type as a factor variable.

Longitude
Coordinates of air quality monitoring stations, to account for spatial dependence.

Random Forests Modelling 260
Our data consisted of three separate tables (afternoon, night and morning data tables) for three models. In each table, there were 11187 rows with unique combinations of fire, AQS and date. For each fire, there could be multiple active dates and each fire could be associated with more than one AQS (i.e. it was within 150 km of multiple AQS). Our data had 48 different AQS and 1429 different days with at least one active fire near an AQS. There were 1883 different combinations of fire and day (we refer to these combinations as "fire-days") consisting of 727 fire-days that had VIIRS hotspots and a fire history record and 265 1156 fire-days that had only VIIRS hotspots. The fire-days from solely VIIRS hotspots were on average smaller than the firedays with a matching NPWS fire history record (209 ha vs 854 ha respectively). 1182 fire-days were from HRBs (mean daily active fire area = 254 ha) and 701 were from WFs (mean daily active fire area = 802 ha). Each fire was observed at a minimum of one AQS, with a mean of six AQS and a maximum of 35 AQS associated with a single fire.
We trained a random forest model using the "ranger" package in R (Wright and Ziegler, 2017). Random forests are robust and 270 efficient machine learning algorithms that involve fitting and averaging of randomized decision trees and have been applied to a range of environmental research problems including fire and emissions (Biau and Scornet, 2016;Hu et al., 2017;Shah et al., 2022). We chose random forests due to several advantages that include high accuracy, fast computation times, easy implementation, robustness and greater interpretability (compared to "black-box" methods) via simple methods to extract variable importance and partial dependence (Rodriguez-Galiano et al., 2015;Biau and Scornet, 2016;Wright and Ziegler, 275 2017).
We split each of our datasets into training (75 %) and test (25 %) sets for analysis, stratified by fire type so that an even proportion of HRBs and WFs appeared within each of the sets. We trained the models using the training set data and used outof-bag (OOB) predictions vs observations for model accuracy checks and we predicted to the test set to calculate test set accuracy statistics. Our accuracy statistics were the correlation coefficient (r), normalised mean error (NME) and normalised 280 mean bias (NMB), as recommended by Emery et al. (2017) for assessing model performance. We ran three different models, one for each analysis period: 1) afternoon mean PM2.5, 2) night mean PM 2.5 , and 3) morning mean PM 2.5 . Predictor variables were the weather variables in Table 1 sampled at both the AQS centroid and fire centroid, distance, daily active fire area, month and AQS coordinates. As highly correlated variables can introduce bias into random forests variable importance calculations (Strobl et al., 2008), we removed variables from analysis where the Pearson correlation was above 0.8: MSLP at 285 AQS and wind speed 850 hPa at AQS were excluded, each of which was correlated with the version sampled at the fire.
We assessed the variable "permutation" importance using the ranger package. Permutation importance is derived from a process where reduction in model accuracy on OOB predictions is calculated after randomly shuffling values for each variable, calculated for all trees and variables (Wright et al., 2016). We assessed predictor variable effects using partial dependence plots calculated in the "pdp" package in R (Greenwell, 2017), and by creating prediction plots where PM2.5 was predicted with 290 all variables held at mean values except two variables of interest, which were each assigned three different levels to illustrate their effects. We also conducted a short descriptive analysis, using satellite images and hourly PM 2.5 of large outliers in the models to understand potential reasons for inaccurate predictions. This is included in Appendix A.

Limitations
There are several limitations to our methods that should be considered when interpreting the results. Our process to identify 295 active fires from VIIRS hotspots excluded hotspots that were outside the 150 km AQS buffer, even if they were part of a fire that straddled the buffer edge. There may be occasions where smoke from hotspots, and entire fires, from > 150 km reached an AQS and influenced PM 2.5, e.g. large WFs during the 2019-20 "Black Summer". The effect of such fires was not captured in our methods.
We set a minimum fire activity threshold of five hotspots (day or night). This may mean that days recorded as having only one 300 fire may have had other smaller fires in the area that may have produced smoke that affected PM 2.5 . Relying on VIIRS had the advantage of being able to better detect when a fire was active, but our process may not have captured all fires on any given day due to cloud cover impeding VIIRS hotspot detection. This may be a form of bias in our analysis as the cloudiest days were selected against. Additionally, VIIRS SNPP hotspots are acquired early afternoon and early morning, meaning that the total burnt area on a day is not measured, only the active area at the time of VIIRS acquisitions. Fire area, or the number of 305 fires, may have been underestimated if clouds were impeding hotspot detection. Our decision to analyse only days with one fire, to better understand distance and direction variables, means that there is a selection bias against the most active fire days (i.e. days with multiple fires). This may include the worst WF days, where multiple fires were more likely to ignite, particularly during 2019-2020. For days that are most suitable for HRBs, authorities are more likely to ignite multiple HRBs. Such days, which could include the worst pollution events, were not included in our analysis but were the subject of separate research 310 (Storey and Price, 2022).
Note that in our VIIRS hotspots clustering process, we used a buffer of 2.5 km to provide a broad "search" area in which to group hotspots: any hotspots within 5 km of each other (two buffer widths) or less would be grouped. This may have meant that on some occasions, separate small fires were grouped. However, we deemed it reasonable to treat these as one fire for our purposes given the similar location meant smoke would be travelling along the same general bearing towards an AQS, which 315 was important for the direct wind variable (Table 1). For example, two fires 5 km apart would have a ~3 degree difference in bearing to an AQS 100 km away (~5 degrees at 50 km). Smaller or larger buffers may have produced different results. Note that if more than one hotspot cluster intersected the same NPWS fire history polygon, we also treated these as the same fire.

Variable summaries 320
Plots of the distribution of PM 2.5 and predictor variables are shown in Figure 3. PM 2.5 was skewed toward low values (afternoon, night, morning mean = 8.1, 10.7, 10 µgm -3 ), with occasional very smoky periods (afternoon, night, morning maximum = 294.2, 394.8, 506.2 µgm -3 ). Most fires were between 75 and 150 km from AQS and only 20 % of fires had their closest AQS within 50 km. Daily active fire area derived from VIIRS hotspots was heavily skewed toward lower values (mean = 458 ha, 95 th centile = 1175 ha). The maximum fire area was 31800 ha, < 1 % of fires (all WF) were over 10000 ha and 94 % 325 were less than 1000 ha.
Afternoon conditions were generally hotter, less humid and had higher PBLH at both fire and AQS locations than nights and mornings. Between WF and HRB, WF afternoons were hotter, drier and had higher PBLH (Fig. 3). MSLP was similar between afternoon, night and morning, but skewed lower for WF compared to HRB. The wind direction variables were clustered around zero, indicating that most of the time wind at the fire and at the AQS was not moving smoke directly from the fire to the AQS 330 ( Fig. 3). For example, only 5 % of rows in the afternoon data indicated that wind sampled at the AQS was coming directly from the fire for at least 3 of the 6 hours. For wind sampled at the fire, this figure was 11 %.  Figure 4 shows the 20 highest mean PM 2.5 values for each six-hour period for HRBs and WFs. The top PM 2.5 values were much greater for WFs than for HRBs in the afternoon, night and morning (~150 to 200 µgm -3 greater for each). >= 80 % of the top 20 PM 2.5 values for WF for afternoon, morning and night were associated with the 2019-2020 wildfires in NSW, many with the Gosper's Mountain wildfire in the Blue Mountains (Boer et al., 2020).

Highest PM2.5 days 340
The top seven afternoon peaks for WF were > 100 µgm -3 (max= 294 µgm -3 ) but only two of the afternoon HRB peaks were > 345 100 µgm -3 . In the night and morning, there were fewer values > 100 µgm -3 , but larger maximums were recorded for HRB and WF for each period (compared to the afternoon). For each rank position, WF values were greater than HRB values, except in the night model where from positions 3 to 20, the HRB values were higher. More information, including satellite images, weather plots and descriptions, on the conditions surrounding the worst PM 2.5 events for each time period for HRBs and WFs is included in Appendix A. 350

Model results
Daily active fire area, PBLH (fire and AQS), temperature and RH at the fire were among the most important variables in the three models (Fig. 5). Some variables were among the most important in only one or two of the models: wind speed at the fire was the fourth and fifth most important in the night and afternoon models but ninth most important in the morning model. The direct wind variables, distance to fire, AQS coordinates, MSLP, month and fire type were all of moderate to lower permutation 360 importance in each model. Partial dependence plots (Fig. 6) indicated that for all models, there was a sharp increase in predicted PM 2.5 when the AQS-tofire distance was below ~20 km, with the morning model displaying the sharpest rise in PM 2.5 as the distance decreased. This effect is despite distance being of middle to lower permutation importance (Fig. 5). Partial plots indicated PM 2.5 increased as fire area increased, particularly in the 0 to 2500 ha range, which is where most training observations were situated (Fig. 3). 365 There was a very large PM 2.5 increase above 10000 ha in the morning and afternoon models, although there is uncertainty here due to a small number of training observations > 10000 ha (Fig. 3). The shape of the PBLH effect differed for each model between the fire PBLH and AQS PBLH. At the AQS, there was a strong negative effect of PBLH (lower PBLH = higher PM 2.5 ), particularly in the night and morning models < 500 m. At the fire, each model had peak PM 2.5 at low and high values of PBLH. For the night and morning models, PM 2.5 peaked when fire PBLH was < ~200 m, with a smaller rise > ~800 m. For 370 the afternoon model, the largest peak was when fire PBLH was high (> ~1500 m), with a smaller rise when < ~500 m. For RH at the fire, predicted PM 2.5 below ~50 % RH was much higher than when RH was above 50 % in the morning and night models.
For wind speed, effects varied between the fire and AQS and with the time period: lower wind speed at the AQS was associated with higher PM 2.5 in all models but at the fire low and high (particularly for the night model) wind speeds were associated with higher PM 2.5 . 375 We calculated model accuracy statistics for the training set (OOB predictions) and the independent test sets and for HRB and WF subsets of each. From the combined statistics, Pearson correlations between predictions and observations (r) for training and test sets ranged from 0.67 to 0.83 (Table 2, Fig. 7). For the statistics by fire type, r was higher for WF than for HRB. For WF, r was 0.7 to 0.88 on the training and test sets. For HRBs, r was 0.59 to 0.69 on the training and test sets. NME for all combinations of training/test set and fire type ranged between 33 % and 39 %, with the lowest NME for the WF subset from 380 the afternoon model (~33 % for training and test set). The NBE indicated that generally there was a slight over-prediction bias that ranged from ~1 % to ~2 %, with a maximum of 6.95 % for WF for the night model test set. The night model had underprediction bias for HRBs on the test set (Table 2, Fig. 7).
The models had large under-predictions for the largest PM2.5 values and a few large over-predictions (Fig. 7). NBE calculated on data that included only where observed PM 2.5 was >= 20 µgm -3 was -30.9 % (training) and -32.8 % (test) for the afternoon 385 were 15 occasions where the model under-predicted in the test set by at least 30 µgm -3 (12 were HRB). The maximum overprediction was by 57 µgm -3 . The morning model had 14 under-predictions on the test set by at least 30 µgm -3 , with the largest 390 under-prediction by 175 µgm -3 for a 2019-2020 WF, although the model correctly predicted this morning as having the highest PM 2.5 in the test set (observed=390 µgm -3 , predicted=215 µgm -3 ). There were 3 over-predictions by at least 30 µgm -3 .
We explored the influence of distance and some selected variables with a series of prediction plots (Fig. 8). PM 2.5 was predicted to increase substantially with decreasing distance within the first 20 km of the fire in all combinations of area, PBLH, RH and temperature in Figure 8. Beyond ~30 km there was minimal to no effect of distance, except in the morning model with a large 395 fire area (Fig. 8a). The effect of temperature at the fire differed between models, such that as temperate increased from 10 to 25 C, PM 2.5 was predicted to decrease in the morning model but increase in the afternoon model. The plots also suggest there is generally a small difference between predicted mean PM 2.5 for WF and HRB for each model once the other predictors including fire area are controlled for.

Discussion 430
Using empirical fire and air quality monitoring station data, we identified important drivers of particulate pollution associated with individual forest fires. The results are important in the context of our first research aim, which was to improve understanding of the fire and weather conditions that influence smoke dispersal and PM2.5 levels. In our models, daily active fire area, PBLH, temperature, relative humidity and wind speed were all important drivers of PM 2.5 from individual fires. The importance of these variables at the fire or at the AQS varied between models. Distance to fire generally had low permutation 435 importance, possibly due to the low number of AQS in the 0 to 50 km range (Fig. 3, Fig. 6). However, partial plots and prediction plots indicated a large influence on model predictions. For example, partial and prediction plots suggested that within 20 km of a fire, PM 2.5 levels rose steeply with decreasing distance. The effect of distance > 50 km was negligible in most cases, suggesting other factors are more important drivers at such distances, although under certain conditions there could be raised PM 2.5 at long distances, such as with higher fire area in the morning model (Fig. 8). Based on Reisen et al. (2018), a 440 1000 ha prescribed burn will emit 160 tonnes of PM 2.5 , enough to fill to exceedance level a cylinder capped by a planetary boundary layer of 500 m to a radius of 64 km. This means there are sufficient particulates available for a distance effect to occur should the weather conditions suit. Other authors have found similar variables to be important in modelling PM 2.5 , including fire size and distance when PM 2.5 was measured within ~10 km of HRBs (Pearce et al., 2012;Price and Forehead, 2021). PBLH was also a consistent predictor of PM 2.5 levels at multiple stations in Sydney during HRB days (Di Virgilio et 445 al., 2018). However, studies such as these have modelled PM 2.5 over smaller scales than we did here or did not attempt to link individual fires to PM 2.5 records. Our data included PM 2.5 measurements up to 150 km from a fire and we built PM 2.5 models using a much larger dataset of fires and PM 2.5 records, which here were from pre-installed permanent AQS. Therefore, the results from our study are more applicable to the individual fire and PM 2.5 relationship across large geographical areas than other studies. 450 Our models suggest the area potentially affected by PM 2.5 from fires is larger than in Price and Forehead (2021), where raised PM 2.5 levels were mostly modelled to be within 5 km of HRBs. Here, our models suggested that raised PM 2.5 levels mostly occurred within 20 km of a fire. Our dataset includes a larger set of fires and includes WFs, which are likely to produce smoke that travels further. In some individual cases in our raw data, fires caused high PM 2.5 levels > 100 km away (e.g. Appendix A Fig A3). Although relatively sparse, analysis using the more remote AQS network is more suited to detecting these longer-455 range effects than when monitors are placed only close to a fire.
Our second aim was to develop predictive models of PM2.5 output from individual forest fires, as a complement to physical models, to improve warnings. There was some success here: r on the test sets indicated moderate to good agreement between predictions and observations: 0.78, 0.70 and 0.83 for the afternoon, night and morning models respectively. The models fit better on the WF portion of the test data (r 0.76 to 0.88) than for HRBs (r 0.65 to 0.69). The better results for WF suggest the 460 models may be more applicable to WFs, e.g. for the issuance of pollution warnings due to WF smoke. An important consideration for using the models for prediction is their accuracy on the largest PM2.5 observations. Events with very high PM 2.5 have the largest health impacts and are therefore the most important to predict, for example, to correctly issue warnings or defer HRBs due to high pollution risk. Our results suggest that, while some predictions for the largest PM 2.5 observations were relatively accurate, the models did not consistently predict larger PM 2.5 events, so may not be suitable as an operational 465 prediction tool without further development.
There are several possible reasons for the biggest outliers and limited accuracy. The AQS network is relatively sparse, being concentrated in greater Sydney, making the distance between any fire and AQS usually large. The mean distance to the closest AQS for each fire-day was 88 km (10 th centile = 31 km). This may partly explain why we did not detect wind direction effects. Price et al. (2012) also did not find significant effects of wind direction when modelling PM 2.5 in relation to MODIS hotspots 470 at similarly broad scales around Sydney and Perth. In contrast, two empirical studies that did detect clear wind direction effects from HRBs, Pearce et al. (2012) and Price and Forehead (2021), placed PM 2.5 monitors close to HRBs, mostly within ~10 km.
The large distances in our data mean smoke was subject to broader weather circulation patterns before reaching an AQS, such as shown in Appendix A. This could create a varying lagged pollution effect that we did not completely account for in our modelling, because smoke may take different amounts of time to reach an AQS depending on circulation patterns. Although 475 we did not focus exclusively on coastal areas, many AQS were in coastal areas, so may have been affected by complex wind patterns. The Sydney basin, for example, can be affected by westerly terrain-related drainage flows, sea breezes and their interaction (Jiang et al., 2017). Differences between land and sea temperatures can influence local wind patterns in coastal areas, creating situations where pollutants emitted overnight or in the morning and blown out to sea are recirculated back over (or near) the source area with a developing sea breeze (Yimin and Lyons, 2003;Levy et al., 2008). Such effects were not 480 accounted for in our study but have been the focus of other research that has used recirculation metrics (Di Bernardino et al., 2022;Wang et al., 2022).
The large distances and sparse network in our data also means that there was a low chance of any particular AQS being downwind from a fire. This is indicated by the wind direction variables being clustered closer to zero (i.e. smoke not blowing from fire to AQS, see Fig. 3) and in cases such as Appendix A Fig. A3, where only two from > 20 AQS detected the smoke 485 from a WF. It may therefore be that the models were mostly optimising for non-smoke-related PM2.5, so it is not surprising that peak events are under-predicted. Our approach is promising, however, more data capturing individual fires burning near monitoring stations is likely required to produce better models. More data could be gathered from the same AQS for another analysis in the future, or by increasing the density of PM 2.5 monitors, either through installing more permanent AQS or via a short-term project that installs a network of temporary AQS in a selected fire-prone area (e.g. Blue Mountains) in times of 490 high-expected fire activity. Some of the variables had interesting non-linear effects. For example, wind speed at the fire during the afternoon was associated with high PM 2.5 both when wind speed was < ~7 km h -1 and > ~15 km h -1 (Fig. 6). Such relationships are due to complex factors. For example, it may be that low wind speeds increase PM 2.5 because previously emitted smoke is more likely to linger, whereas high wind speeds mean that fires are more intense and produce more smoke and particulates. In other words, low 495 wind speed increases smoke concentration at the receiver and high wind speed increases smoke production. The low wind speed effect may be more associated with HRBs, which are conducted in calm weather, and the high wind speed effect associated with WFs. Similar non-linear relationships also exist for other variables, to varying degrees, including PBLH, RH, temperature and MSLP (Fig. 6). Some variables differed in their effects substantially between the fire and AQS. For example, afternoon PBLH at the fire showed increases in PM2.5 at low and high levels, but at the AQS it was only low PBLH that 500 increased PM 2.5 . The PBLH effect at the fire may be similar to the wind effect: low PBLH traps smoke but high PBLH is associated with more active fire behaviour and greater smoke production. Note that there is uncertainty about the strength and directions of the effects at the extremes of the predictor variables, given the lower proportion of observations for model training, as indicated in Figure 6.
Our models predict only small differences between PM 2.5 depending on the fire type variable (HRB or WF), which also had 505 low permutation importance in all three models. Likely, the weather variables and fire area variables included in our model captured most of the differences between HRBs and WFs (e.g. WF on average are larger and burn in hotter windier weather), making the fire type variable mostly redundant in the models. In this case, the models suggest that after accounting for weather and fire size, there are no clear differences in WFs and HRBs in terms of PM 2.5 output. However, other studies have indicated that fundamental differences may exist as WFs inject smoke higher into the atmosphere and consume more fuel per hectare 510 than HRBs (Price et al., , 2018Volkova et al., 2014), thus WF and HRB differences need more investigation.
Our third aim was to make inferences about potential changes in HRB protocols that could reduce PM2.5 impacts. The models indicate the potential combinations of environmental and fire conditions where PM 2.5 is likely to be higher and fire managers must carefully consider whether to undertake HRBs due to PM 2.5 pollution risk. For example, a large HRB < 20 km from a town where PBLH < 300 m during the night and morning (at both fire and receiver site) and < 800 m during the 515 afternoon. When HRBs are > 50 km from a town, a high PM 2.5 impact is much less likely, although certainly still possible (Appendix A). In addition, the HRB area should be a strong consideration as PM 2.5 is predicted to increase as daily active fire area increases between 0 and 2500 ha, although there is uncertainty at larger fire areas because few fires in our data were > 2500 ha (most were < 1000 ha). Note that our fire areas may be an underestimate of total HRB size, as these areas are calculated from VIIRS hotspots, thus based on active fire area at VIIRS overpass times (early afternoon and early mornings), not the total 520 area burnt in a day.
While the models indicate that certain combinations of weather increase PM2.5, this must be weighed with the fact that aspects of HRB implementation cannot always be changed. For example, HRBs are already conducted within the narrow set of weather conditions that allow for ignition and controllable fire spread and often need to be conducted close to populations to have the greatest house protection effect (Clarke et al., 2019). Due to the complex effects and lower predictive accuracy for HRBs, it is 525 difficult to make precise predictions from the models for individual fires. A more detailed model would be required to identify the weather conditions that would allow an HRB to be safely conducted and also for PM2.5 to be low. An assessment that combines predictions from our model of lower-risk PM 2.5 days with a model that predicts the occurrence of within-prescription HRB burning days (Clarke et al., 2019) may be useful to assess the number of overlapping days, i.e. HRB days with low PM 2.5 risk. The effects of different burning strategies, such as breaking a large burn up into multiple blocks, are unknown and could 530 potentially worsen PM 2.5 . Here we did not assess different strategies, and only analysed cases where one fire was burning at a time, not when multiple fires were burning around the same AQS at once. This is a significant limitation of the study, as the smokiest HRB days likely occur when multiple fires are burning at once and/or fires burn for longer periods. Price and Forehead (2021) also suggested overnight burning may have led to the largest PM 2.5 exceedances that they recorded using lowcost monitors near HRBs. Pearce et al. (2012) found burn duration to be an important predictor during their work also 535 monitoring PM 2.5 close to HRBs. The effect of total fire load in a region, i.e. total area of all fires, and regional weather conditions was the subject of separate research (Storey and Price, 2022).

Conclusion
Understanding how individual fires, both wildfires and hazard reduction burns, influence ambient PM 2.5 concentrations is important to allow for proper risk analysis, burn scheduling and issuance of warnings. Our models provide important insights 540 into the influence of weather and fire variables on PM 2.5 concentration from individual fires. We found that daily active fire area, PBLH, temperature and RH all have strong influences, with the effects of the variables varying depending on whether it is measured at the fire site or the receiver location (here, the AQS). The models improve our understanding and may have a place during operational predictions. However, accuracy is similar to existing models, so could be used as a complement.
Further development to improve accuracy would benefit the operational deployment of the models, particularly given the lower 545 correlations between observations and predictions for HRBs. However, our approach is promising and would likely produce better models with a larger set of data, where more cases of single fires near AQS could be found. Increasing the density of PM2.5 monitors (permanent or temporary during fire seasons) would also provide better data to improve the resulting models.
Producing broader scale models of regional level PM 2.5 from regional level fire and weather may be a useful alternative approach for producing operational models. 550 555 560 565

Appendix A 570
This appendix contains case studies of large PM 2.5 exceedance events present in the data used for modelling in the main text.
The purpose is to detail specific events and highlight factors that may have influenced PM 2.5 patterns across the different AQS.
The appendix is organised as seven panel figures of seven different events that each have images and a description. The events selected are the six highest mean 6-hour values from the combinations of fire type (WF, HRB) and period (afternoon, evening, morning), and also the second highest value for afternoon WF, which is included to highlight interesting coastal wind 575 behaviour. Note that the values used in modelling are from AQS data for which only one fire was active within 150 km of the AQS for that day. Higher values were recorded on days with multiple fires, but these are not analysed in this paper. Each figure contains: • Panel (a) in each figure has a background Himawari 8 satellite image for one single hour (time in black text at top) during the relevant time period, with the fire centroid also indicated by an orange circle and general fire area in blue 580 polygon. The background image is overlaid with wind speed (red numbers and red arrow length) and wind direction (red arrow direction) from Bureau of Meteorology weather stations and PM2.5 recorded at all AQS within the image extent at that hour (black circles and white text, larger PM 2.5 value means large circle). The AQS with the highest mean six-hour value is indicated by a red star (same AQS as general location map in panel b). AQS that had multiple fires nearby are not included. Note one extra Himawari image is included for WF night to aid in the 585 description (panel e). Himawari images are provided by Japan Aerospace Exploration Agency (JAXA) and were downloaded from the JAXA P-Tree System (https://www.eorc.jaxa.jp/ptree/terms.html). • Panel (b) in each figure is a map of the general fire location, represented by an orange circle around the fire centroid, with circles representing AQS locations coloured by their mean PM2.5 value (µgm -3 ) for that six-hour period. The highest station values are indicated by the red text and red star. 590 • Panels (c) and (d) in each figure are 10 m and 700 hPa gridded wind speed and direction for the same hour as the Himawari image, sampled from ERA5 gridded reanalysis data. Black arrows indicate wind speed and direction, with longer/larger arrows indicating higher wind speed. The orange fire circle is also in these images for reference. The black solid line is the Australian coastline. The results contain modified Copernicus Climate Change Service information 2021. Neither the European Commission nor ECMWF is responsible for any use that may be made of the Copernicus information or data it contains.

Competing interests
The authors declare that they have no conflict of interest.