© G. Ou et al., Published by EDP Sciences 2023

1 Introduction

Recently, as the demand for crude oil has increased, the crude oil quality has become obviously lower [1,2]. An S Zorb unit can convert high-sulfur gasoline into clean gasoline, which is of importance [3,4]. The S Zorb desulfurization reactor is the core part of an S Zorb unit, wherein gasoline is desulfurized via the gasoline adsorption desulfurization technology based on the principle of adsorption. Desulfurization is achieved via the selective adsorption of sulfur atoms in sulfur-containing compounds by adsorbents [5]. After desulfurization, the sulfur content of ultra-low sulfur gasoline products is less than 10 ppm, which greatly improves the cleanliness of automotive gasoline [6]. All the technological measures adopted in this process technology, such as energy saving, consumption reduction, pollution prevention, whole process control, waste recycling, and resource regeneration, embody the concept of recycling economy: reduction, reuse, and recycling. Currently, more than 60 S Zorb units are being used in China, with more under construction. During the operation of an S Zorb desulfurization unit, the hydrogen-oil mixed gas in the reactor enters the bubble cap from the lifting pipe and impinges the surfacing layer surface of the distribution plate with adsorbent particles. After the surface of the surfacing layer is eroded, the base material of the distribution plate will be directly eroded, which will cause unplanned shutdown of the unit and seriously affect the safe operation of the unit. To improve the safe operation period of the reactor distribution plate under severe conditions of high temperature and pressure and improve the product quality and economic benefit of enterprises, the erosion characteristics of the S Zorb reactor distribution plate need to be studied and an optimized modification scheme needs to be proposed for the existing S Zorb reactor distribution plate.

Since the state of solid particles in an S Zorb reactor is like that in a fluidized bed, a simulation method for a fluidized bed is required, and consequently, a particle loading method that is suitable for this study is presented. The main method for studying the scouring characteristics of a distribution plate in an S Zorb reactor is the calculation of the gas–solid two-phase flow erosion in a fluidized bed [7,8]. Scholars worldwide have conducted considerable numerical simulation and experimental research on the gas–solid two-phase flow in a fluidized bed, providing the flow characteristics and motion rules of gas–solid two-phase flows and presenting the theoretical basis for the design, optimization, and operation of a multi-phase flow reactor [9,10]. However, specific research on the inside of an S Zorb reactor under high temperatures and pressures still needs to be conducted. Alobaid et al. [11] used a DEM model to study the flow characteristics of particles in a quasi-two-dimensional spouted and circulating fluidized bed. They found that the extended CFD-DEM model can accurately predict the particle motion and pressure gradient in the fluidized bed, but the calculation cost is high and an effective program design is required. Saidi et al. [12] studied the gas–solid characteristics in a fluidized bed and found that the wall thickness affects the flow distribution. For a particle phase, the simulation results using a traditional DPM model significantly differ from the actual erosion results. Therefore, this study uses the EDEM software to perform the CFD-DEM coupling simulations.

Gas–solid erosion prediction is an important method for studying the erosion characteristics of an S Zorb reactor distribution plate; thus, the model accuracy directly affects the prediction results. Ahmed et al. [13] fully considered the change of solid particle properties using discrete phase and realizable k-ε models to simulate the erosion of a gas–solid two-phase flow, and their simulation results are consistent with the published data. Zhou et al. [14] used the numerical simulation method to study the influence of particle shape and turbulence intensity of a gas–solid two-phase flow on elbow erosion. By comparing the erosion rates calculated using different classical models, they calibrated a calculation model based on CFD-DEM simulations and verified the effectiveness of the selected model via experiments. Varga et al. [15] employed CFD-DEM numerical simulation and experiments to predict the erosion of a feed pipe in a gas–solid two-phase flow and found that material erosion failure generally occurs at large particle-wall collision angles. Additionally, they indicated that the Hashish empirical erosion model does not consider the critical impact energy and failure phenomenon of the material. Madhusuden et al. [16] stated that the particle-tracking algorithm must properly consider the turbulent diffusion. The use of the random walk model will significantly influence the turbulent diffusion, improving the accuracy and consistency of the CFD erosion prediction. Darihaki et al. [17] considered the effects of particle mass and secondary impact and found that a greater secondary impact is required in the fluid to achieve true erosion predictions. Fan et al. [18] studied the relation between particle aggregation force and radial and axial mixing in a two-dimensional fluidized bed using the DEM model and introduced non-uniform air distribution. They greatly improved the mixing process by changing the motion of air bubbles and helped reduce the particle aggregation. Arkawazi et al. [19] developed a model by combining a moving grid with CFD-DEM coupling to solve the fluid-structure interaction of a free body in a particle-fluid system. Previous studies have modified the gas–solid erosion model according to different research objectives, significantly improving the simulation accuracy.

In this study, starting from the erosion risk of an S Zorb reactor distribution plate, the erosion of the reactor distribution plate is analyzed through numerical simulations and the failure cases are verified. Moreover, the gas–solid two-phase erosion state in an S Zorb reactor distribution plate is studied, and the erosion failure mechanism of the S Zorb reactor is analyzed and summarized. The characterization method is established. Additionally, the high-risk area is predicted, and different parameter optimization schemes are presented through the orthogonal test to improve the operation life of an S Zorb reactor distribution plate.

2 Erosion mechanism of S Zorb reactor distribution plate

2.1 S Zorb desulfurization reaction process

The service process of an S Zorb reactor distribution plate is the S Zorb desulfurization reaction process. The process flow is shown in Figure 1. Sulfur-containing gasoline from the catalytic unit or tank area is filtered into the raw material buffer tank (D-101), mixed with circulating hydrogen after pressurization by the reaction feed pump (P-101), and then enters the adsorption feed heat exchanger (E-101A∼F) via two ways to exchange heat with the top part of the reactor (R-101). After the heat exchange, the mixed hydrogen material is heated in the feed heating furnace (F-101). After reaching the predetermined temperature, the material enters the bottom of the desulfurization reactor and undergoes the adsorption desulfurization reaction. The nominal scale is 900000 tons per year, and the annual operation time is 8400 h.

Fig. 1 S Zorb desulfurization reaction process diagram.

2.2 Modeling and operating conditions of the distribution plate

Figure 2a displays the schematic of the structure of the distribution plate. The distribution plate mainly comprises the lifting pipe, bubble cap, surfacing layer, and the base of the distribution plate. The bubble caps in the distribution plate are arranged in a formation, and the center distance L of the adjacent bubble cap is 152 mm. Figure 2b shows the bubble cap and the magnified diagram of the bubble cap unit. The inner diameter of the bubble cap d₁ is 92 mm, and the inner diameter of the lifting pipe d₂ is 48 mm. The distance from the top of the bubble cap to the surfacing layer h₁ is 133 mm, the shortest distance from the bubble cap to the surfacing layer h₂ is 13 mm, and the distance from the lifting pipe to the bubble cap h₃ is 60 mm. In the bubble cap expansion, a U-shaped outlet and an n-shaped outlet form a tooth; L_U is the width of U-shaped outlet and L_n is the width of n-shaped outlet. Hydrogen-oil mixed gas enters the bubble cap through the lifting pipe and changes the flow path of gas in the bubble cap. Subsequently, the gas can evenly disperse from under the bubble cap, make full contact with sorbent particles, and achieve the best adsorption effect. Since each bubble cap has the same structure and same distance between adjacent bubbles, four adjacent 1/4 bubbles are selected as one bubble cap unit instead of the distribution plate as a whole to perform the simulation calculation. This not only considers the interaction of flow fields between adjacent bubbles and caps but also improves the calculation efficiency.

Table 1 presents the physical properties and operating conditions of the reaction media in the reactor distribution plate.

Fig. 2 (a) Structural schematic of the distribution disc; (b) bubble cap unit and bubble cap expansion diagram.

2.3 Formation mechanism and erosion characteristics

The material of the surfacing layer on the distribution plate is E347 and the base material is 14Cr1MoR (H), both of which are plastic materials. The main components of the sorbent are shown in Table 2. Driven by the drag force of hydrogen-oil mixed gas, the sorbent particles impact the surfacing layer at different speeds and angles and produce micro-cutting erosion. Some particles that impact the surfacing layer rebound and then impact it again, causing secondary erosion, which seriously affects the base material. The failure mechanism diagram is shown in Figure 3.

To express the erosion law of the distribution plate better, the following characterization parameters are introduced herein:

(1) The maximum erosion rate: The relative erosion rate of the material is used to quantitatively characterize the erosion performance of the material. It represents the erosion resistance of the distribution plate under the impact of sorbent particles. The complex erosion mechanism involves many factors. Since the erosion mechanism of the distribution plate in the S Zorb reactor is unclear, therefore, we compare Finnie [20], E/CRC [21] and Oka [22] erosion models.

The Finnie erosion model is as follows:

(1)

(2)

Here, E_r¹ is the wall erosion rate (m³) of the Finnie model, K is a dimensionless constant depending on the particle and wall properties, f(α) is a dimensionless function of the impact angle, m is the mass of particles, u_p is the impact velocity, and p_w is the wall flow stress.

The E/CRC erosion model is defined as follows:

(3)

(4)

where E_r² is the wall erosion rate (kg/kg) of the E/CRC erosion model, and F_s is the coefficient determined by the particle shape, which is 0.2, 0.53 and 1 for circular, semi-circular and serrated particles respectively; K, a, b, C₁, C₂, C₃, C₄ and n are constants.

The Oka erosion model is defined as follows:

(5)

(6)

where E_r³ is the wall erosion rate of Oka model (m³ kg⁻¹); H_v is the wall Vickers hardness; D_ref is the reference diameter; U_ref is the reference speed; A, b, k₁, k₂, k₃, n₁ and n₂ are constants that depend on the properties of particles and walls.

In order to test the compatibility of the classical erosion model with the S Zorb distribution plate erosion prediction calculation, the erosion calculation is calibrated by comparing the erosion of the distribution plate with Finnie model, Oka model and E/CRC model under the same working conditions.

In order to unify, erosion rate is expressed in dimensionless form, which can eliminate the influence of uncertain constant in erosion model. The dimensionless erosion rate is defined as the ratio of the erosion rate of each area to the maximum erosion rate of the model at this time. Figure 4 shows the comparison of dimensionless erosion rates of three erosion models. It can be seen from the figure that the erosion rates calculated by the three models are in good agreement, and the change trend tends to be consistent. Under the three models, the standard deviation (SD) of dimensionless erosion rate of Finnie model, Oka model and E/CRC model are 0.2472, 0.2483 and 0.2580 respectively. The results show that the deviation between Finnie model and the mean value is the smallest. Therefore, Finnie erosion model is used to study the erosion of the distribution plate in the reactor.

(2) Impact rebound model: The impact-rebound recovery coefficient function proposed by Forder et al. [23] is used to modify the momentum changes of different particles before and after collision with the wall. It describes the impact angle with the particle θ. The rebound recovery coefficient of the particles varies in both the normal and tangential directions. It can be used to modify the impact-rebound coefficients of the particles and walls for steady materials and to solve the gas–solid two-phase flow erosion accurately. The formulae and relevant correction coefficients are as follows:

(7)

(8)

Proportion of high erosion area: The proportion of high erosion area in the reactor distribution plate is determined by the structure and load of the distribution plate bubble cap. The proportion of high erosion area in the reactor distribution plate will change with the flow rate at the outlet of the bubble cap. Additionally, it can be calculated as follows:

(9)

In equation (9), R is the proportion of high erosion area; S_h is the main erosion area; and S is the total area under one bubble cap unit.

For an S Zorb reactor unit, the movement of sorbent particles in the distribution plate involves the interaction of different flow fields, as shown in Figure 5. The particles are under upward pressure inside the reactor, and the hydrogen-oil mixed gas in the pipe rises through the bubble cap and drives the sorbent upward at high speeds. Eddy currents also exist outside the bubble cap.

Fig. 3 Failure mechanism diagram of distribution plate.

Fig. 4 Comparison of dimensionless erosion rates of three erosion models.

Fig. 5 Flow field diagram of reactor.

3 Numerical prediction method and result analysis

3.1 Flow control equation

Due to the presence of particle phases, the motion of the fluid phase is governed by the volume averaged Navier-Stokes equation of the compressible fluid, which can be written as

Continuity equation:

(10)

Momentum equation:

(11)

Here, α_f is the volume fraction of the fluid, ρ_f is the density of the fluid, u_f is the velocity of the fluid, and τ_f is the stress tensor of the fluid phase. R_f,p represents the momentum exchange with the particle phase, which is calculated for each unit.

3.2 Selection of turbulence model

Realizable k-ε model:

(12)

(13)

Here,

(14)

(15)

(16)

G_k represents the amount of turbulent energy k generated, G_b represents the turbulent energy caused by the lift, Y_M represents the effect of fluctuating turbulence on the total dissipation rate, and C₂ and C_1ε are constants. Additionally, σ_κ and σ_ε are turbulent Prandtl numbers of k and ε, respectively; S_k and S_ε are user-defined source terms, and η is the expansion parameter.

3.3 Particle motion model selection

In Cartesian coordinates, the force balance equation of particles is

(17)

Here, F_D (u − u_p) is the unit mass drag force of the particles:

(18)

(19)

where u is the fluid phase velocity, u_p is the particle velocity, μ is the fluid dynamic viscosity, ρ is the fluid density, ρ_p is the particle density, and R_e is the particle Reynolds number.

3.4 Grid division and boundary conditions

Based on the structure of the distribution plate, a geometrical model of a combination unit in the reactor distribution plate is established using SolidWorks. Each combination unit comprises four 1/4 bubble caps. The fluid domain extracted from the established geometric model by ANSYS is shown in Figure 6a. The non-structural grid partition method is used to complete the grid partition of the fluid domain, and the structure is shown in Figure 6b. To improve the calculation accuracy and speed, the number of grids is determined to be about 550,000 by the verification of grid independence.

The EDEM-Fluent coupling simulation method is used to calculate the gas–solid two-phase erosion of the reactor distribution plate during plant operation. Fluent software is used to solve the fluid model and realizable k-ε model is used for the turbulence model. The coupled algorithm based on pressure correction and second-order upwind discrete formula is used. The speed inlet and pressure outlet are varied depending on the load conditions. The inlet speeds used herein are 12, 14, 16, 17, 18, and 19 m/s, and the density of the continuous phase medium is 330 kg/m². For sorbent particles, EDEM software is used to solve the problem. The average particle size is 0.4 mm, and the mass flow is 0.6 kg/s. The particle accumulation height is 0.4 m, and the initial particle accumulation time is 15 s.

Fig. 6 (a) Distribution plate fluid domain geometric model; (b) meshing.

3.5 Prediction of erosion characteristics

Figure 7 displays a cloud diagram of the flow rate distribution under different load conditions. Under 90% load, the maximum velocity at the U- and n-shaped outlets is 5.32 and 4.99 m/s, respectively. Under 95% load, the maximum U- and n-shaped outlet speeds increase to 6.30 and 5.92 m/s, respectively. Under 100% load, the maximum U- and n-shaped outlet speeds increase to 7.20 and 6.75 m/s, respectively. Under 105% load, the maximum U- and n-shaped outlet speeds increase to 7.67 and 7.25 m/s, respectively. Under 110% load, the maximum U- and n-shaped outlet speeds increase to 8.15 and 7.72 m/s, respectively. Under 115% load, the maximum U- and n-shaped outlet speeds increase to 8.61 and 8.18 m/s, respectively. The figure shows that under the influence of the structure of bubble cap, a high flow rate area will form at the outlet of the bubble cap, which will result in a high impact velocity of the fluid on the surfacing layer. Compared to the n-shaped outlet, in the U-shaped outlet, the flow rate is higher and the impact on the lower surfacing layer is more severe.

Figure 8 displays the erosion location distribution of the surfacing layer obtained via numerical simulation. The erosion is mainly located in the annular area between 35 and 69 mm below the bubble cap and the corresponding area at the n-shaped outlet of the bubble cap. An obvious erosion area also forms between the adjacent bubble caps. No obvious erosion occurs in the middle diamond area.

Figure 9 presents a cloud diagram of the erosion rate distribution under different loads. At 90%, 95%, 100%, 105%, 110%, and 115% load, the maximum erosion rate of the surfacing layer is 5.947 × 10⁻⁸, 7.944 × 10⁻⁸, 1.066 × 10⁻⁷, 1.186 × 10⁻⁷, 1.253 × 10⁻⁷, and 1.522 × 10⁻⁷kg/(m²· s), respectively. At 90% and 95% loads, erosion mainly occurs in the annular area between 35 and 70 mm below the bubble cap. No obvious erosion occurs between adjacent bubble caps. After 100% load, more obvious erosion occurs between adjacent bubble caps, and the diamond area without obvious erosion significantly decreases in the middle.

Figure 10a shows the variation of the n- and U-shaped outlet flow rates with the load. Figure 10b shows the variation of the flow rate at the n-shaped outlet and the proportion of high erosion area with the load. Figure 10c shows the variation of the proportion of high erosion area of the surfacing layer and the highest erosion rate of the surfacing layer with load conditions. Figure 10d displays the variation of the flow rate at the n-shaped outlet of the bubble cap and the maximum erosion rate of the surfacing layer with load conditions. The figure shows that the flow rate at the n-shaped outlet of the bubble cap is higher than that at the U-shaped outlet under any load, and the flow rate at both locations is the same as that at the load. The change of the maximum erosion rate and proportion of high erosion area to the surfacing layer with load is consistent with the change rule of the flow rate at the bubble cap outlet. Clearly, the local high erosion of the surfacing layer is due to the high flow rate at the n-shaped outlet of the bubble cap, which drives the sorbent particles to impact the surfacing layer at high speeds and yields a high erosion area.

Fig. 7 Cloud diagram of the flow rate distribution under different loads.

Fig. 8 Distribution of the erosion location of distribution plate.

Fig. 9 Cloud diagram of the erosion rate distribution under different loads.

Fig. 10 Result data processing diagram.

3.6 Failure verification

Comparison of the actual failure of the distribution plate in Figure 11 with the simulated failure position shows that the actual failure is in the 35–70 mm annular area with the radius of the bubble cap center, while the simulated result is in the 35–69 mm annular area. The error is within 6%. Therefore, the simulation results are roughly in line with the actual situation and are reliable.

Fig. 11 Failure verification diagram.

4 Erosion resistance parameter optimization based on orthogonal test

4.1 Orthogonal test table design

Based on the erosion characteristics of the surfacing layer in the gas–solid two-phase flow of the reactor distribution plate and based on the reduction of the erosion of the surfacing layer, an orthogonal test table with three factors and three levels is proposed. The orthogonal test method is adopted to predict the erosion, and the optimum parameters are selected. The boundary conditions are set as above.

The factors influencing the erosion rate of the reactor distribution plate are hydrogen-oil ratio, number of bubble cap teeth, and the width-narrow ratio at the outlet of the bubble cap (l_n/l_U). The L₉ (3³) orthogonal table is used for series simulation in this test. Hydrogen-oil ratio is selected according to actual conditions. Number of bubble cap teeth and outlet width-narrow ratio are selected according to the principle of symmetry between the bubble cap and outlet of the bubble cap. Table 3 presents the factors and levels of the orthogonal test. In this simulation, the maximum erosion rate and the proportion of high erosion area are selected as the evaluation indices of the orthogonal test. Factors A, B, and C represent the hydrogen-oil ratio, number of bubble cap teeth, and l_n/l_U, respectively.

4.2 Analysis of orthogonal test results

The evaluation index selected in this orthogonal test is the maximum erosion rate and the proportion of high erosion area. A total of nine groups of simulations are required for three-factor and three-level orthogonal tests. The experimental results are shown in Table 4. Figure 12 more intuitively displays the results of the nine orthogonal tests. Figure 12 and Table 4 show that the maximum erosion rate of the reactor distribution plate and the proportion of high erosion area are affected by the process parameters. Figure 13 displays the cloud diagram of the erosion rate distribution of each group of cases. Among them, when the hydrogen-oil ratio is 0.3, number of bubble cap teeth is 12, and l_n/l_U is 2, the maximum erosion rate and the proportion of high erosion area are the smallest. When the hydrogen-oil ratio is 0.25, number of bubble cap teeth is 12, and l_n/l_U is 0.5, the maximum erosion rate and the high erosion area are the largest.

The order of influence of each factor on the experimental results can be obtained via calculations.

(20)

(21)

(22)

(23)

Here, K_xm represents the average index value of factor _X at level m; each index is marked as y_mn. R_ox and R_1x represent the maximum and minimum values of K_xm at all levels of the same factor, respectively, and R_x is the extreme difference corresponding to each factor.

The calculation results are shown in Tables 5 and 6. The magnitude of the range R represents the importance of the corresponding factors. R₁ denotes the influence of this factor on the erosion rate, and R₂ denotes the influence of this factor on the proportion of high erosion area. The R₁ values corresponding to the hydrogen-oil ratio, the number of bubble cap teeth, and the width-to-width ratio of bubble cap outlet are 2.19, 1.29, and 3.36, respectively; the corresponding R₂ values are 2.8, 3.5, and 3.4 respectively. Thus, the number of bubble cap teeth (factor C) is the main influencing factor, followed by the hydrogen-oil ratio (factor A) and the width-narrow ratio of the bubble cap outlet (factor B), for the maximum erosion rate of the reactor distribution plate. For the proportion of high erosion area in the reactor distribution plate, the width-narrow ratio at the outlet of the bubble cap (factor B) is the main influencing factor, followed by the number of teeth of the bubble cap (factor C) and the hydrogen-oil ratio (factor A). The optimum level combination that is intuitively obtained from the nine groups of simulation is A₃B₃C₃, which is also the optimum level combination through range analysis.

Through the orthogonal test, the following results are obtained:

• The factors influencing the maximum erosion rate of the reactor distribution plate in ascending order are the width- narrow ratio of the bubble cap outlet, the hydrogen-oil ratio, and the number of bubble cap teeth.

• The factors influencing the proportion of high erosion area in the reactor distribution plate in ascending order are the number of bubble cap teeth, the width-narrow ratio of the bubble cap outlet, and the hydrogen-oil ratio.

• The optimum operating conditions of the reactor distribution plate are hydrogen-oil ratio of 0.3, bubble cap teeth number of 12, and bubble cap outlet width-narrow ratio of 2.

Fig. 12 Maximum erosion rate and proportion of high erosion area simulated by the nine groups.

Fig. 13 Cloud diagram of the erosion rate distribution of the nine groups.

5 Conclusions

The erosion of reactor distribution plate is caused by micro cutting erosion caused by adsorbent particles impacting the surfacing layer surface under the bubble cap at high speed driven by high speed hydrogen oil mixed gas. Under the outlet of the bubble cap, a high flow rate is present in the annular area with the radius of the center of the bubble cap between 35 and 69 mm. Particularly, the flow rate was higher under the n-shaped outlet than the U-shaped outlet.

A model of the S Zorb desulfurization reactor was built, and the erosion characteristics were characterized by the maximum erosion rate and the proportion of high erosion area. The relation between the highest erosion rate and the proportion of high erosion area with the erosion characteristics under different conditions were obtained through numerical simulations. The erosion rate in the corresponding erosion area below the n-shaped outlet was higher than that below the U-shaped outlet. The change rule of the highest erosion rate and the proportion of high erosion area to the surfacing layer with load was consistent with the trend in the flow rate at the bubble cap outlet.

Considering the relation between the maximum erosion rate and the proportion of high erosion area, the parameters were optimized. By constructing the orthogonal test parameter table, the optimum operating parameters were obtained as follows: hydrogen-oil ratio of 0.3, bubble cap teeth number of 12, and bubble cap outlet width-narrow ratio of 2. Compared to the existing parameters, the maximum erosion rate and high erosion area were reduced by 29.9% and 17.3%, respectively.

Funding information

This work is supported by the National Natural Science Foundation of China (Nos. U1909216, 51876194, and 52176048).

References

All Tables

All Figures

	Fig. 1 S Zorb desulfurization reaction process diagram.
In the text

	Fig. 2 (a) Structural schematic of the distribution disc; (b) bubble cap unit and bubble cap expansion diagram.
In the text

	Fig. 3 Failure mechanism diagram of distribution plate.
In the text

	Fig. 4 Comparison of dimensionless erosion rates of three erosion models.
In the text

	Fig. 5 Flow field diagram of reactor.
In the text

	Fig. 6 (a) Distribution plate fluid domain geometric model; (b) meshing.
In the text

	Fig. 7 Cloud diagram of the flow rate distribution under different loads.
In the text

	Fig. 8 Distribution of the erosion location of distribution plate.
In the text

	Fig. 9 Cloud diagram of the erosion rate distribution under different loads.
In the text

	Fig. 10 Result data processing diagram.
In the text

	Fig. 11 Failure verification diagram.
In the text

	Fig. 12 Maximum erosion rate and proportion of high erosion area simulated by the nine groups.
In the text

	Fig. 13 Cloud diagram of the erosion rate distribution of the nine groups.
In the text