The radius of the umbrella cloud helps characterize large explosive volcanic eruptions

Eruption source parameters (in particular erupted volume and column height) are used by volcanologists to inform volcanic hazard assessments and to classify explosive volcanic eruptions. Estimations of source parameters are associated with large uncertainties due to various factors, including complex tephra sedimentation patterns from gravitationally spreading umbrella clouds. We modify an advection-diffusion model to investigate this effect. Using this model, source parameters for the climactic phase of the 2450 BP eruption of Pululagua, Ecuador, are different with respect to previous estimates (erupted mass: 1.5–5 × 1011 kg, umbrella cloud radius: 10–14 km, plume height: 20–30 km). We suggest large explosive eruptions are better classified by volume and umbrella cloud radius instead of volume or column height alone. Volume and umbrella cloud radius can be successfully estimated from deposit data using one numerical model when direct observations (e.g., satellite images) are not available. Total erupted volume and umbrella cloud radius can inform estimates of eruption characteristics and volcanic explosivity more reliably than column height or volume alone, according to simulations of the eruption of Pululagua, Ecuador, 2,450 years ago.

T he main evidence for past volcanic eruptions' size and intensity comes from the interpretation of their tephra deposits. Explosive eruptions can be compared by volume of tephra erupted, which varies by many orders of magnitude and creates impacts on local to continental scales. High mass flow rates, associated with high volcanic plumes, cause tephra to disperse widely, resulting in gradual thinning of deposits and gradual fining in granulometry of deposits with distance from the volcanic vent. Volcanic hazard assessments place a premium on the interpretation of deposits because they inform us about the nature of the past activity and help forecast future eruptions. The better we can reconstruct the nature of past eruptions from deposit data, the better we can anticipate potential future volcanic hazards.
The challenge is to estimate eruption source parameters (ESPs) from the geological interpretation of deposits. ESPs include erupted volume and mass, plume height, total grain-size distribution, mass flow rate, and eruption duration. Eruption volume is commonly estimated using alternative statistical models (e.g., exponential, power-law, or Weibull distribution) to describe the thinning of a tephra deposit with distance from the vent [1][2][3][4][5][6][7] . Plume height can be estimated from the distribution of the largest lithic or pumice clasts in the deposit and plume dynamics [8][9][10] , and is typically used to derive the mass flow rate [11][12][13] . However, statistical methods used to estimate eruption volume and column height are sensitive to deposit erosion 14 . In particular, older deposits are often eroded, reworked, or sparsely sampled, leading to uncertainty and bias in volume and column height estimates 15,16 , which may result in misclassification of older eruptions 17 .
Numerical models of tephra dispersal and sedimentation attempt to address potential bias using the advection-diffusion equation to estimate ESPs, by matching observed deposit features with numerical model output [18][19][20][21][22] . Using inversion techniques, deposit data (i.e., mass per unit area, thickness, local grain-size distribution) are used to estimate the erupted mass, plume height, and total grain-size distribution 5,18,19,[23][24][25][26] . One advantage of these models is that they can better estimate ESPs with uncertainty quantification 25,27 .
Success in estimating ESPs is evaluated based on how well the model, statistical or numerical, fits the observed data 18,[23][24][25][26][27][28][29][30] . However, model assumptions, again either statistical or numerical, also may lead to biased estimates of ESPs. In particular, many numerical models assume that tephra is released from either a point, a vertical line atop the volcano or as diffusion along a vertical line [18][19][20][21]30,31 and that the deposit thins monotonically from the vent, with the exception of secondary maxima associated with ash aggregation or local topography and low atmosphere wind fields 7,32-37 . Large explosive eruption plumes deviate significantly from these simplified plume geometries by producing laterally spreading umbrella clouds around the level of neutral buoyancy 13,[38][39][40][41][42][43][44] . A laterally spreading cloud transports a large volume of tephra rapidly away from the volcano in all directions, reaching radii of 10-100 s of kilometers and markedly changing the distribution of tephra on the ground.
If the lateral spreading of the umbrella cloud is not taken into account, other ESPs estimated with the numerical model must compensate, resulting in biased estimates of their values. For example, advection-diffusion models may overestimate eruption column height from deposits if the gravitational spreading of umbrella clouds is not accurately described in the model 26,28,45 . Diffusion itself, represented in most models by a diffusion coefficient, may be overestimated to compensate for the rapidly spreading umbrella cloud 26 . We address this issue by modifying the advection-diffusion algorithm of tephra sedimentation in Tephra2 19,23 , to model particle release from a disk source rather than from a vertical source. The disk model infers the geometry of the laterally spreading cloud and refines the estimate of eruption volume.
Here we show how to estimate the radius of the umbrella cloud from tephra deposits, thus improving our estimation of other ESPs, using data from the climactic phase deposit of the 2450 BP eruption of Pululagua (Ecuador) as a case study. We find that an erupted mass between 1.5-5 × 10 11 kg, umbrella cloud radius between 10 and 14 km, and eruption column height of 20-30 km can describe the thinning of the Pululagua deposit with distance from the vent. This modeling effort suggests that observed variations in the thickness of tephra deposits of large eruptions are better modeled using a disk source and that the radius of the umbrella cloud, rather than the height of the erupting column, is the primary factor controlling the sedimentation. A tephra advection-diffusion-sedimentation model using a disk source can be successfully used to estimate eruption volume, plume height, and umbrella cloud radius, reducing parameter compensation, and reducing uncertainty.
This outcome has potential implications for how volcanic eruptions are classified and used to compile hazard scenarios. The Volcanic Explosivity Index (VEI) 46,47 is a commonly applied scale that uses volume to classify eruptions on a binned quasilogarithmic scale (VEI 0-8) and provides indications of the plume height. We propose that the VEI scale be updated to include the umbrella cloud radius as a metric to characterize large explosive eruptions, acknowledging that any single ESP cannot completely categorize the nature of large explosive eruptions.

Results
Estimation of the optimal ESPs for Pululagua. The 2450 BP caldera-forming event of Pululagua was a Plinian eruption of dacitic composition, that, allegedly, occurred in a negligible wind field (Supplementary Note 1 26,48 ). Previous studies of the stratigraphic sequence of this eruption 26,48 indicate the activity started with small phreatomagmatic explosions followed by three Plinian phases overlain by a final fine ash layer. Although the isopach maps of some eruptive phases suggest deposition under a slight north-westerly wind, the near-circular isopach map of the climactic phase indicates deposition in still atmosphere. This conclusion is also supported by the deposition of the uniform fine ash layer at the end of the eruption. The fine ash has a longer settling time, which would have made it susceptible to dispersion in a disturbed atmosphere. We assume the previous interpretations of the eruptive conditions of the climactic phase were correct and model the tephra data set collected by Volentik et al. 26 using a disk source (Fig. 1) and invert for optimal ESPs using a grid search, which is appropriate given the limited number of model parameters.
Our simulations with a disk source yield a range of ESPs that explain thickness variation of the deposit: a mass estimate of 1.5-5 × 10 11 kg (i.e., a bulk volume of 0.15-0.5 km 3 , for a deposit density of 1000 kg m −3 26 ), umbrella cloud radius of 10-14 km and an eruption column height ranging between 20-30 km (Fig. 2a). We use the chi-square cost function for the goodness of fit for which a value closer to 0 indicates a better agreement between the field measurements and modeled data. We find that a mass of 2.5 × 10 11 kg (0.25 km 3 ), an umbrella cloud radius of 10 km, an eruption column height of 25 km, and a diffusion coefficient of 9500 m 2 s −1 , give a best fit for the distribution of data and yield a VEI 4 classification for the eruption. This range of erupted mass and volume are consistent with estimates by other statistical and numerical methods (Supplementary Note 1 and Supplementary Table 1) 26 .
Assessment of the model's fit. We evaluate model performance using all data from measured stratigraphic sections of the deposit (Fig. 2b). The regression model shows a coefficient of determination of 0.71 and a coefficient of correlation of 0.84. Ideally, the modeled data should fit the observed data perfectly and yield a slope of 1; the best-fit line slope in our case, 0.68 (95% confidence interval 0.58-0.79) suggests that the model tends to slightly underestimate observed thicknesses, especially for a few large values. Wilcoxon sign and signed-rank tests 49 are used to assess the frequency with which the residuals are negative or positive (i.e., over-or under-estimations of the thickness of the deposit) and the impact of the outliers on the model fit 18,28 . For the sign H is the tephra release height (disk height), C is the center of the disk, and R the disk radius, corresponding to the radius at which wind velocity exceeds radial spreading velocity. The different models of tephra fallout from different source geometries are represented by dots of different colors. Note, also, the relative dispersion of particles from their respective source geometries (not to scale). The inset figure represents a 2D view from above of the disk source, resembling the umbrella cloud. , and umbrella cloud radius (R) that can explain the mapped deposit of the climactic phase of the Pululagua 2045 BP eruption. The red dashed line shows the best-fit model. Thickness of the deposit was obtained from mass/area using a bulk deposit density of 1000 kg m −3 26 . The χ 2 -criterion for the goodness of fit is 9300 and 3600 for the maximum and minimum of the estimated range of ESPs, respectively, while for the best fit is~1800 (a value closer to 0 indicating a better agreement between observed and modeled data). b The relationship between the observed and modeled tephra deposit thickness. The dashed line represents a correlation of 1 and the solid line represents the best fit between the observed and modeled data. Both root mean square error (RMSE) and the coefficient of determination (r 2 ) indicate the quality of fit between the modeled and observed thicknesses. test, we have 46 negatives out of 72 residuals (observed − modeled), which indicates that the model tends to overestimate the deposit thickness. For the signed-rank test, we set a null hypothesis that there is no difference between the observed accumulation and the modeled one 18,28 . For n = 72 samples, we calculate the sum ranked of the positives T + = 1197. From this sum, we can estimate the expected value of E(T + ) at 1314 and the variance Var(T + ) = 181, Z = −0.64 (i.e., p > 0.26 for a one-tailed hypothesis), higher than the significance level of 0.05 (5% confidence). Thus, our null hypothesis cannot be rejected, and we consider that the over-estimation of the deposit is caused by outliers in the thickest part of the deposit. We conclude that although our model has a tendency to underestimate thickness proximal to the vent, it reasonably reproduces the overall thickness variation of the deposit. This is encouraging considering that over time several processes (i.e., reworking, compaction, erosion of the fines) have affected the thickness of the Pululagua deposit.
Duration for the umbrella cloud emplacement. Several explosive eruptions in recent decades produced large umbrella clouds that spread quickly across vast distances. Satellite measurements at Pinatubo, 1991, showed that the umbrella cloud reached a 280 km diameter within the first hour of the eruption 50 yielding a spreading rate of~4.6 km min −1 . Although our model does not account for the dynamics of lateral spreading clouds, we can estimate how long it took for the estimated umbrella cloud of Pululagua to reach the model radius. Using models for the rate of a spreading umbrella cloud 40,41,51 and the previously estimated mass flow rate for Pululagua 26 , we calculate that the umbrella cloud formed in 135-335 s after the onset of the eruption. This yields a 2.5-4.4 km min −1 spreading rate for a 10-14 km radius umbrella cloud. High spreading rates of radial (or asymmetrical) clouds facilitate the transport of coarser tephra particles to larger distances resulting in coarser-grained deposits over vast areas 13,52 ; the bulk of particles ejected during a volcanic eruption will fall from this umbrella region 13,36 .

Discussion
Modeling tephra sedimentation using a disk source geometry provides a better estimation of the ESPs without resorting to the use of highly implausible physical parameters (e.g., very high eruption columns; very high diffusion coefficients). Although still simplified in source geometry, our model results are in agreement with the umbrella cloud model, whereby the bulk of tephra sedimentation occurs in regions directly under the leading edge of the cloud 10 with sedimentation occurring past this point largely due to atmospheric diffusion. When compared with other source geometries our model provides an overall better fit of the deposit than models without a disk source. We compare the disk source with a point and vertical line sources in Fig. 3. For this comparison, we keep identical input parameters across the three simulations (i.e., erupted mass, total grain-size distribution, diffusion coefficient, and particle release height (Supplementary Table 2)). We find that both point and line sources tend to overestimate deposit thickness near the vent and show exponential decay of thickness with distance from the vent (Fig. 3). In the Pululagua deposit, there is a significant variation in thickness near the vent in different outcrops (>20 cm) and a slower decay in deposit thickness than expected for exponential thinning. Line and point sources can model this change in deposit thickness with distance from the vent, but only if the model compensates for the narrow source by greatly increasing the diffusion coefficient parameter, thus, smoothing the modeled deposit 26,34 .
The diffusion coefficient is an empirical value used to characterize complex processes occurring in a convective plume and during atmospheric diffusion. Often, the diffusion coefficient in advection-diffusion models is used to increase the dispersal of fine ash in simulations 19,26,40 . Our simulations with a disk source show that high diffusion coefficients (e.g., >20,000 m 2 s −1 ) underestimate the rate of thinning with distance and produce a flat deposit, whereas lower values (e.g., <10,000 m 2 s −1 ) fit the observed thinning (Fig. 4a, b). By introducing the disk source, the observed deposit thickness variation is explained without resorting to a very high diffusion coefficient (e.g., >10 4 m 2 s −1 ) 53 .
Modeling the disk radius with a low diffusion coefficient is useful because it allows us to directly estimate the umbrella cloud radius, defined as the radius the cloud reaches before advection is dominated by wind velocity 13,52 . Our model of the Pululagua deposit constrains disk radius within a few kilometers of uncertainty. An umbrella cloud radius of 10-14 km gives the lowest relative error between the observed and modeled deposit thicknesses (Fig. 2a). We find the model fit is quite sensitive to disk radius; significantly poorer fit is achieved with disk radii that are smaller or larger than this range (Fig. 4a).
Interestingly, introducing the disk radius parameter decreases the model sensitivity to eruption column height. A higher eruption column allows for a longer settling time for finer particles, increasing the distribution of tephra. Our simulations show that a considerable range of eruption column heights (i.e., 10-35 km) equally explain deposit thickness within the uncertainty of the observations, especially at distances >15 km from the vent (Fig. 4c), where the risk to infrastructure and people is likely greater.
We note that sometimes very high volcanic eruption columns (40-60 km, e.g., the Ilopango, Toba, Campanian, and Taupo eruptions) are invoked to explain the observed thinning of very large eruption deposits [53][54][55][56] . High eruption columns are required by these models to increase the total fall time of tephra and to explain the slow thinning of the deposit. The properties of the atmosphere at these very high altitudes are, however, inconsistent with the particle settling characteristics. Tephra fallout models assume that particles reach their settling velocities as a function of atmospheric density. The atmosphere at very high altitudes has such a low density that particles fall from the high altitude into the denser lower atmosphere before they actually reach their settling velocities. This means that increasing the eruption column height is not actually an effective way to increase tephra dispersion very far from the vent. Instead, we suggest that a disk source model for sedimentation is more consistent with the physics of large explosive eruptions as shown by models for laterally spreading umbrella clouds [40][41][42] . Previous models for lateral spreading clouds showed the implications the umbrella has for the accountability of tephra fallout at larger distances from the vent in short time periods 40,41 .
Using the Ash3D numerical model, Mastin et al. 40 noted that for large explosive eruptions umbrella clouds control the dispersal of  Fig. 3, only radius, R, changes between curves. The χ 2 -criterion for the goodness of fit for R = 10 km is 1800, whereas for R = 5 km χ 2 is 2200. b The thinning of the deposit using different diffusion coefficients K; M, H, and R values are the same as in Fig. 3. c Change in the deposit based on different column (disk) heights: M, K, and R are the same as in Fig. 3, only H is varied.
10-100 100-200 200-500 500-1000 >1000 a Adapted from Newhall and Self 47 . b Plumes of transient eruptions will also develop gravitationally spreading clouds at the neutral buoyancy level, although of limited radii due to lower mass discharge rates. tephra in simulations and that simulations are less sensitive to eruption column height. New models to estimate eruption column height from maximum clast dispersal also indicate that the distance traveled by a clast could be explained by a lower lateral spreading cloud instead of a higher sub-vertical plume 8,9 .
Here we confirm, using field data and a revised numerical model, that observed variations in the thickness of tephra deposits of large eruptions are better explained using a disk source. Our analysis of the Pululagua deposit indicates that the radius of the umbrella cloud, rather than the height of the erupting column, is the primary factor controlling the sedimentation. Furthermore, the diffusion coefficient becomes less important in simulations; the thinning of the deposit is driven mainly by the spreading plume with only the distal portion of the cloud being significantly advected by the wind. Altogether, a tephra advection-diffusion-sedimentation model using a disk source can be successfully used to estimate eruption volume, plume height, and umbrella cloud radius and better account for uncertainties in these ESPs.
Our findings have implications for the way volcanologists evaluate large explosive eruptions and offer improvements to volcanic hazard assessments. We suggest for example that the use of the most common eruption magnitude scale, the VEI 46 , can be improved by coupling erupted volume and umbrella cloud radius as metrics to characterize and classify large eruptions (e.g., VEI > 4), rather than relying only on eruption column height to classify the intensity of such events. An additional category dedicated to umbrella clouds can be inserted in the VEI scale (Table 1) using radii ranges based on umbrella clouds observed in near real-time satellite imagery data, and values estimated from field deposits of past eruptions (Table 2). Table 1 shows a tentative update to the VEI scale based on a limited number of observed eruptions and existing published estimations of umbrella clouds of past VEI 7 eruptions. This proposed rubric is semi-logarithmic for critical umbrella cloud radii 10-1000 km as umbrella radius increases more slowly than eruption volume. It is admittedly complicated by factors such as substantial co-pyroclastic density currents plumes from relatively small-volume eruptions (e.g., 1990 Redoubt eruption 57 ), and tentative due to the current sparsity of plume radii estimates from tephra deposit data.
Based on the Pululagua analysis, umbrella cloud radii can be approximated from field deposits with lower uncertainty than eruption column height, and eruption volume can be estimated simultaneously using a single numerical model. Ultimately, the method presented here can be used to further constrain the relationship between volume and umbrella cloud radii of large eruptions. Subsequently, numerical inversions can be conducted to estimate the ESPs of a sufficient number of large past eruptions to update the VEI scale. At the same time, real-time satellite imagery can be used to extract valuable information on the extension of umbrella clouds of modern eruptions. An updated VEI scale that includes an umbrella cloud radius as a metric for large eruptions would better inform models for tephra sedimentation and make these more robust tools for future hazard assessment.

Methods
Advection-diffusion-sedimentation model (ADS). We modify the Tephra2 algorithm 18,19 , an Eulerian model that uses the advection-diffusion-sedimentation (ADS) equation 3,14,18,21,58 to estimate the ground accumulation of particles released from a source above the volcano. The source describes the volcanic plume (i.e., source term) that in sedimentation models comprises a series of empirically derived parameters-total erupted mass (M T ), total grain-size distribution (TGSD), column height (H). Typically, the source term assumes simplified geometries such a point or a line source, but other geometries have been tested (e.g., horizontal line, disk 59 ). We developed a Python algorithm in which we implement a disk-shaped source term to account for large volcanic umbrella clouds observed in nature. The umbrella cloud is inserted as a disk and does not account for the laterally spreading plume dynamics (i.e., not modeling the ash transportation and diffusion within the plume). Instead, we assume the cloud is already emplaced (i.e., the cloud spreading velocity is equal or lower than surrounding wind velocity 44 ) and the modeling of particles' sedimentation starts with the particles falling under gravity from the height at the base of the cloud (H in Fig. 1). We refer the reader to Connor et al. 23 and Bonadonna et al. 19 for a more detailed description of the equations used in the Tephra2 code, including particle settling velocities and particle mass fraction calculations, which remain unchanged in the Python model. Here we describe the implementation of the umbrella cloud source term.
We assume a well-mixed radially spreading cloud of uniform thickness and density 13,44,60,61 with a total mass of: where R is the radius of the umbrella cloud and η the mass fraction of tephra in the cloud 60 . Next, we discretize the disk in a number of equally spaced grid cells (Fig. 1, inset). The distribution of particle classes in the cloud is uniform and based on the values published by Volentik et al. 26 . The mass of tephra in each grid cell (m cj ) encompassed within the disk radius will be: where N cj is the total number of cells describing the disk.
Once the geometry of the source and the distribution of mass are assigned, the code calculates an analytical solution of the ADS equation by integrating for different particle sizes (ϕ) released from the cells describing the umbrella cloud (c j ). We now calculate the accumulation of particles at one point on the ground using: where S is the particle settling velocity, K is the atmospheric diffusion coefficient, and u, v wind field components (direction and speed). In the end, we integrate over a range of particles sizes released from different grid cells and calculate the total mass accumulation of tephra on the ground: Due to the high complexity of the natural processes governing tephra dispersal and deposition, this model is simplified through several assumptions: the surrounding topography is uniform (i.e., flat); particles are released from specific points in the designed disks and fall according to their own settling velocity, which varies with atmospheric density 62 , through an atmosphere with a constant diffusion and no vertical wind.
Tephra transport in umbrella cloud models. We used the density-driven cloud model 13,40,41,51,57 to see if the umbrella cloud radius obtained from our model can be explained given the estimated mass flow rates reported for Pululagua 26 . We implement in code the following equations to estimate the radius of an umbrella cloud and time of emplacement: where λ is an empirical constant evaluated at 0.2 41,51 , N is Brunt-Väisälä frequency estimated at 0.02 40,41,51 . The volumetric flow rate in the umbrella cloud is represented by q and can be estimated using: where M is the mass discharge rate and k efficiency of air entrainment 40,41,51 . We used C = 0.5 × 10 10 m 3 kg −3/4 s −7/8 , a proportionality constant for eruptions in tropical regions 41,51 .

Data availability
The field data used in simulations are available to download at https://github.com/ robertoon/umbrella_cloud_model/tree/main/data, or by request to the corresponding author.

Code availability
The Python code developed to simulate sedimentation from an umbrella cloud and instructions for its use are available to download at https://github.com/robertoon/ umbrella_cloud_model, or by request to the corresponding author.