An uncertainty-based protocol for the setup and measurement of soot–black carbon emissions from gas flares using sky-LOSA

Abstract. 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 sky-LOSA 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,
non-linear impacts on achievable measurement uncertainties. To address this
challenge, this paper presents a prescriptive measurement protocol and
associated open-source software tool that simplify acquisition of sky-LOSA
data in the field. Leveraging a comprehensive Monte Carlo-based general
uncertainty analysis (GUA) to predict measurement uncertainties over the
entire breadth of possible measurement conditions, general heuristics are
identified to guide a sky-LOSA user toward optimal data collection. These
are further extended in the open-source 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 sky-LOSA 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 sky-LOSA technique and ultimately the availability of field
measurement data for flare black carbon emissions.



Introduction
Gas flaring is a routine practice in the oil and gas industry in which producers and refiners burn excess or unwanted gases in open-atmosphere flames, typically from vertical pipe stacks. Flaring is generally preferable to the venting of gases to atmosphere because it reduces carbon dioxide (CO2)-equivalent emissions; however, flaring still emits potent climate-forcing 25 pollutants directly to atmosphere (Allen and Torres, 2011;Johnson et al., 2001Johnson et al., , 2011Johnson et al., , 2013McDaniel, 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 CO2 (Bond et al., 2013;Jacobson, https://doi.org/10.5194/amt-2020-255 Preprint. Discussion started: 24 August 2020 c Author(s) 2020. CC BY 4.0 License. 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(Elvidge et al., , 2009(Elvidge et al., , 2015, 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., 2006). To attributeand, hence, report and regulatesoot/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 10 factors are crude single-valued parameters that link emitted soot mass to volume/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 real-world flares under field conditions, to develop and validate accurate flare soot/BC 15 emission factor models.
At present, there are only two published methods for the quantitative measurement of soot emissions from individual infield flares. One technique employs aircraft-based sampling of a flare plume (Gvakharia et al., 2017;Weyant et al., 2016), where measurements of soot, methane, and CO2 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 ground-based remote optical 20 measurement called sky-LOSA (line-of-sight attenuation using skylight; Conrad and Johnson, 2017;Johnson et al., 2010Johnson et al., , 2011Johnson et al., , 2013. Sky-LOSA quantifies time-resolved 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 off-site compositional analysis, which enables the direct calculation of soot yield for a targeted flare. To date, sky-LOSA has been deployed on 11 field measurement campaigns in Uzbekistan, Mexico, Ecuador, and Canada, providing 28 measurements of 25 soot emissions from 17 unique flares (Conrad and Johnson, 2017;Johnson et al., 2011Johnson et al., , 2013. The key component of a sky-LOSA measurement is the quantification of plume soot loading using image data, via analysis 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 to permit mass emission rate calculation. Uncertainties in sky-LOSA-calculated emission rate are computed under a Monte Carlo (MC) framework and are 30 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 sky-LOSA camera relative to the horizon and sun. Consequently, a sky-LOSA 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 https://doi.org /10.5194/amt-2020-255 Preprint. Discussion started: 24 August 2020 c Author(s) 2020. CC BY 4.0 License.
technique accessible to end-users, enabling an increase in flare soot emissions data, a standardized data acquisition protocol for sky-LOSA is required.
The objective of this work is to complete a general uncertainty analysis (GUA) for the sky-LOSA measurement technique that provides uncertainty-based guidance to an end-user regarding the setup of equipment and acquisition of sky-LOSA data through an accompanying open-source software tool. A summary of sky-LOSA theory, referring to derivations in previous 5 works, is first provided in Section 2 of this manuscript. The GUA methodology is summarized in section 3, including special provisions necessary to reduce the computational burden of the MC-based approach (section 3.1). Representative results from the MC GUA are shown in section 4.1 and general heuristics for the acquisition of sky-LOSA data, including new observations based on MC GUA results are summarized in sections 4.1.1 and 4.1.2. To provide case-by-case guidance, a new open-source software tool to calculate sky-LOSA measurement uncertainty is introduced in section 4.2. Finally, in section 4.3, the software 10 tool is used in a case study that analyses 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 sky-LOSA camera positioning and pointing to minimize measurement uncertainties, ultimately contributing to the standardization of the sky-LOSA measurement technique.

Sky-LOSA measurement 15
The sky-LOSA theory was first validated by Johnson et al., (2010) and is summarized in full by Johnson et al. (2013).
Development of the theory begins with Fig. 1, which shows an example sky-LOSA image for computation of time-resolved soot emission rate from a soot-laden flare plume in the Montney formation of Alberta, Canada. A highly linear, grayscale, scientific-CMOS camera (e.g., pco edge 5.5) is used to obtain upwards of ten minutes of high-speed image data of the flare and turbulent, soot-laden, atmospheric plume. Pseudo 16-bit images are acquired at framerates of 25−50 Hz with a narrow 20 mid-visible bandpass filter (531 ± 20 nm) to yield a scene of spectrally integrated light intensity. Fig. 1 is an example control surface ( ) of specified radius ( [m]), through which the instantaneous mass emission rate of soot (̇ [g s −1 ]) may be computed. For an arbitrary control surface in three dimensions, the instantaneous mass flux of soot through the surface is:
(2), sky-LOSA thus requires knowledge of three items to compute the emission rate: spatial scaling of the image plane to accurately quantify , the velocity field of the plume within the image plane (yielding 5 ( , )), 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 sky-LOSA 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 framerate and exposure, the two-dimensional plume velocity field over the image plane is estimated via image correlation/particle image velocimetry, using a third party software suite such as LaVision DaVis 8.4 that 10 includes a means of uncertainty quantification (Wieneke, 2015). Finally, the novel enabling aspect of sky-LOSA 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 sky-LOSA-measurements of time-resolved soot emission rate are available in previous works (Conrad and Johnson, 2017;Johnson et al., 2011Johnson et al., , 2013). 15 Figure 2 shows an example positioning/pointing of the sky-LOSA camera and an optical axis/LOS within the surrounding skydome. 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 sky-LOSA algorithm. Secondly, the sky is modelled as a diffuse, 20 polarized source concomitant with atmospheric scattering of solar radiation. The distribution of skylight intensity ( ( , , , , ) [W m −2 sr −1 ]) and the incident intensity along the LOS ( ( , , , ) [W m −2 sr −1 ]) are considered using the standard models of the Commission Internationale de l'Eclairage (CIE, 2003), where the index ∈ {1 … 15} indicates a specific CIE sky "type". Finally, ground-level normal solar irradiance ( [W m −2 ]) is estimated using in-field, image-based measurements of the sun (Johnson et al., 2013) or modelled using the CIE models as described in Sect. 3.1.3. With this 25 radiative transfer model, sky-LOSA-quantification of soot mass column density proceeds with the radiative transfer equation where is the measured "transmitted" intensity at the camera [W m −2 sr −1 ], ( ) is the mass-normalized extinction crosssection of soot [m 2 g −1 ] that is a function of eight soot properties represented by the vector , ( , ) is a local source radiant intensity per unit mass [W sr −1 g −1 ] along the measurement path, and the sky-LOSA camera is located at = . Optical properties of soot are computed from the fundamental properties in using Rayleigh−Debye−Gans theory for polydisperse fractal aggregates (RDG−PFA; e.g., Sorensen, 2001).
Using mean value theorem, a path-averaged source radiant intensity ( ̅ ( ) [W sr −1 g −1 ]) can be introduced to simplify the RTE: where ( ) ( ) ⁄ = 1 has also been introduced to the source term. The resulting ratio in front of the second integral 5 represents a path-averaged source intensity ( ̅ [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 = / is the observed transmittance of the plume [−] and = ̅ / is a term that corrects for brightening of the 10 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 ( ( , )) is composed of thermal emission and inscattering components. For sky-LOSA however, where measurements are performed in the mid-visible spectrum and plume temperatures are near ambient, thermal emission negligibly contributes to the source term of the RTE . By contrast, diffuse skylight and direct solar radiation can significantly augment the RTE via inscattering into the optical axis; hence, for sky-LOSA, 15 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 multiple scattering events may occur during light's transmission into the camera. A recent simulation effort, however, has shown that complex multiple scattering effects can be accurately modelled to permit quantification of the inscattering correction ( ) . This approach requires calculation 20 of the inscattering correction using a single-scattering assumption (1SA, subscript "1"): where ( ) is the single-scattering albedo of the polydisperse soot population [−], ( ; ) is the scattering phase function of soot [sr −1 ], and angles and 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 ( 1. ) and sun ( 1. ) components. Following Conrad et al. (2020), the 1SAestimated inscattering correction permits calculation of the idealized transmittancei.e., * = * ( , 1 ( , , , , )).
In the sky-LOSA algorithm, uncertainty in 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 sky-LOSAcomputed soot mass column density is sensitive to , , , , and . Of these five variables, only the pointing of the 10 camera may be controlled by the user, through the selection of angles and . Critically, these user-selected angles are known to influence sky-LOSA measurement uncertainty and, in extreme cases, can even preclude computation of soot emission rate.
For end-users of sky-LOSA, 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 15 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 sky-LOSA-computed 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 sky-LOSA-computation of soot mass column density ( ′ ) as a function of user-selectable ( and ) and uncontrollable parameters ( , , and ) under generalized conditions. These data permit the development of broad 20 heuristics and, ultimately, an easy-to-use software tool to provide specific case-by-case constraints on camera position and pointing.

General uncertainty analysis methodology
The goal of the present general uncertainty analysis (GUA) is to guide a sky-LOSA user in choosing a sky-LOSA camera position/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 post-processing analysis. To provide generalized recommendations, the GUA quantifies measurement uncertainty in soot mass column density for a selected camera 5 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 GUA and special considerations to make this work tractable.

Expansion of the scattering phase function
For a given (modelled) skylight intensity distribution, measured/modelled solar irradiance, camera pointing, and set of soot 10 properties, the 1SA-estimated inscattering correction ( 1 ) can be directly calculated via Eq. (7). One significant challenge, and currently the time-limiting computation in sky-LOSA processing, is the calculation of the skylight component ( 1, ) via numerical integration over three dimensions: , , and , 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, ( , )) of soot, which includes a computationally burdensome hypergeometric series in its solution (Sorensen, 2001). For the present GUA, where the 15 inscattering correction must be computed over a five-dimensional domain ( , , , , ), 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 MC-randomized soot properties , this procedure allows the SPF to be represented as: termed the l th -order soot coefficient) for the set of soot properties ( ) computed via: which can be accurately and efficiently computed via Gauss-Legendre quadrature (Schuster, 2004). Since the soot coefficient decreases towards zero as the order approaches infinity, the infinite series expansion of the SPF can be truncated at a sufficiently large index ( ) with negligible error.

Sky and sun inscattering
Introduction of the Fourier-Legendre-expanded SPF into the sky and sun components of the inscattering correction ( 1, and 1, , Eq. (7)) yields: where Ψ , ( , , , ) and Ψ , ( , , , ) 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 5 1, , which vastly reduces computational burden. Furthermore, with this formulation, the soot coefficients (functions of ) and sky/sun coefficients (functions of , , , and ) do not share any independent variables and can therefore be independently pre-computed.

Solar irradiance
As noted in Sect. 2, ground-level solar normal irradiance -( , ), usually measured via solar images taken in the field -10 can be modelled using the CIE skylight models. To accomplish this, the incident intensity-normalized solar normal irradiance is expanded: where ℎ ( , ) is the diffuse horizontal irradiance, calculable for the CIE models via numerical integration of: which is independent of the value of .
The ratio of solar normal to diffuse horizontal irradiance is complex to quantify in a general sense as it is a function of 15 atmospheric composition. However, for the purposes of the present GUA it is modelled as follows: The numerator of the righthand side is the ground-level solar horizontal irradiance, calculated as the product of the extraterrestrial (ET, subscript "o") solar horizontal irradiance ( ℎ, ( , )) with an exponential representing attenuation through https://doi.org/10.5194/amt-2020-255 Preprint. Discussion started: 24 August 2020 c Author(s) 2020. CC BY 4.0 License.
the atmosphere. In computation of the latter, ( ) is the relative air mass quantifying the amount of air in the atmosphere at the solar zenith angle relative to the vertical direction, * ( ) is the ideal extinction for a clean atmosphere at a given relative air mass, and ( ) is the model-dependent turbidity factor representing the number of clean atmospheres required to represent observations. In the denominator, the cosine is included to transform the numerator into the required ground-level solar normal irradiance. In the present work, typical values of the turbidity factor ( ( )) and irradiance ratio ℎ ( , )/ ℎ, ( , ) are 5 taken from Darula and Kittler (2002) and Kittler et al. (2012), while their recommended formulations for ( ) (Kasten and Young, 1989) and * ( ) (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 ( ( )) 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), ( ) ∼ (0.75,1.25) and for models with an obstructed sun (types 10 1−6), ( ) ∼ (0,1.25). These prior distributions of ( ) were based on observations by Watanabe et al. (2016), who studied the "clearness index" (ground-level horizontal normalized by ET horizontal irradiance) over five 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 variance-equivalent 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-15 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).

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 20 required to represent the soot SPF according to the procedure of Schuster (2004) was typically ( ) = 76 (median). In the most extreme case however, representing a strongly forward scattering particle (corresponding to large soot aggregate size), ( ) reached 698, suggesting that pre-computation 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 increasingly cumbersome as the order gets larger. Importantly though, like the 25 soot coefficients, the magnitude of the sky coefficients approaches zero as the order approaches infinity; therefore, the product of Φ ( ) and Ψ , ( , , , ) more-rapidly decreases in magnitude as increases than either component alone. This randomized sets of ( , , , , ), the median relative difference was just 2.3 × 10 −7 . This implies that it is acceptable to impose that Ψ , ( , , , ) = 0, ∀ > 200.

Sky model groups
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 5 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 sky-LOSAcomputed soot mass column density through use of a single selected CIE sky model in the MC method. To permit capture of CIE sky model error in the GUA, like skies were compiled into sky model "groups" that have similar properties but differing model coefficients and, hence, directional variabilitythe derived sky groups ( ∈ { … }) are summarized in Table 1. By 10 randomly selecting the sky group's component sky models under the MC framework, uncertainty through use of the CIE models is propagated into sky-LOSA-computation of soot mass column density.
Sky group A corresponds to overcast and partly cloudy conditions with an obscured sun where typical turbidity factors of the component skies ( ( )) are high. Sky group 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 ET solar radiation (high 15 ℎ ( , )/ ℎ, ( , )). Sky groups C and D capture the low-turbidity clear CIE sky models. Sky group C includes all five clear sky models, while sky group D excludes CIE sky type 12, the lowest turbidity (i.e., cleanest atmosphere) model. Sky group 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 precludes CIE sky type 12 as a reasonable model, such that sky group D can be used for the case of clear skies in heavily industrial locales. 20 Table 2 summarizes the independent, pre-computed, and randomized 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
In a standard MC analysis, the above procedure would be iterated upon 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 employedspecifically, combined multiple Latin hypercube sampling (CM-LHS) summarized by Nakayama (2011). This variance reduction technique has been used previously for sky-5 LOSA (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 1,000 Latin hypercube sampled data) CM-LHS MC draws were completed. The GUA MC approach permitted pre-computation of the soot coefficients, single-scattering albedo, and mass-normalized extinction cross-section for pre-drawn CM-LHS MC-randomized sets of soot properties ( ). Parallel pre-computation of the sky and sun coefficients was performed for each of the 15 CIE standard skylight intensity distributions and four sky groups over the angles , , and in increments 10 of 2°, and for 18 observed transmittances from 0.25−0.99. This amounted to execution of the CM-LHS MC analysis for almost 66 × 10 6 unique sky-LOSA conditions, permitting derivation of uncertainty statistics over the five independent variables listed in Table 2.

Relative uncertainty metric
To enable an objective comparison of sky-LOSA uncertainty as a function of the independent variables, a parameter describing 15 the relative uncertainty of MC-computed soot mass column density was required. Since relative uncertainty is a comparison of data variability to central tendency, a natural metric is the coefficient of variation (CV), which is the ratio of the standard deviation to the mean of a distribution or dataset. These moment-based statistics, however, are sensitive to outliers (e.g., Devore and Berk, 2012) or long-tailed/skewed distributions, causing some to suggest order statistics as more-robust metrics for variability and central tendency (e.g., Arachchige et al., 2019). For sky-LOSA, MC-and time-averaged soot mass emission 20 rates (moment-based statistic) have historically been reported alongside MC-computed 95% confidence intervals (CIs, order statistic). Therefore, for consistency with past and future sky-LOSA measurements, a CV-estimator based on these statistics was used in this work. For variability, the width of the 95% CI scaled by that of a normal standard distribution (≈ 3.92) was used as an estimator of the standard deviation, while the mean (̅) was employed as the measure of central tendency. The CVestimator for soot mass column density (or any MC-computed data ) is thus: 25 where subscript "95" signifies use of the 95% CI for variability and ℱ −1 ( ) is the q th quantile of data vector . 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, . This is because of solar radiation's influence on 1, and 10 1, . However, the rate at which relative uncertainty changes as a function of is sky model-(as well as solar zenith-and plume transmittance-) dependent. This is partly a consequence of how 1, varies with , but is mostly due to variability in 1, with . For example, the inscattering magnitude for sky group A is effectively due to 1, alone due to high turbidity that strongly attenuates direct solar radiation while, for sky group C, the sun's inscattering contribution ( 1, ) is significant due to the low turbidity of the sky. Thus, the results for sky group A largely represent the effect of on uncertainty through 15 1, while the results for sky group C include an additional (and the most extreme) effect through 1, . These observations imply that a constraint on is necessary, and that a stricter constraint is required for lower turbidity skies.
The mechanism for the decrease in relative uncertainty with increasing stems from the optical characteristics of soot particulate. Figure 4 shows statistics of the "energy distribution function" ( = ( ) ( , )/4 ) discussed by Conrad et al., (2020), which describes how soot particulate directionally scatters light on an energy-basis. The median and 95 of the 20 EDF as a function of angle are shown for the range of the MC-sampled soot properties ( ) used in the GUA. This is useful for the present discussion since the EDF contains all soot-dependent variables in the computation of 1, (see Eq. (7)). The figure shows that the EDF, and its influence on 1, , is much larger and much more uncertain as decreases. The influence of these statistics on column density uncertainty depends on the relative intensity of sunlightspecifically ( , ) ( , , , ) ⁄ in Eq. (7). If this value is small (highly turbid skies), then variability in the EDF is less important, 25 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 ground-level intensity of the sun, the uncertainty of 1, becomes minimal at > 90° and is within 1% of this minimum for > 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 30 to the camera inclination angle ( ). Referring to Fig. 3 Uncertainties for sky groups A and B tend to decrease as the camera inclines, while uncertainties for sky groups C and D increase. The severity of these trends increases with plume transmittance (effectively a reduction in the measured signal) and, 5 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 clear sky 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, low-turbidity skies, intensity gradients are large such that small changes in camera pointing can yield significant changes in inscattered light. By contrast, for higher-turbidity skies, gradients are dampened, and the effect of camera pointing is small. 10 Figure 5 also shows that the 95 of soot mass column density for sky group C tends to upper bound that of sky group D.

This observation holds under most combinations of the MC independent variables and is a result of the somewhat extreme
nature of CIE sky model 12 (lowest turbidity / cleanest atmosphere), which is excluded in sky group D. When considered in sky group 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 largerhence, inclusion of CIE sky model 12 typically increases relative 15 uncertainty. This implies that sky group C reliably imposes the largest constraints on camera positioning and can therefore be used to conservatively locate sky-LOSA equipment under clear skies. However, if in a dense industrial area, where the clearest sky model is not relevant, the less-constraining sky group D can be instead used as noted in Sect. 3.1.5.

Camera pointing heuristics
Upon arrival at a measurement facility, the sky-LOSA user's first task is to determine the position of the sky-LOSA camera 20 for data acquisition. This important decision can be made by considering viable camera pointings from GUA MC data through constraints on and . That is, for a user-identified sky group and plume transmittance and given the position of the sun, the sky-LOSA camera's pointing can be constrained based on a desired threshold in a relative uncertainty, 95 . Table 3  One additional consideration in the pointing of the sky-LOSA camera is the direction of plume propagation. Under quiescent conditions, buoyancy-driven 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 sky-LOSA camera, turbulent plume structures of differing 30 vorticity become overlapped from the camera's perspective. Therefore, it is best to position the sky-LOSA camera such that it points orthogonally to the wind direction, which minimizes out-of-image plane motion of the plume and yields the best data for velocimetry calculations. This should be viewed as a weak constraint on sky-LOSA data acquisition however, since the https://doi.org /10.5194/amt-2020-255 Preprint. Discussion started: 24 August 2020 c Author(s) 2020. CC BY 4.0 License. effect of uncertainty in estimated velocity on mass emission rate is generally negligible compared to that of column density uncertainty.

Further camera heuristics
Following selection of a permissible sky-LOSA camera position, the imaging optics must be chosen. Prime (fixed focal length) lenses are employed in the sky-LOSA technique to avoid ambiguity in optical magnification and, hence, spatial scaling 5 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 unsteadyi.e., moving with the windit is suggested to keep the flare flame approximately one-quarter of the smallest image dimension. By contrast, if the flame is steady, a flame length of approximately one-third of the smallest image dimension should be targeted. 10 For the sky-LOSA camera used by the authors (minimum sensor dimension of ~14 mm), a good rule-of-thumb is that the appropriate focal length ( , [mm]) will be on the order of: 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 full-open (typically f-number ≤ 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 sky-LOSA data for the desired duration; it is recommended, however, that a minimum 20 dataset of ten minutes be obtained to permit good convergence of the time-averaged soot emission rate.

Open-source software tool for simpler sky-LOSA setup -SetupSkyLOSA
While the camera pointing heuristics presented in Table 3 can be used to ensure that 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 25 better approach would be to use the wealth of computed GUA data to provide camera position and pointing constraints on a case-by-case basis. This is made possible using a new open-source software tool, SetupSkyLOSA (Conrad, 2020), that was developed as part of this work using the presented GUA MC data. This MATLAB-based application enables a sky-LOSA user to probe statistics of soot mass column density (such as 95 ) for their specific measurement conditions. The key output https://doi.org/10.5194/amt-2020-255 Preprint. Discussion started: 24 August 2020 c Author(s) 2020. CC BY 4.0 License. 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 flowchart describing the SetupSkyLOSA software's main procedure. For a user-inputted location and time, the software first determines the current position of the sun using an integrated solar position calculatora MATLAB implementation of the National Renewable Energy Laboratory's (NREL's) Solar Position Algorithm (SPA) (Reda and 5 Andreas, 2008). The SPA returns the solar zenith ( ) and absolute bearing of the sun ( , 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 ( ) of the most appropriate CIE sky model/group, an estimate of the observed plume transmittance at the sky-LOSA measurement wavelength ( ), and the desired statistic of soot mass column density ( ). With these inputs, SetupSkyLOSA then loads the GUA MC data for the selected sky model/group and interpolates for the desired statistic using 10 the current solar zenith and estimated plume transmittance.
At this point, the software has computed ( , ) for the user's current set of independent variables ( , , ). However, rather than plotting ( , )i.e., using the relative bearingthe software uses the known absolute bearing of the sun ( , computed by the SPA) to plot ( , ). 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 15 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 sky-LOSA 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 − domain, based on sensor dimensions, employed optics, and the 20 pointing of the centre of the image. This helps a user ensure that the entirety of the image frameincluding the eventual control surface used for emission rate calculationhas reasonable levels of measurement uncertainty, given a user-selected lens of known focal length. To support ideal velocity calculation, the software can also overlay camera pointings ( ) 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% out-of-image plane motion (cos −1 (0.05)). Two 25 additional utilities are not shown in Fig. 6. The "maximizer" utility computes the maximum of a chosen relative uncertainty statistic over a user-defined period. This tool allows a user to seek camera pointings that yield satisfactory uncertainties as the sun moves during the anticipated duration of the sky-LOSA measurement. The "positioner" utility takes the plotted relative uncertainty statistic as a function of camera pointing and provides region(s) where the sky-LOSA camera may be positioned relative to the flare stack given the flare stack height, maximum horizontal stand-off distance from the flare, and relative 30 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.

Case study -Atasta facility
The utility of the novel software tool, SetupSkyLOSA, is shown in this section via a case study of a sky-LOSA 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 5 to the national market. As shown in Fig. 7a For this case study, a sky-LOSA measurement of soot emissions from the central flare stack at the Atasta station was 10 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: 1. The sky-LOSA measurements occur on May 13, 2021, which is the date of that year that the sun most closely reaches the solar zenith ( = 0°). Assumption #2 implies that sky-LOSA data acquisition may occur under skies represented by any of sky groups A−D.

Predicted
Furthermore, assumption #3 suggests that the soot-laden flare plume propagates vertically from the flare stack and, therefore, the constraint on camera position with respect to wind direction is unimportant. The sky-LOSA user wishes to obtain sky-LOSA data with minimal measurement uncertainty, while also avoiding re-location of the sky-LOSA camera throughout the day, if possible. 25 Given the known GPS coordinates of the flare stack, measurement date, and the approximate plume transmittance, SetupSkyLOSA can be used to constrain sky-LOSA camera pointing for any sky condition and time of the day based on the 95 of soot mass column density. Since the user wishes to avoid re-location of the sky-LOSA 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 groups; results are shown in Fig. 8a−8d for sky groups A−D. Noting 30 the differing colour scales in the four figures, there is significant variability in sky-LOSA uncertainties for each of the sky groups, as in Fig. 3. For further context, the path of the sun over the measurement date is overlaid in Fig. 8a−8d in addition to a contour of 95 = 16.5%, which is within 1% of the best attainable uncertainty for these conditions. https://doi.org /10.5194/amt-2020-255 Preprint. Discussion started: 24 August 2020 c Author(s) 2020. CC BY 4.0 License.
Using the uncertainty data in Fig. 8a−8d, the "positioner" utility of SetupSkyLOSA can be employed to highlight where the sky-LOSA camera may be positioned relative to the flare stack. This was performed for each of the sky groups 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 groups. Permissible camera positions for sky group B exist beyond a small region near the stack tip, while those for sky group A are within an annular 5 region surrounding the flare stacksince the lower limit on the camera inclination angle in Fig. 8a imposes a maximum permissible stand-off distance. Sky group D contains two permissible regionsone to the South and one to the North of the flare stackwhile the most-constrained sky group C has one relatively small region to the North of the flare stack. Recalling assumption #2, that predicted sky conditions were uncertain, the sky-LOSA user should ideally position the camera at the intersection of the sky group-dependent permissible regions. This small area is outlined in black in the figure, is ~136 m due 10 North of the flare stack and is 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 sky-LOSA 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 15 uncertainty data enable optimal sky-LOSA camera positioning and pointing and also represent a first-order estimate of soot emission rate uncertainties that are computed in post-processing. Together with the additional utilities and general camera heuristics, this software tool permits a sky-LOSA user to obtain optimal sky-LOSA data that minimizes measurement uncertainties under generalized conditions.

Conclusions 20
A comprehensive Monte Carlo-based general uncertainty analysis (GUA) has been used to develop heuristics constraining the pointing/positioning of sky-LOSA 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 sky-LOSA. 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 25 parameters, the presented results provide generalized guidance to greatly simplify acquisition of optimal sky-LOSA data in the context of complex, non-linear measurement uncertainties. These are further extended in the open-source 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, software-displayed soot mass column density statistics provide the user with a first-order estimate of measurement uncertainty in soot/black carbon emission rate that otherwise is 30 only computable during post-processing. 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 open-source application can https://doi.org /10.5194/amt-2020-255 Preprint. Discussion started: 24 August 2020 c Author(s) 2020. CC BY 4.0 License.
hopefully facilitate broader use of the sky-LOSA technique and ultimately help increase the knowledge base of soot/black carbon emissions from gas flares.

Data availability
The developed software tool and associated data are available online as open-source and build distributions (Conrad, 2020).   ; right linear axis) of the "energy distribution function" ( = ( ) ( , )/ ) that dictates the fraction of incident light energy scattered through angle by soot. The magnitude and uncertainty in the EDF are much larger in the forward scattering direction (small ). As the scattering direction 10 exceeds ~90° however, the relative uncertainty in the EDF approaches a constant minimum, implying that the relative uncertainty in , is minimized for >°. The red line shows that relative uncertainty is within 1% of the minimum if >°.