Global Rates of Subaerial Volcanism on Earth

Knowledge of the global rates of volcanism is fundamental for modeling the Earth, as those rates closely relate to plate tectonics, crustal growth, mantle dynamics, atmospheric evolution, climate change, and virtually any aspect of the global Earth dynamics. In spite of their huge relevance, the global rates of volcanism have remained unknown, hidden within data that appeared disordered, largely fragmented and incomplete, reflecting poor preservation of small eruptions in the geological record, rareness of large eruptions, and distributions far from normal. Here we describe and validate a model that reproduces global volcanism to high statistical significance, and that is so simple to comfortably fit on a t-shirt. We use the model to compute the expected rates of global terrestrial volcanism over time windows from 1 to 100,000 years, and validate it by comparing with observations back to a few million years. Notably, the model can be tested against independent observations collected in the near future, a feature which is relatively uncommon among global models of Solid Earth dynamics.


INTRODUCTION
The global rates of volcanism, or the rates at which magma is transferred from the Earth interior to the surface, represent one fundamental aspect of the dynamics of our planet. In fact, those rates relate to some of the most relevant large scale geophysical processes on Earth, such as mantle dynamics (Huang et al., 1997), Earth crust formation (Ito and Clift, 1998), plate tectonics (Xu et al., 2009), and climate change (Robock, 2000). Quantitative estimates of the rate of volcanism are made difficult by the apparently disordered occurrence of volcanic eruptions, as it emerges from observations and field reconstructions at the basis of volcanic databases. In fact, volcanic eruptions appear on Earth at times and with scale and violence which have largely escaped any attempts of parameterization or understanding. Global subaerial volcanic eruption databases such as the Smithsonian GVP (Global Volcanism Program, 2013: https://volcano.si.edu) and LaMEVE (Crosweller et al., 2012: www.bgs.ac.uk/vogripa) are major data sources to understand and model volcanic activity on Earth. However, the information they bring is highly heterogeneous as a consequence of the enormous differences in the frequency and potential for geological preservation of eruptions spanning from highly frequent, small effusions to rare, globally impacting catastrophic events (Crosweller et al., 2012;Brown et al., 2014;Kiyosugi et al., 2015).
Discovering that over the global scale, eruption inter-event times are exponentially distributed (Papale 2018), or equivalently, that global volcanic eruptions are Poisson distributed, allowed a robust definition of catalogue completeness for each eruption magnitude in terms of VEI (Volcanic Explosivity Index, Newhall and Self, 1982). The exponential distribution of eruption inter-event times helps understand the difficulties in obtaining a global picture from the dataset. In fact, the exponential distribution is the maximum entropy probability distribution for a random variate (with the random variate being greater than zero and over a given range of variability). The implication is that the occurrence of volcanic eruptions appears as much disordered as it could be. When summed up with effective holes and gaps in the record, and with size-dependent deterioration of the information with age, the distribution appears as a messy sequence of data where variably long periods with high eruption rate alternate without any apparent order to periods with low or very low rate, and to more or less long periods of stasis ( Figure 1). That apparent disorder turned into high order once the exponential distribution of return times was recognized, also providing the theoretical basis for a global assessment of volcanic hazard (Papale and Marzocchi, 2019).
While volcanic eruptions appear on Earth at times reflecting an underlying Poisson distribution, their size distributes as a power law over a large portion of its spectrum (Papale et al., 2021). Small eruptions with volume less than order 10 2 Mm 3 , and giant eruption above 10 3 km 3 , depart from the power law, with the former depicting a log-normal distribution, and the latter displaying a more rapid fall in their frequency than for a power law, as for so-called tapered power law distributions (e.g., Kagan and Schoenberg, 2001;Vere-Jones et al., 2001;Geist and Parsons, 2014).
The above describes the global time-size distribution of volcanic eruptions on Earth. The mathematical counterpart is comparably simple (Papale, 2018;Papale et al., 2021): In Eqs 1-3, S represents the complementary cumulative density function (also called "survivor" function), which expresses the probability that the distributed quantity (β i in Eq. 1 and V in Eqs 2, 3) be larger than any considered value. Eq. 1 describes the exponential distribution of eruption interevent times β i around an average corresponding to the inverse of the rate parameter λ i (the subscript i refers to any grouping, such as the one in VEI eruption classes, or no grouping at all if i is taken to refer to the whole class of subaerial terrestrial eruptions). Eruptions distributed in time according to Eq. 1 are associated with volumes distributed according to Eqs 2, 3. Eq. 2 describes the initial log-normal distribution with mean μ and standard deviation σ, and Eq. 3 describes the subsequent tapered power law distribution, with V 0 being the intercept of the power law with the unit frequency axis, V C the corner volume beyond which the frequency decay exceeds that of the power law, and k the power law exponent. All of the above quantities have been numerically determined, and their distributions are reported in Papale (2018) for Eq. 1, and Papale et al. (2021) for Eqs 2, 3.
In the following we show that such a remarkable simplicity embeds the globally observed time-size distribution of subaerial volcanism on Earth, including the complexities introduced above (and partly visible from Figure 1): variably long intervals with average eruption rates significantly above or beyond the longterm average, irregularly alternating and separated by randomly distributed inactive periods with apparently random length. The time-size distribution appears even more disordered, with increasingly large eruptions emerging here and there without any apparent regularity, up to rare cataclysmic events which sometimes emerge as clustered, sometimes are isolated, several tens of thousands years distant from each other.

METHODS
Our objective is twofold: 1) to prove that Eqs 1-3 embed the complexities of the global time-size distribution of terrestrial subaerial volcanism, and 2) to use the model constituted by that set of equations to determine the global rates of volcanism. To achieve such objectives we adopt a Monte Carlo approach whereby Eqs 1-3 are employed to simulate 1,000 different but statistically equivalent, equally possible eruption histories of the Earth, each one 100,000 years long. The memoryless property of Poisson distributions allows us to analyze the entire sequence over mutually independent intervals of any time length. The global rates are analyzed in terms of 1) number of eruptions of any given VEI scale and 2) volume of discharged magma. Extensive discussion of the advantages and limitations in employing different indicators of eruption size, such as VEI, volume, or mass, and on the meaning of individual eruptions in the database (and consequently in our analysis), can be found in Papale (2018) and Papale et al. (2021). Exponential distributions have been used by previous authors (one notable example is Mendoza-Rosas and De La Cruz Reyna, 2008) to model eruption occurrence at various spatial scales. The exponential distribution of inter-event times, or equivalently the Poisson distribution of eruption events, characterizes many other natural and man-made phenomena as diverse as the number of meteorites greater than a given size that strike the Earth in a given period of time; the number of calls to a telephone service, or the number of alpha particles hitting a Geiger detector, in a particular time interval; and many others (Johnson et al., 2005). The fundamental characteristic of all such phenomena is that the occurrence of each individual event is independent from the time since the previous one, which is known as the "memoryless" property (Balakrishnan and Basu, 1995).
The memoryless property of Poisson-distributed events is utilized in this work, as described in the following point-bypoint description of the employed procedure. 1) We adopt a Monte Carlo approach to simulate the occurrence of eruptions of any size in a 100,000 years time frame. Return times are sampled from Eq. 1 through standard inverse transform sampling: Homogeneous variates u 0;1 are generated with the ran2 package (Press et al., 1992), which guarantees genuinely random outcomes up to computer floating point, a characteristic which is relevant in the present case where >4 billion eruptions are generated (see below). Exponentially distributed random inter-event times β for each VEI class 0-8 are then obtained from Eq. 4. Ordering all events according to their time of realization results in one statistical outcome for the eruptive history of the Earth over a 100,000 years long time frame.
2) The operations at point 1 are repeated 1,000 times, obtaining 1,000 possible volcanic eruption histories each 100,000 years long. Because of the memoryless property, the simulated individual eruptions are independent from each other. Accordingly, the obtained histories can be equally regarded as referring to a 1 × 100,000,000 years long time frame, or to any other combination ensuring constant product of number times years. Because the shape of a Poisson distribution depends on the time window, the distributions associated with each different combination differ from each other. In order to ensure consistency and properly analyze the variability of the results with the length of the observational time window, for each different time window of length Δt we consider a number N of intervals such that N · Δt t TOT , with t TOT 10 8 years (that is, 1,000 periods each 100,000 years long, or 10,000 each 10,000 years long, up to 100,000,000 periods each of length 1 year).
3) We use Eqs 2, 3 to assign a volume to each simulated eruptive event. In order to account for the uncertainty in the definition of the continuous volume distribution (Papale et al., 2021), the 1,000 eruption histories of the Earth, each 100,000 years long, have been simulated by varying the parameters of the continuous volume distribution according to their respective variability and interdependence (Papale et al., 2021). For each eruption we first produce a random outcome from a homogeneous u 0;1 variate and compare it with the value of the cumulative distribution at V min , Eqs 2, 3, to determine whether the corresponding volume lies on the log-normal or power law portions of the distribution. Then, we use standard inverse transform sampling to relate the homogeneous variate to the non-homogeneous volume distribution. This is done by numerically solving the two equations below, for the lognormal and the power law sections of the distribution, respectively: The distribution from Eqs 2, 3, obtained by solving Eqs 5, 6, provides the constrained probability (that is, the probability given the occurrence of an eruption) of observing a given eruption volume, from the smallest lava flows to the largest VEI 8 supereruptions. In order to use it in the present analysis, another property of Poisson distributions is employed, namely, that the sum of Poisson events is also Poisson distributed with rate parameter given by the sum of the rate parameters of the individual Poisson events. In other words, the sequence of eruptions with assigned VEI, obtained at points 1 and 2 above, can also be regarded as a sequence of generic eruptions with volume assigned following the procedure at point 3 above (note that the distributions of inter-event times and volumes are entirely disjointed).
The procedure above allows us to separate the assignment of VEI from the assignment of volume to each individual eruption. In fact, especially for low to intermediate 0-5 VEI eruptions, there is no univocal relationship emerging from the database, between assigned eruption VEI and range of erupted volume (Papale et al., 2021). At the same time, the obtained global Earth histories in terms of sequence of eruptions each characterized by its VEI (points 1-2 above) or by its volume (point 3 above) are entirely consistent with each other.

Number Rates of Subaerial Volcanism
Through the procedure at points 1-3 above we obtain a synthetic dataset constituted by more than four billion individual eruptions in total, described in terms of either the eruption VEI or the Frontiers in Earth Science | www.frontiersin.org June 2022 | Volume 10 | Article 922160 eruption volume, and consistent with the reconstructed overall distribution of volcanic eruptions on Earth in terms of inter-event times, relative and absolute probability of each VEI class, and distribution of eruption volumes. Number rates are entirely related to Eq. 1, and they are obtained through Eq. 4. Table 1 and Figures 2A,B report in numerical and graphical form, respectively, the distribution of the number rates of subaerial volcanism on Earth for observational time windows from 1 to 100,000 years. The distributions are normal (a consequence of the central limit theorem) with mean (and median) at 40.5 eruptions per year, and 1st and 99th percentiles over a 1-year window of 26 and 56 eruptions per year, respectively.
The capability of Eq. 1 to provide accurate forecasts of the number of observed eruptions of any VEI size is illustrated in Figure 3 for three different time windows of 1, 10, and 100 years. The comparison involves data for each VEI class up to corresponding catalogue completeness (Papale, 2018). For catalogue completeness extending back by N years, there are N − Δt + 1 different observational time windows of length Δt; in other words, Figure 3 shows the model performance if we were sitting on one random year up to catalogue completeness and using the model to forecast the number of eruptions of any given VEI (or any combination of VEI) expected in the next Δt years. The nice correspondence between the model forecasts and the FIGURE 2 | Computed distribution of the annual eruption number (A,B) and volume (C,D) rates. The left panels show the cumulative density function, while the right panels are probability density functions obtained by fitting the cumulative distributions. For the non-normal volume rate distributions in panels (C,D), the fit is made to a skewed log-normal distribution by means of a method of moments estimator (Ghorbanzadeh et al., 2017) with computation of the Owen function as in Patefield (2000).  (Table 1 and Figures 1A,B). Clearly, the very last years have not been very productive in terms of number of eruptions. By comparison, there were 38 eruptions in 2018, and an average of 38.6 eruptions per year in the five years from 2014 to 2018, in between the 25th and 50th percentiles and close to the computed long-term average of 40.5 eruptions per year. The average rate over 100,000 years for the ten example cases in the figure is in the 3-4.5 km 3 /yr range. The long-term average rate when considering all of the simulated eruption sequences turns out to be 3.92 km 3 /yr (Table 2 and Figures  2C,D), two to four times larger than existing seldom rough estimates in the range 1-2 km 3 /yr (Nakamura and Furomoto, 1974;Fujii, 1975;Crisp, 1984). Figures 2C,D show that the volume distribution progressively changes from close to normal for long time windows, to markedly right-skewed for the shortest analyzed time windows of 1 and 10 years. Zooming into progressively shorter time periods ( Figures  4B-F) reveals more features of the time distribution, as smaller eruptions progressively emerge. At the 10,000 years scale (panel B), VEI 7 eruptions are visible as clear steps in the cumulative trend. Similar steps become more and more frequent by progressively zooming in, down to the shortest 10 and 1 year scales (panels E-F) where the smallest eruptions become visible and the variability in the observed average volume production rates is the highest. At this temporal scale, which corresponds to that of our observations, the complexities of the true distribution clearly emerge. Entropy maximization due to exponentially distributed inter-event times, combined with log-normal or power law distribution of erupted volumes, efficiently mix up resulting in apparently disordered trends, similar to the observed ones, which in fact result entirely by the extreme simplicity of the present model as it is shown in Eqs 1-3. The volume distributions produced by the present model are compared to the observations in Figure 5. Only eruptions with VEI ≥4 can be considered for this comparison, as the databases do not report the volumes for eruptions with VEI 3 or less. Catalogue completeness for VEI ≥4 eruptions extends back to year 1840 CE (Papale, 2018), while the youngest of such eruptions with associated volume in the databases (up to year 2018) dates back to 2014 CE (https://volcano.si.edu). Therefore, we consider the forecasts we would make on any year from 1840 CE projecting into the future up to year 2014 CE, and compare with the corresponding observations. Note that such forecasts are more correctly referred to  catalogue completeness belong to the same distribution as from the model given by Eqs 1-3. Figure 6 shows the yearly eruption volume rate per volume size class. The mean value (thick solid line) does not depend on the length of the observational time window, as it is explained above. The median as well as any other percentile, instead, depends on the observational time window. Panel (A) adopts a linear scale for the yearly eruption volume rate, whereas panel (B) shows it on a logarithmic scale that better highlights the distribution of both the marginal values and the percentiles. The figure shows that most (>50%) of the total volume erupted from global subaerial volcanism is discharged by eruptions with size corresponding to 0.1-1 km 3 . That size interval, therefore, combines with frequency in the most effective way when compared to any other size interval. Summed up with eruptions in the range 1-10 km 3 , this class of medium-sized eruptions makes up nearly three quarters of the global terrestrial subaerial eruption rate. In spite of their order-of-magnitude higher frequency of appearance, small eruption with individual volumes <0.1 km 3 make up less than 17%; and in spite of their order-of-magnitude larger individual volumes, less frequent large to cataclysmic eruptions do not achieve 10% of the total production rate. Volcanic super-eruptions by themselves (individual eruption volume ≥1000 km 3 ) contribute only a negligible proportion (<0.4%) of the global subaerial production rate on Earth.

Volume Rates of Subaerial Volcanism
The model provides a rare opportunity to be formally validated (Marzocchi and Jordan, 2014, for the analysis of   Figure 7 shows such forecasts for the total erupted volume from global subaerial terrestrial volcanism, in terms of probability of exceedance (i.e., frequency of observations exceeding any given volume) over time windows of 1, 10, and 50 years. The curves in Figure 7 provide a synthetic picture of forecasts which are also relevant as robust inputs, constraints or tests for other models of global Earth dynamics.

DISCUSSION
The present model provides a statistical representation of the rates of subaerial volcanism on Earth, consistent with the available data in terms of number and volume distributions. The model is extremely simple. Its foundation rests in the observation that over the global scale, eruption inter-event times are exponentially distributed (Papale, 2018), or equivalently, that the eruptions are Poisson distributed. Eq. 1 expresses that characteristic, and is by itself sufficient to statistically simulate the subaerial eruption history of the Earth in terms of time-VEI distribution. Combined with the global eruption volume distribution at Eqs 2, 3, the model describes the global time vs. erupted volume distribution. With just that little work, the apparently disordered trends observed for the subaerial volcanism on Earth are statistically reproduced: largely variable eruption frequencies, seemingly random, irregularly long quiescent periods, complex succession of events with extreme size variability, and seldom appearance of large eruptions up to globally impacting ones, become fully logical and statistically reproducible; and the rates of subaerial volcanism emerge from the statistical analysis of a large number of possible global volcanic histories, all of them consistent with the observations, and in a statistical sense, all of them equally plausible. Importantly, the number and volume rates presented here are statistically testable, through comparison with additional data on past eruptions as they are constantly produced, as well as against future observations. Those rates can be employed to figure out if currently observed global activities fall close to averages or depict instead periods of significantly lower or higher activities, now quantitatively defined; and can be employed as an input or constraint for other Earth system models, like global tectonics, mantle dynamics, climate change, and others. Previous estimates of global rates of subaerial volcanism on Earth are sparse. Nakamura (1974) and Fujii (1975) suggested 1.85 and 2.2 km 3 /yr, respectively, as the sum of volcanism from subduction areas, intracontinental and intra-oceanic hot spots. Crisp (1984) provided a more systematic analysis, although still limited by the data available at the time, concluding that much of the uncertainty associated with volume estimates was by then impossible to quantify. Crisp's summary suggested 0.65-0.85 km 3 /yr as the sum of subduction-related and intraplate volcanism. That would increase to 0.8-1 km 3 /yr, by including the contribution from the estimated average melt ascent rates beneath Iceland (Eksinchol et al., 2019) and assuming a 5:1 ratio between intrusive and extrusive magmatism, typical of oceanic ridge environments (Crisp, 1984).
Compared with the mean yearly sub-aerial extrusion rate of nearly 4 km 3 /yr determined here, the above early attempts appear significant underestimations. Huybers and Langmuir (2009) and Kutterolf et al. (2013) describe significant variations in eruption rates during deglaciation periods, particularly after Termination I when the former authors describe a global eruption rate increase by a factor of 2-6. By contrast, Watt et al. (2013) find that the rate of arc volcanism, which they rank as making up about 90% of present-day sub-aerial volcanism, did not change at statistically significant levels across deglaciation. They conclude that if there was an increase, that was globally limited to non-arc volcanism and it involved at most a twofold increase, and potentially much less than that. As a matter of fact, the time distributions of volcanic eruptions with VEI 7 and 8, characterized by substantial catalogue completeness extending back across a number of ice ages, do not show fluctuations larger than those expected by constant rate Poisson processes (Papale, 2018). The correspondence between computed and observed eruption number rates in Figure 3 suggests a similar conclusion. Thus, the present results tend to agree with Watt et al. (2013), at least to the extent to which global rate changes across interglacial periods are recorded by the largest explosive eruptions.
The global estimates obtained here, together with their accompanying uncertainties, provide a sound reference for other global models of Earth system dynamics, and for further estimates of relevant global quantities. The example above relates to possible existence and extent of feedback mechanisms between deglaciation, sub-crustal melting, and volcanism. One other relevant aspect is the estimate of global volcanic emissions, particularly carbon dioxide which is relevant for global climate change projections. Current estimates of global volcanic CO 2 fluxes account for sources such as fumaroles, diffuse emissions, and persistent volcanic plumes, either non- eruptive or associated with low-level, VEI 0 to 1 eruptions (Burton et al., 2013;Werner et al., 2019). On the contrary, the amount of carbon dioxide discharged during explosive eruptions is poorly known, largely because of too high proximal risks for accurate measurements combined with limits in remote determinations as a consequence of high atmospheric ash loading during eruptions and large CO 2 background concentration in the atmosphere. In their recent review of carbon dioxide emissions from subaerial volcanic regions, Werner et al. (2019) suggest that "explosive" emissions, that is, CO 2 released during eruptive activity not associated with persistent plumes, may be only a small fraction of the total emissions, and estimated it to be, likely, in the range 0.6-7 Mt/yr. More than 90% of the low VEI eruptions in the 0-1 range, which contribute to current estimates of carbon dioxide emissions from volcanoes, are associated with individual eruption volumes below 100 Mm 3 (Papale et al., 2021). The distributions in Figure 6 suggest therefore that such low VEI eruptions may capture only a minor proportion of the global magma discharges. Quantification of the missing CO 2 flux associated with nonpersistent eruptions requires an estimate of the CO 2 contents in the corresponding magmas. Such contents are unfortunately poorly known, largely because glass inclusions within crystals, which trap volatiles dissolved in magmas at the time of crystallization, are poorly informative of the total content of carbon dioxide due to its extremely low solubility and overwhelming partition in the gas phase (Dixon and Stolper, 1995;Papale, 2005). Additionally, glass inclusion studies are mostly available for basaltic magmas, as the corresponding mafic minerals, first of all olivine, are effective in preserving the trapped volatiles over long times, while sialic minerals are not (Métrich and Wallace, 2008). In basaltic magmas, reconstructed total CO 2 contents can often be as large as several wt% (Papale, 2005;Moretti et al., 2018, and references therein). If the mean rate of magma production from medium to large eruptions (individual volumes >100 Mm 3 ) from Figure 6 is compared with the estimated (Werner et al., 2019) 0.6-7 Mt/yr CO 2 flux from "explosive" emissions, then the average proportion of CO 2 in the discharged magma turns out to be as little as 0.014-0.17 wt % (an average density of 1000 kg/m 3 for the explosive volcano deposits is used to convert from volumes to masses). Differently, global mass balance from estimated total carbon dioxide released during recent, well studied explosive eruptions such as the 1980 Mt. St. Helens and 1991 Piñatubo explosive eruptions (10 and 50 Mt CO 2 , respectively, Werner et al., 2019, compared to 10 11.8 and 10 13.1 kg of erupted magma, Crosweller et al., 2012) suggest that the total CO 2 content in magmas was at least 1.6 and 0.4 wt%, respectively. Thus, when compared to eruption volume rate estimates in this work, the current estimates of "explosive" CO 2 emissions appear too low as they imply exceedingly low CO 2 contents in magmas. If 1 wt% total (dissolved plus exsolved, as for the Mt. St. Helens and Piñatubo cases above) CO 2 is roughly representative as an average for magmas discharged during eruptions with VEI 2 or larger, then Table 2 and Figure 6 show that the corresponding mean flux of carbon dioxide would be close to 40 Mt/yr; over 50 years there would be a 90% confidence to observe an average flux in the range 29-52 Mt/yr, with a 1% probability to observe less than 26 or more than 72 Mt/yr, and a 1‰ probability to observe less than 24 or more than 180 Mt/yr. The above CO 2 flux from non-persistent eruptions would be a significant addition to the current estimates of total flux of volcanic CO 2 , with the most recent ones being in the range 180-225 Mt/yr (Fischer and Aiuppa, 2020) including persistent volcanic plumes, fumaroles and diffuse emissions, and dissolution in volcanic lakes and aquifers. More confident estimates of the contribution by non-persistent volcanic activities to the total carbon dioxide release require better knowledge of total carbon dioxide in magmas, which is unfortunately still poorly known.

CONCLUSION
The examples above are illustrative of how modeling global volcanic rates can impact modeling and understanding of global Earth dynamics. A more comprehensive analysis of such impacts is beyond the aims of this work, which is instead dedicated to presenting the accomplishment of a study initiated years ago by considering first the global distribution of volcanism in terms of discrete VEI classes (Papale, 2018), then prosecuted by analyzing the continuous volume distribution (Papale et al., 2021), finally completed here by quantifying the global rates of subaerial volcanism on Earth. We stress here the beauty deriving from the simplicity in the formulation of the present model, that in just three elementary equations reproduces with satisfactory accuracy one relevant aspect of the large-scale Earth dynamics represented by the time-size distribution of the global subaerial volcanism on Earth. This would never be possible without the availability of global eruption databases, the importance of which cannot be overstated. The relevance of databases is destined to increase further in the coming years together with the explosion of machine learning and artificial intelligence techniques, and we foresee and encourage significant resources dedicated to further developments of databases in volcano and Earth system science.
Our analysis refers to volcanic eruptions as they are identified and categorized in the available databases, and to the global scale represented by the whole Earth system (excluding eruptions on the ocean bottom and flood basalts). At such a scale subaerial volcanic eruptions emerge as Poisson events, separated by exponentially distributed inter-event times (Eq. 1). The Poisson distribution implies that at the global scale volcanic eruptions are independent from each other. Such a global scale behavior does not necessarily apply at the individual volcano scale, where the occurrence of an eruption can be related to previous ones. Hence, the statistical distribution of eruptions for an individual volcano may be different from Poisson. Observations at individual volcano scale are normally too few to discern unambiguously a distribution. The individual volcano scale needs further dedicated investigation, and we recommend that the arguments, modeling, findings and conclusions presented here are limited to the global scale to which our study refers.
Some aspects of the model require further understanding. The initial log-normal portion of the overall distribution may either reflect some fundamental difference in the origin and dynamics of small vs. large eruptions, or it may be an artifact due to grouping of small individual eruptions into larger ones covering longer periods of time, as it is normal practice in database construction (Global Volcanism Program, 2013;Papale et al., 2021). The origin of the power law section of the distribution may relate to highly nonlinear dynamics and to a sort of critical behavior which has been described for other natural systems (including faulting and earthquake generation, Fisher et al., 1997;Huang et al., 1998), resulting in the emergence of self-similarity and scale invariance (e.g., Bak et al., 1988a,b;Marković and Gros, 2014). These aspects need specific investigation to be addressed. Similarly, the relationships with the behaviors in limited geographical or geodynamic regions down to individual volcano scale, which are relevant for local hazard forecasts, require further dedicated investigation.

AUTHOR CONTRIBUTIONS
PP conceived the work, developed the theory, made the calculations, and wrote the manuscript. DG checked the calculations and supported in writing the manuscript. WM provided theoretical supervision and supported in writing the manuscript.

FUNDINGS
DG was supported by an EPOS-IT grant.