Interactive Aerosol Feedbacks on Photolysis Rates in the GEM- MACH v2.4 Air Quality Model in Canadian Urban and Industrial Areas

The photolysis module in Environment and Climate Change Canada’s on-line chemical transport model GEMMACH (GEM: Global Environmental Multi-scale – MACH: Modelling Air quality and Chemistry) was improved, to make use of the on-line size and composition-resolved representation of atmospheric aerosols and relative humidity in GEM-MACH, 10 to account for aerosol attenuation of radiation in the photolysis calculation. We coupled both the GEM-MACH aerosol module and the MESSy-JVAL (Modular Earth Sub-Model System) photolysis module, through the use of the on-line aerosol modeled data and a new Mie lookup table for the model-generated extinction efficiency, absorption and scattering cross sections of each aerosol type. The new algorithm applies a lensing correction factor to the black carbon absorption efficiency (core-shell parameterization) and calculates the scattering and absorption optical depth and asymmetry factor of black carbon, sea-salt, 15 dust, and other internally mixed components. We carried out a series of simulations with the improved version of MESSy-JVAL and wildfire emission inputs from the Canadian Forest Fire Emissions Prediction System (CFFEPS) for two months, compared the model aerosol optical depth (AOD) output to the previous version of MESSy-JVAL, satellite data, ground-based measurements and re-analysis products, and evaluated the effects of AOD calculations and the interactive aerosol feedback on the performance of the GEM-MACH 20 model. The comparison of the improved version of MESSy-JVAL with the previous version showed significant improvements in the model performance with the implementation of the new photolysis module, and with adopting the online interactive aerosol concentrations in GEM-MACH. Incorporating these changes to the model resulted in an increase in the correlation coefficient from 0.17 to 0.37 between the GEM-MACH model AOD one-month hourly output and AERONET (Aerosol Robotic Network) measurements across all the North American sites. Comparisons of the updated mod el AOD with 25 AERONET measurements for selected Canadian urban and industrial sites specifically , showed better correlation coefficients for urban AERONET sites, and for stations located further south in the domain for both simulation periods (June and January 2018). The predicted monthly averaged AOD using the improved photolysis module followed the spatial patterns of MERRA2 re-analysis (Modern-Era Retrospective analysis for Research and Applications Version 2), with an overall under-prediction of AOD over the common domain for both seasons. Our study also suggests that the domain-wide impact of direct and indirect 30 https://doi.org/10.5194/gmd-2021-172 Preprint. Discussion started: 8 July 2021 c © Author(s) 2021. CC BY 4.0 License.

calculated the absorption amplification (ratio of absorption by core-shell black carbon to pure black carbon with the same carbon mass) for a wide range of core-shell thickness, using an implementation of the Bohren and Huffman (1983) Mie scattering algorithm in MatLab (Mätzler, 2002) at 550 nm. They identified five distinct geometric regimes for different black 65 carbon (core) and shell sizes, and calculated the best fit for the absorption amplification for each individual regime , Table 2).
To date, the estimates of AOD in atmospheric models have been based on one or a combination of different mixing states of aerosols. The variation in the resulting aerosol optical properties from the atmospheric models is associated with the assumptions regarding the methods used in AOD calculations, aerosol mixing states, density, refractive index and hygroscopic 70 growth, with the most important factor being the choice of the mixing states of aerosols (Curci et al., 2015). The latter accounts for 30 to 35% of the uncertainty in estimation of AOD and single scattering albedo (Curci et al., 2015). Other studies, e.g., Barnaba et al. (2010), found different spatial patterns in AOD versus surface particulate matter, highlighting the sensitivity of calculated AOD to aerosol vertical profile rather than the aerosol surface concentrations.
The radiative transfer module in chemical transport models contains parameterizations for extinction efficiency (the sum of 75 scattering and absorption efficiencies), single scattering albedo (the ratio of scattering efficiency to total extinction efficiency) and asymmetry factor (the angular direction of the scattered radiation by particles or gases) for each particle type, and calculates scattering and absorption coefficients (a measure of photon scattering and absorption by particles) to predict the radiative state of the atmosphere. AOD is calculated by integrating the extinction of the solar beam due to aerosols over the vertical column.
These optical effects of aerosols may also influence the shorter wavelengths associated with atmospheric gas photolysis, 80 influencing atmospheric reactivity. These processes may be harmonized in an on-line chemical transport model, such as Monitoring of Protection Visual Environments) observations (Gan et al., 2014a;Gan et al., 2015). The aerosol light attenuation method in their model was based on Mie and core-shell scattering (Gan et al., 2014b), and the model used online aerosol feedback on radiation and photolysis (sulfate, nitrate, ammonium, dust and organic aerosols) (Gan et al., 2015). Although their simulations showed the overall observed trends of AOD from SURFRAD, the magnitude of simulated AOD was the effects of different mixing states of black carbon (volume-averaged, core-shell, and externally mixed, as well as the Maxwell-Garnet mixing rule, in which black carbon is assumed to be present in randomly distributed inclusions) on aerosol scattering and absorption properties, for wavelengths between 250 and 700 nm, using an off-line approach of the Aerosol Simulation Program (ASP v2.1, Alvarado et al., 2015). ASP is a single-box aerosol model, with modules to calculate aerosol thermodynamics, gas-to-aerosol mass transfer (condensation/evaporation), coagulation of aerosols, and aerosol optical 120 properties (Alvarado et al., 2016). Using the instruments of the NASA Langley Aerosol Research Group (LARGE; Anderson et al., 1998), they showed that the use of a core-shell mixing state for black carbon, especially for fresh biomass burning episodes, led to the overestimation of aerosol absorption by 29% to 35%, with insignificant dependence on the wavelength, while an external mixture assumption led to the underestimation of aerosol absorption, with a strong dependence on wavelength. Their collected observations suggested using an externally mixed black carbon for the fresh smoke observations, 125 and an internally mixed core-shell approach for the aged Arctic haze and the anthropogenic pollution. Their implementation of a variable mixing state resulted in an average overestimation of aerosol absorption of 10% at 470 nm, 17% at 532 nm and 23% at 660 nm. The mixing state of aerosols has a key impact on radiative transfer, with black carbon's ability to absorb significant amounts of incoming short wavelength light and re-emit this energy as longer wavelengths, resulting in its identification as a short-term 130 climate forcer (IPCC, 2018). However, the mixing state of black carbon, and the impact of that mixing state on the radiative properties of atmospheric aerosols, varies in the literature. Tombette et al. (2008) suggested that the mixing state of black carbon presents an insignificant effect on aerosol optical thickness (AOT) calculations (RMSE difference < 10 -4 ). Liu et al. (2017) recommended using an absorption enhancement in order to account for optical lensing associated with biomass-burning emissions, and no-absorption enhancement for fresh traffic sources. In another study, over 10 European AERONET (Aerosol 135 Robotic Network) sites, Péré et al. (2009) found the mean modelled core-shell single scattering albedo (SSA) provided the closest match to the corresponding measurements, with the spatio-temporal correlation coefficient of 0.51 (compared to 0.04 and 0.35 for the internally homogeneous and externally mixed particles respectively), and therefore is the best representation for simulating anthropogenic and/or biomass burning emissions.
In the work which follows, we make use of the Global Environmental Multiscale -Modelling Air-quality and Chemistry 140 (GEM-MACH) model. The atmospheric chemistry module in GEM-MACH has been included as an on-line component of the core weather-forecast model (GEM) (Côté et al., 1998a, b;Girard et al., 2014;Charron et al., 2012), and consists of airquality processes, including computationally-efficient ADOM-II (Acid Deposition and Oxidant Model, version 2) gas-phase chemistry mechanism with 47 species (Lurmann et al., 1986;Fung et al., 1991), aqueous phase and heterogeneous chemistry, wet and dry deposition, aerosol-cloud processes, and aerosol microphysics (Gong et al., 2003a, b;Moran et al., 2010;Makar 145 et al., 2015a, b). The model's aerosol distribution is based on 12 particle size bins. The aerosol species in GEM-MACH consist of eight components within each size bin: ammonium, sulfate, nitrate, sea salt, crustal material, black carbon, primary organic aerosol (POA) and secondary organic aerosol (SOA). The aerosol and microphysical parameterizations include particle nucleation, condensation and coagulation (Gong et al., 2003a, b), gas and particle dry deposition (Zhang et al., 2001;Makar et al., 2018), cloud processing and in-cloud aqueous-phase chemistry (Gong et al., 2006), direct and indirect feedback options 150 (Makar et al., 2015a, b;Gong et al., 2015) and equilibrium inorganic gas-aerosol partitioning (Makar et al., 2003). The model can be used with either one-way coupling (meteorology drives the chemistry) or with two-way coupling (which enables the model-generated aerosols' impact on radiative transferaerosol direct effect, and on radiative transfer via cloud formationaerosol indirect effect; Makar et al., 2015 a,b). However, the default configuration of GEM-MACH's photolysis calculations makes use of an a priori lookup table as a function of solar zenith angle and height. Here, we update this module, examine the 155 effects of photolysis on the aerosol feedbacks, and show their relative importance to model performance. Table 1 represents the different options for aerosol optical calculations in the current version of GEM-MACH. The original, pre-calculated, clear-sky J-value lookup table in GEM-MACH is a function of solar zenith angle and height. The photolysis rates are calculated using the input data of Peterson (1976) and the radiative transfer model of Dave (1972), with cross-sections and quantum yields taken from DeMore et al. (1988) (Kelly et al., 2012). The model uses the online cloud fraction and liquid 160 water content from the GEM model to scale the pre-calculated clear-sky J-values, based on an algorithm by Chang et al. (1987).
The size distribution and number density profile of aerosols used in the lookup table generation were based on Braslau and https://doi.org/10.5194/gmd-2021-172 Preprint. Dave (1973), and the refractive index of the aerosols were assumed to be independent of height and wavelength, with a single uniform value of 1.5-0.01i (Yamamoto and Tanaka, 1972  The calculations of aerosol optical depth, single-scattering albedo, and asymmetry factor used by the GEM weather forecast 170 model, also used in the offline version of the radiative transfer module in GEM-MACH (Li and Barker, 2005) with no-feedback configuration, are based on a climatology produced by , with specified aerosol loading for continents and oceans, a maximum value at the Equator and a decrease towards the poles. The solar absorption properties are only https://doi.org/10.5194/gmd-2021-172 Preprint. Discussion started: 8 July 2021 c Author(s) 2021. CC BY 4.0 License. assumed to be affected by aerosols in the clear-sky atmosphere (Markovic et al., 2008). In contrast, the on-line version of the radiative transfer module in the feedback-enabled version of GEM-MACH makes use of a homogeneous aerosol mixture Mie 175 scattering (Bohren and Huffman, 1983) code for meteorological radiative transfer calculations, and the model-generated aerosols in the feedback mode (Makar et al., 2015 a, b) are used for radiative transfer calculations. interface of Jöckel et al. (2005). The original formulation of photolysis rates calculations was developed by Landgraf and 180 Crutzen (1998) and has been since adopted in several atmospheric models (Sander et al., 2014). The actinic flux calculations include the effects of aerosols and clouds on photolysis rates. The optical data for cloud scattering and absorption are adopted from Slingo (1989) (Sander et al., 2011a). The previous off-line application of MESSy-JVAL in GEM-MACH (Chen et al., 2020), made use of a climatology of aerosol concentration (a constant aerosol vertical profile above land and a different vertical profile over water). Cloud radiative properties and cloud fraction were calculated online at runtime (Chen et al., 2020). The 185 focus of this study is to improve the representation of aerosol optical properties in MESSy-JVAL, and to determine their impact on model performance.
In the following sections, we describe the methods we used to improve the representation of aerosol optical properties in conclusions of this study are presented in Section 4.

Overview
In order to improve the photolysis module and calculate aerosol feedbacks more accurately, we developed a new Mie lookup table for aerosol optical properties which was accessed within the improved version of MESSy-JVAL. An initial lookup table  200 was calculated using the refractive indices of six representative aerosol components and the Mie scattering code within the VECTOR model (McLinden et al., 2002). Using the hygroscopic growth factor of each aerosol type, we calculated the dry size parameter (Section 2.2, equation1 of the initial lookup table, which was then interpolated into GEM-MACH dry size parameter to calculate the final lookup table as a function of GEM-MACH particle size bins and wavelengths. The calculated aerosol composition and size in GEM-MACH was used as input for photolysis rate calculations, and a hybrid aerosol mixing state was 205 assumed for size bins containing different concentrations of black carbon to calculate the optical properties of four aerosol categories: black carbon, sea-salt, dust and other internally mixed components. Currently, there is a size-resolved (sectional) representation of atmospheric aerosol particles in GEM-MACH, which may be used for determining the impact of aerosol feedbacks on radiation and photolysis. When GEM-MACH's aerosol feedbacks are enabled, the aerosol direct effect makes use of a Mie scattering approach and the assumption of a binary mixture between 210 dry aerosols and water for the complex refractive index (Makar et al., 2015a). For photolysis calculations as carried out here, we ma de use of GEM-MACH's calculated aerosol composition and size as input for photolysis rate calculations. The new algorithm uses the eight dry chemical components of GEM-MACH aerosol (sulfate, nitrate, ammonium, primary organic matter, secondary organic matter, black carbon, dust and sea-salt) and reads in the data from the updated lookup table (see sections 2.2, 2.3, following) to calculate the scattering and absorption optical depth and asymmetry factor of black carbon, sea-215 salt, dust and other internally mixed components. We used the black carbon particle size in GEM -MACH as an indicator of the mixing state of the internally mixed particles. For each GEM-MACH size bin, the mass fraction of black carbon to the total mass of all the other particles in that bin was calculated. Since the externally-mixed black carbon particles are most common near combustion sources, if a GEM-MACH particle size bin contains a black carbon mass fraction that is near or larger than that of a typical combustion emission particle (the black carbon mass fraction of more than 40% (Stevens and Dastoor, 2019), 220 the particle bin was considered to be mostly black carbon, no absorption amplification factor was applied to black carbon photolysis rate calculations, and an external particle mixing state for that size bin was assumed (i.e., black carbon as a separate particle, and a volume-averaged homogeneously mixture of ammonium sulfate, ammonium nitrate, primary and secondary organics). If the black carbon mass fraction for a particle size bin was less than 40%, we assumed that black carbon is coated with other internal particle components, and a core-shell configuration was used for that size bin; the black carbon forms a 225 core and other internal particle components (primary and secondary organic matter, ammonium sulfate, and ammonium nitrate), form a shell. In this case, the black carbon is mixed with other condensed or coagulated components, the bin is more aged, and we apply a lensing correction factor to the black carbon absorption efficiency recommended by .
It should be noted that sea salt and dust were not included in the assumptions of the internally mixed particles, and in both cases mentioned above, those aerosols were considered as separate particles. These underlying assumptions were used for the 230 calculations of the aerosol optical properties for four independent aerosol groups: sea-salt, dust, black carbon and the internally mixed particles (See Fig. 1). We describe our methodology in more detail in the following sections.

Developing a new aerosol optical properties look-up table for MESSy
In order to update the aerosol effects in the MESSy photolysis module in GEM-MACH, we calculated a new lookup using the Mie scattering code within the VECTOR model (McLinden et al., 2002) for extinction efficiency, single scatter albedo and 235 asymmetry factor for six aerosol types, which within the lookup table are treated as pure-composition particles of sea-salt, black carbon, dust, ammonium sulfate, ammonium nitrate and organic carbon. The initial version of the new lookup table was calculated for single components, each with its own water uptake properties derived from the literature. One dimension of the https://doi.org/10.5194/gmd-2021-172 Preprint. Discussion started: 8 July 2021 c Author(s) 2021. CC BY 4.0 License. table is the wet particle radius range for a logarithmically expanding set of cut-sizes of the aerosol distribution of GEM-MACH (0-10, 10-20, 20-40, 40-80, 80-160, 160-320, 320-640, 640-1280, 1280-2560, 2560, 240 five wavelengths (200, 300, 400, 600, 1000 nm) and seven different relative humidity levels (0%, 50%, 70%, 80%, 90%, 95%, 98%, 99%) for all aerosols except dust, which was assumed to have no water uptake. The effective wet particle radii were calculated based on a power-law distribution (Hansen and Travis, 1974) of each particle range. The relative humidity determines the water fraction for all other aerosol types. The selection of wider spacing at lower relative humidity and longer wavelength was due to the growing dependency of optical properties to increasing relative humidity and decreasing 245 wavelengths. We assumed a flat distribution of aerosol radii within each bin size. The complex refractive indices of the watersoluble inorganic aerosols, namely ammonium sulfate, ammonium nitrate and sea salt were calculated using the FORTRAN software developed by Andrew Lacis (https://www.giss.nasa.gov/staff/mmishchenko/ftpcode/lacis/lacis_refrac.rhwmri.f) which has been used in many recent studies (Jeong, 2020;Bozzo et al., 2020;Escribano et al., 2014;Schuster et al., 2009). In this model, all aerosols are treated as homogeneous mixtures. Parametric formulas are derived for the changes in the real part 250 of the refractive index, specific density, size, and water activity as functions of mass fraction. These formulas are used to interpolate spectrally between the refractive indices of dry ammonium sulfate, sea salt, and ammonium nitrate and those of water. For ammonium sulfate, the complex refractive index is from . The water activity for this aerosol type is based on Tang and Munkelwitz (1994). For ammonium nitrate, the real and imaginary parts of the refractive index are from Tang et al. (1981) and Gosse et al. (1997), and the water uptake is from Tang (1996). The complex refractive index of sea salt 255 is based on Shettle and Fenn (1979), which relies on Dorsey (1940).
For dust, no water uptake has been assumed, and thus no dependency on the relative humidity. The complex refractive index of dust is assumed to be independent of particle size and is obtained from the VECTOR model's library of refractive indices with the real part of the dust refractive index varying between 1.55 and 1.57 and the imaginary component increasing monotonically from 0.004i at 1000 nm to 0.025i at 200 nm. For organic carbon, the complex refractive index is wavelength-260 dependent for all relative humidities and was extracted from Fast-JX, the photolysis mechanism used in GEOS-CHEM (Goddard Earth Observing System (GEOS) of the NASA Global Modeling and Assimilation Office). Fast-JX v7.0 calculates aerosol extinction efficiencies at five wavelengths, 200, 300, 400, 600 and 999nm. The file jv_spec_mie.dat (http://rain.ucis.dal.ca/ctm/CHEM_INPUTS/FAST_JX/v2020-02/jv_spec_mie.dat) in Fast-JX v7.0 contains the quantum yields and aerosol cross sections for photolysis. Real and imaginary parts of the refractive index of organic carbon at 98% 265 were obtained by linear interpolation between the indices at RHs of 95 and 99%. For black carbon, the dependence on relative humidity was based on hygroscopic growth factors (HGF; the ratio of the wet particle radius to the dry particle radius) obtained from Lei et al. (2014) (Table 2). The RH-dependent refractive index (nRH) is a HGF 3 -weighted mean of refractive indices of pure black carbon (nBC) and pure water (nwater):  Since the initial version of the new lookup table was calculated based on wet particle sizes, we used the hygroscopic growth 280 factor of each aerosol type (Laskina et al., 2015, Ting et al., 2017, Zamora et al., 2013

285
Using the dry radius size and the wavelengths in the initial lookup table (200, 300, 400, 600, 1000 nm), the initial size parameter was calculated using: where α i is the initial size parameter, r i is the dry radius of the particle, and λ i is the wavelength in the initial lookup table.
Using the same formula, the size parameter for the final lookup table was calculated based on GEM-MACH average dry 290 particle size bins (7.5 nm, 15. nm, 30. nm, 60. nm, 120. nm, 240. nm, 480. nm, 960. nm, 1920. nm, 3840. nm, 7680. nm and 15360. nm) and wavelengths (205 nm, 287 nm, 302 nm, 309 nm, 320 nm, 370 nm and 580 nm). The final lookup table components as a function of particle size, were calculated by linear interpolation of GEM-MACH size parameter from the size parameter values in the initial version of the new lookup table (see Fig. 1). The initial and interpolated extinction efficiency of ammonium sulfate versus dry aerosol size parameter is illustrated in Fig. 2. The final interpolated lookup table of extinction 295 efficiency, single scatter albedo and asymmetry factor was used to calculate the absorption and scattering cross section and asymmetry factor of seven pure aerosol types (sea-salt, black carbon, dust, ammonium sulfate, ammonium nitrate, primary organic carbon and secondary organic carbon) in GEM-MACH. The data in the lookup table is sorted by increasing size parameter values (based on the dry aerosol sizes) for each aerosol type and relative humidity, i.e., the optical properties in the lookup table depend on the water uptake for the given relative humidity. 300

Updating the photolysis module in the GEM-MACH in-line chemical transport model:
After calculating the new lookup table, we modified the MESSy-JVAL code to use the new lookup table, along with the 305 calculated aerosol composition and size by GEM-MACH, as input for photolysis rate calculations. The new updated code uses the eight dry chemical components of GEM-MACH aerosol feedback (sulfate, nitrate, ammonium, black carbon, primary organic matter, secondary organic matter, dust and sea-salt) and the data in the updated lookup table to calculate scattering and absorption optical depth and asymmetry factor of black carbon, sea-salt, dust and other internally mixed components. The volume concentration of each aerosol type (m³ of each aerosol per m³ of air) was calculated using the GEM-MACH predicted 310 mass concentration of each aerosol type (µg kg -1 ). The number concentration of each aerosol (number of aerosols per 1 m³ of air) was calculated by dividing the volume concentration of the aerosol by the volume of each size bin. We used the predicted mass concentration of nitrate in GEM-MACH, and the molecular weight of ammonium nitrate to calculate the mass concentration of ammonium nitrate. Finally, to conserve the mass of ammonium, and since ammonium sulfate, ammonium bisulfate and letovicite have very similar refraction indices, the remaining mass of ammonium and the mass concentration of 315 sulfate were used to calculate the mass concentration of ammonium sulfate.
In order to implement the core-shell parameterization where the mass fraction of black carbon is less than 40% in a particle bin in GEM-MACH, we calculated the number of particles with a black carbon core (NBC), and the mass concentration of black carbon (MBC) for those specific size bins. Using the two values, the dry black carbon core size was calculated as follows: where r BC is the dry radius of the black carbon core of a particle and ρ BC is the density of black carbon ("void-free black carbon core density of 1.8 g cm -3 (McMeeking et al., 2010)"). The total wet radius of a core-shell particle was calculated using the volume-weighted hygroscopic growth factor of the components in the core-shell particle and the total wet radius of the core (black carbon) and shell (secondary organics, ammonium sulfate and ammonium nitrate) in that particle. This information was used to calculate the size parameter of the black carbon core (α bc ) and the total particle size parameter (α total ). The absorption 325 amplification factor for the case of black carbon core-shell, was calculated using the core-shell parameterization by . with the observationally-constrained maximum threshold of 1.93 . As described below, these parameters are used to provide an amplification factor based on previous core-shell Mie scattering calculations carried out by . Sea-salt, dust and black carbon (when its mass fraction was > 40%) aerosols were assumed to be externally mixed at all times 330 -their effective scattering coefficient(sca c ), absorption coefficient(abs c ) and asymmetry factors (asym c ) of each of these aerosols were calculated using the elements of the lookup table and GEM-MACH's predicted concentrations for these aerosol species: where subscript c denotes each aerosol category, for each size bin (i) in GEM-MACH (i = 1 to 12), Q ext is the extinction efficiency of each aerosol type, N i is the number concentration of each aerosol type (cm −3 ), r i is the radius of bin i (cm), ssa i is single scatter albedo of each aerosol type, ag is the asymmetry factor, and mf i is the fraction of the mass concentration of each aerosol type to the total mass concentration of all particles (the asymmetry factor was weighted by mass fraction for pure 340 particles). In the case where the mass fraction of black carbon for a particle size bin was less than 40%, the effective scattering and absorption coefficient of black carbon was calculated using: where r BC is the radius of black carbon core particle and is the absorption amplification factor by , 345 based on the black carbon core and shell size parameters. For the fourth aerosol category in each GEM-MACH size bin (volume-averaged internally mixed particles), we calculated the volume fraction(vf) of each component (ammonium sulfate, ammonium nitrate, primary and secondary organics) to the total volume of the internally mixed particles: vf AS = vi (AS) vi (IM) , vf AN = vi (AN) vi (IM) , vf OC = vi (OC) vi (IM) , vf PC = vi (PC) vi (IM) (10) 350 where vi is the volume concentration of each aerosol in the internally mixed particles , IM stands for the internally mixed particles, AS is ammonium sulfate, AN is ammonium nitrate, OC is secondary organic matter, and PC is primary organic matter. Equations 11 to 13 were used to calculate the effective scattering and absorption coefficients and asymmetry factor of internally mixed particles: (11) 355 where the indices (j = 1,2,3,4) correspond to the aerosol type inside the internally mixed particles (ammonium sulfate, ammonium nitrate, secondary organic carbon, primary organic carbon). Note that the calculations using equations 4 to 8 and 11 to 13 were done for each horizontal grid-point, vertical level, wavelength and relative humidity. 360 We performed a linear interpolation of the relative humidity in the lookup table for all the aerosol types, with the exception of dust aerosol (which had no water uptake), to calculate the scattering and absorption coefficients at each given relative humidity in GEM-MACH. Scattering and absorption optical depth for each model layer and aerosol category were calculated by the following formula: In the above summations, c represents each aerosol category. Integrating the total scattering and absorption optical depth through the entire column at each grid-point and for each wavelength (λ), gives the total modelled AOD at that wavelength.
This information was used to calculate J-values for photolysis reactions depending on the attenuation of the radiation stream 375 by particles.

Simulation Setup and Emissions
For this study, we used GEM-MACH v2.4, with 12-bin average size distribution of particles, and online aerosol feedbacks between weather and air-quality (Makar et al., 2015a, b). The model domain covers most of continental Canada and United States with a horizontal grid-spacing of 10 km, 80 hybrid vertical levels with a model top at 0.1 hPa, 15-minute operator 380 splitting time step and a 1-hour output time step (Fig. 3). The meteorological initial and boundary conditions for our 10 km horizontal resolution simulations were from the operational Regional Deterministic Prediction System, ECCC's operational
For ground-based measurements, we used quality assured AERONET Sun photometer measurements of AOD at 500 nm for North American sites (Fig. 4). The Sun triplet measurements are performed every 15 minutes for older model 4 instruments and every 3 minutes for newer model 5 and CE318-T instruments (Giles et al., 2019). The AERONET data used in this study is cloud screened according to the AERONET V3 algorithm described in Giles et el. (2019). We used the Angstrom Exponents 430 provided by the sun photometer within AERONET data to evaluate the model AOD at 580 nm. The Angstrom exponent (AE) represents the wavelength ( ) dependency of AOD, provides basic information about the size distribution of aerosols and is calculated by the following formula: To calculate AOD at 580 nm, we used a variation of the above formula: where AOD1 and AOD2 are the aerosol optical depth at 500 and 580 nm, and λ 1 and λ2 are 500 nm and 580 nm wavelengths respectively. We used (440 nm -675 nm) Angstrom Exponent to obtain the observed AOD at 580 nm, which was used to compare with GEM-MACH output AOD calculated as described in section 2.3.

Simulation Plan
We  Table 3. It should be noted that the spatial and temporal resolutions are the same (as described in section 2.4) for all the simulations.

Comparison of base and improved MESSy-JVAL
In this section, the results of two simulations with the GEM-MACH model for the month of June, 2018 is compared to AERONET sun photometer measurements of AOD at 500 nm at four North American sites (shown in Fig. 4): simulation Sb 455 with the previous version on MESSy-JVAL (base) and simulation S1 with the improved photolysis module (see Table 3). The

Evaluation of the improved MESSY-JVAL
In this section we evaluate the GEM-MACH output with the improved photolysis module against the observations and assess the effects of (1) AOD calculations (versus an assumed aerosol optical depth of zero) and (2) interactive aerosol feedback with 475 the GEM model, on the resulting calculated photolysis rates (see Table 3).

Comparison with observations
As described in section 2.6, we performed six simulations for the months of January and June 2018, using GEM-MACH model with the updated lookup table and photolysis module. We compared the simulated GEM -MACH AOD with AERONET sun photometer measurements of AOD at 580 nm (converted from AOD at 500 nm using (440 nm -675 nm) Angstrom Exponent) 480 for the entire simulation period (January and June 2018) and for all North American sites in Fig. 4.   As illustrated in Fig. 6, the sample size for Toronto is larger than Fort McKay for both seasons. During the sampling period in summer of 2018, there was an instrument malfunction at Fort McKay from June 19 to August 20, which led to having a smaller 495 number of data points at this site compared to Toronto. As shown in these time series, the maximum value of the modeled AOD was underestimated compared to AERONET data for both sites and both simulation periods. For example, the maximum AERONET AOD for the month of June in Fort McKay and Toronto were measured ⁓0.27 and ⁓0.42 respectively, whereas the simulated GEM-MACH AOD were ⁓0.19 and ⁓0.29 for those sites.

Sensitivity test to AOD calculations and aerosol feedback
In this section we evaluate the impacts of AOD calculations on photolysis rates calculations with the improved photolysis module in the GEM-MACH model. We also assess the effects of the model-generated aerosols on the radiative transfer and cloud formation processes in GEM-MACH, through the comparisons of the optical properties of aerosols in the 'feedback' mode (Makar et al., 2015a, b) vs 'no-feedback' mode, in which the feedback mechanisms (interactions between meteorology 535 and chemistry) have been disabled and the model uses climatological parameterizations for the aerosol direct and indirect effects. Figure 10 shows the monthly average percentage difference in photolysis rate coefficients (J(O 1 D) and J(NO2)) with and without AOD calculations in the photolysis module in GEM-MACH model. The top two panels are the simulation outputs with online aerosol feedback between weather and air-quality in the model. We used the output from simulations S1 and S2 to 540 calculate the percentage difference in photolysis rates in Fig. 10(a) and (b), and the output from simulations W1 and W2 in  Fig. 10(e) and (f), there is not a significant change in monthly-averaged photolysis rates (-0.1% to 0.1%) with no online aerosol feedback on weather in the model (Fig. 10 (e) and (f)). The increase in the photolysis rates differences, and the less organized structure of these changes over the domain in simulations with the online aerosol feedbacks, is due to the presence of the direct and indirect effects from meteorological changes such as changes in cloud patterns, and amplified chemistry perturbations due to weather feedback. 550 Figure 11 is the monthly average percentage difference in AOD, J(O 1 D) and J(NO2) in June 2018 and at the lowest model level with and without interactive aerosol feedback on weather in the model (simulations S1 and S3). This figure shows the effects of the interactive online aerosol feedback on the output of the photolysis module. Note that in both simulations we used AOD calculations in the photolysis module. As shown in this figure, AOD changes from -30% to 30% and J-values from -40% to 40% with and without the GEM-MACH predicted aerosol concentrations in optical properties and photolysis rates calculations. 555 The increase of AOD and the resulting decrease of J-values over the continent, might be due to the increase in the cloud droplet density, and/or in-cloud formation of aerosols in the simulations with interactive aerosol feedback in the model. The results from Fig. 10 and Fig. 11 suggest that the impacts of aerosol feedbacks as parameterized in the model are considerably greater than the impacts of the AOD calculations in the model.  . The more-organized changes in photolysis rates in the no-feedback case, is due to including only the direct effect of photolysis on aerosols, whereas, the effect of clouds and in-cloud formation of aerosols in the feedback case leads to more variability in the photolysis rate changes between simulations S1 and S2.      shows two areas of maxima; one is directly over the forest fire plume, similar to the hotspot on GEM -MACH plot (Fig. 14 (c)) and the weaker hotspot over the same area in MODIS-Terra plot (Fig. 14(a)). The second maximum value in MERRA-2 ( Fig.  600 14(b)) is located downwind of the forest fire plume over the Athabasca Oil Sands region of northeastern Alberta, which is not detected in MODIS-Terra plot (Fig. 14(a)), and is more intensified, compared to GEM-MACH values. One possible reason that GEM-MACH did not show high values of AOD over the second hotspot could be the coarse horizontal resolution of the model for these simulations (10 km × 10 km). Another explanation could be potential deficiencies in Oil Sands emissions or aerosol processes in this simulation. 605  Figure 15 shows the aerosol concentrations and photolysis rates cross sections over the black solid line in Fig. 15(c). As shown 610 in Fig. 15(a), O3 is impacted by titration below the PBL (Planetary Boundary Layer), and there is a low concentration of ozone right above the fire plume. Higher concentrations of O3 can be seen downwind below the PBL. As the fire plume continues to dilute with distance downwind, the NOx concentration in the boundary layer (not shown) decreases and eventually reaches an optimal concentration range for efficient ozone production. This is illustrated in Fig. 15(a) as ozone increases to near 80 ppbv about 170 km downwind from the fire. Figure 15(b) shows the depletion of hydroxyl radical (OH) in the fire plume below the 615 PBL and a maxima above the boundary layer, where there is a high concentration of O3 and low concentration of VOC (Volatile Organic Compound) to deplete OH. The increase in OH downwind from the fire plume is in response to the O3 increase, reaching a maximum value of ⁓2.3×10 -4 µg kg -1 (⁓7.5 ×10 6 molecules cm -3 ) in the upper boundary layer ⁓700-800 hPa. The increase in OH concentration results in delayed production of secondary aerosol components. The PM2.5 predicted in the forest fire plume (Fig. 15(c)) reached a maximum value of 65 µg m -3 near the surface and remained above 10 µg m -3 up to 800 hPa. 620 As shown in Fig. 15(d), the fire plume was predicted to penetrate above the boundary layer height due to the black carbon plume mixing up to 600 hPa. The black carbon concentration decay with distance illustrates the extent of dilution of directly emitted PM2.5 components as it mixes horizontally. As shown in Fig. 15(e), the concentration of SOA increases downwind from the fire plume. Similarly, the concentration of nitrate (Fig. 15(f)) increases downwind between 900 hPa and the boundary (a) (b) (c) layer due to secondary production and the long-lived nature of nitrate. The attenuation of radiation by the fresh fire plume is 625 illustrated by the decrease in J(NO2) below the PBL from 1.26×10 -2 s -1 at ⁓700 hPa to 4.52×10 -3 s -1 at the surface and J(O 1 D) from 1.83×10 -5 s -1 at ⁓700 hPa to 5.48×10 -6 s -1 at the surface ( Fig. 15(g) and (h)). The photolysis rates decrease with distance is due to the attenuation of radiation by the directly emitted fire PM2.5 components.
In order to evaluate the effects of aerosol optical properties in aerosol concentrations and photolysis rates in Lac La Loche fire event, we calculated the percentage difference in aerosol concentrations and photolysis rates with and without AOD 630 calculations (simulations S1 and S2). Figure 16 shows the percentage difference of aerosol concentrations and photolysis rates cross sections with and without AOD calculations over the Lac La Loche forest fire event on June 25, 2018 (black solid line in Fig. 15(c)). Figure 16(a) shows the reduction of O3 in the fresh fire plume below the PBL due to the decrease in J(O 1 D), reaching a difference close to 0.9% near the surface. The OH radical concentration difference ( Fig. 16(b)) responds to the O3 change and decreases by up to 3% at the surface due to the AOD feedback on the photolysis rates. This affects the rate of 635 production of secondary aerosol components such as nitrate ( Fig. 16(f)) and secondary organic material (Fig. 16(e)), although there is a small decrease of up to 7% in nitrate concentration in the fresh plume. There is a reduction in J(NO2) and J(O 1 D) ( Fig. 16(g) and (h)) in the fire plume below the boundary layer. As shown in Fig. 16(c) and (d), after including AOD calculations in the model, the concentrations of PM2.5 and black carbon decrease by 5% in the fire plume from the surface to 900 hPa and increase by 5% above 900 hPa and below the boundary layer. 640 Figure 17 is the percentage difference in aerosol concentrations and photolysis rates, with and without AOD calculations and no aerosol feedbacks in the model (simulations S3 and S4). As shown in Fig. 17(c), (d) and (e), without the interactive aerosol feedbacks in the model, there is an insignificant difference (-0.2% to 0.2%) in PM2.5, black carbon and SOA concentrations compared to the difference in concentrations from simulations S1 and S2 (-50% to 50% with aerosol feedbacks). Similarly, there is -1% to 1% difference in O3 concentrations ( Fig. 17(a)) from simulations S3 and S4. The variability in the production 645 of OH ( Fig. 17(b)) by the photolysis of O3 is more considerable (-5% to 5%). The decrease in the OH concentration with primary-aerosol-component photolysis attenuation results in lower secondary aerosol production downwind, which in turn, slowly counters the primary aerosol attenuation and increases the photolysis rates (maximum of < 1%). The increased photolysis difference results in a positive difference between the ozone and OH concentrations. The predicted OH increases downwind reaching a maximum of 1%. 650

Summary and Conclusions 665
A new lookup table for aerosol optical properties based on a Mie scattering code was calculated and adopted within an improved version of the photolysis module in the GEM-MACH in-line chemical transport model, by interpolating the optical properties of aerosols into GEM-MACH wavelengths and size bins. The modified version of the photolysis module makes use of online interactive aerosol feedback for radiative transfer calculations. Additionally, for particle size bins with black carbon mass fraction of less than 40%, a lensing correction factor to the black carbon absorption efficiency based on   and Fort McKay in winter, whereas there is an under-prediction of AOD in Toronto and Egbert for both seasons. In addition, we calculated a higher correlation coefficient between the model and measurements for urban AERONET sites for both simulation periods (e.g., Toronto vs. Fort McKay). Further investigations on the aerosol processes and emissions at these sites is needed to assess the effects of different climatological and meteorological conditions on photolysis rate calculations.
The sensitivity test to aerosol feedback demonstrates the effects of the model predicted aerosol distribution in the modified 690 photolysis rate calculations. As shown in this study, there was up to ±40% change in monthly averaged photolysis rate calculations with and without online aerosol feedback in the model, whereas with no aerosol feedback in the model this change is very small (-0.1% to 0.1%) between the runs with and without AOD calculations. The sensitivity tests to AOD calculations show a monthly average change of ±10% in photolysis rate coefficients over the North American domain, while as shown in the forest fire case study this number can be as high as ±50% in the fire plume. This study also showed the impact of aerosol 695 feedbacks as parameterized in the model to be considerably greater than the effects of the AOD calculations (by a factor of 3 to 4) on the photolysis rates over the entire domain. Furthermore, this study highlights the importance of model simulations of AOD where satellite observations are obscured by cloud cover. As shown in the case study, the simulated AOD over Lac La Loche forest fire of June 25, 2018, indicates a maximum value directly over the fire plume, while the same hotspot is not detected with the same intensity by MODIS-Terra, possibly due to cloud obscuring surface in satellite retrieval. 700 Further investigation of the effects of the improved photolysis module with a nested configuration of GEM-MACH 10-km domain to a 2.5 km Athabasca Oil Sands domain is needed for more detailed comparisons of model output with observations under the Oil Sands Monitoring Program, 2018 aircraft campaign. This new model capability will also enable us to use model predicted AOD on regional scales over the boreal forest to assess future improvements in biogenic VOC emissions and SOA processes in GEM-MACH, by comparing with clear-sky satellite-derived AOD. June 2021). The executable for GEM-MACH is obtained by providing the chemistry library to GEM when generating its executable. Many of the emissions data used in our model are available online at ECCC web page at http://donnees.ec.gc.ca/data/air/monitor/source-emissions-monitoring-oil-sands-region/source-emissions-oil-sands-715 region/?lang=en (ECCC, 2018) and more recent updates may be obtained by contacting Junhua Zhang or Mike Moran (junhua.zhang@ec.gc.ca; mike.moran@ec.gc.ca). The model output is available upon request to Craig Stroud (craig.stroud@ec.gc.ca). The AERONET version 3 AOD datasets are publicly available from the AERONET website (https://aeronet.gsfc.nasa.gov/, last access: 5 February 2020). Terra MODIS level-3 Atmospheric Daily Global Product of AOD, MERRA-2 Aerosol Optical Depth Analysis V5.12.4, and MERRA-2 Aerosol Scattering AOT were obtained from the 720 publicly accessible Goddard Earth Sciences Data and Information Services Centre platform (GIOVANNI: https://giovanni.gsfc.nasa.gov/giovanni/, last access: 1 October 2020).