Development of Aerosol Optical Properties for Improving the MESSy Photolysis Module in the GEM-MACH v2.4 Air Quality Model and Application for Calculating Photolysis Rates in a Biomass Burning Plume

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 10 use of the on-line size and composition-resolved representation of atmospheric aerosols and relative humidity in GEM-MACH, 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 modelled 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 15 parameterization) and calculates the scattering and absorption optical depth and asymmetry factor of black carbon, sea-salt, 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, 20 and evaluated the effects of AOD calculations and the interactive aerosol feedback on the performance of the GEM-MACH 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 25 Robotic Network) measurements across all the North American sites. Comparisons of the updated model AOD with 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 30

showed the model AOD values ranging from a minimum of ⁓0.06 to a maximum of ⁓0.09, while SURFRAD observations had a minimum of ⁓0.125 and maximum of ⁓0.175. Air Resources Board (CARB)"), was conducted over California, one week prior to ARCTAS-B. These evaluations examined 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 120 thermodynamics, gas-to-aerosol mass transfer (condensation/evaporation), coagulation of aerosols, and aerosol optical 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 125 wavelength. Their collected observations suggested using an externally mixed black carbon for the fresh smoke observations, 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. of incoming short wavelength light and re-emit this energy as longer wavelengths, resulting in its identification as a short-term 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 135 emissions, and no-absorption enhancement for fresh traffic sources. In another study, over 10 European AERONET 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. 140 In the work which follows, we make use of the Global Environmental Multiscale -Modelling Air-quality and Chemistry (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, 145 wet and dry deposition, aerosol-cloud processes, and aerosol microphysics (Gong et al., 2003a, b;Moran et al., 2010;Makar 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 150 et al., 2018), cloud processing and in-cloud aqueous-phase chemistry (Gong et al., 2006), direct and indirect feedback options (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 transfer -aerosol 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 155 makes use of an a priori lookup table as a function of solar zenith angle and height. Here, we update this module, examine the 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 160 and quantum yields taken from DeMore et al. (1988) (Kelly et al., 2012). The model uses the online cloud fraction and liquid 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 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 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 Toon and Pollack (1976), with specified aerosol loading for continents 170 and oceans, a maximum value at the Equator and a decrease towards the poles. The solar absorption properties are only 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 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. 175 The new photolysis module in GEM-MACH (MESSy-JVAL) is based on JVAL14-MESSy of Sander et al. (2014). The photolysis module JVAL was adapted by Sander et al. as a module inside the Modular Earth Sub-model System (MESSy) interface of Jöckel et al. (2005). The original formulation of photolysis rates calculations was developed by Landgraf and 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 180 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 focus of this study is to improve the representation of aerosol optical properties in MESSy-JVAL, and to determine their impact on model performance. 185 In the following sections, we describe the methods we used to improve the representation of aerosol optical properties in MESSy-JVAL, followed by the evaluations of the improved photolysis module and the limitations of the model in photolysis rate calculations. Section 2 provides a description of the model configuration, simulation setup and the observations used to evaluate the model output, calculations of the new lookup table for aerosol optical properties, and a description of the revised photolysis rate calculations in GEM-MACH. The comparison of the updated MESSy-JVAL with the base off-line version is 190 presented in section 3, followed by the results of different evaluations of the improved photolysis module: comparisons with observations, evaluations of the impacts of AOD calculations and model-generated aerosols on the photolysis rates calculations, and a case study of photolysis calculations under high pollutant flux emissions conditions. Summary and 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 was calculated using the refractive indices of seven 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 200 parameter (Section 2.2, equation 1) of the initial lookup table, which was then interpolated into GEM-MACH dry size parameter. The final lookup table which was implemented into GEM-MACH is a function of GEM-MACH dry size parameters 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 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. 205 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 210 are enabled, the aerosol direct effect makes use of a Mie scattering approach and the assumption of a binary mixture between dry aerosols and water for the complex refractive index (Makar et al., 2015a). For photolysis calculations as carried out here, we made 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 215 sections 2.2, 2.3, following) to calculate the scattering and absorption optical depth and asymmetry factor of black carbon, seasalt, 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 220 that of a typical combustion emission particle ("the black carbon mass fraction of more than 40% (Stevens and Dastoor, 2019), 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 225 with other internal particle components, and a core-shell configuration was used for that size bin; the black carbon forms a 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 230 cases mentioned above, those aerosols were considered as separate particles. These underlying assumptions were used for the 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 table 235 using the Mie scattering code within the VECTOR model (McLinden et al., 2002) for extinction efficiency, single scatter albedo and asymmetry factor for seven aerosol types, which within the lookup table are treated as pure-composition particles of sea-salt, black carbon, dust, ammonium sulfate, ammonium nitrate, organic carbon, and pure water. 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 table is the wet particle radius range for a logarithmically expanding set of cut-sizes of the 240 aerosol distribution of GEM-MACH (0-10, 10-20, 20-40, 40-80, 80-160, 160-320, 320-640, 640-1280, 1280-2560, 2560 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 245 relative humidity and longer wavelength was due to the growing dependency of optical properties to increasing relative humidity and decreasing wavelengths. We assumed a flat distribution of aerosol radii within each bin size. The complex refractive indices of the water-soluble 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 250 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 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 Toon et al. (1976). The water activity for this aerosol type is based on Tang and Munkelwitz 255 (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 is based on Shettle and Fenn (1979), which relies on Dorsey (1940).

Aerosol
Density ( Table2: Density and the refractive index of each aerosol type. For dust, no water uptake has been assumed, and thus no dependency on the relative humidity. The complex refractive index 260 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 wavelengthdependent 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 265 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% 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 270 from Lei et al. (2014) (Table 3). The RH-dependent refractive index (nRH) is a HGF 3 -weighted mean of refractive indices of pure black carbon (nBC) and pure water (nwater): with water having a spectrally constant refractive index of 1.333 (with no imaginary component over the relevant wavelength range) and pure black carbon having a spectrally constant refractive index of 1.75 -0.44i (Kou, 1996). The optical properties 275 of the primary organics were calculated using the optical properties of secondary organics and pure water droplets and the hygroscopic growth factor of primary organics. The refractive index of each aerosol type can be found in Table2. In order to employ the new lookup table in GEM-MACH, we developed a stand-alone FORTRAN code to interpolate the optical properties of aerosols into GEM-MACH wavelengths and size bins. The output of this off-line program was used as the new lookup table for the improved version of the photolysis code discussed in this study. Here we provide a brief description of 280 the methods we used in our calculations.
Since the initial version of the new lookup table was calculated based on wet particle sizes, we used the hygroscopic growth factor of each aerosol type (Laskina et al., 2015, Ting et al., 2017, Zamora et al., 2013   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 α is the initial size parameter, r is the dry radius of the particle, and λ is the wavelength in the initial lookup table. 290 Using the same formula, the size parameter for the final lookup table was calculated based on GEM-MACH average dry 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 295 ammonium sulfate versus dry aerosol size parameter is illustrated in Fig. 2. The final interpolated lookup table of extinction 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 300 lookup table depend on the water uptake for the given relative humidity.  volume concentration of each aerosol type (m³ of each aerosol per m³ of air) was calculated using the GEM-MACH predicted mass concentration (µg kg -1 ), and the density (kg m -3 ) of each aerosol type (see Table2). 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 315 ammonium sulfate, ammonium bisulfate and letovicite have very similar refraction indices, the remaining mass of ammonium and the mass concentration of 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: 320 where r is the dry radius of the black carbon core of a particle and ρ 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 325 used to calculate the size parameter of the black carbon core (α ) and the total particle size parameter (α ). The absorption 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 . 330 Sea-salt, dust and black carbon (when its mass fraction was > 40%) aerosols were assumed to be externally mixed at all times -their effective scattering coefficient(sca ), absorption coefficient(abs ) and asymmetry factors (asym ) 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 is the extinction efficiency of each aerosol type, N is the number concentration of each aerosol type (cm ), r is the radius of bin i (cm), ssa is single scatter albedo of each aerosol type, ag is the asymmetry factor, and mf is the fraction of the mass concentration of 340 each aerosol type to the total mass concentration of all particles (the asymmetry factor was weighted by mass fraction for pure 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 is the radius of black carbon core particle and is the absorption amplification factor by , 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: 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: 355 where the indices (j = 1,2,3,4) correspond to each 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 360 11 to 13 were done for each horizontal grid-point, vertical level, wavelength and relative humidity.
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: 365 where τ and τ are the scattering and absorption optical depth of each aerosol category and dz is the vertical level thickness. The total scattering and absorption optical depth at each vertical level for each aerosol category (black carbon, seasalt, dust and internally mixed particles) was calculated using: 370 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. 375 This information was used to calculate J-values for photolysis reactions depending on the attenuation of the radiation stream 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 380 States with a horizontal grid-spacing of 10 km, 80 hybrid vertical levels with a model top at 0.1 hPa, 15-minute operator 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  (Lee et al., 2002). The model has been integrated into the ECCC FireWork air quality forecast system (Pavlovic et al., 2016) and has been operational since May 2019. The CFFEPS v2.03 code is available from: https://doi.org/10.5281/zenodo.2579383 . The version of CFFEPS used in this work is described elsewhere (Makar et al., 2021). 405 km at nadir were used in this study (obtained from https://www.avl.class.noaa.gov/saa/products/welcome). The MERRA-2 atmospheric reanalysis assimilation system produced by NASA's Global Modeling and Assimilation Office (GMAO) is the second and updated version of the original MERRA atmospheric reanalysis (Rienecker et al., 2011, Gelaro et al., 2017.
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 435 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 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: 440 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 445 compare with GEM-MACH output AOD calculated as described in section 2.3.

Simulation Plan
We  Table 4. The spatial and temporal resolutions are the same (as described in section 2.4) for all the simulations.
Note that in the "no-feedback" simulations, aerosol optical (and cloud condensation nucleation) properties come from default climatological properties used in the GEM weather forecast model (Makar et al, 2015a, b). That is, our "no-feedback" simulation is not a "no aerosol" simulation -rather, the "no-feedback" simulation makes use of spatially invariant, "typical" optical properties of our weather forecast model. Having "no AOD" in this case means that model aerosol AOD had not been 460 calculated since it is not used in the feedback code. Further, the aerosol feedback on meteorology code developed in a previous study for the GEM radiative transfer scheme (Makar et al, 2015a, b). In the feedback mode, GEM-MACH uses a separate Mie calculation for generating a lookup-table online, using the particle mass size distribution predicted internally by MACH at each time step, but with a single typical complex refractive index. In contrast, our photolysis rate calculations and the AODs we calculate in the current work, are decoupled from this feedback portion of the model. Our next step in this work is to further 465 modify the GEM radiation code to include the particle chemical composition on-line, and our hybrid particle mixing state assumptions as an alternative to the existing AOD calculation within the feedback portion of the code.

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 with the previous version on MESSy-JVAL (base) and simulation S1 with the improved photolysis module (see Table 4). The 475 simulated GEM-MACH AOD output was compared with AERONET sun photometer measurements of AOD at 580 nm for all North American sites in Fig. 4. As mentioned in section 2.5, we used (440 nm -675 nm) Angstrom Exponent to obtain the AERONET AOD at 580 nm. The statistical analysis of the output from the two GEM-MACH simulations (base and improved) and AERONET AOD at 580 nm for June 2018 at four Canadian sites, Egbert, Fort McKay, Saturna Island and Toronto, are shown in Fig. 5. As shown in these plots, the normalized mean bias ( = ∑ ( ) ∑ ( ) × 100), which represents the mean 480 paired differences between the model and measurements normalized by the mean measurements, is ranging within ±13% for the improved version and 0%-150% for the base version. The NMB calculations from the improved version show an overprediction of AOD in Saturna Island and an under-prediction of AOD in Toronto, Fort McKay and Egbert, whereas the base version shows a significant over-prediction of AOD for all four sites. The root mean square error (RMSE) is significantly smaller in the improved version, with less variability around the mean as shown in the standard deviation (σ) plots. The 485 correlation coefficient plots show better results with the improved version of MESSy-JVAL for all four sites. We calculated a correlation relation of R=0.17 for the base run and R=0.37 for the improved version for all North American AERONET sites.
The base MESSY module uses a climatology for aerosol number density with one fixed vertical profile for grid cells over-land and another fixed vertical profile over water (see Table 1). The uniformity of the fixed profiles does not account for the real atmospheric variability resulting in the larger differences between model and observed aerosol optical depths for the base 490 MESSY version.

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 the GEM model, on the resulting calculated photolysis rates (see Table 4).

3.2.1
Comparison with observations 500 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) for the entire simulation period (January and June 2018) and for all North American sites in Fig. 4. Figure 6 illustrates the time series of AERONET and GEM-MACH hourly AOD output at 580 nm for Fort McKay and Toronto 505 AERONET sites. The output from simulations S1 (with AOD calculations and online aerosol feedbacks, June 2018) and W1 with AOD calculations and online aerosol feedbacks, January 2018) were used to plot these time series.  As shown in this figure, 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 number of data points at this site compared to Toronto. As shown in these time series, the maximum value of the modelled AOD was underestimated compared to AERONET data for both sites and both simulation periods. For example, the 520 maximum AERONET AOD for the month of June in Fort McKay and Toronto were measured ⁓0.27 and ⁓0.42 respectively, whereas the GEM-MACH AOD (simulations S1 and W1) were ⁓0.19 and ⁓0.29 for those sites.

(d) (c)
on both GEM-MACH and MERRA-2 plots. As shown in Fig. 9(d), there is an AOD hotspot over the central Alberta in January, which was not captured by GEM-MACH ( Fig. 9(b)). The location of the hotspot in central Alberta in winter season in the MERRA-2 product is coincident with the location of coal-fired power plants and oil refinery emissions of SO2. This could be the result of different SO2 emissions used in MERRA-2 and GEM-MACH models, or an over-prediction error when sulfate aerosol is retrieved over a snow covered surface. Over the common domain, there is an underestimation of monthly average 555 AOD in the GEM-MACH model compared to MERRA-2 for both seasons.

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' 560 mode (Makar et al., 2015a, b) vs 'no-feedback' mode, in which the feedback mechanisms (interactions between meteorology 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 565 with online aerosol feedback between weather and air-quality in the model. We used the output from simulations S1 and S2 to calculate the percentage difference in photolysis rates in Fig. 10(a) and (b), and the output from simulations W1 and W2 in calculations. The output from simulations S3 and S4 were used in Fig. 10(e) and (f) , meaning there is no online aerosol feedback in these simulations. As shown 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 575 cloud patterns, and amplified chemistry perturbations due to weather feedback. The hotspots of greatest difference in Fig. 10 (f) reflect the direct effects on aerosol with no weather feedback in these simulations. Figure 11 is the monthly average percentage difference in AOD, number mixing ratio of clouds in the boundary layer at three heights, J(O 1 D) at the lowest model level, and J(NO2) at the lowest model level for June 2018 with and without interactive aerosol feedback on weather in the model (simulations S1 and S3). This figure shows the effects of the interactive online 580 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. The decrease of Jvalues 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. As we can see in Fig. 11(b), (c) and (d), the moist marine air 585 masses have the greatest impact on J-values change likely due to changes in low clouds. The areas impacted are mostly over the ocean or along ocean boundaries (northeast US/Canada, Northwest Territories, north Pacific). Over land, in convective air masses, there is less impact on J-values averaged over a region. These changes in humidity impact low clouds which, in turn, impact the radiative transfer and photolysis rates. 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. 590  Fig. 12 and  Fig. 12(a), (b)) . In the cases without the online 605 aerosol feedback in the model (simulations S3 and S4), the standard deviation for both stations is considerably smaller (0.55 for Toronto and 0.21 for Fort McKay) and the values of percentage difference in J(O 1 D) are more clustered around the mean ( Fig. 12(c), (d)). 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. 610       (Fig. 14(b) and 14(c)) show two areas of maxima; one is directly over the forest fire plume, similar to the hotspot on GEM-MACH plot (Fig. 14(c)), and a weaker hotspot over the north-east of the major forest fire plume, which is more 630 intensified compared to the GEM-MACH secondary hotspot. The aging of the major fire plume downwind is evident in all three plots. The maximum GEM-MACH AOD (0.625) is underestimated compared to the MAIAC (maximum of 3) and VIIRS (maximum of 1.7). One possible explanation for the underestimation in both primary and secondary hotspots, could be the potential deficiencies in modeled biomass burning emissions or aerosol processes in this simulation. Curci et al (2015) suggest that AOD under-predictions may be a common problem for current air-quality models. Other studies e.g., Pan et al. (2020) 635 andJohnson et al. (2016) show underestimations of modeled AOD over biomass burning areas, which could be related to the shortcomings in biomass burning processes in air-quality models.  Loche. As shown 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 650 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 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 655 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. 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 660 SOA increases downwind from the fire plume. Similarly, the concentration of nitrate (Fig. 15(f)) increases downwind between by the fresh fire plume is 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. 665 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 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), 670 reaching a difference close to 1% to 6% from the surface up to the mid-boundary layer. The OH radical concentration difference ( Fig. 16(b)) responds to the O3 change and decreases by 10% to 20% from the surface up to the boundary level, due to the AOD feedback on the photolysis rates. This affects the rate of 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. 675 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. The changes in black carbon concentration (±5%) in Fig. 16(d) at the surface up to the mid-boundary layer are mostly due to aerosol feedbacks on meteorology. Figure 17 is the percentage difference in aerosol concentrations and photolysis rates, with and without AOD calculations and 680 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 of OH (Fig. 17(b)) by the photolysis of O3 is more considerable (-5% to 5%). The decrease in the OH concentration with 685 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%.

(g) (h)
A new lookup table for aerosol optical properties based on a Mie scattering code was calculated and adopted within an 705 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  core 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 730 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 735 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.
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 740 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.