1 2-D tomography of volcanic CO 2 from scanning hard target differential absorption LIDAR : The case of Solfatara , Campi Flegrei ( Italy )

2-D tomography of volcanic CO2 from scanning hard target differential absorption LIDAR: The case of Solfatara, Campi Flegrei (Italy) Manuel Queiβer, Domenico Granieri, Mike Burton School of Earth, Atmospheric and Environmental Sciences, University of Manchester, Oxford Road, Manchester M13 9PL, 5 UK Istituto Nazionale di Geofisica e Vulcanologia (INGV), Sezione di Pisa, 50126 Pisa, Italy Correspondence to: Manuel Queiβer (manuel.queisser@manchester.ac.uk), Tel.: +44(0)161 2750778, Fax.: +44(0)161 306 9361


Introduction
Subaerial volcanoes emit a variety of gaseous species, dominated by water vapor and CO 2 , and aerosols.Originating from exsolution processes that may take place deep in the crust due to the low solubility of CO 2 in magmas, volcanic CO 2 is a powerful tracer for magmatic recharge and ascent processes (Burton et al., 2013;Frezzotti and Touret, 2014;Chiodini et al., 2015;La Spina et al., 2015).Measuring volcanic CO 2 emission rates is therefore also a feasible pathway towards improved forecasting of volcanic activity, such as seismicity or eruptions (Petrazzuoli et al., 1999;Carapezza et al., 2004;Aiuppa et al., 2011).Unfortunately, magmatic CO 2 is released not only actively via vents such as the volcano mouth but also diffusively via soil or flank degassing (Baubron et al., 1991;Hards, 2005;Chiodini et al., 2007).In addition, in most cases the volcanic CO 2 signal is modest compared with ambient concentrations (Burton et al., 2013) and quickly diluted into the atmosphere.A common approach to determine the magmatic CO 2 flux is based on a gridded sampling of the CO 2 distribution in the volcanic plume itself (Gerlach et al., 1997;Lewicki et al., 2005;Diaz et al., 2010;Lee et al., 2016) from which 2-D CO 2 concentration maps are retrieved by secondary data processing, such as statistical methods (Lewicki et al., 2005;McGee et al., 2008) and dispersion modeling (Aiuppa et al., 2013;Granieri et al., 2014).Integrating the CO 2 concentrations over the cross-sectional plume area and multiplying the result with the transport speed perpendicular to the cross section yields CO 2 fluxes.The in situ method has two drawbacks.Firstly,  (427495,4519949), P4 (427507, 4519967) and P5 (427520, 4519986).Also shown are the respective range vectors (rays) for all five scans and the numbered locations of the LI-COR measurements.(d) Photo taken during the scan at P5 looking towards east.The largest clouds of condensed water aerosol appeared near the main vents (Bocca Nuova, BN, and Bocca Grande, BG) on the left.The CO 2 DIAL, visible in the lower right corner, comprised of the tripod carrying the telescope (with transmitter unit) and the main unit (red box).
it may be dangerous to perform in situ measurements from within the volcanic plume (e.g., due to toxic gases or low visibility near the crater mouth).Secondly, in situ methods allow for a very accurate estimation of CO 2 concentration, but only in the vicinity of the measurement point, potentially missing significant contributions from in between the measurement points.
Remote sensing techniques (see Platt et al., 2015, for overview of state of the art), notably active remote sensing platforms, including differential absorption lidar (DIAL) and spectrometers (Menzies and Chahine, 1974;Weibring et al., 1998;Koch et al., 2004;Kameyama et al., 2009), acquire columns of range-resolved (Sakaizawa et al., 2009;Aiuppa et al., 2015) or column-averaged (Amediek et al., 2008;Kameyama et al., 2009) CO 2 concentrations.They provide a powerful tool to overcome the aforementioned drawbacks of in situ measurement techniques by offering a faster, safer and more comprehensive acquisition (spatial coverage yields inclusive CO 2 concentration profiles).Moreover, there is no need for receivers or retroreflectors at the opposite end of the measurement column, which not only increases flexibility and timeliness of the acquisition but also is crucial for some measurements, including airborne or spaceborne acquisitions.
Active remote sensing platforms based on hard-target DIAL (topographic target DIAL) using continuous wave lasers can be made compact, rugged and portable, which is desirable for platform-independent measurement of atmospheric CO 2 , be it ground based or airborne (Sakaizawa et al., 2013;Queißer et al., 2015a).However, the drawback compared with "traditional", pulsed DIAL is that no rangeresolved CO 2 concentrations are measured, but column densities (in m −2 ) or, as in this work, path length concentration products are (called "path amount" hereafter, in ppm m).Hence, during a scanning measurement, no 2-D concentration maps are obtained, but 1-D profiles of path amounts, that is, values of path amounts versus scanning angle.To obtain gas fluxes 1-D profiles suffice.However, there are good reasons to obtain gas distribution information (e.g., to localize flux hot spots), which is particularly important for heterogeneous CO 2 emissions, including sites with both diffuse and vented degassing.It would be very desirable, and this was the motivation of this work, to use the 1-D profiles of path amounts to reconstruct a 2-D map that contains the geometry of the anomalous CO 2 release, leaving aside precise CO 2 mixing ratios.Note that tomographic reconstructions of volcanic gas plumes have already been performed, however, for SO 2 and using passive remote sensing techniques (Kazahaya et al., 2008;Wright et al., 2008;Johansson et al., 2009).
The study was focusing on a zone of mixed diffuse and vented degassing of magmatic CO 2 within the Solfatara crater (Italy) reported previously (e.g., Bagnato et al., 2014).Solfatara is a fumarolic field and part of the active volcanic area of Campi Flegrei (CF; Fig. 1).CF is a nested caldera erbium doped fiber amplifier.The beam pick is a glass wedge and reflects part of the transmitted light back into an integrating sphere (depicted as diffuser).ADC: analog-to-digital converter; DAC: digital-to-analog converter.The CO 2 cell is used to calibrate the seed laser wavelengths.To minimize hard-target and turbulence-related speckle noise the collimator used had a relatively high divergence of 1.7 mrad, while the telescope field of view was 1.5 mrad.For mechanical reasons the optical band pass filter was mounted before the collimating lens.The change in transmission spectrum can be neglected.
resulting from two large collapses, the last one ∼ 15 ka ago (Scarpati et al., 1993).CF is in direct vicinity to the metropolis of Naples and thus a direct threat to millions of residents.Thanks to its accessibility and strong CO 2 degassing Solfatara provides almost a model like volcano, a natural laboratory, to test new sensing approaches.However, it is part of one of the most dangerous volcanic zones in the world, showing ground uplift coupled with seismic activity with magma degassing likely having a significant role in triggering unrest (Chiodini at al., 2010).Solfatara therefore merits particular monitoring efforts and any new results on observables, may they stem from well-tried or new methods, are of direct importance to understand the fate of this active volcanic system.

Measuring 1-D profiles of CO 2 path amounts
The CO 2 DIAL (Fig. 2) is an active remote sensing platform based on the differential absorption lidar principle (Koch et al., 2004;Amediek et al., 2008).It is a further development of the portable instrument described in Queißer et al. (2015a, b).Specifics about the instrument and how CO 2 path amounts were retrieved are explained in detail in Queißer et al. (2016).Here, only a brief overview is given.By taking the ratio of the optical powers associated with the received signals for the wavelengths coinciding with an absorption line of CO 2 and the wavelength at the line edge, λ ON and λ OFF , respectively, one arrives at which is the differential optical depth τ , the key quantity to be measured.N CO 2 is the CO 2 number density R is the range, i.e., the distance between the instrument and the hard target, σ is the difference between the molecular absorption cross sections of CO 2 associated with λ ON and λ OFF , P (λ) is the received ("science") and P (λ) ref is the transmitted optical power ("reference").The latter is measured as a reference to normalize fluctuations of the transmitted power.The normalized optical power in Eq. ( 1) is referred to as grand ratio (GR): The two distributed feedback fiber seed lasers emit at λ ON = 1572.992nm and λ OFF = 1573.173nm (Rothman et al., 2013).To be able to easily reject background noise (such as solar background), lock-in detection is used.Consequently, both seed laser beams (for λ ON and λ OFF ) are amplitude modulated using two LiNbO 3 electro-optical modulators (EOM) at slightly different sine tones near 5 kHz and simultaneously amplified by an erbium doped fiber amplifier (EDFA) before being simultaneously transmitted.The transmitted optical power can be adjusted between ∼ 80 mW and a maximum of 1.5 W. A glass wedge scatters a fraction of the transmitted light into an integrating sphere where the reference detector is mounted, which measures P (λ ON ) and P (λ OFF ).The transmitted light is diffusively backscattered by a hard target, which can be any surface located up to ∼ 2000 m away from the instrument, and is received by a 200 mm diameter Schmidt-Cassegrain telescope with a focal length of 1950 mm and field of view of ∼ 1.5 mrad.Typically the received optical power is a couple of nW at a bandwidth integrated noise of ∼ 1 pW (root mean squared noise equivalent power).The analog-to-digital converter (ADC) operates at 250 kSamples s −1 and has a resolution of 16 bit.The integration time per scan angle was set to 4000 EOM modulation periods, which corresponds to data chunks of length 784 ms (integration time) for both science and reference channel.Each of these four chunks of data is demodulated using a digital lock-in routine following Dobler et al. (2013).After the lock-in operation one arrives at four DC signals, associated with the optical powers P (λ ON ), P (λ OFF ), P (λ ON ) ref and P (λ OFF ) ref .τ is calculated using the right-hand side of Eq. ( 1), after taking the mean of each of the four signals.However, to retrieve CO 2 path amounts the left-hand side of Eq. ( 1) was not used.Instead, to account for the instrumental offset of τ , prior to scanning the volcanic plume, values of τ were acquired for different R in the ambient atmosphere (Queißer et al., 2016).The points were used to fit a calibration curve.The ordinate at R = 0 gave the instrumental offset.The calibration curve was used to convert the measured in-plume τ to CO 2 path amounts Y CO 2 (in ppm m).Column-averaged CO 2 mixing ratios X CO 2 , av (in ppm) were obtained by dividing path amounts by the hard-target range R. Note that X CO 2 , av was mainly retrieved for display purposes rather than for data processing.The range was measured by an onboard range finder (DLEM, Jenoptik, Germany), based on a 1550 nm lidar with pulse energy of 500 µJ and accuracy < 1 m.By pivoting the receiver-transmitter unit using a step motor, values for Y CO 2 per heading were attained and hence 1-D profiles.
The precision of the column-averaged CO 2 mixing ratio was evaluated as with the signal-to-noise-ratio (S/N) of where GR and σ GR are the mean and standard deviation of the grand ratio, respectively.They were estimated from time series acquired at fixed angles in between the scans at CF.The S/N accounts for all noise sources occurring during acquisition, including instrumental noise, non-stationary baseline drift, solar background noise, atmospheric noise (mostly air turbulence) and perturbation by aerosol scattering (e.g., condensed water vapor).The second term depicts the relative range uncertainty (standard deviation of ranges σ R over mean of ranges R ), which is typically ∼ 1 m.The relative uncertainty due to hard-target speckle was estimated as (MacKerrow et al., 1997) where D is the spot diameter on the hard target (in m) and ξ is the dimension of the telescope field of view (in m) on the hard target.Typical values were 25 and 42 cm, respectively.

Reconstructing a 2-D map of CO 2 mixing ratios
The goal is to obtain CO 2 mixing ratios (X CO 2 , in ppm) at a given point (x, y).The region of interest (area bounding the scans) was divided into grid cells with length x (in x direction) and y (in y direction).Within a given grid cell the associated mixing ratio is uniform.X CO 2 was inferred from the measured Y CO 2 using an inverse technique following Pedone et al. (2014a).The CO 2 path amount is associated with the product of a range segment and the (uniform) CO 2 mixing ratio X CO 2 along that range segment.For a given beam path and heading angle (hereafter also referred to as ray) and for p grid cells traversed by the ray the corresponding governing equation can be written as where r i depicts the length of the ray segment in grid cell i p i=1 r i = R .X CO 2 ,i is the (unknown) CO 2 mixing ratio within grid cell i (in ppm).Including all rays, one arrives at a system of linear equations, which can be written as where L is a m × n matrix, called geometry matrix.n depicts the number of model grid cells.The elements of L are the length of all m paths for n grid cells.For instance, the element l m=21,n=6 depicts the length of ray 21 in grid cell number 6.To be able to apply Eq. ( 7) to the measured Y CO 2 , the associated rays were mapped onto (x,y) coordinates using x j,k = r j,k cos ϕ j and y j,k = r j,k sin ϕ j .ϕ j is the cumulative heading angle of the j th ray determined from the scanning angular velocity and the time interval between the rays retrieved from the time stamps of the data.For a given ray j , r j,k is the kth path length increment, such that N is chosen such that adjacent x j and y j were no more than ∼ 1 m apart.Since the ϕ i are not absolute angles, but rather relative to the first ray, the x j,k and y j,k were rotated to match the world coordinate system, i.e., rotated such that the first ray's (x, y) coordinates match the real world ray trajectory.
Table 1.Inversion data input file showing (x, y) coordinates with associated measured Y CO 2 and X CO 2 , av = Y CO 2 /R j of the first ray (j = 1) and x j =1,k , y j =1,k for the first five and last four path length increments.N = 100.This was performed using the well-known rotation matrix with a rotation angle, which is the angle between the first ray and the x axis.The Y CO 2 versus the resulting (x, y) coordinates formed the inversion input file (see Table 1 for an example).From the data in that input file an algorithm built the geometry matrix L by using the (x, y) coordinates and computing the length of each ray in each model grid.c is a n × 1 matrix containing uniform X CO 2 per grid cell and is the desired quantity to be inverted.a is a m × 1 matrix containing the measured (observed) Y CO 2 for each ray, retrieved from the inversion input file.For simplicity, n x = n y , where n x and n y are the number of grid cells in x and y direction, respectively.Thus, n = n 2 x .Since m > n the system is overdetermined.This means L in Eq. ( 7) cannot be inverted to arrive at c, but c has to be approximated.This is done here using a least square fitting procedure.To solve Eq. ( 7) for c a least square solver, the MATLAB LSQR routine was used.The algorithm iteratively seeks values for c, which minimize the normalized misfit a − Lc .Therefore, c represents a model with a maximized likelihood of explaining the observed data a.By reshaping c into the measurement 2-D grid a 2-D map was obtained.Ordinary Kriging interpolation (Oliver, 1996) may then be applied to obtain a higher resolution estimate of the map.

CO 2 flux retrieval
From the inverted 2-D map of X CO 2 the CO 2 flux (in kg s −1 ) was computed as where X CO 2 ,pl are the inverted, background-corrected CO 2 mixing ratios computed as where X CO 2 ,bg = 380 ppm is the background CO 2 mixing ratio at Solfatara measured in situ from an air sample, which was collected far from the degassing area, in accordance with Chiodini et al. (2011).The global average CO 2 level is now ∼ 400 ppm but this level may be quite different from one location to another for various reasons (latitude, weather conditions, seasonal effects, anthropogenic contributions).Inside and around the Solfatara there are many green areas (e.g., the nearby crater of Astroni), where the photosynthetic activity by plants likely causes the decrease in CO 2 levels compared to the more urbanized surroundings.u is the magnitude of the component of the plume transport speed perpendicular to the scanned cross section (in m s −1 ).Since the scans were performed along the horizontal plane u refers to the vertical component of the plume transport speed.N air is the number density of air (in m −3 ), computed using meteorological data (pressure, temperature, humidity) acquired by a portable meteorological station close to the instrument.M CO 2 is the molar mass of CO 2 (in kg mol −1 ) and N A is Avogadro's constant (in mol −1 ).
The plume transport speed was evaluated from digital video footage acquired during the measurement (Queißer et al., 2016), employing a video analysis program (Tracker from Open Source Physics).Condensed water vapor aerosol emitted by various vents in the region of interest was assumed to propagate with the same velocity as the volcanic CO 2 .At a given video frame a pixel was fixed and the calibrated propagated distance (in pixels) was measured as the video proceeded.Since the frame rate of the video was known (30 frames per second), the speed by which the tracked point and hence a parcel of gas was transported could be estimated.
The relative error of the CO 2 flux was estimated as where u is the absolute uncertainty of the plume speed, X CO 2 ,pl is the absolute error of the CO 2 mixing ratio (measurement precision) at a given point within the integrated area and σ X CO 2 is an estimate of the error of the inverted model.Figure 3 summarizes all main steps that were involved in the methods to retrieve the flux.

Results
The experiment took place on 4 March 2016 inside the crater of Solfatara (Fig. 1) and focused on the diffuse CO 2 release alongside the Solfatara crater edge, located south of the main vents Bocca Nuova (BN) and Bocca Grande (BG), although they were included in the scans.Elevated CO 2 mixing ratios, up to 1500 ppm at places, could be affirmed by It is assumed that during the complete acquisition the CO 2 distribution did not change ("frozen plume").For each scan and for each heading differential optical depths τ have been retrieved and converted into Y CO 2 (and X CO 2 , av ), as detailed in the method section.The resulting 1-D profiles of CO 2 path amounts are shown in Fig. 3. Numerous wiggles indicate widespread and heterogeneous degassing activity, suggesting diffuse degassing or CO 2 advected by local wind eddies.In addition, there are symmetric features, such as around 26 • in Fig. 3a, which appeared in scans carried out prior to the experiment and the day before, thus suggesting vented degassing activity.The angular scanning velocity was 2.1 mrad s −1 , associated with an angular resolution of 1.65 mrad, which corresponds to a lateral resolution of around 24 cm between points in Fig. 3.
Before inverting X CO 2 , the Y CO 2 and the associated rays were mapped onto a (x, y) grid as detailed in the method section.The coordinate system was chosen such that the instrument positions of all five scans were located on the y axis (Fig. 1c).Table 1 shows an excerpt of the result of this procedure.Figure 5 shows the plot of X CO 2 , av versus x j,k and y j,k (shown in Table 1).It represents a semi-quantitative map, indicating where high CO 2 concentrations are likely to be expected and thus contains useful a priori information for the inversion.
The LSQR algorithm was tested using a synthetic, realistic scenario.Synthetic data Y CO 2 were generated from a true model comprising of known X CO 2 at each grid point using Eq. ( 7) and the real geometry matrix L, which contained the actual instrument positions and measured ranges.X CO 2 of the true model started at 380 ppm at grid 1 and increased by 60 ppm per increase in grid number (Fig. 6a).The inversion result is shown in Fig. 6b and c.By running the inversion with varying number of grid cells and random arrangements of X CO 2 the viable number of grid cells was found to be between n = 4 up to 49 without considerable loss of capabil-ity to recover the true X CO 2 .Beyond that recovery became poorer.In addition, various inversions with random distribution of random X CO 2 (drawn from a normal distribution with mean and standard deviations from the data presented in Fig. 5) were carried out that all yielded good recovery of X CO 2 .For the real data, however, already for n > 16 the inversion yielded unreasonable high or low X CO 2 , which can be described as oscillating X CO 2 .The difference between the synthetic and the real data vector a is that the former exhibits a higher variability of Y CO 2 .This is illustrated in Fig. 7a, which shows the components of a for one of the synthetic model realizations and the real data.The real data show high frequency components.One issue is the model discretization.In reality the scanned area was far from consisting of 16 grids with homogeneous X CO 2 only, but instead was a complex distribution of an unknown number of CO 2 patches, which leads to that spiky appearance of the real data vector.Further on, the more the "frozen" plume assumption holds the more the "humps" in Fig. 7a should resemble each other, since the x axis on Fig. 7a and b scales with the time of acquisition.For the synthetic case, there are five somewhat similar, but distorted, humps (Fig. 7b) due to the different viewing angle of each scan and differences in angle spans.This is, however, partially also the case for the real data vector (Fig. 7a), more so for the first four scans and less so for the last one.The spikes, i.e., the high frequency components, however, differ a lot, suggesting short-lived fluctuations on a small spatial scale that do not reappear between subsequent scans.This high frequency "noise" is due to the fact that actual fluctuations of the X CO 2 occurred during the measurement.They happen on a faster timescale than the measurements of the five profiles (which took 142 min) and the scans themselves, violating our "frozen" plume assumption, at least for small patches of CO 2 .This leads to inconsistencies between the system of linear equations (there is no solution that satisfies all equations simultaneously) and hence problems to fit the spiky data.
Adding these high frequency components to the synthetic data by perturbing them using random perturbations drawn from normal distributed Y CO 2 (Fig. 7b) and performing inversions gave an outcome very similar to the real data inversion.The inversion struggled to recover the true model X CO 2 (Fig. 7c), yielding oscillations.To improve the real data inversion result synthetic inversions were performed after smoothing the perturbed synthetic data vector a (Fig. 7b).Low-pass filtering the data, however, removes peaks and therefore falsifies the data, which the iteration then tries to fit, yielding poorly recovered X CO 2 .Smoothing had therefore to be applied with care.As a result, moderately smoothing the data vector a using a running average with window size 14 (rays) yielded reasonable recovery of the true model (Fig. 7d).However, performing synthetic tests with over 500 random arrangements of the grid cells occasionally led to oscillations in the model, meaning that the true X CO 2 was over-and undershot again.It is important to note that this  1 for first ray), which yields column-averaged mixing ratios X CO 2 , av .Also shown are the instrument positions (squares on y axis) starting with P1 at y = 20 m.For display purposes the data have been regridded on a regular grid of 90×90 points using natural interpolation.One would expect high anomalous CO 2 mixing ratios near the main vents (BN and BG near x = 120 m, y = 140 m) and the southern part of the area.Low anomalous CO 2 mixing ratios are to be expected in the northwestern part.Note that due to the abundance of data some data points were masking each other.They were thus averaged, leading to a maximum mixing ratio lower than actually observed (e.g., in Fig. 4b).
occurred for grid cells at the model edges.This is the area where ray coverage is the poorest.It appeared that by increasing the number of grid cells n this issue worsened.The inverse problem is then still overdetermined in the sense of ray coverage (for instance, for n = 64 on average there are still 10 rays per grid cell).However, this is not the case for the edge grid cells, as some grid cells are not traversed by some rays, which is aggravated with increasing n.Another contribution to the inconsistency of the system of linear equations was given by rays that are not traversing the edge grid cells due to the topography of the hard target (slope).These two issues are illustrated in Fig. 7e.Moreover, albeit being a largely overdetermined (in fact mix-determined) problem, the rays are close to being parallel as the acquisition was onesided and the instrument positions were not widely spread out.At best, the relatively large number of rays with similar angles provides redundancy of information.At worst, instead of adding extra information the lack of angle diversity may have contributed to the inconsistency of the problem, making the solution more unstable, meaning that a small change in data caused large variations in the model.Increasing the model grid number n increased inconsistency.
A possible explanation for the limited spatial resolution for the real data inversion is therefore the following: the problem of recovering an unknown and dynamic distribution of model parameters, in this case X CO 2 , is inherently ill posed (Zhang et al., 2016), which implies non-uniqueness of the solution.This non-uniqueness is caused, for instance, by the unknown number of CO 2 patches or sources and their heterogenic and anisotropic distribution (Hamilton et al., 2014).In addition, a lack of observational data (mainly "frozen plume" assumption and angle diversity) caused inconsistency of the problem, which led to both unrealistic model parameters ("oscillations") and a further amplified non-uniqueness of the model.It is apparent that increasing the grid number of the model would aggravate inconsistency and thus model error and non-uniqueness.
A more promising measure therefore was to omit parts of the rays, which is effectively smoothing the data vector without falsifying it, so that inconsistencies in Eq. ( 7) are reduced.Best results were obtained by selecting only every third ray and the corresponding Y CO 2 as input to the inversion.
As for the synthetic tests, a constant X CO 2 , the mean of the raw data (Fig. 5), was used as a starting model for inverting the real data.As a consequence of the discussion above the maximum feasible number of grid cells was n = 16 for a robust inversion.The resulting grid length was x = 38 m and y = 33 m.The inverted model is shown in Fig. 8a.To increase spatial resolution the inverted model was interpolated onto a grid with grid spacing x/8 and y/8 using ordinary Kriging interpolation.The result is shown in Fig. 8b.The zones without ray coverage were excluded from the 2-D map.Overlaying the resulting 2-D map of CO 2 mixing ratios with the map of Solfatara reveals a zone of increased anomalous CO 2 degassing activity along the southeastern edge of Solfatara, which is in reasonable agreement with in situ data from the LI-COR CO 2 analyzer (Fig. 8c).
The resulting 2-D map of CO 2 mixing ratios was used to compute the CO 2 flux.The vertical component of the plume transport speed was estimated to be 1.1 ± 0.2 m s −1 .The plume speed uncertainty was retrieved from the standard deviation of various plume speeds retrieved from different tracks carried out across the plume (Queißer et al., 2016).Equation ( 10) was employed to estimate the flux uncertainty.Since the synthetic test inversions behave very similar to the real data inversion, to obtain a meaningful model uncertainty of the real data inversion 1000 synthetic inversions with random true models, drawn from a normal distribution using mean and standard deviation of the data presented in Fig. 5, were carried out.Computing the absolute difference between true and recovered X CO 2 at each model grid point yielded 16 × 1000 values with a mean of 373 ppm.This was adopted as a measure for the model standard deviation of the real data inversion in Eq. ( 10), assuming it was constant throughout the model.Following the measurement error, the model error gave the most dominant uncertainty, associated to relative errors of the inverted X CO 2 between 34 and 90 %.A constant X CO 2 ,pl = max( X CO 2 , av ), i.e., the maximum error of all five scans, was considered for the measurement error.Using Eq. ( 8) the resulting CO 2 flux was computed as 17.0 ± 9.2 kg s −1 (±1 SD) or 1465 ± 798 tons day −1 .6) should yield almost identical path amounts for both rays.However, this is not the case since, due to the topography, their lengths differ.This adds inconsistency to Eq. ( 7).If as a remedy the number of model grid cells would be increased, for instance such that the grid spacing is halved (dotted lines), this would lead to grid cells without ray coverage, as illustrated here, and hence erroneous model values.

Discussion
Retrieving CO 2 path amounts using calibration based on mixing ratios as done here is not optimal as the mixing ratios depend on temperature, pressure and humidity of the air.It was a solution that was adopted for practical reasons.Compared with the precision of the instrument the resulting bias was, however, low and could be neglected.Typically, the air temperature varied between 11 and 15 • C and relative humidity varied between 55 and 65 %.The air pressure was steady around 974 hPa.In a simplified estimation, this corresponds to a bias of ∼ 1.5 % of the total column-averaged mixing ratio, which would correspond to around 9 ppm, assuming typical measured path amounts and ranges.To put this into context, the error of the column-averaged mixing ratio was typically ∼ 110 ppm (Fig. 3) and is dominated by the instrument S/N.The associated increase in flux uncertainty would be less than 1 %.In the future, however, a better retrieval scheme will be adopted.This is important as we strive to increase in the measurement precision of the CO 2 DIAL to sense more subtle CO 2 emissions than at Solfatara and a bias related to uncertainty of the number density becomes relevant.Note that a bias of 20 ppm due to uncertainties in background mixing ratio would also entail a small decrease in measurement precision by ∼ 4 % and an increase in flux error by roughly 2 %.
The retrieved 2-D map (Fig. 8c) indicates an elongated zone of intense anomalous degassing along the eastern edge of the Solfatara crater.Encouragingly, this is a persistent feature in different inversions performed with different number of beam paths and smoothing parameters of input data and underpins that it is real.Previous measurements sampling the Solfatara area with accumulation chambers yielded an increased anomalous CO 2 degassing activity in the corresponding area too (Granieri et al., 2010;Tassi et al., 2013;Bagnato et al., 2014).The retrieved elongated zone of anomalous CO 2 degassing likely encompasses at least two major vents (Fig. 8c).The locations of the peaks in CO 2 mixing ratio in Fig. 8c fairly agree with the 1-D input data.For instance, the peak near the center of the crater corresponds to the peak www.atmos-meas-tech.net/9/5721/2016/Atmos.Meas.Tech., 9, 5721-5734, 2016 To project the inverted model onto the map the function "overlay image" of Google Earth software has been used using altitude 0 relative to instrument altitude above ground.Also shown are X CO 2 from in situ measurements (measurement points 3 to 10) using the LI-COR CO 2 analyzer.Note that the in situ values had been acquired a day before the scans and thus serve as an approximate reference only.Moreover, since measured column-integrated mixing ratios were below 1500 ppm, there is no X CO 2 of 1500 ppm expected in the inversion result.

26
• in the first scan in Fig. 3a.The second scan (Fig. 3b) indicates a rather abrupt decrease in X CO 2 , av at 28 • , in line with the edge of the zone of elevated CO 2 mixing ratios at the crater center (left part of Fig. 8c).This central degassing feature is coherent with results of recent campaigns (Granieri et al., 2010;Tassi et al., 2013;Bagnato et al., 2014).The increase in X CO 2 , av near 9 • in Fig. 3d matches the position of the local peak in X CO 2 between in situ points 7 and 8 in Fig. 6c.Provided sufficient rays, which was the case for the zones away from the edges of the 2-D map, disagreement between the peaks in the 1-D data (Fig. 3) and those in the 2-D map (Fig. 6) are likely due to physical fluctuations in CO 2 concentration ("frozen" plume approximation), as evinced in the result section.During acquisition one could visually identify at least five small vents emitting water vapor and therefore most likely also CO 2 .Though not recovered due to the limited spatial and temporal resolution of the inversion, this advocates that there are in fact separate vents south of the main vents, near the edge of the Solfatara crater.
Retrieved X CO 2 peak near 1300 ppm (2 m above ground), in line with the in situ LI-COR data, although not spatially matching them in places.This can be explained by the fact that the in situ values were acquired the day before so that local wind and thus dispersion patterns were different.Nevertheless, both the LI-COR in situ data and the inversion result indicate high X CO 2 near the main vents and along the crater edge.Near the main vents the highest CO 2 mixing ratios in the 2-D map are located ca. 17 m west of BN.In fact, the whole zone of high X CO 2 is shifted ∼ 17 m northwest from where one would expect it.Since the predominant wind direction at the time of acquisition was around 300 • , to first order one would expect the CO 2 to disperse rather towards southeast, along the crater edge.The main vent area was at the edge of the scanned area.Note that the relative inversion residual a−Lc a was 0.18, which means on average 18 % of X CO 2 ,pl is unexplained by the model in Fig. 8a.Note that due to the non-uniqueness of the problem a lower misfit does not imply that the inverted model is closer to the "true" model.For instance, increasing the grid number decreased the misfit but yielded models with XCO 2 much larger than the maximum or lower than the minimum in Fig. 4.Although being close to the expected distribution, the mismatch is therefore likely due to a non-unique realization of the model.The latter is related to poor ray coverage for grid 16 (northeastern most), which is partly corresponding to the zone containing the main vents, and the error resulting from the low model discretization (in reality X CO 2 was certainly not constant over 30 m).Note that the acquisition focused on the zone south of the main vents.Possibly, but less likely, CO 2 was advected slightly towards west due to dispersive mechanism related to local wind eddies decoupled from the main wind direction.These dispersive mechanisms take place in any case and make a distinction between CO 2 from the main vents and the surrounding diffuse degassing challenging.For that reason, in future acquisitions at that site the region of interest shall be scanned from instrument positions aligned along a half circle around the zone rather than using a "flat" scan geometry as chosen here.
For a comparison, CO 2 fluxes were computed directly from the 1-D profiles (Fig. 4), that is, similar to Eq. ( 8) but using path amounts, ignoring any heterogeneity in the CO 2 distribution.Given the uncertainty, the average flux of all five scans (1055 ± 389 tons day −1 ) is compatible with the result obtained from the 2-D map (1465 ± 798 tons day −1 ).
However, both the CO 2 flux from the 2-D map and from the 1-D profiles are higher than fluxes previously estimated.Pedone et al. (2014a) report a CO 2 flux of only ∼ 300 tons day −1 in early 2013, but it focused on the area around the Solfatara main vents, that is, 8000 m 2 .In this study the area considered for flux computation was over 21 000 m 2 .Moreover, the average degassing rate at Solfatara has been increasing by ∼ 9 % each year over the past 10 years or so (Chiodini et al., 2010;d'Auria, 2015).The magnitude of the retrieved flux of this work equals the total CO 2 flux of diffuse degassing structure (which encompasses the Solfatara crater, the Pisciarelli fumaroles and the area around Solfatara crater) reported by Granieri et al. (2003), however, 13 years prior to this study, when emission rates were 1.09 13 (or 3) times lower.A very rough extrapolation of the 1465 tons day −1 obtained here to the whole of Solfatara crater, that is, by assuming a constant flux per unit area (Queißer et al., 2016), would yield a flux between 2400 and 4700 tons day −1 .This makes 1465 tons day −1 a coherent figure.Extrapolating the flux for the main vent area of 300 tons day −1 from 2013 would yield a flux of 390 tons day −1 in early 2016.Integrating CO 2 mixing ratios of the area around the main vents only (bounded to the south by in situ point 6; Fig. 8c) yields a flux of 364 ± 206 tons day −1 of CO 2 , in excellent agreement with the extrapolated flux.Finally, we note that to our knowledge all former studies except one (Pedone et al., 2014a) inferred X CO 2 and hence CO 2 fluxes from a grid of point measurements.If the area associated with lower CO 2 concentration was larger than the area with very high CO 2 concentration, point sampling may have missed degassing activity in between the measurement points and so tended to yield lower flux values.
All five scans were performed one-sided, i.e., from a single half space, as often the case in geophysical tomography problems (e.g., Hobro et al., 2003).As evinced above, this lack of angle diversity appears to have contributed to inconsistency of the inverse problem, leading to unrealistic models, which aggravated with increasing model grid number (non-uniqueness).Non-uniqueness was, however, counterbalanced by the abundance of hard data available, which facilitated the rejection of unlikely models, increasing our confidence to have arrived at a model close to the "true" model.A coarse grid spacing (small n) was necessary to exclude non-unique solutions with high confidence.As a consequence the inverted model is missing small-scale features.Nevertheless, given the fair agreement with the hard data, the inverted 2-D model (Fig. 6c) is quantitatively sound and outlines the geometry of the diffuse degassing probed at Solfatara.This case therefore enabled us to demonstrate that one may obtain useful tomographic results from one-sided scanning of a degassing feature.We expect an enlarged angle diversity to be the most feasible means to enhance consistency of the problem as studies using tunable diode lasers have shown (Pedone et al., 2014a, b) and potentially this will increase the maximum number of grids usable for stable inversion, boosting 2-D resolution.The goal of this study was not to design and demonstrate a universal tomography method for generic gas emission scenarios, but instead to extract, with the data at our disposal, useful information for the case of Solfatara.Although the results suggest this method to be useful, it is not clear whether it will be applicable to other degassing manifestations with few or no a priori constraints.This is particularly true for cases with an anisotropic or heterogenic distribution and unknown number of gas emitters that pose a highly ill-posed inverse problem.Future measurements of this type are envisaged, including sites with wellconstrained degassing characteristics, using a wider variety of viewing angles, which will allow a better assessment as to which extent this method is useful for one-sided tomography of highly non-isotropic CO 2 plumes in nature.This will help to assess the general applicability of this method, that is, increase confidence in the method, for sites where, unlike for the present case, few or no a priori constraints are available.Another point of improvement, especially for ill-posed problems, is related to finding meaningful data covariance matrices that could be used to weight (regularize) the fitting procedure (Snieder and Trampert, 1999)

Conclusions
As magmatic CO 2 degassing rates are tracers for the dynamics and chemistry of the magma plumbing system beneath Campi Flegrei and at volcanic areas in general, a comprehensive quantification of magmatic CO 2 degassing strength is of interest for volcanology and of vital importance for civil protection.
Scanning hard-target DIAL measurements has been performed at Solfatara crater (Campi Flegrei, Italy), which allowed an inclusive measurement of CO 2 amounts in the form of 1-D profiles of CO 2 path amounts.From the 1-D profiles a 2-D map of CO 2 mixing ratios X CO 2 has been reconstructed, outlining the main CO 2 distribution.Primarily the one-sided acquisition geometry of the input 1-D profiles and assumptions about the emission dynamics ("frozen" plume) increased inconsistency of the inverse problem and limited the feasible spatial resolution of the model.
However, the resulting map of X CO 2 was in satisfactory agreement with independent measurements of CO 2 mixing ratio.Therefore, the 2-D map was deemed to be useful to retrieve the CO 2 flux.The fluxes retrieved using that map is compatible with previous results.The 1-D profiles have been acquired from a single half space (one-sided), which indicates this tomography method to be potentially beneficial for scanning strongly non-isotropic CO 2 distributions, such as diffuse emissions, that can be viewed from limited angles only, sites difficult to access (e.g., volcanic craters), or airborne acquisitions.The authors expect a more adapted acquisition geometry to enhance resolution of the 2-D map, which would allow a more accurate gas flux estimation.However, more tests are necessary to assess general applicability and increase confidence in this method.

Data availability
The data as well as the MATLAB script used for the tomographic inversion are available upon request from the Research Data Repository of the University of Manchester (UK).

Figure 1 .
Figure 1.Geography and measurement geometry.(a) Location of the Solfatara crater as part of the volcanic area of Campi Flegrei, near Naples (Italy).(b) Nadir photo of Solfatara crater.The rectangle contains the region of interest.(c) Zoom of area outlined by the rectangle depicting the five instrument positions P1 to P5 with the following UTM coordinates: P1 (427476, 4519921), P2 (427485, 4519935), P3 (427495, 4519949), P4 (427507, 4519967) and P5 (427520, 4519986).Also shown are the respective range vectors (rays) for all five scans and the numbered locations of the LI-COR measurements.(d) Photo taken during the scan at P5 looking towards east.The largest clouds of condensed water aerosol appeared near the main vents (Bocca Nuova, BN, and Bocca Grande, BG) on the left.The CO 2 DIAL, visible in the lower right corner, comprised of the tripod carrying the telescope (with transmitter unit) and the main unit (red box).

Figure 2 .
Figure2.Scheme of the CO 2 DIAL as used for this experiment.EOM: electro-optical modulator; DLEM: range finder module; EDFA: erbium doped fiber amplifier.The beam pick is a glass wedge and reflects part of the transmitted light back into an integrating sphere (depicted as diffuser).ADC: analog-to-digital converter; DAC: digital-to-analog converter.The CO 2 cell is used to calibrate the seed laser wavelengths.To minimize hard-target and turbulence-related speckle noise the collimator used had a relatively high divergence of 1.7 mrad, while the telescope field of view was 1.5 mrad.For mechanical reasons the optical band pass filter was mounted before the collimating lens.The change in transmission spectrum can be neglected.

Figure 3 .
Figure 3. Scheme summarizing the main data processing steps involved in this paper.

Figure 5 .
Figure5.Contour plot of Y CO 2 after mapping them onto (x, y) coordinates (which is the input to the inversion), divided by the range for all 627 beam paths (same as shown in Table1for first ray), which yields column-averaged mixing ratios X CO 2 , av .Also shown are the instrument positions (squares on y axis) starting with P1 at y = 20 m.For display purposes the data have been regridded on a regular grid of 90×90 points using natural interpolation.One would expect high anomalous CO 2 mixing ratios near the main vents (BN and BG near x = 120 m, y = 140 m) and the southern part of the area.Low anomalous CO 2 mixing ratios are to be expected in the northwestern part.Note that due to the abundance of data some data points were masking each other.They were thus averaged, leading to a maximum mixing ratio lower than actually observed (e.g., in Fig.4b).

Figure 6 .
Figure 6.Synthetic inversion result with n = 16 grid cells.(a) True model used to generate synthetic column-averaged Y CO 2 .Each grid cell is identified by a grid number.The dotted line outlines the ray coverage.The instrument positions are indicated.(b) Inverted model.(c) True and inverted X CO 2 versus grid cell.The inverted X CO 2 for grid 13 is off since the ray coverage associated with that area was poor.

Figure 7 .
Figure 7. Synthetic inversion tests with perturbed data.(a) Data vector a (Eq.7), showing the real and synthetic data (used for inversion result in Fig. 6) versus ray number.(b) Synthetic data vector, perturbed synthetic data vector and smoothed perturbed synthetic data vector.(c) Inversion result using the perturbed synthetic data from (b).(d) Inversion result using the smoothed perturbed synthetic data.As the true model values were drawn from a random distribution, the model arrangement differs between (c) and (d).(e) Scheme of a single grid cell at the model edge facing the instrument position with two subsequent rays that have traversed the same grid cells before reaching the displayed grid.As they are in the same grid and nearly parallel, Eq. (6) should yield almost identical path amounts for both rays.However, this is not the case since, due to the topography, their lengths differ.This adds inconsistency to Eq. (7).If as a remedy the number of model grid cells would be increased, for instance such that the grid spacing is halved (dotted lines), this would lead to grid cells without ray coverage, as illustrated here, and hence erroneous model values.

Figure 8 .
Figure 8. Retrieved 2-D model of X CO 2 .(a) Inverted model of X CO 2 .(b) Inverted X CO 2 after ordinary Kriging interpolation (interpolated grid size is 3.75 by 4.38 m).The ray coverage is depicted by the dotted line.The straight lines outline the grid.(c) X CO 2 superposed onto nadir photo of Solfatara for those areas covered by the rays.To project the inverted model onto the map the function "overlay image" of Google Earth software has been used using altitude 0 relative to instrument altitude above ground.Also shown are X CO 2 from in situ measurements (measurement points 3 to 10) using the LI-COR CO 2 analyzer.Note that the in situ values had been acquired a day before the scans and thus serve as an approximate reference only.Moreover, since measured column-integrated mixing ratios were below 1500 ppm, there is no X CO 2 of 1500 ppm expected in the inversion result.