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

Abstract

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 reflectivity 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.

1. 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 sediments are limited in that the properties and sedimentary environment of volcanic ash particles change depending on weather conditions, such as wind and rain, during and after sedimentation [9], such that sediment data are often unsuitable for remote sensing and numerical forecasting of ash fall distribution and transportation, which change from moment to moment.

Recent studies have attempted to clarify the physical properties of ash particles at high temporal resolution using meteorological optical disdrometers, which were developed to measure the physical characteristics of precipitation particles [10,11,12] including particle size distribution (PSD), shape, falling velocity, and canting. Optical disdrometers have also been used to study microphysics in precipitation clouds, to develop quantitative precipitation estimates, and to provide ground truth data for satellite remote sensing of precipitation. Because the optical disdrometer is a non-contact sensor, it can measure falling volcanic ash particles without altering local environmental conditions. A certain type of optical disdrometer can acquire the shape of a falling single particle using its fast scanning frequency, which is typically a few thousands per second.

Free-fall experiments conducted in a large experimental facility measured the falling velocity, shape, and canting angles of ash particles using a two-dimensional video disdrometer (2DVD) [10], which was developed to observe precipitation particles [13,14]. Another type of optical disdrometer, Parsivel² [15,16], was installed on Sakurajima, a volcano in Kagoshima Prefecture, Japan to automatically measure falling ash particles [11]. This instrument measured temporal changes in the PSD and fall velocity of tephra from a total of 76 eruptions of Sakurajima over a 2-year period. An analysis of these records showed that temporal changes in the Parsivel² PSD data were synchronized with the crustal movements associated with eruptions. Based on the success mentioned above, automated measurements of ash particles are ongoing with 13 Parsivel² instruments that have been arranged in all directions centering on the Minami-dake summit of Sakurajima [12].

The goal of the present study was to find characteristics of volcanic ash PSD. One of the unique features of the present study is an analysis of PSD data measured by disdrometers every 1 min. The results obtained will be useful for the radar monitoring of ash falls and improvement of volcanic ash transport and diffusion models.

The structure of the present paper is as follows. In Section 2, we describe the principles by which ash particles are measured by Parsivel², the data analysis, and processing methods. In Section 3, we describe the functional forms and integrated parameters of ash particle PSD, which provide a theoretical basis for quantitative ash fall estimation and analysis of volcanic ash cloud microphysics. Our results are described and discussed in Section 4 and Section 5, respectively. The phenomenon of fragmentation, which is an important eruption mechanism based on the physical and chemical characteristics of volcanic ash particles, is beyond the scope of this study, and will not be discussed.

2. Measurements of Volcanic Ash Particles

2.1. Parsivel²

In this study, we used the Parsivel² (PARticle SIze and VELocity) optical disdrometer; its main specifications are listed in

Table 1. Main specifications of the Parsivel² disdrometer, according to the instruction manual (OTT, 2011).

Category	Subcategory	Specification
Sensor	Type	Laser diode
	Wavelength	650 nm
	Output rating	0.2 mW
	Laser class	1 (IEC/EN 60825-1: 2014)
Sampling	Light strip (W × D)	30 × 1 mm (54 cm2)
	Area (W × D)	180 × 30 mm (54 cm2)
	Interval	10 s to 60 min
Output	Particle size	32 classes: 0.2, …, 8 mm (liquid particles) 32 classes: 0.2, …, 25 mm (solid particles)
	Fall velocity	32 classes: 0.2, …, 20 m s−1
	Precipitation intensity (accuracy)	0.001, …, 1200 mm h−1 (±5% for liquid, ±20% for solid)
	Precipitation type	8 types (drizzle, drizzle/rain, rain, snow, mixed rain/snow, snow grains, sleet, hail)
	Radar reflectivity factor	−9.999, …, 99.999 dBZ ± 20%
	Kinetic energy	0, …, 999.999 J/(m2 h)
	Visibility in precipitation	0, …, 20,000 m

[17]. Parsivel² 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⁻¹) 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.

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² was evaluated by comparing measured data with those from an impact disdrometer observing raindrops [16,21]. The results showed that Parsivel² has improved the accuracy of rainfall and PSD measurements of particles 0.5 and 4 mm in size. The accuracy of Parsivel² is acceptable for estimating the fall velocity of particles < 1 mm in diameter. Parsivel² was found to underestimate fall velocity near 1 m s⁻¹, but this trend decreased as particle size increased.

2.2. Data Collection

The ash particle data analyzed in this study were collected by the Parsivel² 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² database for six volcanic eruptions that occurred in 2018 (

Table 2. List of Sakurajima volcanic eruptions analyzed in the present study.

Case	Eruption Onset (LST)	Cloud Top Height (m)	Movement	Disdrometer Code
1	08:08, 15 May 2018	2500	Top	NABE
2	19:48, 22 May 2018	Unknown	Unknown	HIKP, HART
3	08:01, 30 May 2018	2500	Eastward	KURP
4	11:35, 10 June 2018	3500	Top	HART
5	07:19, 16 June 2018	4700	Westward	SBTT
6	15:38, 16 July 2018	4600	Westward	AKAM

), 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² 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² was within the cumulative reflectivity distribution area were selected for analysis.

2.3. Data Processing

Various physical quantities of particles falling through the atmosphere were calculated from PSD data measured by Parsivel² [26]. The PSD of falling particles is calculated as follows:

$$N(D_i)={\displaystyle \sum _{j=1}^{nd}\frac{C_{ij}}{A_i\cdot \Delta t\cdot V_j\cdot \Delta D_i}};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}N(D_i)\hspace{0.17em}[{\mathrm{mm}}^{-1}{\mathrm{m}}^{-3}],\hspace{0.17em}{D}_i[\mathrm{mm}]$$

(1)

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:

$$A_{\mathrm{i}}=\left[180\times (30-0.5D_i)\right]\times {10}^{-6}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}[{\mathrm{m}}^2]$$

(2)

The mean fall velocity of particles with a mean diameter D_i is calculated as follows:

$$V(D_{\mathrm{i}})={\displaystyle \sum _{j=1}^{nv}{V}_j\cdot C_{ij}}/{\displaystyle \sum _{j=1}^{nv}{C}_{ij}}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}[\mathrm{m}\hspace{0.17em}{\mathrm{s}}^{-1}]$$

(3)

where nv is the number of velocity bins (32 for this instrument). The ash fall rate R_A was calculated as follows:

$$R_{\mathrm{A}}=\frac{\pi {\rho }_{\mathrm{p}}}{6}{\displaystyle \sum _{i=1}^{nd}(V(D_i)\cdot N(D_i)\cdot }{D}_i^3);\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{R}_{\mathrm{A}}\hspace{0.17em}[\mathrm{kg}\hspace{0.17em}{\mathrm{m}}^{-2}{\mathrm{s}}^{-1}]$$

(4)

where ρ_p is the solid density of ash particles. We assumed that ρ_p = 2.5 g cm⁻³ (= 2.5 × 10³ kg m⁻³). From Equation (1),

$$R_{\mathrm{A}}=\frac{\pi }{6}\hspace{0.17em}{10}^{-9}{\rho }_{\mathrm{p}}{\displaystyle \sum _{i=1}^{nd}{\displaystyle \sum _{j=1}^{nv}{D}_i^3\frac{C_{ij}}{A_i\hspace{0.17em}\Delta t}}};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{R}_{\mathrm{A}}\hspace{0.17em}[\mathrm{kg}\hspace{0.17em}{\mathrm{m}}^{-2}{\mathrm{s}}^{-1}]$$

(5)

If R_A is expressed in [mm h⁻¹],

$$R_{\mathrm{A}}=6\pi \hspace{0.17em}{10}^{-4}\frac{{\rho }_{\mathrm{p}}}{{\rho }_{\mathrm{b}}}{\displaystyle \sum _{i=1}^{nd}{\displaystyle \sum _{j=1}^{nv}{D}_i^3\frac{C_{ij}}{A_i\hspace{0.17em}\Delta t}}};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{R}_{\mathrm{A}}\hspace{0.17em}[\mathrm{mm}\hspace{0.17em}{\mathrm{h}}^{-1}]$$

(6)

where ρ_b is the bulk density of ash deposits. We assumed that ρ_b = 1.25 g cm⁻³ (= 1.25 × 10³ kg m⁻³). 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:

$$C_{\mathrm{A}}=\frac{\pi }{6}{10}^{-6}{\rho }_{\mathrm{p}}{\displaystyle \sum _{i=1}^{nd}{\displaystyle \sum _{j=1}^{nv}{D}_i^3\frac{C_{ij}}{A_i\hspace{0.17em}\Delta t\hspace{0.17em}{V}_j}}};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{C}_{\mathrm{A}}\hspace{0.17em}[\mathrm{kg}\hspace{0.17em}{\mathrm{m}}^{-3}]$$

(7)

and the reflectivity factor Z is calculated as follows:

$$Z={\displaystyle \sum _{i=1}^{nd}{\displaystyle \sum _{j=1}^{nv}{D}_i^6}\frac{C_{ij}}{A_i\hspace{0.17em}\Delta t\hspace{0.17em}{V}_j};\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}Z\hspace{0.17em}[{\mathrm{mm}}^6{\mathrm{m}}^{-3}]}$$

(8)

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_{\mathrm{m}}=\frac{{\displaystyle \sum _{i=1}^{nd}N(D_i)D_i^4\Delta D_i}}{{\displaystyle \sum _{i=1}^{nd}N(D_i)D_i^3\Delta D_i}}; D_{\mathrm{m}} [\mathrm{mm}]$$

(9)

D_m is closely related to the median volume particle diameter D₀ (Section 3.2 and Section 3.4). The normalized PSD intercept parameter N_w is calculated as follows:

$$N_w=\frac{4^4}{\pi {\rho }_{\mathrm{p}}}\left[\frac{C_{\mathrm{A}}}{D_{\mathrm{m}}^4}\right]; N_{\mathrm{w}} [{\mathrm{mm}}^{-1}{\mathrm{m}}^{-3}]$$

(10)

The basic parameters of volcanic ash particles described above were directly calculated from PSD data measured by Parsivel².

3. Functional Representation of Falling Ash Particles

3.1. 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.

3.2. 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]:

$$N(D)=N_0{D}^{\mu }\exp (-\mathsf{\Lambda }D)$$

(11)

where N(D) is the number of ash particles per unit diameter increment within a unit volume, N₀ (mm^−1−μ m⁻³) is the intercept parameter characterizing PSD, Λ (mm⁻¹) is the slope parameter, and μ is the shape parameter. D is expressed in mm, and μ, which is normally dimensionless, takes the units of N₀. Thus, N₀ and μ are not completely independent, such that N₀ is affected by the variation of μ and N₀ has no physical meaning [29]. Equation (11) can be expressed using the median volume diameter D₀, which has a physical meaning and is used in radar meteorology instead of Λ. Thus,

$$N(D)=N_0{D}^{\mu }\exp [-G(D/D_0)]$$

(12)

When D_max/D₀ ≥ 2.5, and μ ≥ −3, the following approximation holds, with an accuracy of ≤0.5% [27]: G = ΛD₀ = 3.67 + μ. The median volume diameter D₀ 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₀ is cumbersome, the mass-weighted mean D_m is often used instead of D₀. The relationship between D_m and D₀ is as follows [27]:

$$D_{\mathrm{m}}=(4+\mu )/(3.82+\mu )D_0$$

(13)

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. List of coefficients and exponents of ash particle fall velocity formulae proposed in previous studies. V_t (D) = a D^b, where V_t is expressed in m s⁻¹ and D is expressed in mm.

a	b	Note	Reference
5.90	0.530	0.25 ≤ D < 4 mm (free fall experiments, Parsivel2)	Miwa et al. [31]
5.96	0.530	0.25 ≤ D < 4 mm (free fall experiments, 2DVD)	Suh et al. [10]
4.04	0.830	0.25 ≤ D < 2 mm (Sakurajima, Parsivel2, 26,075 particles)	Kozono et al. [32]
3.96	0.748	0.25 ≤ D < 4 mm (Sakurajima, Parsivel2, 79,170 particles)	Kozono [33]
3.18	0.728	0.25 ≤ D < 4 mm (Sakurajima, Parsivel2, 63,237 particles)	Present study
3.14	0.817	0.25 ≤ D < 2 mm (Sakurajima, Parsivel2, 63,237 particles)	Present study
3.42	0.371	0.2 ≤ D ≤ 6.0 mm (Mt. St. Helens, 19 March 1982)	Harris and Rose [34]
5.02	0.335	0.22 ≤ D ≤ 1.3 mm (Mt. St. Helens, 19 March 1982)	Harris and Rose [34]
6.87	1.0	0.02 ≤ D ≤ 0.38 mm (Mt. St. Helens, 18 May 1980)	Harris and Rose [34]
4.33	0.437	0.02 ≤ D ≤ 6.0 mm (Mt. St. Helens, all data)	Present study (average)

). The value 3.67 is typically used for raindrops. From Equation (13), when μ ≥ 0, the difference between D_m and D₀ 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:

$$N(D)=N_0{D}^{\mu }\exp [-(4+\mu )(D/D_{\mathrm{m}})]$$

(14)

The intercept parameter Λ can be expressed as follows:

$$\mathsf{\Lambda }=(4+\mu )/D_{\mathrm{m}}$$

(15)

As a special case of the gamma distribution, μ = 0,

$$N(D)=N_0\exp (-\mathsf{\Lambda }D)=N_0\exp [-G(D/D_0)]$$

(16)

or

$$N(D)=N_0\exp [-4(D/D_{\mathrm{m}})]$$

(17)

This PSD model is a well-known exponential PSD function [30].

3.3. 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]:

$$N(D)=N_{w}F(D/D_{\mathrm{m}})=N_{w}F(X)$$

(18)

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:

$$M_n={\displaystyle {\int }_0^{\infty }{D}^{n}N(D)dD}=\frac{N_0\mathsf{\Gamma }(n+1+\mu )}{{\mathsf{\Lambda }}^{n+\mu +1}}$$

(19)

D_m and N_w are expressed as follows, in terms of M_n:

$$D_{\mathrm{m}}=M_4/M_3=(4+\mu )/\mathsf{\Lambda }$$

(20)

$$N_{\mathrm{w}}=\frac{4^4}{\mathsf{\Gamma }(4)}\frac{M_3^5}{M_4^4}=N_0{D}_{\mathrm{m}}^{\mu }\frac{\mathsf{\Gamma }(4+\mu )}{\mathsf{\Gamma }(4)}\frac{4^4}{{(4+\mu )}^{4+\mu }}$$

(21)

From Equations (14), (18), (20), and (21),

$$F(X)=\frac{\mathsf{\Gamma }(4)}{4^4}\frac{{(4+\mu )}^{4+\mu }}{\mathsf{\Gamma }(4+\mu )}{X}^{\mu }\exp [-(4+\mu )X]$$

(22)

Therefore, the normalized gamma PSD can be represented as follows:

$$\frac{N(D)}{N_w}=\frac{\mathsf{\Gamma }(4)}{4^4}\frac{{(4+\mu )}^{4+\mu }}{\mathsf{\Gamma }(4+\mu )}{\left(\frac{D}{D_{\mathrm{m}}}\right)}^{\mu }\exp [-(4+\mu )\left(\frac{D}{D_{\mathrm{m}}}\right)]$$

(23)

Importantly, the normalized gamma PSD described in Equation (23) is represented by three independent parameters, N_w, D_m, and μ.

3.4. Calculation of PSD Parameters

Three parameters of the gamma PSD, N₀ (or N_w), Λ (D₀ or D_m), and μ can be calculated using the nth moment M_n of the PSD [41,42], as follows:

$$M_n={\displaystyle {\int }_0^{\infty }{D}^{n}N(D)dD}$$

(24)

M_n can be calculated by the following equation using PSDs measured by Parsivel²:

$$M_n={\displaystyle \sum _{i=1}^{32}{D}_i^n}{N}_i(D_i)$$

(25)

For convenience of computation, we introduce a variable η, which is defined as the combination of M₂, M₄, and M₆, as follows:

$$\eta =\frac{M_4^2}{M_2{M}_6}=\frac{{\left[\mathsf{\Gamma }(5+\mu )\right]}^2}{\mathsf{\Gamma }(3+\mu )\mathsf{\Gamma }(7+\mu )}=\frac{(3+\mu )(4+\mu )}{(5+\mu )(6+\mu )}.$$

(26)

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:

$$\mu =\frac{(7-11\eta )-{({\eta }^2+14\eta +1)}^{1/2}}{2(\eta -1)}.$$

(27)

If μ is known, the slope parameter Λ of the gamma profile is given by [42]:

$$\mathsf{\Lambda }={\left[\frac{M_2\mathsf{\Gamma }(5+\mu )}{M_4\mathsf{\Gamma }(3+\mu )}\right]}^{1/2}={\left[\frac{M_2(4+\mu )(3+\mu )}{M_4}\right]}^{1/2}$$

(28)

where $\mathsf{\Gamma }(x)=(x-1)\mathsf{\Gamma }(x-1)$. N₀ is calculated as follows:

$$N_0=\frac{M_n{\mathsf{\Lambda }}^{\mu +n+1}}{\mathsf{\Gamma }(\mu +n+1)}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}(n=2,4,\hspace{0.17em}\mathrm{or}\hspace{0.17em}6)$$

(29)

where N₀ 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₀ for each case is about 1%. D₀ can be computed as follows if μ and Λ are obtained:

$$D_0=\frac{3.82+\mu }{\mathsf{\Lambda }}$$

(30)

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.

3.5. 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:

$$N_{\mathrm{T}}={\displaystyle {\int }_0^{D_{\max }}N(D)dD}=m_0$$

(31)

$$C_{\mathrm{A}}=\frac{{\rho }_w\pi }{6}{\displaystyle {\int }_0^{D_{\max }}{D}^{3}N(D)dD}=\frac{{\rho }_w\pi }{6}{m}_3$$

(32)

$$\begin{array}{l}{R}_{\mathrm{A}}=\frac{{\rho }_w\pi }{6}{\displaystyle {\int }_0^{D_{\max }}{D}^3{v}_{\mathrm{t}}(D)N(D)dD}\\ \hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}=\frac{a{\rho }_w\pi }{6}{\displaystyle {\int }_0^{D_{\max }}{D}^{3+b}N(D)dD=\frac{{\rho }_w\pi }{6}{m}_{3+b}}\end{array}$$

(33)

where v(D) is the falling velocity of volcanic ash particles of diameter D, calculated as follows:

$$v_{\mathrm{t}}(D)=a{D}^b$$

(34)

$$Z={\displaystyle {\int }_0^{D_{\max }}{D}^{6}N(D)dD}=m_6$$

(35)

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:

$$N_{\mathrm{T}}=m_0=\frac{N_0}{{\mathsf{\Lambda }}^{\mu +1}}\gamma \left(\mu +1,\alpha \frac{D_{\max }}{D_{\mathrm{m}}}\right)$$

(36)

$$C_{\mathrm{A}}=\frac{{\rho }_w\pi }{6}{m}_3=\frac{{\rho }_w\pi }{6}\frac{N_0}{{\mathsf{\Lambda }}^{\mu +4}}\gamma \left(\mu +4,\alpha \frac{D_{\max }}{D_{\mathrm{m}}}\right)$$

(37)

$$R_{\mathrm{A}}=\frac{{\rho }_w\pi }{6}{m}_{3+b}=\frac{{\rho }_w\pi }{6}\frac{N_0}{{\mathsf{\Lambda }}^{\mu +b+4}}\gamma \left(\mu +b+4,\alpha \frac{D_{\max }}{D_{\mathrm{m}}}\right)$$

(38)

$$Z=m_6=\frac{N_0}{{\mathsf{\Lambda }}^{\mu +7}}\gamma \left(\mu +7,\alpha \frac{D_{\max }}{D_{\mathrm{m}}}\right)$$

(39)

The gamma PSD parameters obtained by the full moment method in Section 3.4 may then be substituted into these equations.

4. Results

4.1. 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², 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⁻¹, with larger particles occurring less frequently but with greater variance. This variation may be caused by Parsivel² 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:

$$v_{\mathrm{t}}=3.18D^{0.728};\hspace{0.17em}{v}_{\mathrm{t}}\hspace{0.17em}[{\mathrm{ms}}^{-1}],\hspace{0.17em}D\hspace{0.17em}[\mathrm{mm}]$$

(40)

For volcanic ash particles (D ≤ 2 mm), the fall velocity equation is:

$$v_{\mathrm{t}}=3.14D^{0.817};\hspace{0.17em}{v}_{\mathrm{t}}\hspace{0.17em}[{\mathrm{ms}}^{-1}],\hspace{0.17em}D\hspace{0.17em}[\mathrm{mm}]$$

(41)

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², 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⁻³).

4.2. PSD Parameters

4.2.1. Temporal Changes in PSD

Changes in the PSDs of the six selected eruptions over time, as eruption clouds passed over the Parsivel² device, are shown in Figure 5. The maximum R_A (11.7 mm h⁻¹) was observed during eruption 4. The number of particles observed per unit volume (N_T) ranged from 10 to 10³ m⁻³. R_A and log₁₀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₁₀N_T and log₁₀ 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² 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² 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² 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.

4.2.2. 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⁰ to 10⁴). 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.

To clarify the statistical characteristics of the PSDs of volcanic ash particles, we investigated the frequency distribution and probability density function of four parameters (log₁₀N₀, Λ, D₀, and μ) of gamma PSDs, and of the additional PSD parameters D_max and N_T observed in all eruptions. The results are shown in Figure 8. The six parameters (log₁₀N₀, D₀, Λ, μ, D_max, and N_T) showed modes (standard deviation [SD]) of 8.8 (9.6), 0.65 (0.29), 13.4 (25.0), 8.3 (13.7), 0.97 (0.679), and 182 (380), respectively. It interesting that each of three parameters (log₁₀N₀, Λ, μ) 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. Statistical values of gamma, normalized, and integrated particle size distribution (PSD) parameters. ρ_b = 1.25 × 10³ kg m⁻³ and ρ_p = 2.50 × 10³ kg m⁻³ were assumed for the estimation of R_A.

Parameter	Unit	Mode	Med	Mean	SD	Max	Skew	D10	D90
log10 N0	mm−1 m−3	8.83	9.64	15.8	9.64	47.9	0.581	4.79	28.2
L	mm−1	13.4	27.2	33.5	25.0	118	0.742	6.05	65.4
D0	mm	0.65	0.69	0.79	0.292	2.13	1.88	0.54	1.21
m	–	8.3	17.1	18.9	13.7	67.3	0.661	3.10	36.8
Dmax	mm	0.97	1.19	1.42	0.679	3.75	1.42	0.81	2.36
log10Nw	mm−1 m−3	4.60	4.42	4.15	0.704	5.15	−0.875	3.07	4.87
Dm	mm	0.853	1.19	0.853	0.350	2.36	1.83	0.559	139
NT	m−3	182	252	381	380	1788	1.75	26	998
RA	mm h−1	0.52	1.08	2.01	2.17	11.7	1.64	0.23	5.43
CA	g m−3	0.0931	0.180	0.225	0.189	0.845	1.21	0.0323	0.478
10log10Z	dBZ	17.1	18.2	19.3	7.43	36.9	0.087	9.28	28.8

.

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.

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. Summary of relationships among gamma PSD parameters. The numbers in the brackets are the correlation coefficient and the root-mean-square error, respectively.

	log10N0	Λ	μ	Dmax
log10N0	-	log10N0 = 0.383 Λ + 2.942 (0.995, 0.91)	log10N0 = 0.684 μ + 2.89 (0.973, 2.211)	log10N0 = 19.67 Dmax −1.606 (0.920, 3.784)
Λ	Λ = 2.585 log10N0 − 7.308 (0.995, 2.361)	-	Λ = 1.751 μ + 0.469 (0.960, 7.005)	Λ = 42.23 Dmax−2.046 (0.937, 8.76)
μ	μ = 1.386 log10N0 − 3.02 (0.974, 3.148)	μ = 0.527 Λ + 1.222 (0.960, 3.842)	-	μ = 23.87 Dmax−1.713 (0.853, 7.178)
Dmax	Dmax = 5.863 log10N0 − 0.595 (0.922, 0.263)	Dmax = 4.862 Λ−0.418 (0.937, 0.238)	Dmax = 3.018 μ−0.314 (0.779, 0.393)	-

. 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₀ and Λ are estimated by the N₀–μ and Λ–μ relationships, respectively.

4.2.3. Normalized PSD

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², 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₁₀N_w) used to normalize the PSD. The shape parameter μ is shown in Figure 8d. The modes (SD) of D_m and log₁₀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).

4.2.4. 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:

$$R_{\mathrm{A}}/N_{\mathrm{w}}=2.54\times {10}^{-4}{D}_{\mathrm{m}}{}^{4.97}$$

(42)

By rearranging Equation (42), we obtained:

$$R_{\mathrm{A}}=2.54\times {10}^{-4}\hspace{0.17em}{D}_{\mathrm{m}}{}^{4.97}{N}_{\mathrm{w}}$$

(43)

The power law equations for the upper and lower 95% confidence limits of Equation (43) are as follows:

$$R_{\mathrm{A}}=2.72\times {10}^{-4}\hspace{0.17em}{D}_{\mathrm{m}}{}^{5.07}{N}_{\mathrm{w}};\hspace{0.17em}\hspace{0.17em}95\%\hspace{0.17em}\mathrm{upper}\hspace{0.17em}\mathrm{limit}$$

(44)

$$R_{\mathrm{A}}=2.36\times {10}^{-4}\hspace{0.17em}{D}_{\mathrm{m}}{}^{4.87}{N}_{\mathrm{w}};\hspace{0.17em}\hspace{0.17em}95\%\hspace{0.17em}\mathrm{lower}\hspace{0.17em}\mathrm{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].

4.3. Conventional Relationship for Quantitative Ash Fall Estimation

Figure 13 shows the frequency and probability density distributions of the three integrated PSD parameters (R_A, C_A, and 10log₁₀Z). In this study, the modal values (SD) of R_A, C_A, and 10log₁₀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. R_A–Z and C_A–Z relationships obtained from regression analyses of PSD data measured by Parsivel². Note that R_A (mm h⁻¹), C_A (g m⁻³), Z (mm⁶ m⁻³). R: correlation coefficient, RMSE: root mean square error. ρ_b = 1.25 × 10³ kg m⁻³ and ρ_p = 2.50 × 10³ kg m⁻³ were assumed for the estimation of R_A.

Case	Site	Sampling Period (LST)	(min)	RA–Z Relationship (R, RMSE)	CA–Z Relationship (R, RMSE)	Range
1	NABE	08:22–08:35, 15 May 2018	14	RA = 8.62 × 10−2 Z 0.809 (0.963, 0.300)	CA = 1.51 × 10−2 Z 0.680 (0.912, 0.046)	0 ≤ Z ≤ 1 × 102
2	HART	20:03–20:11, 22 May 2018	9	RA = 5.44 × 10−2 Z 0.595 (0.910, 0.907)	CA = 0.629 × 10−2 Z 0.569 (0.950, 0.058)	0 ≤ Z ≤ 2 × 103
	HIKP	20:00–20:14, 22 May 2018	12	RA = 3.03 × 10−2 Z 0.826 (0.977, 0.397)	CA = 0.496 × 10−2 Z 0.770 (0.980, 0.043)	0 ≤ Z ≤ 1 × 103
3	KURP	08:09–08:13, 30 May 2018	4	RA = 14.6 × 10−2 Z 0.522 (0.959, 0.613)	CA = 1.72 × 10−2 Z 0.423 (0.903, 0.055)	0 ≤ Z ≤ 1 × 103
4	HART	12:03–12:20, 10 June 2018	12	RA = 9.45 × 10−2 Z 0.568 (0.972, 0.929)	CA = 4.49× 10−2 Z 0.294 (0.899, 0.073)	0 ≤ Z ≤ 5 × 103
5	SBTT	07:35–08:21, 16 June 2018	31	RA = 3.55 × 10−2 Z 0.778 (0.917, 0.339)	CA = 1.20 × 10−2 Z 0.692 (0.960, 0.036)	0 ≤ Z ≤ 4 × 102
6	AKAM	15:53–16:12, 16 July 2018	18	RA = 10.0 × 10−2 Z 0.653 (0.851, 1.123)	CA = 0.553 × 10−2 Z 0.782 (0.913, 0.100)	0 ≤ Z ≤ 7 × 102
All	-	-	88	RA = 22.8 × 10−2 Z 0.436 (0.819, 1.247)	CA = 6.56 × 10−2 Z 0.261 (0.634, 0.147)	0 ≤ Z ≤ 5 × 103

. A scatter plot of R_A and Z values obtained from Parsivel² observations for all eruptions is shown in Figure 14. The R_A–Z relationship was defined using the nonlinear least-squares method, as follows:

$$R_{\mathrm{A}}=22.8\times {10}^{-2}\hspace{0.17em}{Z}^{0.436};\hspace{0.17em}\hspace{0.17em}{R}_{\mathrm{A}}\hspace{0.17em}[{\mathrm{mm} \mathrm{h}}^{-1}],\hspace{0.17em}Z\hspace{0.17em}[{\mathrm{mm}}^6{\mathrm{m}}^{-3}]$$

(46)

The R_A–Z relationships for the upper and lower 95% confidence limits are as follows:

$$R_{\mathrm{A}}=31.0\times {10}^{-2}\hspace{0.17em}{Z}^{0.489};\hspace{0.17em}\hspace{0.17em}95\%\hspace{0.17em}\mathrm{upper}\hspace{0.17em}\mathrm{limit}$$

(47)

$$R_{\mathrm{A}}=14.6\times {10}^{-2}\hspace{0.17em}{Z}^{0.383};\hspace{0.17em}\hspace{0.17em}95\%\hspace{0.17em}\mathrm{lower}\hspace{0.17em}\mathrm{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.

5. Discussion

In radar meteorology, various mathematical formulas for ‘precipitation particles’ have been proposed by researchers. These include exponential [30], gamma [27,28], log-normal [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². 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₁₀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₀, μ, Λ) calculated by the momentum method and non-linear regression analysis. Regression lines obtained were as follows.

Log₁₀ N₀_mom = 0.934 × log₁₀ N₀_reg + 1.47,

μ_mom = 0.904 × μ_reg + 2.14,

Λ_mom = 0.959 × Λ_reg + 2.72.

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₀, μ, 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.

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², 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⁰ to 10⁴, these statistical values confirmed that the gamma PSD model is appropriate for describing volcanic ash PSD.

6. 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² 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₁₀N₀, Λ, 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₁₀N₀, Λ, 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₀ (or Λ) can be estimated from the differential reflectivity Z_DR, which is measured by polarimetric radar, and N₀ is estimated from D₀ 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.