Characteristics of Particle Size Distributions of Falling Volcanic Ash Measured by Optical Disdrometers at the Sakurajima Volcano, Japan

: In the present study, we analyzed the particle size distribution (PSD) of falling volcanic ash particles measured using optical disdrometers during six explosive eruptions of the Sakurajima volcano in Kagoshima Prefecture, Japan. Assuming the gamma PSD model, which is commonly used in radar meteorology, we examined the relationships between each of the gamma PSD parameters (the intercept parameter, the slope parameter, and the shape parameter) calculated by the complete moment method. It was shown that there were good correlations between each of the gamma PSD parameters, which might be one of the characteristics of falling volcanic ash particles. We found from the normalized gamma PSD analysis that the normalized intercept parameter and mass-weighted mean diameter are suitable for estimating the ash fall rate. We also derived empirical power law relationships between pairs of integrated PSD parameters: the ash fall rate, the volcanic ash mass concentration, the reﬂectivity factor, and the total number of ash particles per unit volume. The results of the present study provide essential information for studying microphysical processes in volcanic ash clouds, developing a method for quantitative ash fall estimation using weather radar, and improving ash transport and sedimentation models.


Introduction
The physical properties and chemical components of ash particles emitted by explosive volcanic eruptions represent useful basic data for investigating the mode, scale, and mechanisms of eruptions. Previous studies have investigated the generation mechanism of individual eruptions [1,2] and the fragmentation process during explosive eruptions by estimating the total grain size distribution following sediment surveys [3][4][5].
In recent years, volcanic sediment data analysis has facilitated many types of research, such as the detection, tracking, and prediction of ash fall to prevent volcanic disasters [6][7][8], estimation of the content fraction of particles smaller than 63 µm (which affect aircraft operation) from total grain size distribution [7], and investigation of variation in tephra features with distance from the volcano crater using the 100-year eruption records of volcanoes along the western coast of North America [8].
Although investigating ground sediments deposited following volcanic eruptions is an important method in volcanology research, and for volcanic disaster prevention, sediment data acquired on the ground are inevitably time-integrated, they are usually accumulated hours or days after an eruption. Large historical eruptions have resulted in centuries of accumulated sediment. Conventional geological methods for investigating

Parsivel 2
In this study, we used the Parsivel 2 (PARticle SIze and VELocity) optical disdrometer; its main specifications are listed in Table 1 [17]. Parsivel 2 measures falling particles as follows. A flat, horizontal laser measuring surface (180 mm × 30 mm × 1 mm) is formed between the light-emitting and light-receiving devices. The diameter (0.2-25 mm) and fall speed (up to 20 ms −1 ) of the particles are calculated in real time based on the voltage drop and duration as falling particles pass through the laser measurement surface (Figure 1). In these calculations, it was assumed that particles < 1 mm in diameter are spherical, the axial ratio of particles 1-5 mm in diameter varies linearly from 1 to 0.7, and the axial ratio of particles > 5 mm in diameter is 0.7. These assumptions are based on the original function of Parsivel, which was designed to measure raindrops, and may lead to errors in measurements of volcanic ash particles. At present, it is difficult to estimate these errors because the relationship between the axis ratio and diameter of volcanic ash particles is unclear.  . Schematic diagram illustrating ash particle diameter and fall velocity measurements by Parsivel 2 (Löffler-Mang and Joss, 2000). As a particle passes through the measurement laser beam sheet, voltage reduction ∆V is generated during time interval ∆t at the light-receiving device. The particle diameter is estimated from ∆V, and fall velocity is estimated from ∆t.
The PSD estimation error for Parsivel has been compared with that for other instruments [15,16,[18][19][20]. The accuracy of the first-generation Parsivel and Parsivel 2 was evaluated by comparing measured data with those from an impact disdrometer observing raindrops [16,21]. The results showed that Parsivel 2 has improved the accuracy of rainfall and PSD measurements of particles 0.5 and 4 mm in size. The accuracy of Parsivel 2 is acceptable for estimating the fall velocity of particles < 1 mm in diameter. Parsivel 2 was found to underestimate fall velocity near 1 m s −1 , but this trend decreased as particle size increased.

Data Collection
The ash particle data analyzed in this study were collected by the Parsivel 2 network (Figure 2), which was installed on Sakurajima and operated by the Disaster Prevention Research Institute of Kyoto University. We extracted PSD data, including eruption onset time and cloud top height, from the Parsivel 2 database for six volcanic eruptions that occurred in 2018 (Table 2), according to reports of Sakurajima volcanic eruptions published by the Kagoshima Local Meteorological Observatory of the Japan Meteorological Agency [22,23].  To confirm that an eruption cloud had passed over a Parsivel 2 site, we used X-band multi-parameter radar (XMP), which was installed approximately 11 km from Sakurajima's Minami-dake summit. XMP is operated by the Ministry of Land, Infrastructure and Transport to monitor debris flows on Sakurajima; although it is used to observe precipitation, it has also been shown to capture ash particles falling from eruption clouds [24].
Eruption cases were selected for analysis as follows. First, eruptions in which the eruption cloud top height was >2000 m above the vent were extracted. Then, XMP data for the extracted eruptions were downloaded from the Extended Rainfall Indicator Network (XRAIN) data download system [25] to draw time-cumulative plan position indicator (PPI) images of reflectivity at 2 elevation angles, of 1.7 • and 6 • . Finally, eruptions in which the installation point of Parsivel 2 was within the cumulative reflectivity distribution area were selected for analysis.

Data Processing
Various physical quantities of particles falling through the atmosphere were calculated from PSD data measured by Parsivel 2 [26]. The PSD of falling particles is calculated as follows: where N(D i ) is the number of particles from D i to D i + ∆D i in a unit volume, and D i , ∆D i are the average and bin spacing of the i-th size bin, respectively. C ij is the number of particles measured in the i-th size bin and the j-th speed bin, and nd is the number of size bins (32 in this study). V j is the measured time interval (60 s for this instrument), and A i is the fall rate measured for the j-th speed bin. The effective measurement area of the i-th size bin is calculated as follows: The mean fall velocity of particles with a mean diameter D i is calculated as follows: where nv is the number of velocity bins (32 for this instrument). The ash fall rate R A was calculated as follows: where ρ p is the solid density of ash particles. We assumed that ρ p = 2.5 g cm −3 (= 2.5 × 10 3 kg m −3 ). From Equation (1), If R A is expressed in [mm h −1 ], where ρ b is the bulk density of ash deposits. We assumed that ρ b = 1.25 g cm −3 (= 1.25 × 10 3 kg m −3 ). The volcanic ash mass concentration C A is defined as the amount of ash particles in a unit volume of air, is calculated as follows: and the reflectivity factor Z is calculated as follows: Atmosphere 2021, 12, 601 6 of 24 The mass-weighted diameter D m is a parameter characterizing the normalized PSD, and is calculated as the ratio of the fourth moment of PSD to the third moment, as follows: D m is closely related to the median volume particle diameter D 0 (Sections 3.2 and 3.4). The normalized PSD intercept parameter N w is calculated as follows: The basic parameters of volcanic ash particles described above were directly calculated from PSD data measured by Parsivel 2 .

Radar Meteorological Approach
In radar meteorology, precipitation PSD is indispensable for quantitative estimation of precipitation and hydrometeor classification. Therefore, functions to approximate observed precipitation particles have been proposed in several previous studies. Functional forms of PSD can be used to investigate spatiotemporal variation in the microphysical processes of precipitation clouds, derive theoretical formulae for quantitative precipitation estimates, and simulate precipitation particle scattering.
In this study, we applied radar meteorology concepts and techniques to volcanic ash fall phenomena. In this section, we describe volcanic ash PSD using either gamma or normalized gamma functions. In addition, we show that the integrated PSD parameters such as R A , C A , Z, etc. can be expressed in terms of either gamma or normalized gamma PSD models.

Gamma PSD Model
Raindrop size distribution is commonly described using the gamma PSD model. In this study, we replaced raindrops with falling ash particles in the gamma PSD model, which is expressed as follows [27,28]: where N(D) is the number of ash particles per unit diameter increment within a unit volume, is the intercept parameter characterizing PSD, Λ (mm −1 ) is the slope parameter, and µ is the shape parameter. D is expressed in mm, and µ, which is normally dimensionless, takes the units of N 0 . Thus, N 0 and µ are not completely independent, such that N 0 is affected by the variation of µ and N 0 has no physical meaning [29]. Equation (11) can be expressed using the median volume diameter D 0 , which has a physical meaning and is used in radar meteorology instead of Λ. Thus, When D max /D 0 ≥ 2.5, and µ ≥ −3, the following approximation holds, with an accuracy of ≤0.5% [27]: G = ΛD 0 = 3.67 + µ. The median volume diameter D 0 is defined as the diameter that divides the total volume of all falling ash particles within a unit volume into two equal parts in the falling ash PSD. Since the computation of D 0 is cumbersome, the mass-weighted mean D m is often used instead of D 0 . The relationship between D m and D 0 is as follows [27]: Atmosphere 2021, 12, 601 The value 3.82 in the above equation is attributed to a power law approximation of the fall velocity of ash particles (D ≤ 2 mm; Table 3). The value 3.67 is typically used for raindrops. From Equation (13), when µ ≥ 0, the difference between D m and D 0 is less than 9%, and as µ becomes larger, this difference becomes smaller; when µ = 10, the difference is about 2%. Using D m to represent the gamma PSD, the model becomes: The intercept parameter Λ can be expressed as follows: As a special case of the gamma distribution, µ = 0, or This PSD model is a well-known exponential PSD function [30]. Table 3. List of coefficients and exponents of ash particle fall velocity formulae proposed in previous studies

Normalized Gamma PSD Model
Generally, falling ash PSDs fluctuate depending on C A or R A . To investigate the shape of the ash PSD, which is not affected by C A , normalized (scaled) PSD concepts have been proposed [35][36][37][38][39][40]. The normalized PSD is generally expressed as follows [39]: where N w is the normalized intercept parameter, which normalizes N (D), and D m is the mass-weighted mean diameter, which normalizes D. F is a factor that determines the geometry of the PSD, and X = D/D m . The normalized PSD expressed in Equation (18) is not a functional form of PSD. In this section, we derive a normalized PSD assuming a gamma function for measured PSD. When N(D) is represented by the gamma function, i.e., Equation (12), the nth moment M n of PSD is expressed as follows: Atmosphere 2021, 12, 601 8 of 24 D m and N w are expressed as follows, in terms of M n : From Equations (14), (18), (20), and (21), Therefore, the normalized gamma PSD can be represented as follows: Importantly, the normalized gamma PSD described in Equation (23) is represented by three independent parameters, N w , D m , and µ.

Calculation of PSD Parameters
Three parameters of the gamma PSD, N 0 (or N w ), Λ (D 0 or D m ), and µ can be calculated using the nth moment M n of the PSD [41,42], as follows: M n can be calculated by the following equation using PSDs measured by Parsivel 2 : For convenience of computation, we introduce a variable η, which is defined as the combination of M 2 , M 4 , and M 6 , as follows: such that η is a monotonically increasing function with respect to µ. By solving Equation (26), the shape parameter µ of the gamma PSD model can be obtained from the following equation: If µ is known, the slope parameter Λ of the gamma profile is given by [42]: . N 0 is calculated as follows: where N 0 can be obtained by specifying any value of n if µ and Λ are obtained. Thus, n can be 2, 4, or 6, and the difference in the calculated value of N 0 for each case is about 1%. D 0 can be computed as follows if µ and Λ are obtained: The present paper utilizes the complete moment method to calculate the gamma PSD parameters. However, it may be necessary to mention that the calculation of gamma PSD parameters is more complex than the complete moment method when the lower and upper bounds of the integral are finite in Equation (24), i.e., an incomplete gamma distribution. The effects of D min and D max on the calculated gamma PSD parameters have been examined assuming an incomplete gamma distribution [42]; the effect of D min on the gamma PSD parameters is small, such that D min = 0 may be assumed [43,44]. The present paper use the observed D max to estimate the nth moment of the PSD.

Integrated PSD Parameters
The integrated PSD parameters include the total number of falling ash particles N T , C A , R A , and Z. The integrated parameters are expressed by PSD moments, as follows: where v(D) is the falling velocity of volcanic ash particles of diameter D, calculated as follows: The integrated PSD parameters defined in Equations (31)-(33) are determined by PSD moments calculated from the measured PSD data, and the relationships among integrated PSD parameters can be obtained by regression analyses of the integrated PSD parameters. If the PSD data are not available, then the relationships among integrated PSD parameters are derived by assuming the functional form of the PSD. Assuming the gamma PSD model, the following relationships are derived: The gamma PSD parameters obtained by the full moment method in Section 3.4 may then be substituted into these equations.

Fall Speed of Volcanic Ash Particles
The terminal velocity of falling particles is generally determined according to the balance among gravitational, buoyancy, and aerodynamic drag forces [45]. In the present study, we express the particle fall velocity as a power law function of particle diameter, as shown in Equation (33), which is obtained by applying regression analysis to fall velocity data measured using Parsivel 2 , as shown in Equation (3).
The distribution of the number of ash particles in the parameter space of fall velocity and diameter for the six selected volcanic eruptions is shown in Figure 3. The total number of samples was 63,237, and most of the observed data were for ash particles, i.e., D ≤ 2 mm. Particles with diameters ranging from 0.4 to 0.8 mm were the most frequent. The fall velocities corresponding to these particles were widely distributed, from 0.2 to 4 ms −1 , with larger particles occurring less frequently but with greater variance. This variation may be caused by Parsivel 2 measurement error or particle density variation. The falling velocity of particles increased in proportion to the square of the particle density [45]. Another reason for the large scatter observed in Figure 3 may be the particle shape, which determines the drag force acting on the particle [10]. The flatter the particle shape, the slower the falling velocity. The largest particle observed in the present study was 4 mm in diameter. The power law equation for the fall velocity of all data shown in Figure 3 is: For volcanic ash particles (D ≤ 2 mm), the fall velocity equation is: which is quite similar to Equation (40). Table 3 summarizes the fall velocity equation derived in the present study and those proposed in previous studies. The equations for fall velocity versus diameter listed in Table 3 are shown in Figure 4. Among all curves shown in Figure 4, the formulae derived from experiments of ash particle free fall from a height of 17 m [10,31] give the largest fall velocity for the same diameter, likely because the fall velocity was accelerated by the fingering phenomenon [46,47]. Parsivel 2 , installed on Sakurajima, was used to observe the falling velocity of actual ash particles. [11,32]; the falling velocity of ash particles for a total of 76 eruptions from 2014 to 2016 was described as a power law (Figure 4) derived from observed data [33] and previously developed theoretical equations [11,34,48]. The power law formulae obtained in the present study provided the lowest fall velocity for the same diameter, and smaller than with the theoretical equations [34] in which the particle density was assumed to be 2500 (kg m −3 ).

Temporal Changes in PSD
Changes in the PSDs of the six selected eruptions over time, as eruption clouds passed over the Parsivel 2 device, are shown in Figure 5. The maximum R A (11.7 mm h −1 ) was observed during eruption 4. The number of particles observed per unit volume (N T ) ranged from 10 to 10 3 m −3 . R A and log 10 N T were somewhat correlated. Particles with a maximum diameter of 3.75 mm were observed in eruptions 2b and 3. D max was not necessarily correlated with R A . In all eruptions except eruption 1, large values of D max were observed immediately after the onset of ash fall ( Figure 5). D max and D m were positively correlated, and D max was negatively correlated with both log 10 N T and log 10 N w . These new findings on the PSD of ash particles will be examined in detail in the following sections.
A temporal change in the PSD at a certain point is considered to reflect a spatial change in PSD within the eruption cloud. Therefore, it is necessary to know which part of the eruption cloud was observed by Parsivel 2 when discussing a temporal change in the PSD. The distribution of the time-integrated reflectivity factor in eruption clouds during eruption 6 is shown in Figure 6. We also obtained the trajectories of the centroids of the eruption clouds (data not shown). Calculations of time-integrated radar reflectivity and the trajectories are based on PPI data at an elevation angle of 1.7 • . Parsivel 2 data for ash particles in the central part of eruption cloud as it passed over the observation site are shown in Figure 5. Unfortunately, there is no Parsivel 2 site in the area where the accumulated reflectivity factor reached a maximum, since this area is 2 km from the vent and is therefore restricted.

Gamma PSD Parameters
Temporal changes in the observed PSD profiles for eruption 6 are shown in Figure 7, as a representative example. Two fitting curves of the gamma PSD function are superimposed on each set of observed PSD data; one curve was obtained by the momentum method, and the other curve, by non-linear regression analysis. It should be noted that non-linear regression analysis must be carried out to the logarithmic form of the gamma PSD formula because N(D) distributes over a wide range (from 10 0 to 10 4 ). A good fit was obtained for all samples, except at 15:54 LST, when large particles were observed. The gamma PSD obtained by the momentum method was approximately exponential, while the gamma PSD obtained by non-linear regression analysis exhibited an upward convex shape that fitted small particles. It should be noted that the results obtained by the two methods employed agree quite well for PSDs at 16:00 LST to 16:04 LST, a period during which only particles smaller than 2 mm were observed. The results suggest that the difference in the gamma PSD functional forms obtained by non-linear regression analysis and the momentum method is acceptable. We analyzed all the PSD data obtained in order to quantitatively confirm the results, which are shown in the discussion section of the present paper.  (380), respectively. It interesting that each of three parameters (log 10 N 0 , Λ, µ) had a wider range compared with that of precipitation particles. Other descriptive statistics, such as the maximum, median, and skewness, are summarized in Table 4. Table 4. Statistical values of gamma, normalized, and integrated particle size distribution (PSD) parameters. ρ b = 1.25 × 10 3 kg m −3 and ρ p = 2.50 × 10 3 kg m −3 were assumed for the estimation of R A .  The main objective of radar remote sensing is determining the physical quantities of eruption clouds from measured radar parameters. Normally, the number of physical quantities, which are unknowns, is larger than the number of parameters that can be measured by radar. Therefore, the number of unknowns is reduced by theoretically or empirically defining relationships between unknown parameters. If two parameters are strongly correlated, then the obtained relationship will be incorporated into a weather radar monitoring system for eruption clouds. Theoretical interpretations of such correlations are useful for understanding the microphysical processes of eruption clouds.

Parameter
Clear correlations were detected between gamma PSD parameters in this study (Figure 9), which suggests that the gamma PSD model well describes the observed volcanic ash PSDs. The relationships between gamma PSD parameters are summarized in Table 5. Thus, if one gamma PSD parameter is obtained, the remaining parameters can be estimated using the information listed in Table 5. A scatter plot of D max and µ is also shown in Figure 9; the shape parameter µ approaches 0 (i.e., exponential PSD) when D max >~3, and becomes large (i.e., upward concave PSD) when D max is small. Thus, once D max is given for an initial condition of the PSD, µ is determined by the µ-D max relationship. Next, N 0 and Λ are estimated by the N 0 -µ and Λ-µ relationships, respectively.   Figure 10 is a scatter plot showing PSDs before and after normalization. The normalized PSD parameters (N w and D m ) were calculated directly from the PSD data measured by Parsivel 2 , without assuming a functional form of the PSD. Figure 10a shows the PSD before normalization, which is characterized by wide scattering due to variation in ash fall strength, such that the characteristic distribution could not be observed. After normalization, the PSDs converged to reduce the scatter (Figure 10b), allowing us to observe the form of the distribution. The shape parameter µ was determined by least-squares fitting of the measured N(D)/N w data according to Equation (23). It may be interesting to examine the associations between N w or D m with the type and magnitude of volcanic eruption in a future study, using larger sample sizes. Despite the small sample sizes in the present study, they were useful for determining the statistical characteristics of normalized gamma PSD parameters. Figure 11 shows the frequency distributions and probability density functions of the two parameters (D m and log 10 N w ) used to normalize the PSD. The shape parameter µ is shown in Figure 8d. The modes (SD) of D m and log 10 N w of the normalized gamma PSD parameters were 0.656 (0.350) and 4.60 (0.704), respectively, and the mode (SD) of µ was 8.31 (13.7).

Quantitative Ash Fall Estimation Using Normalized PSD Parameters
Next, we examined correlations between R A and gamma PSD parameters (N w , D m , and µ). Figure 12 shows a scatter plot of R A /N w and D m . Applying nonlinear regression analysis, we obtained the following power law equation: By rearranging Equation (42), we obtained: The power law equations for the upper and lower 95% confidence limits of Equation (43) are as follows: R A = 2.72 × 10 −4 D m 5.07 N w ; 95% upper limit (44) R A = 2.36 × 10 −4 D m 4.87 N w ; 95% lower limit (45) In Figure 12, the correlation coefficient (root-mean-square error) between R A /N w and D m was 0.991 (0.000177), which was higher (lower) than that between R A and Z, as described in the next section. This result suggests that Equation (42) provides an advanced method for estimating R A if N w and D m are estimated with a sufficient degree of accuracy. For precipitation phenomena, two-frequency radar observations based on satellite remote sensing [49] and dual polarization radar observation [50] have been proposed to estimate N w and D m . Similar techniques may be available for ash fall phenomena, and should be explored in a future study.
As physical quantities that represent the microphysical processes within precipitation clouds, logN w and D m have also been used to distinguish precipitation types (stratiform or convective) [51,52]. Similarly, D m and N w are expected to provide new insight into volcanic eruption clouds and their microphysical processes in future studies [53][54][55]. Figure 13 shows the frequency and probability density distributions of the three integrated PSD parameters (R A , C A , and 10log 10 Z). In this study, the modal values (SD) of R A , C A , and 10log 10 Z of the integrated PSD parameters were 0.52 (2.17), 0.0931 (0.189), and 17.1 (7.48), respectively. Next, we investigated the correlations among these integral parameters. In particular, the relationship between R A and Z is crucial as a practical equation for estimating the R A from weather radar observations. The C A -Z relationship is used to estimate the C A of volcanic ash. These relationships were derived from regression analyses of observed data for each eruption, and are summarized in Table 6. A scatter plot of R A and Z values obtained from Parsivel 2 observations for all eruptions is shown in Figure 14. The R A -Z relationship was defined using the nonlinear least-squares method, as follows:

Conventional Relationship for Quantitative Ash Fall Estimation
The R A -Z relationships for the upper and lower 95% confidence limits are as follows: R A = 14.6 × 10 −2 Z 0.383 ; 95% lower limit (48) It must be mentioned that the R and RMSE for all eruptions are not as good as those of an individual eruption because the PSDs are changeable depending on the eruption case, which is similar to the variation of PSD in precipitation. Table 6. R A -Z and C A -Z relationships obtained from regression analyses of PSD data measured by Parsivel 2 . Note that R A (mm h −1 ), C A (g m −3 ), Z (mm 6 m −3 ). R: correlation coefficient, RMSE: root mean square error. ρ b = 1.25 × 10 3 kg m −3 and ρ p = 2.50 × 10 3 kg m −3 were assumed for the estimation of R A .

Discussion
In radar meteorology, various mathematical formulas for 'precipitation particles' have been proposed by researchers. These include exponential [30], gamma [27,28], lognormal [56,57], Poisson [58], and Weible [59]. The present study assumed a gamma PSD model for 'volcanic ash particles', for the following two reasons. First, gamma PS models have been used widely in previous radar meteorology studies, because the model is relatively simple and its error structure has been investigated [43]. Second, we can compare the results of our analysis with those of previous studies using the same PSD model. The differences in gamma PSD parameters between precipitation particles and volcanic ash particles, if any, can be used to discriminate eruption clouds from precipitation clouds, one of the goals of monitoring volcanic ash falls using radar.
However, when the gamma PSD model is applied to volcanic ash particles, one point remains unclear. 'Does the gamma PSD model represent observed volcanic ash PSD?' This question results from the fact that validation studies have been limited due to the lack of observed ash PSD. In this section, we discuss the validity of the gamma PSD model of volcanic ash using PSD data obtained by Parsivel 2 . First, we examined the validity of the momentum method used to calculate the gamma PSD. As we have already shown in Figure 7 in this paper, the gamma N(D) obtained by the momentum method agreed with that obtained by non-linear regression analysis. However, as the number of analyzed PSD samples was limited in Figure 7, we examined all PSD samples obtained in six eruption cases. It should be mentioned again here that non-linear regression analysis must be executed not to N(D) but to log 10 N(D). If we used linear N(D) for the analysis, the obtained gamma PSD regression curve would represent only small particles that have large N(D) values. Figure 15 shows comparisons of gamma PSD parameters (N 0 , µ, Λ) calculated by the momentum method and non-linear regression analysis. Regression lines obtained were as follows.
Log 10 N 0 _mom = 0.934 × log 10 N 0 _reg + 1.47, µ_mom = 0.904 × µ_reg + 2.14, where, the subscripts 'mom' and 'reg' in the variables mean the momentum method and regression analysis, respectively. The correlation coefficient R (the root mean square error RMSE) for N 0 , µ, and Λ were 0.962 (2.23), 0.940 (4.02), and 0.964 (5.38), respectively. We can confirm from these statistical values that the gamma PSD parameters estimated by the moment method agreed with those obtained by non-linear regression analysis.
Atmosphere 2021, 12, x FOR PEER REVIEW 21 of Next, we evaluated the validity of the gamma PSD model itself, which is in tu based on the complete momentum method. The result is shown in Figure 16a. Compar N(D)_g based on the gamma model with N (D)_o observed by a Parsivel 2 , we obtained Next, we evaluated the validity of the gamma PSD model itself, which is in turn based on the complete momentum method. The result is shown in Figure 16a. Comparing N(D)_g based on the gamma model with N (D)_o observed by a Parsivel 2 , we obtained N (D)_g = 0.917 × N (D)_o + 0.175. The R and RMSE were 0.825 and 0.554, respectively. According to Figure 16b, most of the residuals of N (D)_g were within ±10% of their values. Given that the value of N (D) ranged from 10 0 to 10 4 , these statistical values confirmed that the gamma PSD model is appropriate for describing volcanic ash PSD. Next, we evaluated the validity of the gamma PSD model itself, which is in turn based on the complete momentum method. The result is shown in Figure 16a. Comparing N(D)_g based on the gamma model with N (D)_o observed by a Parsivel 2 , we obtained N (D)_g = 0.917 × N (D)_o + 0.175. The R and RMSE were 0.825 and 0.554, respectively According to Figure 16b, most of the residuals of N (D)_g were within ±10% of thei values. Given that the value of N (D) ranged from 10 0 to 10 4 , these statistical value confirmed that the gamma PSD model is appropriate for describing volcanic ash PSD.

Summary
In the present study, we examined the characteristics of the PSD of volcanic ash particles and the relationships among PSD parameters. We used a total of 166 PSD samples collected by Parsivel 2 during six explosive eruptions of the Sakurajima volcano in 2018. The PSD data were analyzed using methods developed for use in radar meteorology.
We proposed gamma and normalized PSD models to describe the PSDs of volcanic ash particles. The observed PSD data were well-described by both models. Strong correlations among the gamma PSD parameters (log 10 N 0 , Λ, and µ) were found and relationships between the gamma PSD parameters were derived. Interestingly, the µ-R max relationship showed that µ changes inversely with R max . When D max was large, µ approached zero; i.e., the PSD became exponential. These results could be applied to set up initial PSD conditions for a volcanic ash transport and diffusion model. It must be noted that the range of each of the gamma PSD parameters (log 10 N 0 , Λ, and µ) of ash particles was wider than that of precipitation particles.
The relationships between gamma PSD parameters can be used to solve the inverse problem in radar meteorology, i.e., the retrieval of PSD parameters from weather radar observations. In precipitation studies, D 0 (or Λ) can be estimated from the differential reflectivity Z DR , which is measured by polarimetric radar, and N 0 is estimated from D 0 and observed Z, assuming that µ = 0 [60]. This method is based on the relationship between raindrop shape and diameter, such that flatter raindrops are associated with larger raindrop diameters. Recently, volcanic ash particle shape has been studied using 2DVD, which was designed to measure the shape of a falling raindrop. The classification of ash particle shapes can then be used to identify the volcanic eruption type [61]. To date, a clear relationship between particle shape and diameter has not been found for volcanic ash particles. Further studies of volcanic ash particle shape are required to establish a method for weather radar retrieval of ash particle PSD.
In this study, we also proposed a normalized PSD model, which does not require consideration of the PSD variation associated with ash fall intensity. We found a clear correlation between R A /N w and D m , and derived its power law functional relationship. Once the two parameters N w and D m are retrieved from radar measurements, R A can be estimated from the R A /N w -D m relationship proposed in this study. However, further theoretical and observational studies are needed to retrieve parameters N w and D m ; for example, scattering simulations based on observed PSD could be used to establish theoretical relationships between N w and D m , and polarimetric radar parameters such as Z DR and K DP . Polarimetric radar observations of volcanic ash particles are necessary to validate such theoretical relationships.
Finally, we proposed conventional formulae for estimating the R A and C A . The R A -Z relationship was derived through regression analysis and the C A -Z relationship was calculated from measured PSD data. Although these relationships have large errors associated with instantaneous monitoring of ash fall caused by spatiotemporal variation in PSD, they are convenient and applicable for radar ash fall monitoring. It should be noted that the relationships obtained in this study were derived from a limited number of data. Thus, further analysis using larger PSD datasets is necessary.

Symbol Description Unit A i
Effective measured area of the i-th size bin mm 2 C A Volcanic ash mass concentration kg m −3 C ij Number of particles measured in i-th diameter bin and j-th Velocity bin -D Particle diameter mm D 0 Median volume diameter mm D m Mass-weighted mean diameter mm D max Maximum particle diameter mm