Research article 26 Feb 2021
Research article  26 Feb 2021
An uncertaintybased protocol for the setup and measurement of soot–black carbon emissions from gas flares using skyLOSA
 Energy and Emissions Research Laboratory, Department of Mechanical and Aerospace Engineering, Carleton University, Ottawa, ON K1S 5B6, Canada
 Energy and Emissions Research Laboratory, Department of Mechanical and Aerospace Engineering, Carleton University, Ottawa, ON K1S 5B6, Canada
Correspondence: Matthew R. Johnson (matthew.johnson@carleton.ca)
Hide author detailsCorrespondence: Matthew R. Johnson (matthew.johnson@carleton.ca)
Gas flaring is an important source of atmospheric soot–black carbon, especially in sensitive Arctic regions. However, emissions have traditionally been challenging to measure and remain poorly characterized, confounding international reporting requirements and adding uncertainty to climate models. The skyLOSA optical measurement technique has emerged as a powerful means to quantify flare black carbon emissions in the field, but broader adoption has been hampered by the complexity of its deployment, where decisions during setup in the field can have profound, nonlinear impacts on achievable measurement uncertainties. To address this challenge, this paper presents a prescriptive measurement protocol and associated opensource software tool that simplify acquisition of skyLOSA data in the field. Leveraging a comprehensive Monte Carlobased general uncertainty analysis (GUA) to predict measurement uncertainties over the entire breadth of possible measurement conditions, general heuristics are identified to guide a skyLOSA user toward optimal data collection. These are further extended in the opensource software utility, SetupSkyLOSA, which interprets the GUA results to provide detailed guidance for any specific combination of location, date–time, and flare, plume, and ambient conditions. Finally, a case study of a skyLOSA measurement at an oil and gas facility in Mexico is used to demonstrate the utility of the software tool, where potentially small regions of optimal instrument setup are easily and quickly identified. It is hoped that this work will help increase the accessibility of the skyLOSA technique and ultimately the availability of field measurement data for flare black carbon emissions.
Gas flaring is a routine practice in the oil and gas industry in which producers and refiners burn excess or unwanted gases in openatmosphere flames, typically from vertical pipe stacks. Flaring is generally preferable to the venting of gases to the atmosphere because it reduces carbon dioxide (CO_{2})equivalent emissions; however, flaring still emits potent climateforcing pollutants directly to the atmosphere (Allen and Torres, 2011; Johnson et al., 2001, 2011, 2013; McDaniel, 1983; Pohl et al., 1986). These pollutants have public health implications (e.g., Anenberg et al., 2012) and include unburnt hydrocarbons, volatile organic compounds, and particulate matter (U.S. EPA, 2018). Soot particulate matter (commonly referred to as black carbon, BC) has been suggested by some to be the second most potent climate forcer after CO_{2} (Bond et al., 2013; Jacobson, 2001; Ramanathan and Carmichael, 2008; Sato et al., 2003). Annual flaring is estimated by satellite imagery to be ∼140 billion m^{3} (Elvidge et al., 2007, 2009, 2016), making it one important source of global soot emissions. Although other industrial sectors dominate gas flaring in absolute soot emissions, the locations of flaring activities (particularly in Russia) likely have a disproportionate impact on the sensitive Arctic climate due to efficient transport pathways penetrating the Arctic air mass (e.g., Popovicheva et al., 2017; Stohl et al., 2007).
With the addition of BC to the Gothenburg protocol in 2012 (United Nations Economic Commission for Europe (UNECE), 2012), 34 countries are now legally bound to report, where data are available, BC emissions under UNECE's Convention on LongRange Transboundary Air Pollution, including the European Union, Russia, the United States of America, and Canada. To attribute – and, hence, report and regulate – soot–BC emissions from various sources, emission factors that relate soot–BC emissions to a measure of industrial activity are required. Unfortunately, for gas flaring, commonly employed soot emission factors are crude singlevalued parameters that link emitted soot mass to volume or mass of gas flared regardless of flare design, gas composition, or operating conditions. This contrasts with numerous studies that have observed a significant influence of flare gas composition and flame aerodynamics on soot emissions (Becker and Liang, 1982; Conrad and Johnson, 2017; McEwen and Johnson, 2012) and even soot properties (Conrad and Johnson, 2019; Trivanovic et al., 2020). Further soot yield data are needed, particularly for realworld flares under field conditions, to develop and validate accurate flare soot–BC emission factor models.
At present, there are only two published methods for the quantitative measurement of soot–BC emissions from individual infield flares. One technique employs aircraftbased sampling of a flare plume (Gvakharia et al., 2017; Weyant et al., 2016), where measurements of soot, methane, and CO_{2} concentrations during transects through the plume are used to provide flarespecific estimates of soot yield, using assumed flare gas compositions. The second technique is a groundbased remote optical measurement called skyLOSA (lineofsight attenuation using skylight; Conrad and Johnson, 2017; Johnson et al., 2010, 2011, 2013). SkyLOSA quantifies timeresolved soot mass emission rates through analysis of highspeed image data. Parallel access to flare infrastructure permits simultaneous measurement of flare gas flow rate and gas sample extraction for offsite compositional analysis, which enables the direct calculation of soot yield for a targeted flare. To date, skyLOSA has been deployed on 11 field measurement campaigns in Uzbekistan, Mexico, Ecuador, and Canada, providing 28 measurements of soot emissions from 17 unique flares (Conrad and Johnson, 2017; Johnson et al., 2011, 2013).
The key component of a skyLOSA measurement is the quantification of plume soot loading using image data, via analysis and modelling of radiative transfer through the atmospheric flare plume at the measurement wavelength. For each acquired image, soot mass column density is resolved pixel by pixel over a control surface within the image plane and combined with simultaneously measured velocity (obtained via image correlation velocimetry) to permit mass emission rate calculation. Uncertainties in skyLOSAcalculated emission rate are computed under a Monte Carlo (MC) framework and are dominated by uncertainties that affect computation of soot mass column density. While these uncertainties are influenced by numerous parameters considered within the MC analysis, they are also sensitive to the positioning and pointing of the skyLOSA camera relative to the horizon and sun. Consequently, a skyLOSA user must position the camera according to several constraints, which may be heuristic but can also vary with uncontrollable measurement parameters. To make the measurement technique accessible to endusers, enabling an increase in flare soot emissions data, a standardized data acquisition protocol for skyLOSA is required.
The objective of this work is to complete a general uncertainty analysis (GUA) for the skyLOSA measurement technique that provides uncertaintybased guidance to an enduser regarding the setup of equipment and acquisition of skyLOSA data through an accompanying opensource software tool. A summary of skyLOSA theory, referring to derivations in previous works, is first provided in Sect. 2 of this paper. The GUA methodology is summarized in Sect. 3, including special provisions necessary to reduce the computational burden of the MCbased approach (Sect. 3.1 and Appendix A). Representative results from the MC GUA are shown in Sect. 4.1, and general heuristics for the acquisition of skyLOSA data, including new observations based on MC GUA results, are summarized in Sect. 4.1.1 and 4.1.2. To provide casebycase guidance, a new opensource software tool to calculate skyLOSA measurement uncertainty is introduced in Sect. 4.2. Finally, in Sect. 4.3, the software tool is used in a case study that analyzes optimal camera positioning for flare measurements at a gas refining and transport facility in Campeche, Mexico. This work enables a consistent approach for the selection of skyLOSA camera positioning and pointing to minimize measurement uncertainties, ultimately contributing to the standardization of the skyLOSA measurement technique.
The generalized skyLOSA theory was summarized in full by Johnson et al. (2013) and has been the subject of a variety of validation efforts (Conrad et al., 2020a, b; Johnson et al., 2010). Development of the theory begins with Fig. 1, which shows an example skyLOSA image for computation of timeresolved soot emission rate from a sootladen flare plume in the Montney formation of Alberta, Canada. A highly linear, grayscale, scientificCMOS camera (e.g., pco edge 5.5) is used to obtain upwards of 10 min of highspeed image data of the flare and turbulent, sootladen, atmospheric plume. Pseudo 16bit images are acquired at frame rates of 25–50 Hz with a narrow midvisible bandpass filter (531±20 nm) to yield a scene of spectrally integrated light intensity.
Overlaid in Fig. 1 is an example control surface (C) of specified radius (r [m]), through which the instantaneous mass emission rate of soot (${\dot{m}}_{\mathrm{s}}$ [g s^{−1}]) may be computed. For an arbitrary control surface in three dimensions, the instantaneous mass flux of soot through the surface is
where ρ_{s} is the mass concentration of soot [g m^{−3}], u is the local plume velocity vector [m s^{−1}], n is the unit vector locally normal to the control surface [–], and dA is an infinitesimal area [m^{2}]. For skyLOSA, where threedimensional data along a pixel's line of sight (LOS) are collapsed to two dimensions through projection, an equivalent formulation for the instantaneous mass emission rate is
where ${\mathit{\rho}}_{\mathrm{s}}^{\prime}\left(r,\mathit{\varphi}\right)=\int {\mathit{\rho}}_{\mathrm{s}}\left(r,\mathit{\varphi},x\right)\mathrm{d}x$ is the soot mass column density along a LOS (where the x dimension is orthogonal to the image plane) [g m^{−2}] and ${u}_{n}\left(r,\mathit{\varphi}\right)=\left(\int {\mathit{\rho}}_{\mathrm{s}}\left(r,\mathit{\varphi},x\right)\mathit{u}\left(r,\mathit{\varphi},x\right)\cdot \mathit{n}\left(\mathit{\varphi}\right)\mathrm{d}x/\phantom{\int {\mathit{\rho}}_{\mathrm{s}}\left(r,\mathit{\varphi},x\right)\mathit{u}\left(r,\mathit{\varphi},x\right)\u2022\mathit{n}\left(\mathit{\varphi}\right)\mathrm{d}x{\mathit{\rho}}_{\mathrm{s}}^{\prime}\left(r,\mathit{\varphi}\right)}{\mathit{\rho}}_{\mathrm{s}}^{\prime}\left(r,\mathit{\varphi}\right)\right)$ represents the component of the massconcentrationweighted LOSaveraged velocity of the plume [m s^{−1}] that is normal to the control surface from the camera's perspective (as shown in the figure). Via Eq. (2), skyLOSA thus requires knowledge of three items to compute the emission rate: spatial scaling of the image plane to accurately quantify r, the velocity field of the plume within the image plane (yielding u_{n}(r,ϕ)), and the soot mass column density resolved over the control surface. Spatial scaling of the image is obtained through use of a pinhole analogy for the skyLOSA optics, coupled with a measurement of the distance to the flare stack tip by laser rangefinder (e.g., Laser Technology Inc. TruPulse 360R; Johnson et al., 2013). Given that imaging is performed with a global shutter and at a sufficiently rapid frame rate and exposure, the twodimensional plume velocity field over the image plane is estimated via image correlation velocimetry using a thirdparty software suite such as LaVision DaVis 8.4 that includes a means of uncertainty quantification (Wieneke, 2015). Finally, the novel enabling aspect of skyLOSA is the use of bounded knowledge of soot optical properties from literature data to compute the soot mass column density with accurate uncertainties via radiometric observations and modelling of radiative transfer within the atmosphere and plume. Example skyLOSA measurements of timeresolved soot emission rate are available in previous works (Conrad and Johnson, 2017; Johnson et al., 2011, 2013).
Figure 2 shows an example positioning and pointing of the skyLOSA camera and an optical axis (LOS) within the surrounding sky dome. For a given LOS, a Cartesian coordinate system is defined where the positive x direction is the path that light travels into the camera, the positive z direction is the general direction of plume motion, and the y dimension completes the orthogonal system. To model radiative transfer, there are three boundary conditions that must be considered. Firstly, the ground is treated as a cold, black surface and is thus ignored within the skyLOSA algorithm. Secondly, the sky is modelled as a diffuse, polarized source concomitant with atmospheric scattering of solar radiation. The distribution of skylight intensity ($I\left(\mathit{\alpha},{\mathit{\alpha}}_{\mathrm{s}},Z,{Z}_{\mathrm{s}},a\right)$ [W m^{−2} sr^{−1}]) and the incident intensity along the LOS (${I}^{\mathrm{o}}\left({\mathit{\alpha}}_{\mathrm{s}},\mathit{\beta},{Z}_{\mathrm{s}},a\right)$ [W m^{−2} sr^{−1}]) are considered using the standard models of the Commission Internationale de l'Eclairage (CIE, 2003), where the index a∈{1…15} indicates a specific CIE sky “type”. Finally, groundlevel normal solar irradiance (E_{sn} [W m^{−2}]) is estimated using infield, imagebased measurements of the sun (Johnson et al., 2013) or modelled using the CIE models. With this radiative transfer model, skyLOSA quantification of soot mass column density proceeds with the radiative transfer equation (RTE):
where I^{t} is the measured “transmitted” intensity at the camera [W m^{−2} sr^{−1}], ${\mathit{\sigma}}_{\mathrm{m}}^{\mathrm{e}}\left(\mathit{b}\right)$ is the massnormalized extinction cross section of the polydisperse soot population [m^{2} g^{−1}] that is a function of eight soot properties represented by the vector b, J^{s}(x,b) is a local source radiant intensity per unit mass [W sr^{−1} g^{−1}] along the measurement path, and the skyLOSA camera is located at x=L. Populationaveraged optical properties of soot are computed from the fundamental properties in b using Rayleigh–Debye–Gans theory for polydisperse fractal aggregates (RDGPFA; e.g., Sorensen, 2001). The vector b is composed of the absorption function of soot at the measurement wavelength (E(m_{λ}) [–], where m_{λ} is the spectral complex refractive index of soot [–]), the scatteringtoabsorption function of soot at the measurement wavelength ($F\left({m}_{\mathit{\lambda}}\right)/E\left({m}_{\mathit{\lambda}}\right)$ [–]), the monodisperse primary particle diameter of the soot population (d_{p} [nm]), the geometric mean soot aggregate size (N_{g} [–]), the geometric standard deviation of the lognormal soot aggregate size distribution (σ_{g} [–]), the diameterbased fractal prefactor (k_{f} [–]), the fractal (Hausdorff) dimension (D_{f} [–]), and the material density of soot (ϱ_{s} [g cm^{−3}]). Consistent with laboratory observations of soot structure in the overfire region of flarelike flames (e.g., Köylü and Faeth, 1992), these eight properties are assumed to be spatially and temporally uniform within a single MC draw; however, they are treated as random variables within the skyLOSA MC method (Johnson et al., 2013). The prior probability distributions employed within the MC method inherently link light absorption measurements and computed mass column density and emissions using skyLOSA; as such, in keeping with Andreae and Gelencsér (2006) and Petzold et al. (2013), skyLOSAinferred soot–BC mass might therefore be called “equivalent BC” as is recommended for all lightabsorptionbased diagnostics, especially where absorptionenhancing nonBC material may be present. For the case of flaregenerated soot–BC however, studies in the field (Schwarz et al., 2015; Weyant et al., 2016) have identified that the presence of nonBC aerosols in flare plumes is “not statistically different from zero” (Weyant et al., 2016), which is supported by laboratory observations (e.g., Kazemimanesh et al., 2019). This justifies use of soot property probability distributions derived from literature data of freshly emitted soot particulate by Johnson et al. (2013).
Using the mean value theorem, a pathaveraged source radiant intensity ($\stackrel{\mathrm{\u203e}}{{J}^{\mathrm{s}}}\left(\mathit{b}\right)$ [W sr^{−1} g^{−1}]) can be introduced to simplify the RTE:
where ${\mathit{\sigma}}_{\mathrm{m}}^{\mathrm{e}}\left(\mathit{b}\right)/{\mathit{\sigma}}_{\mathrm{m}}^{\mathrm{e}}\left(\mathit{b}\right)=\mathrm{1}$ has also been introduced to the source term. The resulting ratio in front of the second integral represents a pathaveraged source intensity ($\stackrel{\mathrm{\u203e}}{{I}^{\mathrm{s}}}$ [W m^{−2} sr^{−1}]), and the exponential in the first term corresponds to the transmittance of the plume in the absence of radiative sources, defined as the idealized transmittance τ^{*} (Johnson et al., 2013). It can be shown that the integral in the second term is equal to the complement of the idealized transmittance, permitting solution of the latter from the three intensities in Eq. (4):
where ${\mathit{\tau}}_{\mathrm{obs}}={I}^{\mathrm{t}}/{I}^{\mathrm{o}}$ is the observed transmittance of the plume [–] and $S=\stackrel{\mathrm{\u203e}}{{I}^{\mathrm{s}}}/{I}^{\mathrm{o}}$ is a term that corrects for brightening of the plume by radiative sources [–]. Noting the definition of the idealized transmittance, the column density of soot is simply
In general, the local source radiant intensity (J^{s}(x,b)) is composed of thermal emission and inscattering components. For skyLOSA however, where measurements are performed in the midvisible spectrum and plume temperatures are near ambient, thermal emission negligibly contributes to the source term of the RTE (Conrad et al., 2020b). By contrast, diffuse skylight and direct solar radiation can significantly augment the RTE via inscattering into the optical axis; hence, for skyLOSA, S represents an inscattering correction. Unfortunately, even if exact knowledge of skylight intensity distribution and solar radiation were available, it is not possible to fully account for the effect of inscattering without prior knowledge of the spatial distribution of soot within the plume. This is because multiplescattering events may occur during light's transmission into the camera. A recent simulation effort, however, has shown that complex multiplescattering effects can be accurately modelled to permit quantification of the inscattering correction (S) (Conrad et al., 2020b). This approach requires calculation of the inscattering correction using a singlescattering assumption (1SA, subscript “1”):
where ω(b) is the singlescattering albedo of the polydisperse soot population [–], p(θ;b) is the scattering phase function of soot [sr^{−1}], and angles θ and θ_{s} represent the angles between the LOS and a region of sky or the sun as shown in Fig. 2. The inscattering correction can be parsed into sky (S_{1,sky}) and sun (S_{1,sun}) components. Following Conrad et al. (2020b), the 1SAestimated inscattering correction permits calculation of the idealized transmittance – i.e., ${\mathit{\tau}}^{*}={\mathit{\tau}}^{*}\left({\mathit{\tau}}_{\mathrm{obs}},{S}_{\mathrm{1}}\left({\mathit{\alpha}}_{\mathrm{s}},\mathit{\beta},{Z}_{\mathrm{s}},a,\mathit{b}\right)\right)$. Ultimately, this allows for calculation of the soot mass column density by
According to Eq. (8), skyLOSA computation of soot mass column density is a function of the position of the camera and sun, fieldobserved plume transmittance, skylight intensity distribution and solar irradiance (through a), and soot properties. In the skyLOSA algorithm, soot mass column density is directly calculated by MC analysis over uncertain variables, which include soot properties and other intermediate parameters. This implies that uncertainty in skyLOSAcomputed soot mass column density is sensitive to α_{s}, β, Z_{s}, τ_{obs}, and a. Of these five variables, only the pointing of the camera may be controlled by the user, through the selection of angles α_{s} and β. Generally, overall uncertainty in skyLOSAcomputed soot mass column density is driven by uncertain soot properties through their effect on the massnormalized extinction cross section and the singlescattering estimate of the inscattering correction. Under optimal camera positioning, the latter is negligible relative to the former, and 95 % confidence intervals on skyLOSAcomputed soot mass emission rates are on the order of −26 %$/+\mathrm{36}$ % (Conrad and Johnson, 2017; Johnson et al., 2013). However, the magnitude and uncertainty of the inscattering correction are strongly sensitive to the userselected angles, α_{s} and β, and, in extreme cases, can even preclude accurate computation of the soot emission rate.
For endusers of skyLOSA, the sensitivity of measurement uncertainty to camera pointing necessitates a standardized (and ideally simple) data acquisition protocol to optimize camera position and pointing under general conditions. This would allow a priori setup decisions to minimize or constrain uncertainties within reasonable limits. An acquisition protocol must therefore be constructed using quantitative knowledge of measurement uncertainty in skyLOSAcomputed soot mass column density. Restated in the context of the above theory, the objective of this work is to quantify via a comprehensive general uncertainty analysis the uncertainty in skyLOSA computation of soot mass column density (${\mathit{\rho}}_{\mathrm{s}}^{\prime}$) as a function of userselectable (α_{s} and β) and uncontrollable parameters (Z_{s}, τ_{obs}, and a) under generalized conditions. These data permit the development of broad heuristics and, ultimately, an easytouse software tool to provide specific casebycase constraints on camera position and pointing.
The goal of the present general uncertainty analysis (GUA) is to guide a skyLOSA user in choosing a skyLOSA camera position and pointing to minimize measurement uncertainties. The developed software tool can also be used to give an initial estimate of uncertainties in the measured soot emission rate ahead of a more detailed postprocessing analysis. To provide generalized recommendations, the GUA quantifies measurement uncertainty in soot mass column density for a selected camera pointing and other independent variables via MC analysis over uncertain variables that include all relevant soot properties. This section describes the MC method used in the present GUA including novel updates to the MC approach that are necessary to make this present work tractable. This new methodology is a significant improvement to the skyLOSA algorithm that enables accelerated MC computation of soot column density and, hence, emission rates from skyLOSA image data.
3.1 Updated skyLOSA MC method
3.1.1 Expansion of the scattering phase function
For a given (modelled) skylight intensity distribution, measured or modelled solar irradiance, camera pointing, and set of soot properties, the 1SAestimated inscattering correction (S_{1}) can be directly calculated via Eq. (7). One significant challenge, and currently the timelimiting computation in skyLOSA processing, is the calculation of the skylight component (S_{1,sky}) via numerical integration over three dimensions: α, Z, and N, where the latter represents the aggregate size distribution of the soot population. The complexity of this task is exacerbated by the scattering phase function (SPF, p(θ,b)) of soot, which includes a computationally burdensome hypergeometric series in its solution (Sorensen, 2001). For the present GUA, where the inscattering correction must be computed over a fivedimensional domain (${\mathit{\alpha}}_{\mathrm{s}},\mathit{\beta},{Z}_{\mathrm{s}},a,\mathit{b})$, an alternative, more rapid means of computing the inscattering correction was required to avoid combinatorial explosion.
One such means is through a Fourier–Legendre expansion of the SPF. For an arbitrary set of randomized soot properties b, this procedure allows the SPF to be represented as a weighted sum of Legendre polynomials (P_{l}(cos θ) [–]):
where L(b) is the order at which the expansion is truncated [–] and Φ_{l}(b) is the lthorder Legendre coefficient [sr^{−1}] (hereinafter termed the lthorder soot coefficient) for the set of soot properties (b) computed via
which can be accurately and efficiently computed using Gauss–Legendre quadrature (Schuster, 2004). Since the soot coefficient decreases towards zero as the order l approaches infinity, the infinite series expansion of the SPF can be truncated at a sufficiently large index L(b) with negligible error (refer to Appendix A Sect. A1 for further details).
Introduction of the Fourier–Legendreexpanded SPF into the sky and sun components of the inscattering correction (S_{1,sky} and S_{1,sun}, Eq. 7) yields
where ${\mathrm{\Psi}}_{\text{l,sky}}\left({\mathit{\alpha}}_{\mathrm{s}},\mathit{\beta},{Z}_{\mathrm{s}},a\right)$ and ${\mathrm{\Psi}}_{\text{l,sun}}\left({\mathit{\alpha}}_{\mathrm{s}},\mathit{\beta},{Z}_{\mathrm{s}},a\right)$ are denoted as the sky and sun coefficients, respectively. Importantly, these equations show that use of the expanded SPF removes reference to soot properties from the integral in the computation of S_{1,sky}, which vastly reduces computational burden. Furthermore, with this formulation, the soot coefficients (functions of b) and sky and sun coefficients (functions of α_{s}, β, Z_{s}, and a) do not share any independent variables and can therefore be independently precomputed.
While the incident intensitynormalized solar horizontal irradiance (${E}_{\mathrm{sn}}\left({Z}_{\mathrm{s}},a\right)/{I}^{\mathrm{o}}\left({\mathit{\alpha}}_{\mathrm{s}},\mathit{\beta},{Z}_{\mathrm{s}},a\right)$ [sr]) is typically measured in the field by obtaining neutral densityfiltered images of the sun (Johnson et al., 2013), this parameter must be modelled for the GUA. This was accomplished using the CIE sky models and modeldependent typical turbidity factors and diffusetoextraterrestrial solar horizontal irradiance ratios as further detailed in Appendix A Sect. A2.
3.1.2 Sky model categorization
The standard CIE sky models have found good utility in a variety of fields, from urban planning (e.g., Acosta et al., 2014) to building design (e.g., Wong, 2017); however, the models naturally suffer from directionally dependent error in skylight intensity. This is particularly true for overcast and partly cloudy skies since the models, which are smooth functions, do not capture steep gradients in skylight intensity due to cloud structures. Thus, there is some additional uncertainty in skyLOSAcomputed soot mass column density through use of a single CIE sky model in the MC method. To permit capture of CIE sky model error in the GUA, like skies were sorted into sky “categories” that have similar properties but differing model coefficients and, hence, directional variability. The derived sky categories (a∈{A…D}) are summarized in Table 1 and can be selected by a skyLOSA user based on simple observations in the field such as the visibility of the sun (obstructed vs. unobstructed) and presence of clouds (overcast, partly cloudy, or clear). By randomly selecting a sky category's component sky models under the MC framework, uncertainty through use of the CIE models is propagated into skyLOSA computation of soot mass column density.
Sky category A corresponds to overcast and partly cloudy conditions with an obscured sun. Typical turbidity factors of the component skies (T(a)) are high in these conditions, corresponding to low groundlevel solar irradiance. Sky category B represents partly cloudy conditions with an unobscured sun where turbidity factors of the component skies are moderate and hemispherical skylight is strong relative to extraterrestrial solar radiation. Sky categories C and D capture the lowturbidity clear CIE sky models. Sky category C includes all five clearsky models, while sky category D excludes CIE sky type 12, the lowestturbidity (i.e., cleanest atmosphere) model. Sky category D was defined based on the notion that oil and gas activities can be relatively dense geographically. In this case, field experience suggests that local emissions from industrial infrastructure likely preclude CIE sky type 12 as a reasonable model, such that sky category D can be used for the case of clear skies in heavily industrial locales.
3.2 MC implementation
Table 2 summarizes the independent, precomputed, and random variables required to compute soot mass column density under the GUA MC framework. There are five independent variables that define the pointing of the camera relative to the sun (α_{s}, β, Z_{s}), the observed plume transmittance (τ_{obs}), and the skylight intensity distribution (CIE sky type, a). Each MC draw randomly chooses soot properties (b) and the scalar multiplier to the sun component of inscattering (ξ_{ED}(a), described in Appendix A Sect. A2). For analyses using the defined sky categories A–D, one CIE sky model is randomly obtained (i.e., a_{k}) from the selected sky category. The kth MC estimate of the soot mass column density is then calculated via
where ${\mathit{\tau}}^{*}\left({\mathit{\tau}}_{\mathrm{obs}},{S}_{\mathrm{1},k}\left({\mathit{\alpha}}_{\mathrm{s}},\mathit{\beta},{Z}_{\mathrm{s}},{a}_{k},{\mathit{b}}_{k}\right)\right)$ is deterministically computed while considering multiplescattering effects, as described by Conrad et al. (2020b), and
In a standard MC analysis, the above procedure would be iterated upon K times to yield a collection of soot mass column density estimates from which a posterior distribution of soot mass column density could be computed. To accelerate MC procedures in this work, a MC variance reduction technique was employed – specifically, combined multiple Latin hypercube sampling summarized by Nakayama (2011). This variance reduction technique has been used previously for skyLOSA (Conrad and Johnson, 2017) and was found to reduce computational burden by a factor of 2–3. For the GUA, 5×10^{5} (500 sets of 1000 Latinhypercubesampled data) MC draws were completed. The GUA MC approach permitted precomputation of the soot coefficients, singlescattering albedo, and massnormalized extinction cross section for predrawn random sets of soot properties (b). Parallel precomputation of the sky and sun coefficients was performed for each of the 15 CIE standard skylight intensity distributions and four sky categories over the angles α_{s}, β, and Z_{s} in increments of 2^{∘} and for 18 observed transmittances from 0.25–0.99. This amounted to execution of the MC analysis for almost 66×10^{6} unique skyLOSA conditions, permitting derivation of uncertainty statistics over the five independent variables listed in Table 2.
^{a} The means by which variables are randomly drawn or computed.
^{b} The method by which a variable is considered in the GUA MC method. Given a set of “independent” and “random” variables, precomputed data (${\mathrm{\Psi}}_{\text{l,sky}}\left({\mathit{\alpha}}_{\mathrm{s}},\mathit{\beta},{Z}_{\mathrm{s}},{a}_{k}\right)$, ${\mathrm{\Psi}}_{\text{l,sun}}\left({\mathit{\alpha}}_{\mathrm{s}},\mathit{\beta},{Z}_{\mathrm{s}},{a}_{k}\right)$, ${\mathit{\sigma}}_{\mathrm{m}}^{\mathrm{e}}\left(\mathit{b}\right)$, and Φ_{l}(b)) are obtained and used to calculate soot mass column density.
^{c} Sky model coefficients are independent variables in the analysis of a single skylight intensity model (a∈{1…15}), but, for the derived sky categories (a∈{A…D}), skylight intensity models are randomized following a discrete uniform distribution with support over the skylight intensity models included in the sky category.
^{d} Typical values as listed for each CIE sky model as per the literature (Darula and Kittler, 2002; Kittler et al., 2012).
^{e} Authorselected distributions: 𝒰(0.00,1.25) for CIE sky types 1–6 and 𝒰(0.75,1.25) for CIE sky types 7–15.
^{f} Prior probability distributions for soot properties were derived by Johnson et al. (2013).
To enable an objective comparison of skyLOSA uncertainty as a function of the independent variables, a parameter describing the relative uncertainty of MCcomputed soot mass column density was required. A natural means of representing relative uncertainty is a coefficient of variation (CV)like metric that describes a measure of data variance normalized by a measure of central tendency. For consistency with skyLOSA measurements, a CV estimator based on the mean and 95 % confidence interval was selected for this work. For variability, the width of the 95 % CI was scaled by that of the standard normal distribution (≈3.92) to yield an equivalent standard deviation, while the mean ($\stackrel{\mathrm{\u203e}}{\mathit{x}}$) was employed as the measure of central tendency. The CV estimator for soot mass column density (or any MCcomputed data x) was therefore
where subscript “95” signifies use of the 95 % CI for variability, and ${\mathcal{F}}_{\mathit{x}}^{\mathrm{1}}\left(q\right)$ is the qth quantile of data vector x.
4.1 Representative results
Figure 3 shows relative uncertainty results at different camera pointings for an example skyLOSA measurement scenario of a flare with 90 % observed plume transmittance. The selected solar zenith angle (Z_{s}=32.8^{∘}) represents the annual minimum for the Canadian city of Fort St. John, British Columbia, which is located in the Montney oil and gasproducing formation. CV_{95} data are plotted for sky categories A–D in Figs. 3a–d, respectively, as a function of relative solar azimuth (α_{s} in Fig. 2) and camera inclination (β); the position of the sun in α_{s}–β coordinates (180, 57.2^{∘}) is overlaid in each subfigure. Additionally, contours displaying the solar scattering angle (θ_{s}) are overlaid in Fig. 3d.
There are two observable trends in the data of Fig. 3 that persist through all measurement conditions. Firstly, the relative uncertainty is a strong function of the solar scattering angle, θ_{s}. This is because of solar radiation's influence on S_{1,sky} and S_{1,sun}. However, the rate at which relative uncertainty changes as a function of θ_{s} is dependent on the sky model (as well as solar zenith and plume transmittance). This is partly a consequence of how S_{1,sky} varies with θ_{s} but is mostly due to variability in S_{1,sun} with θ_{s}. For example, the inscattering magnitude for sky category A is effectively due to S_{1,sky} alone due to high turbidity that strongly attenuates direct solar radiation, while, for sky category C, the sun's inscattering contribution (S_{1,sun}) is significant due to the low turbidity of the sky. Thus, the results for sky category A largely represent the effect of θ_{s} on uncertainty through S_{1,sky} while the results for sky category C include an additional (and the most extreme) effect through S_{1,sun}. These observations imply that a constraint on θ_{s} is necessary and that a stricter constraint is required for lowerturbidity skies.
The mechanism for the decrease in relative uncertainty with increasing θ_{s} stems from the optical characteristics of soot particulate. Figure 4 shows statistics of the “energy distribution function” ($\mathrm{EDF}=\mathit{\omega}\left(\mathit{b}\right)p\left(\mathit{\theta},\mathit{b}\right)/\mathrm{4}\mathit{\pi}$) discussed by Conrad et al. (2020b), which describes how soot particulate directionally scatters light on an energy basis. The median and CV_{95} of the EDF as a function of angle θ are shown for the range of the MCsampled soot properties (b) used in the GUA. This is useful for the present discussion since the EDF contains all sootdependent variables in the computation of S_{1,sun} (see Eq. 7). The figure shows that the EDF, and its influence on S_{1,sun}, is much larger and much more uncertain as θ_{s} decreases. The influence of these statistics on column density uncertainty depends on the relative intensity of sunlight – specifically ${E}_{\mathrm{sn}}\left({Z}_{\mathrm{s}},a\right)/{I}^{\mathrm{o}}\left({\mathit{\alpha}}_{\mathrm{s}},\mathit{\beta},{Z}_{\mathrm{s}},a\right)$ in Eq. (7). If this value is small (highly turbid skies), then variability in the EDF is less important, but if this value is large (low turbidity skies), then variability in the EDF can dominate relative uncertainty in soot mass column density. Interestingly, relative uncertainty in the EDF approaches a constant value towards 90^{∘}. This implies that, regardless of the groundlevel intensity of the sun, the uncertainty of S_{1,sun} becomes minimal at ${\mathit{\theta}}_{\mathrm{s}}\stackrel{\mathrm{\u0303}}{>}\mathrm{90}$^{∘} and is within 1 % of this minimum for ${\mathit{\theta}}_{\mathrm{s}}\stackrel{\mathrm{\u0303}}{>}\mathrm{60}$^{∘} as indicated by the red line in Fig. 4.
The second trend in Fig. 3 that is generally seen across all measurement conditions is the sensitivity of relative uncertainty to the camera inclination angle (β). Referring to Fig. 3, much of the observed trend in β can likely be attributed to the effect of θ_{s}, since ${\mathit{\theta}}_{\mathrm{s}}={\mathit{\theta}}_{\mathrm{s}}({\mathit{\alpha}}_{\mathrm{s}},\mathit{\beta},{Z}_{\mathrm{s}})$. However, β still influences measurement uncertainty as can be observed in the region where the θ_{s} effect is small (${\mathit{\theta}}_{\mathrm{s}}\stackrel{\mathrm{\u0303}}{>}\mathrm{60}$^{∘}). Figure 5a through e show example trends in soot mass column density uncertainty as a function of β for plume transmittances of 0.25 to 0.95 and a fixed solar azimuth of α_{s}=60^{∘}. These panels each plot CV_{95} as a function of β, averaged over the range of Z_{s} that ensures θ_{s}>60^{∘}. The trends are different for each sky model. Uncertainties for sky categories A and B tend to decrease as the camera inclines, while uncertainties for sky categories C and D increase. The severity of these trends increases with plume transmittance (effectively a reduction in the measured signal), and, as plume transmittance increases towards unity, local minima–maxima in uncertainties as a function of camera inclination may appear. The differing influence of β between overcast/partly cloudy and clearsky models is largely due to the specific CIE sky models, which dictate the gradient in skylight intensity near the horizon and the sun. For clear, lowturbidity skies, intensity gradients are large such that small changes in camera pointing can yield significant changes in inscattered light. By contrast, for higherturbidity skies, gradients are dampened, and the effect of camera pointing is small.
Figure 5 also shows that the CV_{95} of soot mass column density for sky category C tends to upper bound that of sky category D. This observation holds under most combinations of the MCindependent variables and is a result of the somewhat extreme nature of CIE sky model 12 (lowest turbidity/cleanest atmosphere), which is excluded in sky category D. When considered in sky category C however, CIE sky model 12 tends to increase both the variability and central tendency of soot mass column density, but the relative change in the former is larger – hence, inclusion of CIE sky model 12 typically increases relative uncertainty. This implies that sky category C reliably imposes the largest constraints on camera positioning and can therefore be used to conservatively locate skyLOSA equipment under clear skies. However, if in a dense industrial area, where the clearestsky model is not relevant, the lessconstraining sky category D can instead be used as noted in Sect. 3.1.2.
4.1.1 Camera pointing heuristics
Upon arrival at a measurement facility, the skyLOSA user's first task is to determine the position of the skyLOSA camera for data acquisition. This important decision can be made by considering viable camera pointings from GUA MC data through constraints on β and θ_{s}. That is, for a useridentified sky category and plume transmittance and given the position of the sun, the skyLOSA camera's pointing can be constrained based on a desired threshold in a relative uncertainty, CV_{95}. Table 3 provides an example of such constraints. The table lists bounds on the camera inclination angle (β) and solar scattering angle (θ_{s}) for each sky category, given the plume transmittance: where bounds were computed by determining where CV_{95}≤17 % for all values of ${Z}_{\mathrm{s}}\in \left[\mathrm{0},\mathrm{90}{}^{\circ}\right]$, such that the results are independent of solar zenith.
One additional consideration in the pointing of the skyLOSA camera is the direction of plume propagation. Under quiescent conditions, buoyancydriven flare plumes will propagate vertically away from the flare stack; however, under sufficiently strong crosswinds, the flame and plume can bend over and propagate horizontally, parallel to the wind direction. In this latter case, if the plume propagates towards or away from the skyLOSA camera, turbulent plume structures of differing vorticity become overlapped from the camera's perspective. Therefore, it is best to position the skyLOSA camera such that it points orthogonally to the wind direction, which minimizes outofimage plane motion of the plume and yields the best data for velocimetry calculations. This should be viewed as a weak constraint on skyLOSA data acquisition however, since the effect of uncertainty in estimated velocity on mass emission rate is generally negligible compared to that of column density uncertainty.
4.1.2 Further camera heuristics
Following selection of a permissible skyLOSA camera position, the imaging optics must be chosen. Prime (fixed focal length) lenses are employed in the skyLOSA technique to avoid ambiguity in optical magnification and, hence, spatial scaling of the image. The most appropriate prime lens for the studied flare is one that maintains the entirety of the flare flame well within the image during the data acquisition period. This helps to ensure that a control surface within the image plane that transects the plume and encloses the flame can be derived, as shown in Fig. 1. For a flame that is relatively unsteady – i.e., moving with the wind – it is suggested to keep the flare flame approximately onequarter of the smallest image dimension. By contrast, if the flame is steady, a flame length of approximately onethird of the smallest image dimension should be targeted. For the skyLOSA camera used by the authors (minimum sensor dimension of ∼14 mm), a good rule of thumb is that the appropriate focal length (f, [mm]) will be on the order of
where H is the horizontal standoff distance from the flare stack [m], and L_{f} is the length of the flare flame [m]. However, use of a prime lens necessitates a tradeoff between lens focal length, horizontal standoff distance, and size of the flame within the image.
With an appropriate lens selected, the user must then choose imaging parameters that influence the exposure and focus of the image. The objective is to obtain an image that maximizes the digital signal while minimizing exposure time and ensuring the flame is in focus. In the authors' experience, this can be obtained with a lens aperture close to fullopen (typically fnumber ≤5.6) and an exposure time less than ∼2 ms. Prior to acquiring the image data, the flame and flare stack should be brought into focus. The user can then obtain skyLOSA data for the desired duration; it is recommended, however, that a minimum dataset of 10 min be obtained to permit good convergence of the timeaveraged soot emission rate.
4.2 Opensource software tool for simpler skyLOSA setup – SetupSkyLOSA
While the camera pointing heuristics presented in Table 3 can be used to ensure that CV_{95}≤17 % for the listed plume transmittances regardless of solar zenith angle, this simplified set of constraints is also necessarily overly conservative and excludes specific combinations of inputs that might produce similar or better uncertainties in different scenarios. An even better approach would be to use the wealth of computed GUA data to provide camera position and pointing constraints on a casebycase basis. This is made possible using a new opensource software tool, SetupSkyLOSA (Conrad, 2020), that was developed as part of this work using the presented GUA MC data. This MATLABbased application enables a skyLOSA user to probe statistics of soot mass column density (such as CV_{95}) for their specific measurement conditions. The key output of the software tool is an image of the desired soot mass column density statistic plotted as a function of absolute camera pointing. The software tool is briefly described in this section and employed in a case study in Sect. 4.3.
Figure 6 shows a flow chart describing the SetupSkyLOSA software's main procedure. For a userinputted location and time, the software first determines the current position of the sun using an integrated solar position calculator – a MATLAB implementation of the National Renewable Energy Laboratory's (NREL's) Solar Position Algorithm (SPA) (Reda and Andreas, 2008). The SPA returns the solar zenith (Z_{s}) and absolute bearing of the sun (α_{sN}, where the subscript “N” implies the absolute bearing measured clockwise from true north) at the current time and over the measurement date. The user also inputs the index (a) of the most appropriate CIE sky model or category, an estimate of the observed plume transmittance at the skyLOSA measurement wavelength (τ_{obs}), and the desired statistic of soot mass column density (η). With these inputs, SetupSkyLOSA then loads the GUA MC data for the selected sky model or category and interpolates for the desired statistic using the current solar zenith and estimated plume transmittance.
At this point, the software has computed η(α_{s},β) for the user's current set of independent variables $({Z}_{\mathrm{s}},{\mathit{\tau}}_{\mathrm{obs}},a$). However, rather than plotting η(α_{s},β) – i.e., using the relative bearing – the software uses the known absolute bearing of the sun (α_{sN}, computed by the SPA) to plot η(α_{N},β). That is, the requested statistic is plotted as a function of absolute camera bearing and inclination, which together define the camera pointing. The user can then easily determine an acceptable camera pointing using a laser rangefinder (for inclination) and a compass (corrected to true north), laser rangefinder, or GPS device (for absolute bearing).
SetupSkyLOSA also includes several added utilities to support optimal positioning and pointing of the skyLOSA camera. Firstly, using the same pinhole camera model that enables spatial scaling of the image, the software tool can optionally overlay the approximate extent of the image sensor in the α_{N}−β domain, based on sensor dimensions, employed optics, and the pointing of the centre of the image. This helps a user ensure that the entirety of the image frame – including the eventual control surface used for emission rate calculation – has reasonable levels of measurement uncertainty, given a userselected lens of known focal length. To support ideal velocity calculation, the software can also overlay camera pointings (α_{N}) that are closely orthogonal to the wind. This follows the heuristic discussed in Sect. 4.1.1. The software shows the optimal range of camera pointing as orthogonal to the wind ±18.2^{∘}, which corresponds to 5 % outofimage plane motion (cos ^{−1}(0.95)). Two additional utilities are not shown in Fig. 6. The “maximizer” utility computes the maximum of a chosen relative uncertainty statistic over a userdefined period. This tool allows a user to seek camera pointings that yield satisfactory uncertainties as the sun moves during the anticipated duration of the skyLOSA measurement. The “positioner” utility takes the plotted relative uncertainty statistic as a function of camera pointing and provides region(s) where the skyLOSA camera may be positioned relative to the flare stack given the flare stack height, maximum horizontal standoff distance from the flare, and relative uncertainty threshold. The user can optionally print these permissible regions to a .kml file for use in mapping software like Google Earth. These latter two utilities are employed in the following case study.
4.3 Case study – Atasta facility
The utility of the novel software tool, SetupSkyLOSA, is shown in this section via a case study of a skyLOSA measurement at a real oil and gas facility. The Atasta Gas Processing and Transport Centre (Centro de Proceso y Transporte de Gas Atasta) is a midstream oil and gas facility near Atasta, in the Mexican state of Campeche. The facility is under the jurisdiction of Petróleos Mexicanos and receives sour gas and condensates from the Cantarell offshore oil field for processing and transport to the national market. As shown in Fig. 7a, the Atasta facility is located 35 km west of Ciudad (Cd) del Carmen and approximately 5 km south of the shore of the Bay of Campeche. The facility occupies approximately 1 km^{2}, with most infrastructure in the southeast corner of the site as visible in Fig. 7b. Flaring activities include multiple pitstyle and vertical stack flares, which are in the northwest corner of the site.
For this case study, a skyLOSA measurement of soot emissions from the central flare stack at the Atasta station was considered as indicated in Fig. 7b, the base of which is located at 18^{∘}38^{′}41.46^{′′} N and 92^{∘}10^{′}08.59^{′′} W. The following example measurement details are assumed.

The skyLOSA measurements occur on 13 May 2021, which is the date of that year that the sun most closely reaches the solar zenith (Z_{s}=0^{∘}).

Predicted sky conditions are uncertain and may change between overcast and fully clear conditions throughout the day.

Wind speed is predicted to be low and the flare is strongly buoyant.

The flare stack is 30 m in height, and the horizontal standoff distance of the skyLOSA camera is limited to 250 m or less due to available optics.

The flare is lightly sooting, with an observed transmittance of approximately 90 % (τ_{obs}≈0.90).
Assumption 2 implies that skyLOSA data acquisition may occur under skies represented by any of sky categories A–D. Furthermore, assumption 3 suggests that the sootladen flare plume propagates vertically from the flare stack, and, therefore, the constraint on camera position with respect to wind direction is unimportant. The skyLOSA user wishes to obtain skyLOSA data with minimal measurement uncertainty, while also avoiding relocation of the skyLOSA camera throughout the day, if possible.
Given the known GPS coordinates of the flare stack, measurement date, and approximate plume transmittance, SetupSkyLOSA can be used to constrain skyLOSA camera pointing for any sky condition and time of the day based on the CV_{95} of soot mass column density. Since the user wishes to avoid relocation of the skyLOSA camera, camera position and pointing should be constrained using the relative uncertainty maximized over the day. The “maximizer” utility of SetupSkyLOSA permits this calculation for each of the sky categories; results are shown in Fig. 8a–d for sky categories A–D. Noting the differing colour scales in the four figures, there is significant variability in skyLOSA uncertainties for each of the sky categories, as in Fig. 3. For further context, the path of the sun over the measurement date is overlaid in Fig. 8a–d in addition to a contour of CV_{95}=16.5 %, which is within 1 % of the best attainable uncertainty for these conditions.
Using the uncertainty data in Fig. 8a–d, the positioner utility of SetupSkyLOSA can be employed to highlight where the skyLOSA camera may be positioned relative to the flare stack. This was performed for each of the sky categories using the uncertainty threshold of 16.5 %. Permissible camera positions were output in .kml format by the positioner utility and are overlaid on a map of the Atasta facility in Fig. 8e and are quite different for each of the sky categories. Permissible camera positions for sky category B exist beyond a small region near the stack tip, while those for sky category A are within an annular region surrounding the flare stack – since the lower limit on the camera inclination angle in Fig. 8a imposes a maximum permissible standoff distance. Sky category D contains two permissible regions – one to the south and one to the north of the flare stack – while the mostconstrained sky category C has one relatively small region to the north of the flare stack. Recalling assumption 2 that predicted sky conditions were uncertain, the skyLOSA user should ideally position the camera at the intersection of the skycategorydependent permissible regions. This small area is outlined in black in the figure, ∼136 m due north of the flare stack and just 604 m^{2} in size (∼0.31 % of the 250 m radius region). It is apparent in Fig. 8e that this ideal position intersects a clearing in the treed area where the skyLOSA camera should be positioned for the specific conditions of this case study.
This case study shows the remarkable utility of the SetupSkyLOSA software tool. The tool quickly provides resolved measurement uncertainty data from the GUA that would otherwise require millions of MC analyses to compute. These uncertainty data enable optimal skyLOSA camera positioning and pointing and also represent a firstorder estimate of soot emission rate uncertainties that are computed in postprocessing. Together with the additional utilities and general camera heuristics, this software tool permits a skyLOSA user to obtain optimal skyLOSA data that minimize measurement uncertainties under generalized conditions.
A comprehensive Monte Carlobased general uncertainty analysis (GUA) has been used to develop heuristics constraining the pointing and positioning of skyLOSA equipment for measurement of soot–black carbon emissions from gas flares. The GUA identifies generalized constraints based on predicted measurement uncertainties in soot mass column density, computed using skyLOSA. The results show that equipment setup constraints can be classified based on the conditions of the sky, relative positioning of the sun, and inclination angle of the camera. With additional heuristics on camera optics and imaging parameters, the presented results provide generalized guidance to greatly simplify acquisition of optimal skyLOSA data in the context of complex, nonlinear measurement uncertainties. These are further extended in the opensource software utility, SetupSkyLOSA, which interprets the GUA results to provide detailed guidance for any specific combination of location, date–time, and flare, plume, and ambient conditions. Furthermore, softwaredisplayed soot mass column density statistics provide the user with a firstorder estimate of measurement uncertainty in soot–black carbon emission rate that otherwise is only computable during postprocessing. The case study using SetupSkyLOSA to identify optimal equipment setup at a real oil and gas facility in Mexico demonstrates the utility of this new software tool, which as an opensource application can hopefully facilitate broader use of the skyLOSA technique and ultimately help increase the knowledge base of soot–black carbon emissions from gas flares.
A1 Truncation of the expanded SPF
For the prior probability distributions of soot properties derived by Johnson et al. (2013), the total order of soot coefficients required to represent the soot SPF according to the procedure of Schuster (2004) was typically L(b)=76 (median). In the most extreme case however, representing a strongly forward scattering soot population (corresponding to large soot aggregate size), L(b) reached 698, suggesting that precomputation of the sky and sun coefficients up to Ψ_{698} would be necessary to compute the inscattering correction in the worst case. While calculation of the sun coefficients to this large order is generally trivial, calculation of the sky coefficients becomes overwhelmingly cumbersome as the order increases. Importantly though, like the soot coefficients, the magnitude of the sky coefficients approaches zero as the order approaches infinity; therefore, the product of Φ_{l}(b) and ${\mathrm{\Psi}}_{\text{l,sky}}\left({\mathit{\alpha}}_{\mathrm{s}},\mathit{\beta},{Z}_{\mathrm{s}},\mathit{a}\right)$ more rapidly decreases in magnitude as l increases than either component alone. This permits further truncation of the series for the calculation of S_{1,sky}. Specifically, calculation of S_{1,sky} via Eq. (11) using L(b)=200 was consistently in close agreement with direct numerical integration via Eq. (7) – where, over 10^{6} randomized sets of $\left({\mathit{\alpha}}_{\mathrm{s}},\mathit{\beta},{Z}_{\mathrm{s}},a,\mathit{b}\right)$, the median relative difference was just $\mathrm{2.3}\times {\mathrm{10}}^{\mathrm{7}}$. This implies that it is acceptable to impose that ${\mathrm{\Psi}}_{\text{l,sky}}\left({\mathit{\alpha}}_{\mathrm{s}},\mathit{\beta},{Z}_{\mathrm{s}},a\right)=\mathrm{0},\forall l>\mathrm{200}$.
A2 Incident intensitynormalized solar normal irradiance
As noted in Sect. 2, groundlevel solar normal irradiance – E_{sn}(Z_{s},a), usually measured via solar images taken in the field – can be modelled using the CIE skylight models. To accomplish this, the incident intensitynormalized solar normal irradiance is expanded:
where D_{h}(Z_{s},a) is the diffuse horizontal irradiance, calculable for the CIE models via numerical integration of
which is independent of the value of α_{s}.
The ratio of solar normal to diffuse horizontal irradiance is complex to quantify in a general sense as it is a function of atmospheric composition. However, for the purposes of the present GUA it is modelled as follows:
The numerator of the righthand side is the groundlevel solar horizontal irradiance, calculated as the product of the extraterrestrial (subscript “o”) solar horizontal irradiance (${E}_{\mathrm{sh},o}\left({Z}_{\mathrm{s}},a\right))$ with an exponential representing attenuation through the atmosphere. In computation of the latter, m(Z_{s}) is the relative air mass quantifying the amount of air in the atmosphere at the solar zenith angle relative to the vertical direction, ${{\mathit{\sigma}}^{\mathrm{e}}}^{*}\left(m\right)$ is the ideal extinction for a clean atmosphere at a given relative air mass, and T(a) is the modeldependent turbidity factor that is used to consider realistic atmospheres and describes the number of clean atmospheres required to represent attenuation through the nonideal, polluted atmosphere. In the denominator, the cosine is included to transform the numerator into the required groundlevel solar normal irradiance. In the present work, typical values of the turbidity factor (T(a)) and irradiance ratio ${D}_{\mathrm{h}}\left({Z}_{\mathrm{s}},a\right)/{E}_{\mathrm{sh},o}\left({Z}_{\mathrm{s}},a\right)$ are taken from Darula and Kittler (2002) and Kittler et al. (2012), while their recommended formulations for m(Z_{s}) (Kasten and Young, 1989) and ${{\mathit{\sigma}}^{\mathrm{e}}}^{*}\left(m\right)$ (Navvab et al., 1984) are employed:
To model some amount of unknown uncertainty due to the use of “typical” metrics listed in the literature, an additional randomized variable (ξ_{ED}(a)) was introduced as a scalar multiplier to the sun component of the inscattering correction. For CIE sky models with an unobstructed sun (types 7–15), ${\mathit{\xi}}_{\mathrm{ED}}\left(a\right)\sim \mathcal{U}\left(\mathrm{0.75},\mathrm{1.25}\right)$, and for models with an obstructed sun (types 1–6), ${\mathit{\xi}}_{\mathrm{ED}}\left(a\right)\sim \mathcal{U}\left(\mathrm{0},\mathrm{1.25}\right)$. These prior distributions of ξ_{ED}(a) were based on observations by Watanabe et al. (2016), who studied the “clearness index” (groundlevel horizontal normalized by extraterrestrial horizontal irradiance) over 5 years at 47 observation stations across Japan. They found that the relative variation in the clearness index was approximately 4.3 % for skies with unobscured suns and 35 % for skies with obscured suns; corresponding varianceequivalent uniform distributions would have a range of 15 % and 121 %, respectively. For skies with an unobscured sun, this range was expanded to 50 % (0.75–1.25) to give a conservatively broad prior since measurement uncertainty can be quite sensitive to solar irradiance. By contrast, for skies with an obscured sun, where the solar irradiance has a small contribution to measurement uncertainty, the uniform distribution was only slightly widened to 125 % (0–1.25).
The developed software tool and associated data are available online as opensource and build distributions (https://doi.org/10.5281/zenodo.3908540; Conrad, 2020).
Both authors conceptualized the research and developed the methodology. MRJ was responsible for funding acquisition, project administration, provision of resources, and supervision. BMC curated the data, performed the formal analysis and investigation, developed the software, and produced the original manuscript including visualizations. Both authors reviewed and edited the manuscript throughout the publication process.
The authors declare that they have no conflict of interest.
We are grateful for the support of Michael Layer (project manager, Natural Resources Canada) for championing this and related projects and to Brian Crosland (Natural Resources Canada) for lending computational resources.
This research has been supported by Natural Resources Canada (grant no. CHGHG IETS19103) and the Natural Sciences and Engineering Research Council of Canada (NSERC, grant nos. 479641, 06632, 522658).
This paper was edited by Oleg Dubovik and reviewed by three anonymous referees.
Acosta, I., Navarro, J., and Sendra, J. J.: Lighting design in courtyards: Predictive method of daylight factors under overcast sky conditions, Renew. Energ., 71, 243–254, https://doi.org/10.1016/j.renene.2014.05.020, 2014.
Allen, D. T. and Torres, V. M.: TCEQ 2010 Flare study final report, Texas Commission on Environmental Quality (TCEQ), Austin, Texas, USA, available at: http://www.tceq.texas.gov/assets/public/implementation/air/rules/Flare/2010flarestudy/2010flarestudyfinalreport.pdf (last access: 22 February 2021), 2011.
Andreae, M. O. and Gelencsér, A.: Black carbon or brown carbon? The nature of lightabsorbing carbonaceous aerosols, Atmos. Chem. Phys., 6, 3131–3148, https://doi.org/10.5194/acp631312006, 2006.
Anenberg, S. C., Schwartz, J., Shindell, D., Amann, M., Faluvegi, G., Klimont, Z., JanssensMaenhout, G., Pozzoli, L., Van Dingenen, R., Vignati, E., Emberson, L., Muller, N. Z., West, J. J., Williams, M., Demkine, V., Hicks, W. K., Kuylenstierna, J., Raes, F., and Ramanathan, V.: Global Air Quality and Health Cobenefits of Mitigating NearTerm Climate Change through Methane and Black Carbon Emission Controls, Environ. Health Persp., 120, 831–839, https://doi.org/10.1289/ehp.1104301, 2012.
Becker, H. A. and Liang, D.: Total emission of soot and thermal radiation by free turbulent diffusion flames, Combust. Flame, 44, 305–318, https://doi.org/10.1016/00102180(82)900803, 1982.
Bond, T. C., Doherty, S. J., Fahey, D. W., Forster, P. M., Berntsen, T. K., De Angelo, B. J., Flanner, M. G., Ghan, S., Kärcher, B., Koch, D., Kinne, S., Kondo, Y., Quinn, P. K., Sarofim, M. C., Schultz, M. G., Schulz, M., Venkataraman, C., Zhang, H., Zhang, S., Bellouin, N., Guttikunda, S. K., Hopke, P. K., Jacobson, M. Z., Kaiser, J. W., Klimont, Z., Lohmann, U., Schwarz, J. P., Shindell, D. T., Storelvmo, T., Warren, S. G., and Zender, C. S.: Bounding the role of black carbon in the climate system: A scientific assessment, J. Geophys. Res.Atmos., 118, 5380–5552, https://doi.org/10.1002/jgrd.50171, 2013.
CIE: Spatial Distribution of Daylight – CIE Standard General Sky, Commission Internationale de l'Eclairage, Vienna, Austria, 2003.
Conrad, B. M.: SetupSkyLOSA: A MATLAB Tool to support Acquisition of SkyLOSA Data, Zenodo, https://doi.org/10.5281/zenodo.3908540, 2020.
Conrad, B. M. and Johnson, M. R.: Field measurements of black carbon yields from gas flaring, Environ. Sci. Technol., 51, 1893–1900, https://doi.org/10.1021/acs.est.6b03690, 2017.
Conrad, B. M. and Johnson, M. R.: Mass absorption crosssection of flaregenerated black carbon: Variability, predictive model, and implications, Carbon, 149, 760–771, https://doi.org/10.1016/j.carbon.2019.04.086, 2019.
Conrad, B. M., Thornock, J. N., and Johnson, M. R.: Beam steering effects on remote optical measurements of pollutant emissions in heated plumes and flares, J. Quant. Spectrosc. Ra., 254, 107191, https://doi.org/10.1016/j.jqsrt.2020.107191, 2020a.
Conrad, B. M., Thornock, J. N., and Johnson, M. R.: The effect of multiple scattering on optical measurement of soot emissions in atmospheric plumes, J. Quant. Spectrosc. Ra., 254, 107220, https://doi.org/10.1016/j.jqsrt.2020.107220, 2020b.
Darula, S. and Kittler, R.: CIE general sky standard defining luminance distributions, in eSIM 2002, available at: http://www.ibpsa.org/?page_id=291 (last access: 22 February 2021), International Building Performance Simulation Association, Montreal, Canada, 11–13 September 2002.
Elvidge, C. D., Baugh, K. E., Tuttle, B. T., Howard, A. T., Pack, D. W., Milesi, C., and Erwin, E. H.: A Twelve Year Record of National and Global Gas Flaring Volumes Estimated Using Satellite Data: Final Report to the World Bank, available at: https://eogdata.mines.edu/interest/flare_docs/DMSP_flares_20070530_b.pdf (last access: 22 February 2021), 2007.
Elvidge, C. D., Ziskin, D., Baugh, K. E., Tuttle, B. T., Ghosh, T., Pack, D. W., Erwin, E. H., and Zhizhin, M.: A Fifteen Year Record of Global Natural Gas Flaring Derived from Satellite Data, Energies, 2, 595–622, https://doi.org/10.3390/en20300595, 2009.
Elvidge, C. D., Zhizhin, M., Baugh, K., Hsu, F., and Ghosh, T.: Methods for global survey of natural gas flaring from visible infrared imaging radiometer suite data, Energies, 9, 14, https://doi.org/10.3390/en9010014, 2016.
Gvakharia, A., Kort, E. A., Brandt, A. R., Peischl, J., Ryerson, T. B., Schwarz, J. P., Smith, M. L., and Sweeney, C.: Methane, black carbon, and ethane emissions from natural gas flares in the Bakken Shale, North Dakota, Environ. Sci. Technol., 51, 5317–5325, https://doi.org/10.1021/acs.est.6b05183, 2017.
Jacobson, M. Z.: Strong radiative heating due to the mixing state of black carbon in atmospheric aerosols, Nature, 409, 695–697, https://doi.org/10.1038/35055518, 2001.
Johnson, M. R., Wilson, D. J., and Kostiuk, L. W.: A fuel stripping mechanism for wakestabilized jet diffusion flames in crossflow, Combust. Sci. Technol., 169, 155–174, https://doi.org/10.1080/00102200108907844, 2001.
Johnson, M. R., Devillers, R. W., Yang, C., and Thomson, K. A.: SkyScattered solar radiation based plume transmissivity measurement to quantify soot emissions from flares, Environ. Sci. Technol., 44, 8196–8202, https://doi.org/10.1021/es1024838, 2010.
Johnson, M. R., Devillers, R. W., and Thomson, K. A.: Quantitative field measurement of soot emission from a large gas flare using skyLOSA, Environ. Sci. Technol., 45, 345–350, https://doi.org/10.1021/es102230y, 2011.
Johnson, M. R., Devillers, R. W., and Thomson, K. A.: A generalized skyLOSA method to quantify soot/black carbon emission rates in atmospheric plumes of gas flares, Aerosol Sci. Tech., 47, 1017–1029, https://doi.org/10.1080/02786826.2013.809401, 2013.
Kasten, F. and Young, A. T.: Revised optical air mass tables and approximation formula, Appl. Optics, 28, 4735–4738, https://doi.org/10.1364/AO.28.004735, 1989.
Kazemimanesh, M., Dastanpour, R., Baldelli, A., Moallemi, A., Thomson, K. A., Jefferson, M. A., Johnson, M. R., Rogak, S. N., and Olfert, J. S.: Size, effective density, morphology, and nanostructure of soot particles generated from buoyant turbulent diffusion flames, J. Aerosol Sci., 132, 22–31, https://doi.org/10.1016/j.jaerosci.2019.03.005, 2019.
Kittler, R., Kocifaj, M., and Darula, S.: Daylight Science and Daylighting Technology, Springer, New York, USA, 2012.
Köylü, Ü. Ö. and Faeth, G. M.: Structure of overfire soot in buoyant turbulent diffusion flames at long residence times, Combust. Flame, 89, 140–156, https://doi.org/10.1016/00102180(92)90024J, 1992.
McDaniel, M.: Flare efficiency study, United States Environmental Protection Agency, Research Triangle Park, North Carolina, USA, available at: http://nepis.epa.gov/Exe/ZyPURL.cgi?Dockey=P1003QGZ.txt (last access: 22 February 2021), 1983.
McEwen, J. D. N. and Johnson, M. R.: Black carbon particulate matter emission factors for buoyancydriven associated gas flares, J. Air Waste Manage., 62, 307–321, https://doi.org/10.1080/10473289.2011.650040, 2012.
Nakayama, M. K.: Asymptotically Valid Confidence Intervals for Quantiles and ValuesatRisk When Applying Latin Hypercube Sampling, available at: http://www.iariajournals.org/systems_and_measurements/tocv4n12.html (last access: 22 February 2021), Int. J. Adv. Syst. Meas., 4, 86–94, 2011.
Navvab, M., Karayel, M., Ne'eman, E., and Selkowitz, S.: Analysis of atmospheric turbidity for daylight calculations, Energ. Buildings, 6, 293–303, https://doi.org/10.1016/03787788(84)900616, 1984.
Petzold, A., Ogren, J. A., Fiebig, M., Laj, P., Li, S.M., Baltensperger, U., HolzerPopp, T., Kinne, S., Pappalardo, G., Sugimoto, N., Wehrli, C., Wiedensohler, A., and Zhang, X.Y.: Recommendations for reporting “black carbon” measurements, Atmos. Chem. Phys., 13, 8365–8379, https://doi.org/10.5194/acp1383652013, 2013.
Pohl, J. H., Tichenor, B. A., Lee, J., and Payne, R.: Combustion Efficiency of Flares, Combust. Sci. Technol., 50, 217–231, https://doi.org/10.1080/00102208608923934, 1986.
Popovicheva, O. B., Evangeliou, N., Eleftheriadis, K., Kalogridis, A. C., Sitnikov, N., Eckhardt, S., and Stohl, A.: Black Carbon Sources Constrained by Observations in the Russian High Arctic, Environ. Sci. Technol., 51, 3871–3879, https://doi.org/10.1021/acs.est.6b05832, 2017.
Ramanathan, V. and Carmichael, G.: Global and regional climate changes due to black carbon, Nat. Geosci., 1, 221–227, https://doi.org/10.1038/ngeo156, 2008.
Reda, I. and Andreas, A.: Solar Position Algorithm for Solar Radiation Applications (Revised), Golden, Colorado, USA, 2008.
Sato, M., Hansen, J., Koch, D., Lacis, A., Ruedy, R., Dubovik, O., Holben, B., Chin, M., and Novakov, T.: Global atmospheric black carbon inferred from AERONET, P. Natl. Acad. Sci. USA, 100, 6319–6324, https://doi.org/10.1073/pnas.0731897100, 2003.
Schuster, G. L.: Inferring the Specific Absorption and Concentration of Black Carbon From Aeronet Aerosol Retrievals, PhD thesis, Pennsylvania State University, University Park, PA, USA, 85 pp., 2004.
Schwarz, J. P., Holloway, J. S., Katich, J. M., McKeen, S., Kort, E. A., Smith, M. L., Ryerson, T. B., Sweeney, C., and Peischl, J.: Black carbon emissions from the Bakken oil and gas development region, Environ. Sci. Tech. Let., 2, 281–285, https://doi.org/10.1021/acs.estlett.5b00225, 2015.
Sorensen, C. M.: Light Scattering by Fractal Aggregates: A Review, Aerosol Sci. Tech., 35, 648–687, https://doi.org/10.1080/02786820117868, 2001.
Stohl, A., Berg, T., Burkhart, J. F., Fjǽraa, A. M., Forster, C., Herber, A., Hov, Ø., Lunder, C., McMillan, W. W., Oltmans, S., Shiobara, M., Simpson, D., Solberg, S., Stebel, K., Ström, J., Tørseth, K., Treffeisen, R., Virkkunen, K., and Yttri, K. E.: Arctic smoke – record high air pollution levels in the European Arctic due to agricultural fires in Eastern Europe in spring 2006, Atmos. Chem. Phys., 7, 511–534, https://doi.org/10.5194/acp75112007, 2007.
Trivanovic, U., Sipkens, T. A., Kazemimanesh, M., Baldelli, A., Jefferson, A. M., Conrad, B. M., Johnson, M. R., Corbin, J. C., Olfert, J. S., and Rogak, S. N.: Morphology and size of soot from gas flares as a function of fuel and water addition, Fuel, 279, 118478, https://doi.org/10.1016/j.fuel.2020.118478, 2020.
U.S. EPA: AP 42 – Compilation of Air Pollutant Emission Factors, Volume I, Sect. 13.5, Industrial Flares, Research Triangle Park, North Carolina, USA, available at: https://www3.epa.gov/ttn/chief/ap42/ch13/index.html (last access: 22 February 2021), 2018.
United Nations Economic Commission for Europe: 1999 Protocol to Abate Acidification, Eutrophication and Groundlevel ozone to the Convention on Longrange Transboundary Air Pollution, ECE/EB. AIR/114, available at: https://unece.org/environmentpolicyair/protocolabateacidificationeutrophicationandgroundlevelozone (last access: 22 February 2021), 2012.
Watanabe, T., Takamatsu, T., and Nakajima, T. Y.: Evaluation of Variation in Surface Solar Irradiance and Clustering of Observation Stations in Japan, J. Appl. Meteorol. Clim., 55, 2165–2180, https://doi.org/10.1175/JAMCD150227.1, 2016.
Weyant, C. L., Shepson, P. B., Subramanian, R., Cambaliza, M. O. L. L., Heimburger, A., Mccabe, D., Baum, E., Stirm, B. H., and Bond, T. C.: Black carbon emissions from associated natural gas flaring, Environ. Sci. Technol., 50, 2075–2081, https://doi.org/10.1021/acs.est.5b04712, 2016.
Wieneke, B.: PIV uncertainty quantification from correlation statistics, Meas. Sci. Technol., 26, 074002, https://doi.org/10.1088/09570233/26/7/074002, 2015.
Wong, I. L.: A review of daylighting design and implementation in buildings, Renew. Sust. Energ. Rev., 74, 959–968, https://doi.org/10.1016/j.rser.2017.03.061, 2017.
 Abstract
 Introduction
 SkyLOSA measurement
 General uncertainty analysis methodology
 Results and discussion
 Conclusions
 Appendix A: Implementation details of the updated MC method
 Data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 SkyLOSA measurement
 General uncertainty analysis methodology
 Results and discussion
 Conclusions
 Appendix A: Implementation details of the updated MC method
 Data availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References