Ground-Based Remote Sensing of Volcanic CO2 Fluxes at Solfatara (Italy) - Direct Versus Inverse Bayesian Retrieval

CO2 is the second most abundant volatile species of degassing magma. CO2 fluxes carry information of incredible value, such as periods of volcanic unrest. Ground-based laser remote sensing is a powerful technique to measure CO2 fluxes in a spatially integrated manner, quickly and from a safe distance, but it needs accurate knowledge of the plume speed. The latter is often difficult to estimate, particularly for complex topographies. So, a supplementary or even alternative way of retrieving fluxes would be beneficial. Here, we assess Bayesian inversion as a potential technique for the case of the volcanic crater of Solfatara (Italy), a complex terrain hosting two major CO2 degassing fumarolic vents close to a steep slope. Direct integration of remotely sensed CO2 concentrations of these vents using plume speed derived from optical flow analysis yielded a flux of 717 ± 121 t day−1, in agreement with independent measurements. The flux from Bayesian inversion based on a simple Gaussian plume model was in excellent agreement under certain conditions. In conclusion, Bayesian inversion is a promising retrieval tool for CO2 fluxes, especially in situations where plume speed estimation methods fail, e.g., optical flow for transparent plumes. The results have implications beyond volcanology, including ground-based remote sensing of greenhouse gases and verification of satellite soundings.


Introduction
Magma contains dissolved volatiles, which are released when magma rises and depressurizes. Carbon dioxide (CO 2 ) is the second most abundant volatile species of degassing magma [1,2]. Measuring CO 2 mass emission rates (fluxes) therefore provides a powerful sampling tool that can be used, for example, to constrain magma crystallization processes related to the magma dynamics [3], learn about the geochemical carbon cycle [4] or for monitoring and early warning of volcanic unrest [5,6]. The nested caldera of Campi Flegrei (CF) is currently in a state of unrest [7,8] and is, as far as immediate human casualties are concerned, perhaps the most dangerous volcano on Earth, since it is located right at the edge of the metropolitan area of Naples in Italy (Figure 1a). Approximately at the center of CF lies the crater of Solfatara. Solfatara is a tuff cone associated with hot soils, intense hydrothermal activity and fumaroles discharging water vapor and CO 2 as the two most abundant volatile species [7,9,10]. This area is therefore a crucial sampling spot for quantifying degassing in order to understand the future evolution of CF, immediately impacting civil protection measures.
Remote sensing, notably active, laser-based remote sensing is feasible to acquire spatially integrated measurements of CO 2 concentrations. This technique is therefore ideal to probe extended degassing geometry by means of a forward model. However, inversion has its own flaws. For instance, inverted model parameters are generally non-unique and the forward model may not account for all real-world parameters. All this may lead to erroneous fluxes. This paper compares two approaches with the aim to assess the usefulness of an inverse approach to supplement methods based on Equation (1), using Solfatara as an example case. To that end, the recently developed portable Laser Absorption Remote Sensing Spectrometer (LARSS) was used to acquire CO 2 concentrations at the two main vents of Solfatara, Bocca Nuova (BN) and Bocca Grande (BG), which were then used as input for two different retrieval methods to obtain the CO 2 flux of these vents:

•
A direct integration over the scanned area, following Equation (1), using plume speed retrieved from optical flow analysis. • A Bayesian inversion using a simple Gaussian dispersion forward model.
It was found that, under certain circumstances, Bayesian inversion might be a useful supplement to volcanic CO 2 flux retrieval methods, particularly in situations where plume speed estimation is difficult. This could be the case where image-based methods fail due to a lack of features traceable by digital image processing methods.
The rest of the paper is organized as follows. First, LARRS is briefly outlined, followed by an introduction of the two flux retrieval methods, direct integration and Bayesian inversion. Then, the field experiment is described and the data are presented. Finally, the retrieved CO 2 fluxes are discussed.
Remote Sens. 2018, 10,125 3 of 17 powerful option for plume speed estimation, provided the plume is detectable by the imaging device [2,12,29,30]. Since plume speed retrieval is a challenge in its own right, the central motivation of this paper is as follows: Can the plume speed estimation be circumvented to arrive at the flux? Inverse modeling provides a potential way to do so [22]. It accounts for gas dispersion and environmental factors such as degassing geometry by means of a forward model. However, inversion has its own flaws. For instance, inverted model parameters are generally non-unique and the forward model may not account for all real-world parameters. All this may lead to erroneous fluxes. This paper compares two approaches with the aim to assess the usefulness of an inverse approach to supplement methods based on Equation (1), using Solfatara as an example case. To that end, the recently developed portable Laser Absorption Remote Sensing Spectrometer (LARSS) was used to acquire CO2 concentrations at the two main vents of Solfatara, Bocca Nuova (BN) and Bocca Grande (BG), which were then used as input for two different retrieval methods to obtain the CO2 flux of these vents:

•
A direct integration over the scanned area, following Equation (1), using plume speed retrieved from optical flow analysis. • A Bayesian inversion using a simple Gaussian dispersion forward model.
It was found that, under certain circumstances, Bayesian inversion might be a useful supplement to volcanic CO2 flux retrieval methods, particularly in situations where plume speed estimation is difficult. This could be the case where image-based methods fail due to a lack of features traceable by digital image processing methods.
The rest of the paper is organized as follows. First, LARRS is briefly outlined, followed by an introduction of the two flux retrieval methods, direct integration and Bayesian inversion. Then, the field experiment is described and the data are presented. Finally, the retrieved CO2 fluxes are discussed.

Remote Sensing of CO 2 Concentrations with LARRS
LARRS has been developed in the ERC proof-of-concept project CarbSens with the aim to assess the commercial potential of a turnkey system to remotely measure greenhouse gas fluxes (www.carbsens.com). The instrument and its working principle are detailed elsewhere [13,31] so only a brief overview is given. LARSS consists of a main unit and a transmitter/receiver unit (TX/RX unit, Figure 2a). The latter comprises of the telescope, transmitter and an integrating sphere for power reference measurement. It is portable (mass: 10 kg main unit + 6 kg TX/RX unit), which allows it to be transported easily and set up at any kind of surface, such as house roofs or airplanes and has a range between~100 m and~1500 m. The CO 2 absorption line at 1572.335 nm is sampled at 40 wavelengths by sweeping the emission wavelength of a diode laser. The laser light is amplified, transmitted, backscattered at a topographic target and received by the telescope (Figure 2b). After the detected signal is digitized, the optical transmittance of the telescope's viewing path is deduced for each of the 40 wavelengths ( Figure 2b). A model absorption spectrum is fitted to the 40 measured transmittances, resulting in a best estimate of the path-averaged CO 2 path density (in m −2 ), which is computed as where F is the best fit transmittance and ∆σ the differential absorption cross section (in m 2 ) of CO 2 . The ranges r i (in m) for a given angular position i are measured with a range finder LIDAR aligned with the telescope.
Remote Sens. 2018, 10, 125 4 of 17 start angle; (d) 3D view of the scanning geometry. The apparent bending of the beam paths results from projecting them onto the ground.

Remote Sensing of CO2 Concentrations with LARRS
LARRS has been developed in the ERC proof-of-concept project CarbSens with the aim to assess the commercial potential of a turnkey system to remotely measure greenhouse gas fluxes (www.carbsens.com). The instrument and its working principle are detailed elsewhere [13,31] so only a brief overview is given. LARSS consists of a main unit and a transmitter/receiver unit (TX/RX unit, Figure 2a). The latter comprises of the telescope, transmitter and an integrating sphere for power reference measurement. It is portable (mass: 10 kg main unit + 6 kg TX/RX unit), which allows it to be transported easily and set up at any kind of surface, such as house roofs or airplanes and has a range between ~100 m and ~1500 m. The CO2 absorption line at 1572.335 nm is sampled at 40 wavelengths by sweeping the emission wavelength of a diode laser. The laser light is amplified, transmitted, backscattered at a topographic target and received by the telescope (Figure 2b). After the detected signal is digitized, the optical transmittance of the telescope's viewing path is deduced for each of the 40 wavelengths ( Figure 2b). A model absorption spectrum is fitted to the 40 measured transmittances, resulting in a best estimate of the path-averaged CO2 path density (in m −2 ), which is computed as where is the best fit transmittance and σ Δ the differential absorption cross section (in m 2 ) of CO2. The ranges i r (in m) for a given angular position i are measured with a range finder LIDAR aligned with the telescope.

Flux Retrieval from Direct Integration
Profiles of CO2 concentrations, i.e., CO2 concentrations versus heading angle, are attained by scanning the TX/RX unit across a degassing plume. CO2 fluxes are computed following Equation (1) in discretized form as

Flux Retrieval from Direct Integration
Profiles of CO 2 concentrations, i.e., CO 2 concentrations versus heading angle, are attained by scanning the TX/RX unit across a degassing plume. CO 2 fluxes are computed following Equation (1) in discretized form as where v pl refers to the component of the plume transport speed perpendicular to the plane of the CO 2 concentration profile, i.e., the component perpendicular to the plane of the scan. In this study, this is the vertical component. For simplicity, it is referred to as plume speed hereafter. M CO 2 is the molar mass of CO 2 (in kg mol −1 ) and N A is Avogadro's constant (6.02214086 × 10 23 mol −1 ). ∆β is the constant scan angle increment. N col pl , the background-corrected, or in-plume path density of CO 2 , is retrieved by subtracting the total CO 2 path density by the ambient CO 2 path density measured outside the plume, i.e., N col pl (r i ) = N col (r i ) − N col bg (r i ), where N col (r i ) is the total path density from Equation (2) and N col bg (r i ) depicts the ambient CO 2 path density. Background-corrected in-plume mixing ratios, hereafter referred to as concentrations, are needed for the Bayesian inversion and are retrieved as where N air is the number density of air, computed using meteorological data.

Plume Speed Retrieval Using Optical Flow Analysis
In previous works, we used a video tracking method to estimate the plume speed [12]. The reasonable assumption was made that close to the vent the visible condensed water vapor aerosol of the plume propagated with the same velocity as the volcanic CO 2 . The implementation of a simple optical flow (OF) analysis presented here is a further evolution of this technique, serving the need for a more precise plume velocity and a more rigorous uncertainty estimate of that plume velocity. The two major physical sources of this uncertainty are turbulences and heterogeneities of plume speed across the plume related to different gas temperatures.
Optical flow analysis is becoming increasingly popular in conjunction with UV-sensitive cameras to measure SO 2 fluxes in volcanology [32][33][34] and has recently been used in the visible region for CO 2 flux estimation [5]. OF calculates the displacement of image features between subsequent video frames, yielding a displacement vector field for each analyzed pixel by solving where I is the pixel intensity, which is assumed constant between two frames under the assumption of small displacements in space (∆x horizontal and ∆y vertical, in pixels) and time ∆t.
There is a fundamental difference between applying OF to SO 2 camera footage and the present application. While for SO 2 flux measurement, the camera's imaging chip measures the gas concentration and in parallel tracks plume movement, here the camera's sole purpose is tracking the plume movement, while the light detector to measure the gas concentration is located inside the telescope's optical path ( Figure 2a). Therefore, a single estimate of the predominant plume velocity is required, that is, a velocity v pl representing the entire plume motion. In general, this would be a plume speed vector describing the dominant plume speed and direction (as in [35]). For the present application, however, it is not the magnitude of the plume speed vector that is needed but its component perpendicular to the scanning plane, which here is the vertical plume speed component, i.e., the vertical displacement ∆y with time.
The Farnebäck algorithm [36] implemented in the Python OpenCV computer vision library is used to solve Equation (5) for the pixel displacement vector components ∆x and ∆y for each pixel per frame, yielding a displacement vector field. This allows the spatial variability of the displacements to be assessed and hence the plume speed. In the following, the basic equations to retrieve the plume speed v pl from OF analysis of the plume motion are detailed by considering vertical displacements only. The considerations for horizontal displacements are equivalent. The predominant displacement ∆y j between frames j and j + 1, associated with the jth time step is computed as The first factor coverts pixels into meters. R is the distance (in m) between the camera and the plume, ∆ pix is the effective pixel pitch (in m) of the photo sensor inside the camera and f the focal length (in m) of the camera. In practice, the direction of v pl is not always perpendicular to the plane of the scan (and thus to the camera's line of sight), underestimating the displacement by a factor of sin α y , which is therefore used to correct for that. α y is the angle between the horizontal plane of the scan and the vertical plane of the plume propagation (or horizontal plane of plume propagation for a vertical scan). For each frame pair, the displacements follow a statistical distribution (number of pixels per displacement). The second factor (the sum) depicts the first moment of this distribution, where p(∆y l ) is the relative occurrence (or weight) of the vertical displacements ∆y l in bin l and m is the number of bins.
Optical flow analysis needs two subsequent frames to compute the optical flow velocity as v pl = ∆y /∆t, where ∆t is the magnitude of the time step between the frames. However, fluctuations in displacement due to turbulent plume motion that may result in non-linear displacement with time may not be captured by using two frames only. Figure 3 shows the optical flow field for the main vents BN and BG, computed from only three consecutive frames with ∆t = 0.5 s, which does not appear to be completely laminar at places. Therefore, all available frames are considered, which means the aforementioned analysis is repeated for a set of subsequent frame pairs, leading to a cumulative vertical displacement over time where M is the total number of frames considered. The slope of the linear fit of y (t j ) versus t j = j∆t yields the plume speed, representing a best estimate of a constant plume speed. The uncertainty of the plume speed is derived from the root mean square error (RMSE) of the fit and is a measure of the laminarity of the plume. In the following, techniques applied here to minimize noisy displacements are detailed. To focus the analysis on physical plume motion only and exclude parts of the plume not probed by the laser beam of LARRS, only those displacements inside a region of interest (ROI) are considered ( Figure 3). Figure 4 shows a video frame captured from position A (Figure 1b,c), containing the Solfatara fumarolic main vents BG and BN in the background. The ROI of the foreground (Figure 4a) does not contain any distinctively moving feature. The corresponding histogram of the displacements is thus associated with random displacement (noise) only, which is fairly normally distributed, centered at an expectation value that randomly changes with analyzed frame pair. In contrast, the distribution of the displacements for the ROI of the plume in Figure 4b has two distinct peaks. Both are normally distributed. This indicates that the physical plume displacements are superposed with noise displacements. Note that the bimodal distribution in Figure 4b also appeared if the ROI was completely filled out by the plume. The otherwise monomodal distribution suggests a very directional movement of the plume and that the plume and only the plume would be tracked by an OF algorithm. This moreover suggests that the sub-plumes coming from the main vents can be well described by a single plume moving at a single velocity. Displacements associated with the aforementioned noise are excluded as follows. Before computing the cumulative displacements y (t j ), the mean noise displacement is retrieved from the static image scenery (Figure 4a) as Displacements less than ∆y noise + 3σ noise (excluding 99.73% of the noise displacements), where σ noise is the sample standard deviation of the noise, are excluded from further analysis. As a result, ∆y noise = 0.013 pixels and σ noise = 0.06 pixels, leading to a lower limit vertical displacement threshold of 0.19 pixels. Furthermore, negative vertical displacements are discarded.

Flux Retrieval from Bayesian Inversion
As a method that sidesteps plume speed estimation, Bayesian inversion is used here to solve for the flux of the two main vents (BN, BG). The code, implemented in Python by [37], is freely available on Github. The inversion assumes a point source with known position. The two vents are therefore assumed to be a single vent with an effective vent location and so their combined CO 2 flux is inverted for. This is a reasonable assumption given their proximity and their common plume speed. The inverse problem can be formulated as P(Q|C) ∼ P(Q)P(C|Q), where P(Q|C) is the posterior probability density function (pdf), i.e., the probability of the flux Q given a set of measured CO 2 concentrations = C pl (Equation (4)). P(Q) is the prior pdf, which provides a constraint on the boundaries within which the flux is expected to be. P(C|Q) is called likelihood and is the probability that a flux gives rise to a set of CO 2 concentrations assuming a theoretical model, i.e., a forward model, which connects concentrations with flux. Taking input parameters of measurement geometry and meteorology into account, for each measured CO 2 mixing ratio, assuming a CO 2 flux, the forward model computes a corresponding theoretical value of the CO 2 mixing ratio. These values are then used to derive the likelihood function of the form where N is the number of measured concentrations (number of points in the concentration profile from LARRS), C mod depicts the modeled concentrations and τ the constant spread in the errors. The forward model assumes a Gaussian plume described as where u, v, w represent the three dimensions in space, χ the wind speed, H the height of the vent, and σ v and σ w are the standard deviations of the time-averaged plume concentration in downwind and vertical directions, respectively. They depend on the downwind distance and height, which are modeled using a power law taking into account the atmospheric stability class [37]. The latter is determined by means of the Monin-Obhukov length. The Monin-Obhukov length (in m) specifies the stability of the atmospheric boundary layer. It has been calculated for various meteorological conditions and times of day for the case of Solfatara in [38]. A value of −5.5 m was adopted, which accounts for the hydrodynamic instability of the atmosphere in broad daylight, representing the daytime atmospheric conditions at Solfatara very well. While the model predicts point concentrations, the input data from LARRS (Equation (4)) come in path-averaged concentrations. A subroutine therefore converts point concentrations to line averaged values, taking into account air temperature and pressure. The actual inversion process consists of iteratively maximizing the likelihood, i.e., finding the most likely flux, by performing a Markov Chain Monte Carlo (MCMC) search in parameter space, which in this case is one-dimensional and consists of a single flux value only. Thinning, i.e., considering only every nth sample of the posterior pdf, is applied to reduce autocorrelation between samples.

Field Measurements at Solfatara and the Data
LARRS was placed on a fixed location inside the crater, scanning a sector of 70 • (Figure 1b,c) on 26 May 2017 between 9:56 and 10:01 local time. The step-motor mounted TX/RX unit was pivoted with a velocity of 4 mrad s −1 using 1 s integration time per point, corresponding to a radial section of 40 cm at the target, the latter being the Solfatara crater wall (Figure 1b-d). Figure 5a shows the target height versus scanning angle relative to LARRS. The scanned sector contained the two main fumarolic vents BN and BG. Figure 5b shows the resulting path-averaged mixing ratios computed with Equation (4). Concentration peaks with path-averaged mixing ratios around 1600 ppm are located at~110 • and around~120 • , i.e., further northwards than the location of the main vents. This can be explained by the wind direction, which was 210 • (SW) on average, thus pushing the plumes northeastwards towards the crater edge.
N col bg were measured by a scan upwind, outside of any volcanic plume, using a hill range between 700 and 900 m distance as the target. The corresponding path-averaged CO 2 mixing ratio was found to be 499 ppm. For comparison, two in-situ measurements with a LI-COR analyzer were performed at points near the optical paths of LARSS, yielding CO 2 mixing ratios of 550 ppm and Remote Sens. 2018, 10, 125 9 of 16 560 ppm, respectively. These are remarkably high CO 2 concentrations, given that the wind came from the sea. The proximity of the measurement points to the road and the dense network of roads in that area may well cause these values [39]. Consequently, an ambient CO 2 mixing ratio of 499 ± 61 ppm was considered (Figure 5b).
Uncertainties of path-averaged CO 2 mixing ratios were usually between 2% and 5% (or 10 to 30 ppm), associated with a path-averaged detection limit of~10,000 ppm·m. The main source of uncertainty was the contribution of the instrument itself (baseline drift) and the fitting error. The uncertainty calculations and a detailed description of influences of various error sources are provided in [31]. Meteorological data (temperature, pressure, humidity) were retrieved from the INGV meteorological station at Solfatara. The average data during the scan is shown in Table 1. These data were used as input for the Bayesian inversion. Table 2 summarizes the inversion input data. 26 May 2017 between 9:56 and 10:01 local time. The step-motor mounted TX/RX unit was pivoted with a velocity of 4 mrad s −1 using 1 s integration time per point, corresponding to a radial section of ~40 cm at the target, the latter being the Solfatara crater wall (Figure 1b-d). Figure 5a shows the target height versus scanning angle relative to LARRS. The scanned sector contained the two main fumarolic vents BN and BG. Figure 5b shows the resulting path-averaged mixing ratios computed with Equation (4). Concentration peaks with path-averaged mixing ratios around 1600 ppm are located at ~110° and around ~120°, i.e., further northwards than the location of the main vents. This can be explained by the wind direction, which was 210° (SW) on average, thus pushing the plumes northeastwards towards the crater edge. N col bg were measured by a scan upwind, outside of any volcanic plume, using a hill range between 700 and 900 m distance as the target. The corresponding path-averaged CO2 mixing ratio was found to be 499 ppm. For comparison, two in-situ measurements with a LI-COR analyzer were performed at points near the optical paths of LARSS, yielding CO2 mixing ratios of 550 ppm and 560 ppm, respectively. These are remarkably high CO2 concentrations, given that the wind came from the sea. The proximity of the measurement points to the road and the dense network of roads in that area may well cause these values [39]. Consequently, an ambient CO2 mixing ratio of 499 ± 61 ppm was considered (Figure 5b).
Uncertainties of path-averaged CO2 mixing ratios were usually between 2% and 5% (or 10 to 30 ppm), associated with a path-averaged detection limit of ~10,000 ppm·m. The main source of uncertainty was the contribution of the instrument itself (baseline drift) and the fitting error. The uncertainty calculations and a detailed description of influences of various error sources are provided in [31]. Meteorological data (temperature, pressure, humidity) were retrieved from the INGV meteorological station at Solfatara. The average data during the scan is shown in Table 1. These data were used as input for the Bayesian inversion. Table 2 summarizes the inversion input data.

CO 2 Flux from Direct Integration
For the plume speed estimation, video footage (Video S1) of a length of 32 s with a resolution of 480 × 640 pixels 2 was acquired from the plume using a commercial digital camera (Canon IXUS 850 IS, pixel pitch 1.8 µm, frame rate 15 Hz) placed next to the telescope of LARRS at position A (Figure 1b). The time interval between two frames was chosen small enough to prevent undersampling of turbulent plume motion, which could lead to erroneous displacement magnitudes. A time step of ∆t = 0.5 s was chosen, corresponding to every 8th frame. To estimate α y , video footage was acquired from position B (Figure 1b) with the camera's line of sight almost in line with the prevailing wind. Deviations from the vertical plane were negligibly small and thus sin α y ≈ 1. Although not needed for the vertical plume speed component, α x was computed as the angle of the camera's line of sight and the wind direction, resulting in 79 • . Figure 6 shows the cumulative optical flow field of 32 s of video projected onto the last video frame. Figure 7a shows the cumulative vertical component of the displacement and Figure 7b the horizontal component for comparison. The error of the slope was of the order of 1% (cm s −1 ) and significantly lower than the uncertainty arising from shifting the ROI by a few pixels (a pixel corresponds to~5 cm). As a consequence, the uncertainty was derived by systematically varying the ROI location by a small amount and computing the standard deviation of the resulting plume speeds. The resulting mean slope and error of the linear regression were 0.95 ± 0.05 ms −1 for the vertical displacement. For comparison, the horizontal velocity component from the fit was −1.38 ± 0.24 ms −1 (Figure 7b) and agrees with the average wind speed of 1.4 m s −1 , suggesting that wind momentum dominated the horizontal plume motion, which is consistent with our own observation (Section 3). Note that the wind speed in Table 2 was not measured inside the crater and therefore not expected to be equal. Using Equation (2), the resulting CO 2 flux was computed as 8.3 ± 1.4 kg s −1 (717 ± 121 t day −1 ).

CO2 Flux from Bayesian Inversion
A constant prior pdf was assumed, ranging from 0 to a certain maximum value, in this case 30 kg s −1 , being over six times above the expected flux, so that biasing of the inversion would be avoided. Several tests with changing parameters were performed until reasonable convergence of the MCMC process was achieved. Given the meteorological data (Table 1), no fluxes satisfactorily explained the observed concentrations, that is, the inverted fluxes were orders of magnitude above the anticipated flux magnitude of Section 4.1. This changed once the modeled wind direction was altered, keeping all other parameters fixed. A very good match was achieved assuming a wind direction of 130°. The corresponding inversion result is shown in Figure 8.

CO 2 Flux from Bayesian Inversion
A constant prior pdf was assumed, ranging from 0 to a certain maximum value, in this case 30 kg s −1 , being over six times above the expected flux, so that biasing of the inversion would be avoided. Several tests with changing parameters were performed until reasonable convergence of the MCMC process was achieved. Given the meteorological data (Table 1), no fluxes satisfactorily explained the observed concentrations, that is, the inverted fluxes were orders of magnitude above the anticipated flux magnitude of Section 4.1. This changed once the modeled wind direction was altered, keeping all other parameters fixed. A very good match was achieved assuming a wind direction of 130 • . The corresponding inversion result is shown in Figure 8.

CO 2 Flux from Direct Integration
Thanks to continuous monitoring by the Vesuvius Observatory of Istituto Nazionale di Geofisica e Vulcanologia (INGV), the overall CO 2 flux of the Solfatara crater is well documented [40]. Dedicated measurements of the flux from the main vents only, however, are rare. In late 2012, using a MultiGAS instrument, a combined flux of 2.5 ± 0.7 kg s −1 for BG and BN has been measured [10]. Around the same time, using a path-averaging instrument, the authors of [11] measured a CO 2 flux of 3.5 ± 0.2 kg s −1 (306 ± 20 t day −1 ), which is higher, likely because their flux retrieval method was integrating over all mixing ratios within an extended area around the main vents. This method is therefore more similar to the method used here than the MultiGAS approach. It also illustrates the uncertainty associated with isolating a vented emission within an area of diffuse degassing and degassing from several minor vents under constant dispersion. In other words, CO 2 concentrations at the main vents are likely always influenced by influx from the immediate surroundings, regardless of the flux measurement tool [41]. In fact, trying to isolate the CO 2 concentrations at the main vent using the data in Figure 5b is challenging as it requires knowledge of the ambient CO 2 concentration in direct proximity upwind of the vents. This was attempted here, but it was found that the ambient concentration could be retrieved with a relatively large uncertainty only, which would have led to a large uncertainty of the flux. In the following, therefore, the path-integrated flux of 717 ± 121 t day −1 is compared with the figure obtained by [11]. Clearly, that figure is roughly half of the flux retrieved in the present work. A hint as to why that is is given when comparing the actual mixing ratios measured in 2012 and 2017. In 2012, plume-related values of 825 ppm on average were observed, whereas during this campaign plume-related values were 1300 ppm on average, which could be confirmed by simultaneous measurements around the plume area using the Li-COR analyzer. This corresponds to a~1.6-fold increase in the mixing ratio and a 1.7-to 2-fold increase in the background-corrected (in-plume) CO 2 mixing ratio (accounting for error margins in [11]). Assuming, to first order, that this factor applies to all in-plume mixing ratios at the main vent area entails an extrapolated flux for May 2017 between 497 and 662 t day −1 , in agreement with the flux retrieved here. In conclusion, a flux of 717 ± 121 t day −1 from the main vents is consistent with independent measurements of the flux of the BN/BG vent area.
Since 2012, the average soil CO 2 flux from Solfatara has increased by a factor of~1. 5 [40] in line with the factor of~1.6, which is arguably quite simplified but suggests that the CO 2 flux at the main vents has, between 2012 and 2017, increased at a similar rate (~1.6) to the medium flux (~1.5) across Solfatara.

CO 2 Flux from Bayesian Inversion
Given the actual prevailing wind direction of 210 • (SSW) during the measurement, the modeled plume and hence the modeled CO 2 concentrations (Equation (11)) are expected around the area marked by C (Figure 1b), where there is a steep slope (Figure 1d), which the forward model does not account for because for that location no CO 2 concentrations could be measured. So, the only way for the inversion to match the observed concentrations is to increase the CO 2 flux towards positive infinity. This is illustrated in Figure 9, which shows inverted fluxes for different wind directions. For wind directions between about 20 • and 180 • , fluxes are moving towards infinity. Once the wind blows towards areas away from the slope (between~50 • and~160 • ), i.e., when coming from between the NE and SE, the flux decreases until it eventually approaches the expected flux of around 8 kg s −1 . The best match is achieved for the two wind directions 90 • and 130 • , the latter having the highest confidence ( Figure 9). Clearly, this is only a sensitivity test and a different wind direction would have led to a concentration profile different from Figure 5b. However, the integrated CO 2 amount would have remained constant and so would have the flux. This suggests that the forward model is suitable for situations where the full model domain is covered by measured concentrations. While this may sound trivial, in practice, it is straightforward only for a fairly flat topography. Generally, the wind direction and the measurement geometry have to be consistent, which, due to the topography, was not the case here.
The best estimate of the flux (Figure 8) is moreover quite sensitive to the position of the vent as another sensitivity test shows. We recall that while the code solves for the flux of a point source of known position, here an effective vent position was considered, which was assumed to be that of BN as this yielded a flux value compatible with the value from direct retrieval of Section 4.1. Figure 9 shows the fluxes versus wind direction for an effective vent location at the mean coordinates of BN and BG (halfway between the two). For all wind directions, the inverted fluxes have increased. In order for the modeled concentrations to match the concentrations measured by LARRS, a vent located at the new effective location has to produce a higher flux because this effective vent location is further upwind than BN (Figure 1b). Realistic fluxes are still retrieved, namely for wind directions between 80 • and 130 • , yet their mean values, shown in Figure 9, are all higher than 8 kg s −1 . When moving the effective vent location even further away from BN, to BG, which is only 24 m eastwards from BN (~12% of the model dimension), the inverted fluxes eventually remain above an unrealistic 15 kg s −1 , for all wind directions (not shown in Figure 9). Note that the residuals, that is, the difference between modeled and measured CO 2 concentrations, are equal for all three scenarios.
Although the OF result suggests that the wind speed of 1.4 m s −1 used for the inversion is consistent, it is still interesting to consider an inversion with a different wind speed of say 1 m s −1 . Then, the resulting fluxes would be smaller by~20%, for both scenarios, as qualitatively expected (Equation (11)). Moreover, the temperature measured~1 km away could differ from the actual temperature in the Solfatara crater. Given a wind direction of 130 • , a change in temperature by 1 K would yield a change in inverted flux by only about 5%. A change in atmospheric pressure by 2 hPa (the largest deviation measured at Pozzuoli that day) would entail a change of the inverted fluxes by less than 5%.
These tests indicate the vent location to be the chief source of error. Thus, the vent location should be reasonably well known, or if this is not possible, it should be inverted for in parallel (e.g., [42]). The envelope depicts , the lower and upper 95% credibility interval (dotted lines in the histogram in Figure 8). The blue line depicts the fluxes for an anticipated vent location at BN, which matches the flux from direct integration (dotted line) for angles between ~80° and ~140°. The green line depicts the result of a second such test for an effective vent location halfway between BN and BG with increasing trend, indicated by the arrow, not matching fluxes from direct integration.

Conclusions
A remote sensing spectrometer has been used to remotely probe the fumaroles of Solfatara volcano in Italy, yielding path-averaged CO2 concentrations (concentration profiles). The combined flux from the two main vents Bocca Nuova and Bocca Grande has been computed using plume speed from optical flow analysis and integration of the concentration profiles, yielding 717 ± 121 t day −1 , in line with independent measurements. A Bayesian inversion based on a Gaussian plume model yielded a very similar value, however, only when the wind direction was consistent with the measurement geometry. For the present case, this means the area downwind the plume has to be covered by measurements. While for the given wind direction and terrain this was not possible, for many situations it is.
The inversion results indicate that precise (few m) knowledge of the vent location is important to arrive at accurate fluxes. As in many degassing sites where vent locations are unknown, the vent location could be inverted for in parallel to the corresponding fluxes.
Not only does the consistency of the result solidify confidence in the direct integration method for flux retrieval, including the optical flow analysis applied here, but it is also promising for a Figure 9. Sensitivity test results showing the inverted CO 2 flux for different presumed wind directions. Fluxes for wind directions between 220 • and 360 • are positive infinitive and not shown. The envelope depicts , the lower and upper 95% credibility interval (dotted lines in the histogram in Figure 8). The blue line depicts the fluxes for an anticipated vent location at BN, which matches the flux from direct integration (dotted line) for angles between~80 • and~140 • . The green line depicts the result of a second such test for an effective vent location halfway between BN and BG with increasing trend, indicated by the arrow, not matching fluxes from direct integration.

Conclusions
A remote sensing spectrometer has been used to remotely probe the fumaroles of Solfatara volcano in Italy, yielding path-averaged CO 2 concentrations (concentration profiles). The combined flux from the two main vents Bocca Nuova and Bocca Grande has been computed using plume speed from optical flow analysis and integration of the concentration profiles, yielding 717 ± 121 t day −1 , in line with independent measurements. A Bayesian inversion based on a Gaussian plume model yielded a very similar value, however, only when the wind direction was consistent with the measurement geometry. For the present case, this means the area downwind the plume has to be covered by measurements. While for the given wind direction and terrain this was not possible, for many situations it is.
The inversion results indicate that precise (few m) knowledge of the vent location is important to arrive at accurate fluxes. As in many degassing sites where vent locations are unknown, the vent location could be inverted for in parallel to the corresponding fluxes.
Not only does the consistency of the result solidify confidence in the direct integration method for flux retrieval, including the optical flow analysis applied here, but it is also promising for a possible future application of Bayesian inversion to retrieve volcanic CO 2 fluxes.
In conclusion, depending on the terrain, a Bayesian inversion may be a supplement or even an alternative retrieval method of gas fluxes, in particular for situations where plume speed estimation is impaired or impossible, for instance, for invisible plumes. A Bayesian inversion is relatively simple to implement, which suggests that both direct and inverse retrieval could be used side by side to improve confidence of CO 2 flux measurements from volcanoes and gas flux measurements in general, such as greenhouse gas remote sensing or verification of gas fluxes from satellite remote sensing.