Simulation, image generation, and tomographic reconstruction of idealized volcanic plumes

ABSTRACT This paper presents a volcanic plume simulation and image generation framework, alongside a method for the tomographic reconstruction of volcanic ash plumes. The simulation framework facilitates the generation and processing of imagery analogous to that produced by real-world multi-spectral infrared observations of volcanic emissions. With this required imagery simulated, methods for the tomographic reconstruction of volcanic plumes can be tested. This simulation approach was undertaken due to the lack of suitable real-world multi-spectral and multi-angle imagery, and the difficulties and dangers of arranging multiple ground-based cameras or operating UASs in proximity to an active volcano. The efficacy of the simulation framework is demonstrated through a series of sensitivity analyses, assessing the change in reconstruction accuracy when modifying simulation variables such as the number and distribution of images, spatial resolution, and camera pointing inaccuracy. It is proposed that the results of these analyses can be used to inform the design and optimization of real-world observation campaigns.


Introduction
Volcanic eruptions can pose significant hazards across a wide range of domains. These hazards can extend from the local environment through to airspace thousands of kilometres away. Beyond direct volcanic hazards such as pyroclastic and lava flows, tephra, and lahars, volcanic ash suspended in the atmosphere can pose a longer lasting and more widespread risk (Auker et al. 2013). Ash can cause acute health impacts to humans at tens to hundreds of kilometres distance, and damage to aircraft airframes or engines at up to 1000 km (Durand and Grattan 2001;Guffanti, Casadevall and Budding 2010). Following the Eyjafjallajökull eruption in 2010, the UK Civil Aviation Authority approached aircraft engine manufacturers to determine volcanic ash concentration thresholds for safe aircraft operation. As a result, the Rolls-Royce Safe-to-Fly chart was produced, using evidence from historical aircraft volcanic ash encounters and test flights during the 2010 eruption (Clarkson, Majewicz and Mack 2016). The main result derived from this investigation was a guideline indicating safe operation up to ash concentrations of 2 mg/m 3 .
However, more recent re-evaluation of the incidents used to inform this initial operational limit suggested that the historical events may have occurred at much lower ash concentrations than previously thought. Additionally, it was suggested that a dosagebased approach, incorporating both ash concentration and duration of exposure, should be considered going forward (Clarkson, Majewicz and Mack 2016). This was further reinforced by Rolls-Royce's Volcanic Ash and Aviation Position Statement, stating that all Trent and RB211 engine types operating in ash concentrations of 2 mg/m 3 for up to 2 h, translating to a dose of 14.4 g s/m 3 , would not experience significant reductions in flight safety margins (Rolls-Royce 2017). As such, the measurement and forecasting of volcanic ash concentrations in the atmosphere is vitally important. This is generally performed through the use of remote sensing, usually satellite-based, and atmospheric transportation and dispersion modelling handled by local Volcanic Ash Advisory Centres (Jones et al. 2007; Thomas and Watson 2010).
Remote sensing of volcanic plumes can provide 'snapshots' of the current state of a target plume. This is generally performed by utilizing measurements from a variety of regions of the electromagnetic spectrum. For volcanic ash this is usually imaging with visible, ultraviolet or infrared light (Thomas and Watson 2010). Of particular interest is the use of multi-spectral infrared remote sensing, taking advantage of volcanic ash's unique spectral response in the Thermal Infrared (TIR) region. By measuring the difference in brightness temperature of a target plume at two narrow-band TIR wavelengths (11 μm and 12 μm), multi-spectral analysis facilitates the identification and discrimination of ash from meteorological clouds (Prata 1989) and quantitative retrieval of particle sizes, mass loading, and optical depths (Wen and Rose 1994). Dispersion modelling allows the forecasting of the location and movement of volcanic ash clouds days into the future. These models make use of the remote sensing snapshots and derived eruption source parameters, such as plume height, mass eruption rate, and particle size distributions to inform their predictions (Webley, Stunder and Dean 2009). Under real-world conditions, these parameters can be poorly defined, leading to large uncertainties in longer term forecasting (Mastin et al. 2009).
One proposed method to advance the monitoring of volcanic ash clouds is the use of multi-angle imagery to provide 3-dimensional information about a target plume. The use of a multi-spectral infrared imaging system along with the quantitative retrieval of ash properties would facilitate a tomographic three-dimensional reconstruction of a volcanic ash cloud, including internal ash concentration distributions Bernardo 2009, 2014;Wood et al. 2019). This reconstruction would act as both an invaluable snapshot of ash contaminated airspace and provide a source for more well-defined inputs into dispersion models. Additionally, with three-dimensional ash concentration distributions, it would be possible to accurately route aircraft through areas of lower concentration, taking into account the safe ash dosage limits of each individual aircraft.
However, there is very limited availability of real-world data to develop these threedimensional reconstruction methods. This is, namely, due to the lack of suitable multispectral and multi-angle satellite imagery, along with the significant dangers of operating multiple ground-based cameras or unoccupied aerial systems (UASs) in proximity to an active volcano. Additionally, the geometrical constraints imposed by both local geography and satellite orbits can limit the available viewing angles. As such, this paper details the development of a Blender-based volcanic plume simulation and image generation framework capable of producing the required imagery. The framework facilitated the development and testing of an initial tomographic reconstruction scheme.
A review of previous related work is provided in Section 2, describing the quantitative retrieval methods and previous stereoscopic and three-dimensional reconstructions. Following this is a description of both the Blender simulation framework and the tomographic reconstruction method in Section 3. Reconstruction results of an idealized 'Shepp-Logan' style test case (Gach, Tanase and Boada 2008) are presented in Section 4 as a series of sensitivity analyses, investigating the effects of a range of observation parameters of reconstruction accuracy. The performance of the reconstruction method, along with the efficacy of the simulation environment are discussed in Section 5, with conclusions and future work detailed in Section 6.

Previous work
The use of multi-spectral infrared imaging in the monitoring and analysis of volcanic ash clouds was theorized in 1989 (Prata 1989) and has been in regular use since (Thomas and Watson 2010;Blackett 2017). Methods for the quantitative retrieval of ash properties using similar multi-spectral infrared imagery were presented shortly afterwards, capable of generating per-pixel line-of-sight mass loading (Wen and Rose 1994). Demonstrations of this quantitative retrieval is less common, with many focusing on proof of concepts with prototype ground-based imaging systems Bernardo 2009, 2014).
Similarly, multi-angle imagery from space has previously been used to determine three-dimensional properties of volcanic plumes. However, these retrievals have generally been limited to simpler properties, such as cloud top heights rather than full threedimensional reconstructions. These retrievals are mostly performed using two stereoscopic images generated by satellite imaging payloads, such as the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER), Advanced Along Track Scanning Radiometer (AATSR), and the Sea and Land Surface Temperature Radiometer (SLSTR), with the limited viewing angles preventing more complex retrievals (Pieri and Abrams 2004;Virtanen et al. 2014;Fernandez-Moran et al. 2021).
More recently, the use of multi-angle imagery from ground-based cameras and UASs has been investigated (Wood et al. 2019;Albadra et al. 2020). Of note is a 2017 Volcán de Fuego field campaign, utilizing four nicAIR multi-spectral imaging systems arranged around the volcano and imaging simultaneously (Wood et al. 2019). However, all such work has been limited to the reconstruction of three-dimensional convex hull representations of the target plumes. By making use of both multi-angle and multi-spectral imagery (from space or ground), tomographic methods would allow the three-dimensional reconstruction of the internal structure and concentration distributions of a volcanic ash plume. This method has been demonstrated using visible spectrum satellite imagery from the Multi-angle Imaging SpectroRadiometer (MISR) payload. This remote sensing instrument boasts nine independent sensors, each arranged at a different along-track viewing angle, and has been used to generate two-dimensional tomographic reconstructions of aerosol plumes (Garay, Davis and Diner 2016). However, the visible spectrum imagery from MISR is unsuitable for volcanic ash tomographic reconstruction as it does not facilitate the required quantitative retrieval.
Similarly, Sulphur Dioxide (SO 2 ) emissions, also present during volcanic eruptions as well as in many industrial plumes, are commonly targeted for observation. Line-of-sight retrievals of SO 2 concentration can be acquired using differential optical absorption spectroscopy (DOAS), facilitating the tomographic reconstruction of atmospheric SO 2 (Frins et al. 2014). When considering volcanic SO 2 , the data collection methods share similar limitations to those introduced in Section 1, mainly the difficulties operating near active volcanoes and geometric constraints limiting viewing angles (Casaballe et al. 2020). As such, simulations have been used to aid the development of these methods, though both simulated and real-world reconstructions of SO 2 have generally been limited to twodimensional slices of the target plumes (Casaballe et al. 2017(Casaballe et al. , 2020Valente de Almeida, Matela and Vieira 2020). Whilst these reconstructions are based on DOAS observations of SO 2 , the tomographic methods employed are still relevant to volcanic ash reconstruction as the input data, the per pixel line-of-sight loading, is of a similar form. Generally, the tomographic methods employed to be used with DOAS observations are iterative reconstruction algorithms, usaully based on algebraic reconstruction techniques (ARTs) such as Simultaneous Algebraic Reconstruction Technique (SART) and Maximum Likelihood Expectation Maximization (MLEM) (Valente de Almeida, Matela and Vieira 2020). Other iterative methods, such as least-squares regression, are also employed (Casaballe et al. 2017). More recently, advances in the miniaturization of DOAS devices are facilitating UAS mounted DOAS observations. Such a system would be capable of producing observations along many more viewing angles, allowing full three-dimensional reconstruction of the target plumes (Valente de Almeida, Matela and Vieira 2020).

Plume simulation and image generation
This section describes the development of a volcanic plume simulation environment built in the 3D computer graphics software tool-set 'Blender', facilitating the generation of imagery suitable for tomographic reconstruction. This simulation environment and the imagery it generated was used to aid in the development and testing of the tomographic method detailed below. A simulated approach to the problem of data collection was used due to the scarcity of suitable real-world imagery and the difficulties involved in collecting new imagery through field campaigns. A simulated approach also provided a reliable 'ground-truth', generated by directly exporting the simulated plume from the Blender framework. This ground-truth was then used for comparison and evaluation of the tomographic reconstruction.
The Blender software package (version 2.93, Blender Foundation (2021)) was chosen due to its focus on physically accurate rendering and its built-in plume simulation. Blender's path-traced rendering engine allows the generation of plume imagery with pixel intensities proportional to the line-of-sight mass loading through the plume that are required for tomographic reconstruction. The built-in plume simulation facilities the native generation and simple configuration of semi-realistic plumes to image. Blender also has some heritage in this area, with smoke simulations being successfully used to generate training data for smoke detection algorithms (Zhang et al. 2018;Nguyen, Quach and Pham 2020).

Plume simulation
Blender offers two plume simulation methods: a relatively simple particle system simulation and a more advanced fluid simulation system based on the 'mantaflow' computational fluid dynamics framework. The particle system method is fairly computationally inexpensive and capable of generating semi-realistic plumes, though it does not easily model internal density distributions and thus is not suitable for tomographic reconstruction. The mantaflow-based method uses a fluid dynamics solver to simulate fluid flow within a domain. This flow can then be sampled at divisions within the domain, allowing the extraction of data, such as the fluid density, velocity, or temperature. As such, it is possible to model the internal structure and concentration of a simulated plume and extract this information for use as ground-truth for comparison. This process can then be repeated over a range of time-steps, simulating and extracting a ground-truth of a dynamically evolving plume if desired. In Blender, this mantaflow simulation is described as a 'smoke' simulation and will be referred to as such during this paper.
The fundamentals of this smoke simulation in Blender are described in the Blender Reference Manual and so will not be covered here (Blender Documentation Team 2021). Of note however, are some specific requirements for the smoke density settings. These settings are necessary to generate suitable semi-transparent plumes that represent the images of ash clouds produced via multi-spectral TIR imaging. Valid values for the smoke density for the specific setup utilized in this work were found to be between 0.1 and 3 'Blender units', with the material rendering density at 0.1 'Blender units'. The lower end of these values produces a fairly optically thin plume, with the upper end of these values produces an optically thick, but not yet opaque, plume. Increasing either of these values further would lead to saturation, where the plume becomes optically opaque. At this point, any further increase in the amount of smoke, either through increased density or increased line-of-sight depth, would result in no change to the measured pixel intensities, thus making any further internal structure indistinguishable. Whilst saturated plumes were avoided for during this study, the simulation environment is capable of modelling them by simply increasing these simulation densities further.

Image generation
Once a plume has been simulated, images can be rendered. Again, the fundamentals of this task are detailed in the Reference Manual (Blender Documentation Team 2021). This rendering process and the subsequent post processing detailed below aims to produce final imagery equivalent to that produced during the infrared multi-spectral retrieval process introduced above. Whilst the simulated rendering process currently does not consider proper radiative transfer calculations or other atmospheric effects, the output images do successfully mimic the key features of actual multi-spectral retrievals, namely a per-pixel line-of-sight mass loading. As such, the simulated imagery can be used to test and develop tomographic methods.
The Blender based rendering process produces images with pixel intensities proportional to, but not equal to, the line-of-sight mass loadings that are needed for reconstruction. An additional processing step is thus required before the rendered images can be used for tomographic reconstruction. As such, a transfer function relating the rendered pixel values to smoke column loading, i.e.; the total mass of smoke within the column defined by the pixel's line-of-sight, was required. To generate this transfer function, a Blender scene was set up to render a series of images of a domain with a precisely controlled smoke column loading. This was achieved by simulating a 'column' of smoke and varying both the length of the column and density of the smoke. Figure 1(a) presents a schematic of this Blender scene. An image metric was then measured for each loading value, defined by the current smoke density and column length. Pixel transparency was chosen for this image metric. Figure 1(b) shows the relationship between pixel transparency and mass loading, with pixel transparency reducing non-linearly to near saturation with increasing mass loading.
This mimics the phenomenon predicted in volcanic ash clouds where, once a plume is optically saturated, unambiguous retrieval of microphysical properties is no longer possible (Prata 2001). Figure 1(c) presents a data set, showing the same non-linear effect, extracted from Figure 2 of Wen and Rose's theoretical radiative transfer calculations for volcanic ash plumes (Wen and Rose 1994). For these calculations, Wen and Rose found measured ash mass loading to be dependant on particle density, effective particle radius, and optical depth. Thus, with a constant particle density and a theoretical plume with a uniform effective particle radius, optical depth is directly proportional to ash mass. The specific values plotted in 1c for optical depth and the corresponding brightness temperature at 11 μm were retrieved for a constant particle radius of 2 μm. This shows the same non-linear decrease in the observed value, in this case brightness temperature, with increased mass loading, in this case directly proportional to optical depth.
With the relationship between pixel response and smoke mass of the simulated Blender imagery displaying a similar non-linear trend to the theoretical volcanic ash response, the images generated from the Blender-based image simulation process can be used as substitutes for the quantitative retrieval process used on real-world TIR imagery. An exponential curve of best fit was calculated for the simulated Blender data to be used as a transfer function in the following reconstructions. The units for the calculated mass loading are in Blender's internal unit system, which is consistent with the extracted ground-truth. The output images produced using this processing step share the same form as the quantitative retrievals of actual volcanic ash clouds, namely the perpixel line-of-sight mass loading. The data also mimics the real world non-linear response to increased mass loading apparent in actual volcanic ash clouds. An additional real-world consideration is the fact that at least two images are required for real-world multi-spectral retrieval, one in each of the narrow-band TIR wavelengths. As such, either two adjacent cameras would be required, or one camera with a filtering mechanism taking two separate images. This could in theory introduce either a positional or temporal displacement between the two narrow-band images making up the single processed image. An additional temporal effect is the time synchronization between cameras at different viewpoints distributed around the target volcano. For the positional displacement, it is reasonable to assume that the distance between the cameras would be negligible compared to the distance to the target plume, leading to no significant effect on the output images. For the temporal effects two timescales are important (a) the time it takes to move the filter wheel in any given camera and (b) the accuracy of synchronization between different cameras. Current filter wheel mechanisms can operate fairly rapidly, leading to the different individual wavelength images likely being separated by under 2 s (Thomas and Prata 2018). At typical distances between the vent and the camera (on the order a few kilometres), and the fairly low spatial resolutions possible with TIR cameras, plume movement and structural change is negligible at this scale and on these timescales. Additionally, clock synchronization between cameras at different viewpoints, either from an internet connected device such as a laptop or smartphone, or though GNSS timing signals, allows for images from different cameras to be acquired within a second of one another.

Tomographic reconstruction methods
This section details the reconstruction approach to be applied to the simulated imagery. The approach utilizes a computed tomography (CT) scheme to retrieve the threedimensional internal concentration distribution and structure of the target plume. Computed tomography is most commonly used in medical imaging as X-Ray computed tomography, but it also has other applications in astronomy, microscopy, and atmospheric science (Deans 2007).
The fundamental process of CT relies on the fact that a pixel's line-of-sight loading (column loading) y, passing through a scalar field of voxels with masses x, can be calculated by where S represents a 'system matrix' containing the path length of each line-of-sight in y through each of the voxels in x. The column loadings in y can be calculated from the measured pixel values by applying the transfer function described above. This is the equivalent step to the Infrared qualitative retrieval methods used for real-world imagery. To calculate the system matrix, S, the voxels intersections for each line-of-sight in y must be calculated. To perform this calculation, a ray can be cast through the predefined voxel field for each relevant pixel within a given image. A voxel traversal algorithm is then applied, identifying the intersecting voxels for each ray. The traversal algorithm used in this work utilizes a Digital Differential Analyser (DDA) to march each ray through the field of voxels, identifying intersections. An adapted implementation of Amantides and Woo's 'Fast Voxel Traversal Algorithm for Ray Tracing' was used (Amanatides and Woo 1987;Mena-Chalco 2021). This process of calculating the system matrix can be very computationally expensive. In order to reduce the size of this problem, an initial Space Carving reconstruction stage was employed. This reconstruction stage takes the same input images and produces a convex hull representation of the target plume. This process is relatively computationally inexpensive and limits the number of voxels to be considered for the voxel traversal, and thus tomographic reconstruction, to only those actually in the plume. The Space Carving method has been used for a range of scenarios including industrial plume measurement (Rusch and Harig 2010), meteorological cloud tracking (Peng et al. 2015), and recently convex hull volcanic plume reconstructions (Wood et al. 2019).
By taking this calculated system matrix S, along with the known set of line-of-sight loadings y, the inverse of Equation 1 can be calculated to solve for the unknown voxel values x. There are many ways to solve this general inversion problem for tomographic reconstruction, mainly focused around the iterative reconstruction algorithms mentioned in Section 2. For atmospheric tomography, these methods can be tuned, making use of prior information and assumptions regarding boundary conditions, gas distributions, and spatial derivatives within the target plume. These can be applied as constraints, weightings, and regularization parameters, ideally resulting in smoother and more accurate reconstructions (Casaballe et al. 2020).
However, as the main goal of this paper is the demonstration of a simulation environment to aid the development and optimization of volcanic plume tomographic reconstruction methods, a simple bounded variable least-squares (BVLS) method was chosen as an initial proof of concept (Garay, Davis and Diner 2016). This method allows upper and lower limits to be placed on the voxel values in x, before solving the constrained leastsquares problem to determine an optimum solution to the inversion problem. For this simulated reconstructions, the bounds were based on the known minimum and maximum voxel masses. With the simulation environment in place, more advanced reconstruction methods and additional tuning steps can be evaluated going forward, allowing optimized solutions to be developed.

Idealized 'Shepp-Logan' style test case
For this initial test of the simulation environment and reconstruction method, an idealized test case was used. Inspiration for the test case was taken from the standard 'Shepp-Logan' phantom developed in 1974 for early computed tomography development (Shepp and Logan 1974). This two-dimensional model utilized a set of ellipses of varying sizes and grey-levels (analogous to density) distributed to mimic the geometric and x-ray attenuation properties of a human head. Three-dimensional models have since been developed for testing of three-dimensional reconstruction schemes (Gach, Tanase and Boada 2008). However, as the Shepp-Logan phantom was designed with medical imaging in mind, the boundaries between the ellipses representing different densities are, generally, very sharp. This mimics the boundaries found in the human body, for example, between bone and brain-tissue, but is not necessarily a good representation of a aerosol plume, such as a volcanic ash cloud. Aerosols or gas concentrations, which are far from their sources generally have low spatial gradients, indicating smooth density changes (Price et al. 2001).
As such, the Shepp-Logan phantom was adapted to develop a simple test case suited to the problem of volcanic ash reconstruction. This test case consisted of a large spherical volume with three smaller internal spherical inclusions with differing smoke densities. These spheres acted as smoke sources, the smoke density for each being allowed to diffuse outwards. The diffuse sphere radius was roughly twice that of the source sphere with the density decreasing by around 10% over that range. The locations, radii, and densities of the spheres are presented in Table 1, the density values are given in Blender's internal unit system, denoted as Blender Units (BU). Figure 2 shows this three-dimensional test case used during the initial tests and development of the reconstruction scheme.
With this test case and the corresponding ground-truth, the accuracy of the tomographic reconstruction, the simulated plume can be evaluated for a variety of user generated observation setups. These setups can vary the number of cameras, and thus the number of viewpoints, along with the position of each camera and its optical properties, such as sensor size and resolution or lens focal length. These variables can be chosen to meet the requirements of any individual study. The view direction for each camera is tracked to the centre of the Shepp-Logan plume using Blender's 'Track To' object constraint.
Within the Blender simulation, the specific scale of both the plume test case and the camera distances is not important as the results are normalized as a percentage of the voxel size of the reconstruction. Additionally, the apparent plume size or 'pixels on target' for the final image depends on both the target size and distance. Imaging a target twice as large but twice as far away will produce a similar image. As such, smaller numbers were use during this investigation for modelling simplicity. This assumption would need to be re-evaluated, however, if atmospheric effects were being considered in the image rendering and processing stages.

Evaluation metric
To evaluate the performance of the tomographic reconstructions, an evaluation metric was needed to compare the reconstructed model with the previously extracted ground truth. As perfect knowledge of the ground truth was available at the same voxel resolution as the reconstruction, a direct comparison on a voxel-by-voxel level was possible. The origins of both the ground-truth and the reconstruction were also consistent, negating the 'double penalty problem' that can occur with grid point-based evaluation metric, such as RMS (Wernli et al. 2008). This makes a root mean squared error calculation both accurate and straightforward and so was chosen for this investigation. To calculate the RMS error (RMSE), the difference in mass between each reconstructed voxel and its corresponding ground-truth voxel was calculated as the absolute error. The sum of this absolute error for each voxel was then used to calculate the root mean squared error for the full reconstruction. Equation 2 presents this calculation, with x i equal to the ground truth mass of voxel i, x i equal to the reconstructed mass, and N equal to the number of voxels within the plume. N will depend on both the size of the reconstructed plume and the voxel resolution decided on by the user. This metric was used to compare the performance of multiple reconstructions of the same ground truth when modifying variables such as the number or resolution of the cameras.

RMSE ¼
ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi

Results
The power of the simulated environment is the ability to run many different reconstructions with varied settings to assess performance over a wide range of conditions. For an initial test of the simulation environment, a reconstruction was performed using the idealized Shepp-Logan test case and simply varying the number and position of the cameras to assess the effect on reconstruction accuracy. Following this is a more general sensitivity analysis investigating the effect of changing a number of variables, such as the pitch angle of the cameras, the spatial resolution of the observations, and the effect of uncertainty in the measurement of camera attitude, specifically, the camera yaw angle. This generalized analysis facilitates the comparison with a variety of possible real-world observation geometries. A detailed discussion of each set of results is presented in Section 5. For all runs, a baseline camera setup was used, with a sensor resolution of 644 � 512 pixels, a pixel size of 25 μm, and a lens focal length of 25 mm. The lenses were assumed to be perfect, with no (or perfectly corrected) optical aberration. These values were chosen to be similar to the specifications for portable microbolometer-based TIR imagers that have seen use in the field, such as the nicAIR cameras (Lopez et al. 2015).

Idealized Shepp-Logan test case
As an idealized initial test of simulation environment, this analysis investigated the effect of the number of unique viewpoints on the reconstruction accuracy of the Shepp-Logan test case. Eight runs were performed, with the number of viewpoints ranging from 3 to 10. In each of the runs, the cameras were positioned to be evenly distributed on a circle with a 100 m radius around the plume. Each of the cameras was fixed to the x À y plane with a z-axis coordinate of 0. Figure 3 presents the 3, 4, and 5 viewpoint camera layouts.
Images were then generated and a reconstruction performed for each run. Figure 4 presents the three rendered and processed images for the three camera run. These images have been cropped to only show the central pixels containing the plume. These images, with pixel intensities equal to line-of-sight mass loading through the plume, are then used as inputs to the tomographic reconstruction process.
With the rendered and processed images generated, the tomographic reconstruction of the plume can be calculated, providing estimates for the mass of each voxel within the plume. As this process produces a three-dimensional plume reconstruction, the results are presented as two-dimensional slices of this three-dimensional plume. Figure 5 presents these internal slices of both the ground-truth plume as well as slices for the runs with 3 and 8 viewpoints. Figure 6(a) presents the RMS error of the tomographic reconstruction against the number of viewpoints used for the reconstruction, with all viewpoints being evenly distributed around the plume as in Figure 3. The general trend is an increase in reconstruction accuracy with an increasing number of viewpoints; however, both the 4 and 6 viewpoint runs have significantly worse performance than their immediate neighbours, with 8 and 10 showing the same effect to a lesser degree. The cause of this, a reduction of unique lines-of-sight when viewpoints were placed directly opposite one another, is discussed in greater detail in Section 5.1. To remedy this issue, an additional set of reconstructions were performed with offset camera positions shown in Figure 10. The results of these offset runs are shown in Figure 6b.

Camera pitch angle and spatial resolution
Whilst the analysis of the number of unique viewpoints presented above provides some useful insights into the distribution of viewpoints when planning an observation setup, in a real-world observation it is unlikely to have cameras perfectly level with the plume. Most ground-based volcanic observations will involve a camera positioned on the slopes of the Figure 5. Slices of the ground truth and reconstruction for the idealised Shepp-Logan test case, each column presents a slice at the x coordinates x ¼ À 2, x ¼ 0, and x ¼ 2 respectively. Each row presents a different reconstruction, with (a-c) the ground truth, (d-f) the reconstruction using 3 viewpoints, and (g-i) the reconstruction using 8 viewpoints.
target volcano, some distance from the vent. For many volcanoes, this would leave the camera below the volcano vent, requiring some amount of pitch angle to frame a given plume.
Additionally, for a real-world observation setup the possible positions of the cameras, and thus the distances from camera to the target plume, will be heavily influenced by the local geography. This can have a significant influence on the spatial resolution of the observations. The spatial resolution, or spatial sampling distance (SSD), is a measure of the size of an individual pixel when projected over the distance from the camera location to the target. The SSD can also be controlled by camera and lens properties, such as the sensor pixel size and the lens focal length. As such, multiple of these variables can all be approximated by a single SSD variable, with the current consideration that atmospheric effects of increased viewing distances are not considered. The SSD can be calculated by first calculating the instantaneous field of view (iFoV) of a single sensor element (pixel) using Equation 3, with d as the pixel size and f the lens focal length. This iFoV can then be projected over the distance to the target plume using Equation 4, with h as the distance to the plume.
Thus, a sensitivity analysis investigating the effect of both the camera pitch angle and the SSD can cover a wide range of plausible real-world observation setups. For this sensitivity analysis, the four viewpoint run with offset cameras, introduced in Section 4.1 and presented in Figure 10b, was chosen as a starting point. Each camera was pitched up by a specified angle and then repositioned downward, maintaining the baseline distance of 100 m, to emulate a position on the slopes of a volcano. The pitch angles were ranged from 0° to 25°, with examples of the 0° and the 25° configurations presented in Figure 7. For each pitch angle, the lens focal length was then changed with all other camera properties remaining fixed. Focal lengths were ranged from 550 mm to provide a range of SSDs from 0.05 to 0.5 m at the target distance of 100 m.
To generalize these results, the SSDs were then normalized as a percentage of the reconstructed voxel resolution of 0.3125 m. This reconstructed voxel resolution is the side length of an individual cubic voxel in the final reconstruction, that is, the smallest feature size of the reconstruction. For example, a SSD of 0.1 m would normalize to 32% of the reconstruction resolution of 0.3125 m. This same normalization can then be used on different scales, for example, a SSD of 10 m and a reconstruction resolution of 31.25 m would also normalize to 32%. Figure 8 presents RMS error against the normalized SSD. Each of the pitch angles investigated is presented as a separate line. Generally, the pitch angle has little effect on the reconstruction accuracy. Each SSD line has a fairly consistent RMS error up to roughly 50% of the reconstruction resolution, at which point the reconstruction accuracy starts to diminish, indicating that there are diminishing returns on improving the SSD beyond around 50% of the desired reconstruction resolution. Further discussion is presented in Section 5.2, along with the application of this generalized set of results to a real-world observation setup.

Camera pointing inaccuracy
For a real-world setup, the measurement of camera location and attitude is an important consideration. Over the distances generally required for volcanic observations, even small pointing inaccuracies can lead to large offsets in the assumed and actual location of a pixel when projected onto the target plume. For this analysis, the effect of inaccuracy in camera pointing measurements was investigated, specifically, the measurement of each camera's yaw angle. The baseline camera setup discussed in Section 4 was utilized, with five viewpoints distributed as discussed in Section 5.1. A range of rotations from 0° to 0.75° were applied to the camera's yaw angle after the images had been rendered for each run. The modified camera attitude was then used in the reconstruction, simulating an error in measurement.
The size of this pointing error, equal to the size of the assumed and actual pixel location offsets at the target, was then normalized by calculating this offset distance and taking it as a percentage of the reconstruction resolution. For example, at the 100 m range used in this simulation, a pointing inaccuracy of 0.2° produces an angular offset of roughly 0.35 m which is around 110% the reconstructed voxel resolution of 0.3125 m. Similarly, at a range of 10 km and with a reconstruction resolution of 31.25 m, the same 0.2° error would also be calculated as 110%, even though the actual pixel offset is now roughly 35 m. Figure 9 shows the RMS error plotted against this normalized pointing inaccuracy. The RMS error can be seen to be fairly constant at low pointing errors. For pointing errors between 0% and 50% of the reconstruction resolution, the RMS error is essentially constant. There is a very slight rise in RMS error between 50% and 150%, and then a dramatic rise beyond 150%. This sharp threshold can be used as an effective maximum pointing error for a given observation setup that depends on both the distance to target and the desired reconstruction resolution. Further discussion is presented in Section 5.3, along with the application of this generalized set of results to a realworld observation setup.

Discussion
The goal for the development of the Blender-based simulation environment described above was to provide a tool to aid the development and refinement of tomographic reconstruction methods for volcanic plumes. Additionally, the tool could be used to aid the set-up and optimization of real-world volcanic observations. Blender was chosen as  the basis for the simulation environment as it natively provides both semi-realistic plume simulation and the path-traced rendering required to produce imagery suitable for tomographic reconstruction.
As mentioned in section 3.1.2, the pixel intensities of the simulated images rendered by Blender are proportional but not equal to the required line-of-sight mass loadings. The transfer function presented in Figure 1 is thus used to convert from the rendered pixel values to mass loading, analogous to the real-world quantitative retrieval process starting from measured brightness temperatures. This transfer function, and its similarity to the theoretical non-linear response of volcanic ash clouds presented in Figure 1c is essential for the simulation environment to be useful in aiding real-world observation set-ups. It was vital that the simulated transfer function had a similar non-linear response to the real-world process, with measured values tending to a limit as the plume becomes optically saturated.
With the simulation environment producing responses similar to real-world imaging, the reconstruction accuracy of some real-world scenarios could be evaluated. The sensitivity analyses presented in the results and discussed below provides some useful insights and suggestions into optimizing a real-world observation campaign.

Number and distribution of viewpoints
This analysis, presented in Section 4.1, investigated the effects of the number of unique viewpoints on the tomographic reconstruction accuracy. This is important as it is a significant investment, or sometimes impossible due to local geography, to set up multiple groundbased cameras to provide unique viewpoints for volcanic observations. Optimizing the number of cameras needed is thus worth-wile. This is less of an issue if utilizing the manoeuvrability of a UAS or satellite to provide multiple viewpoints, though optimizing a UAS's flight path or the desired number of images to capture would still prove useful.
The predicted trend for this analysis was that an increasing number of viewpoints would result in an increase in reconstruction accuracy, up to some limit of diminishing returns. As the number of images from unique viewpoints increases, the amount of unique information captured about the internal density distributions within the plume increases. This assumed trend is confirmed in Figure 6a, with diminishing returns being reached at around 8 to 10 viewpoints for this setup.
However, the figure shows some outliers, most notably for the 4 and 6 viewpoint runs, though also at 8 and 10 runs to a lesser degree. It was theorized that these outliers were the result of having some of the viewpoints mirrored, thus reducing the number of truly unique viewpoints. With two viewpoints directly opposite one another, many of the linesof-sight of the two viewpoints will essentially be identical or very similar, just mirrored.
With the viewpoints evenly distributed on a circle around the origin, the runs with an even number of viewpoints end up with half their viewpoints being perfectly mirrored (e.g. four viewpoints would be evenly distributed 90° apart from one another, two at ½�100; 0; 0� and two at ½0; �100; 0�, see Figure 10a). This issue does affect the 8 and 10 viewpoint runs, though by that point there are enough unique views to somewhat compensate. By shifting one viewpoint in each opposing pair, a truly unique view from each of the viewpoints can be guaranteed, see Figure 10b. Figure 6b presents the RMS error for a set of runs in this configuration, with the 4, 6, 8, and 10 viewpoint runs having half their viewpoints shifted so they are not opposite each other. The results for these shifted viewpoint runs are clearly an improvement and follow the general trend more closely.
Whilst this issue should be considered for similar real-world configurations of cameras, it is unlikely that the perfectly mirrored alignment provided by the simulation would be replicated in the field, even if no consideration is taken for it. Due to the difficult geographical conditions around many volcanoes, imposing constraints on ground-based observation locations, a single optimum configuration for the distribution of the cameras would not be much use. Instead, the simulation framework would be used on a case-by-case basis to optimize for a specific volcano or field campaign.

Pitch angle and spatial resolution
This analysis investigated the effects on the reconstruction accuracy of varying both the pitch angle and the spatial resolution of the imagers. These were modified together as both the spatial resolution, or spatial sampling distance, and the camera pitch angle are likely to be proportional to the camera's distance from the target volcano. The SSD is also proportional to camera's sensor pixel size and the lens focal length, meaning a single variable can be used to model a variety of observation configurations, with the current consideration that atmospheric effects of increased viewing distances are not considered. The SSD is also especially important for volcanic observations with infrared imagers as generally these have much lower resolutions than visual spectrum imagers.
For the pitch angle, the predicted trend for varying this property was uncertain. At very extreme and unrealistic pitch angles (e.g. beyond 60° or 70°) it is clear that the viewpoints of each would become less unique, with a 90° pitch angle essentially requiring all cameras to be at a similar location directly below the target. However, it was unclear whether these effects would appear at more realistic angles up to around 25°. Figure 8 does show a small but consistent increase in reconstruction error as pitch angle increased, though this increase is much smaller than the other sources of error investigated in this paper, indicating that other properties such as the distribution of viewpoints and viewing distances are of more concern.
The predicted trend for the change in SSD was similar to the 'number of viewpoints' analysis. As the spatial resolution improves, meaning the SSD decreases, the reconstruction accuracy was expected to improve, up to a limit of diminishing returns where any improvement in spatial resolution would have no further impact on reconstruction accuracy. This trend can be seen to be confirmed in Figure 8. The point at which diminishing returns is reached can be seen to be when the SSD, the size of a pixel when projected from the camera sensor to the target, is roughly 50% that of the reconstructed voxel resolution. It is assumed that with a SSD larger than or close to the desired voxel resolution, it is likely that a single projected pixel will overlap multiple voxels. As such, it is impossible to determine unique values for these voxels, requiring the pixel value to be 'blurred' across them. With a SSD roughly 50% of the voxel resolution, it is guaranteed that each voxel will have unique information.
With this information, an optimum SSD requirement of at least half the desired reconstruction resolution can be derived. The reverse is also possible. If SSD is fixed due to constrained imager properties or locations, the optimum minimum reconstruction resolution can be taken as double this fixed SSD. The effects of camera pitch angles, when constrained to realistic values, can generally be considered negligible compared to other considerations.
Additionally, this generalized set of results can be used to perform an initial comparison with real-world observation setups, for example, the November 2017 field campaign at Volcán de Fuego in Guatemala, presented by Wood et al. (2019). This field campaign utilized four nicAIR cameras distributed fairly evenly around the target volcano, though with no perfectly mirrored viewpoints. The cameras all had pitch angles between 15° and 20°, and SSDs between 5.7 m and 11.5 m. Whilst tomographic reconstruction has yet to be attempted on the data recorded during the campaign, a simpler Space Carving reconstruction was performed, producing a convex hull model of the plume. The reconstruction resolution selected for this process was 25 m, giving 46% as a worst case SSD, normalized by reconstruction resolution. These values can then be compared with the results presented in Figure 8 as the circle marked 'Fuego 2017'. This comparison indicates that, for a tomographic reconstruction using the same dataset, the same reconstruction resolution of 25 m ideal. This result could also be used to indicate that for a future field campaign, the SSD would likely need to be improved if a finer reconstruction resolution was desired, with additional 'bespoke' simulations allowing further optimization.

Camera pointing inaccuracy
This analysis investigated the effects of simulating inaccuracy in the measurement of the camera attitude. In the simulation framework, it is possible to extract perfect attitude information but for real-world observations this is a significant challenge (Wood et al. 2019). Position measurement is also prone to real-world error but with modern GNSS systems it is possible to achieve accuracies of less than a metre (Dabove and Di Pietra 2019). At the distances generally required when observing volcanic ash and the spatial resolutions generally available with TIR imagers over these distances, errors in the measurement of a cameras position on this scale will have a negligible effect on reconstruction. As such, this analysis only focused on the camera pointing errors, specifically, the yaw angle.
The predicted trend for this analysis was a fairly rapid rise in reconstruction inaccuracy as pointing error increased. Figure 9 confirms this trend, though indicates a distinct threshold below which reconstruction accuracy is not significantly effected. As the apparent size of the offset arising from a given pointing error depends on the distance to the target, these results were normalized. The final pointing error was calculated as the offset of assumed vs. actual pixel locations, taken as a percentage of the reconstructed voxel resolution.
From Figure 9, it can be seen that an apparent offset of up to roughly 150% of the reconstruction resolution has only a small effect on reconstruction accuracy. However, beyond this point the error rises sharply. At low pointing errors with offsets below 50%, it is assumed that the small offset in the position of each pixel when projected onto the plume does not dramatically change which voxels the pixel's line-of-sight intersects. Thus, the system matrix S is not dramatically changed, leading to a similar reconstruction accuracy. For offsets up to around 150%, the calculated pixel line-of-sight intersections will start to change compared to the 'deal' case of zero offset, though usually by only a single voxel. However, beyond this point it is assumed that the offset will always result in the pixel's calculated line-of-sight differing from the ideal. The maximum offset a pixel can have and still be centred on the correct voxel would be the diagonal of size of the voxel. For a square sided voxel, this would be equal to ffi ffi ffi 2 p of the voxel resolution, or roughly 140%. Additionally, beyond 150% the line-of-sight intersections will generally be offset by multiple voxels, further worsening the reconstruction accuracy.
With this information, a limit on the maximum viewing angle measurement error can be calculated. For example, an acceptable attitude measurement error can be calculated for the 2017 Fuego field campaign, described in Section 5.2. With a desired reconstruction resolution of 25 m and the normalized acceptable pixel offset at the target of 150%, an actual pixel offset limit of 37.5 m can be defined. This displacement, at a range of roughly 10 km, means an acceptable viewing angle measurement error of around 0.21°, a level of accuracy, generally, only possible with specialist surveying instruments such as a Theodolite (Piwetz et al. 2018;Wood et al. 2019). One potential workaround, however, is the pre-processing of recorded data to reduce this error. An example of this kind of preprocessing is the alignment of the skyline captured in real-world imagery with simulated skylines generated though the use of a high-accuracy digital elevation model of the target (Wood et al. 2019). Other examples could include alignment using observations of a common geographical feature visible to all cameras. Either way, the very strict pointing accuracy requirement for accurate tomographic reconstruction should be considered when planning a real-world field campaign.

Conclusions and future work
This paper has demonstrated the simulation, rendering, and tomographic reconstruction of an idealized volcanic plume. The Blender-based simulation and rendering framework has been described, along with the definition of a transfer function mapping the rendered pixel values onto the column mass loading needed for tomographic reconstruction. A Shepp-Logan inspired test case suitable for volcanic plume reconstruction was also defined and used in subsequent sensitivity analyses.
The tomographic reconstruction method used has been detailed. The line-of-sight, and thus voxel intersections, of each target pixel were calculated using a voxel traversal algorithm. These intersections, along with the column-mass-loadings, are then used in a bounded variable least squares (BVLS) solver to determine the smoke mass in each of the target voxels.
The simulated approach has facilitated the direct comparison of results with a perfect ground truth model extracted separately from the Blender model. This allowed the accurate and precise quantification of the performance of the tomographic reconstruction algorithm. The simulated approach also enabled easy configuration of the reconstruction settings. This allowed multiple sensitivity analyses to be performed, investigating the affect of variables, such as the number and distribution of observation viewpoints, the spatial resolution and pitch angles of the simulated cameras, and the uncertainties in the measurement of camera attitude.
The configurable environment will enable further work including simulations and reconstructions emulating previous real-world observation setups. This will aid in further developing and optimizing the reconstruction process to be used on realworld datasets. The simulation environment can also be used to aid future real-world observation campaigns by optimizing variables such as the distribution and number of cameras before entering the field. Some initial investigations into this kind of optimization are presented. Some further improvements to the model are also suggested, such as the modelling of atmospheric effects in the image generation process and the development and testing of more advanced tomographic reconstruction methods.
It is hoped that this simulation framework will aid the further development of tomographic reconstruction methods for use with volcanic plumes. It is also hoped that the planning and execution of similar real-world attempts to image and reconstruct volcanic plumes can be informed by the results of further sensitivity analyses, facilitating the optimization of variables, such as camera positioning before entering the field.