## Abstract

Typhoon Bolaven caused significant damage with severe rainfall all over South Korea, including Cheju Island, which received more than 250 mm in 2 days in August 2012. It was regarded as the most powerful storm to strike the Korean Peninsula in nearly a decade. The rainfall-rate datasets were obtained from S-band radar operated by the Korea Meteorological Administration to be analyzed and compared with the mesoscale Cloud Resolving Storm Simulator (CReSS) model simulation. Multifractal analysis was conducted to understand the structure of the rainfall rate with height in the typhoon system. The radar rainfall data presented with strong intermittency across scales at lower altitudes (1 and 2 km) and a more homogeneous rainfall field at high altitude (5 km) with two parameters (fractal codimension and multifractality index). The statistical scaling moment function and maximal singularities show clear significant differences between radar and the CReSS model.

## 1. Introduction

Multifractals represent a framework for the analysis and simulation of geophysical fields, such as rainfall, over a wide range of spatiotemporal scales. They were found and introduced in the 1980s following discussions on the scale-invariance properties of the geophysical field (Schertzer and Lovejoy 2011). Multifractals are based on the assumption of the rainfall generated by a multiplicative cascade process, and they distribute intensity structures from large to small scales. For many years, multifractal tools have been commonly used for studying rainfall fields as well as geophysical fields over a wide range of scales (Schertzer and Lovejoy 1987; Gupta and Waymire 1993; Harris et al. 1996; Marsan et al. 1996; De Lima and Grasman 1999; Deidda 2000; Biaou 2004; Macor et al. 2007; De Montera et al. 2009; Tchiguirinskaia et al. 2011; Gires et al. 2013; Hoang et al. 2014; etc.). The concept of multifractal analysis is scale invariant. It uses the method of dividing structures, and the random multiplicative increments in probability distribution are the same at each step of the cascade process, which makes it more fitted to rainfall fields than an ad hoc statistical tool since it is physically based on the agreement of the scale-invariance properties of the Navier–Stokes equation. It can govern the atmospheric behavior and is assumed to be transmitted to the unknown equations for rainfall processes. Chigirinskaya et al. (1994) pointed out that inhomogeneity plays a significant role by increasing the stability of the structures such as typhoons and indicated that the most straightforward framework for considering the extreme nonlinear variability over a wide range of scales is multifractals, and it is a result of an elementary scale-invariant process when the generator of the field reproduces itself from scale to scale.

Furthermore, multifractals are understood as the tool that provides the natural framework for scale-invariant nonlinear dynamics (Tessier et al. 1993) and has become a somewhat standard tool to analyze and simulate meteorological and hydrological data, especially radar data, which have the rare advantage of providing space–time (3D+1) fields. A unified multifractal model of atmospheric dynamics was developed by Schertzer and Lovejoy (1983, 1985) and Lovejoy et al. (1993).

Despite the capacity of dealing with extreme multiscale phenomena with high-quality data, there have been few multifractal studies of typhoons since Chigirinskaya et al. (1994) and Lazarev et al. (1994). These studies relied on time series data obtained from 1D aircraft or balloon trajectories without any radar data. Therefore, to see the detailed structure of rainfall rate with height in the typhoon, the multifractal framework was applied in this study by using radar.

Additionally, with the selected case of Typhoon Bolaven, most of the studies were focused on the Fujiwhara effect between Typhoons Bolaven and Tembin. Although the typhoon itself caused much damage and the scale was large, very few case studies were conducted with different observation equipment such as rain gauges, satellite, radar, etc. Also, the different model simulations were performed to investigate the structure and the characteristics of the typhoon, but the studies were focused on the prediction of the track of the typhoon directions or the intensities. Figure 1 shows the track of the typhoon and the accumulated rainfall amount.

Looking at previous studies about Typhoon Bolaven, Origuchi et al. (2013) explained that updrafts of the middle eyewall were stronger than those of other eyewalls, and the middle eyewall had abundant water substances below 6 km above ground level (AGL) in cloud-resolving ensemble model experiments. The regions of the strong tangential velocity existed at the outer edges of the eyewalls. Though the outer eyewall was not compared with the central eyewall, it had a similar structure. Downdrafts existed between the eyewalls. Although the updrafts and the liquid water substances in the lower layer appeared at the 10-km radius, the formation of the innermost eyewall was insufficient by the simulation of the Japan Meteorological Agency (JMA) nonhydrostatic model (JMANHM). Ryu and Lee (2012) showed the reliability of the radar rainfall-rate dataset obtained by the Korea Meteorological Administration (KMA) during this event by comparing with the satellite *Chullian* dataset. Also, Liu et al. (2015) presented a framework of Fujiwhara effect between Typhoons Bolaven and Tembin by using remote sensing imagery and image processing techniques, which explains how Bolaven was intensified and changed its direction by the impact of combined interactions. Sun et al. (2015) mentioned the torrential rain during Bolaven was caused by increased cyclonic vorticity. The warm–moist air mass from the southeast of the rainfall area caused the strong ascending motion upward, while the conditional symmetric instability was an important instability mechanism for the torrential rain enhancement.

There are several references that describe the development mechanisms and the structure of the convection bands that cause heavy rainfall over the Korean Peninsula, on the basis of case studies and numerical modeling (Sun and Lee 2002; Shin and Lee 2005; Cho and Lee 2006; Hong and Lee 2009; Yu and Lee 2010; Choi et al. 2011). By performing the high-resolution numerical experiment for Typhoon Bolaven, the maintenance mechanism of a long-lived concentric eyewall is explained by the lack of dissipation of the inner eyewall and the constancy of the large radius of the outer eyewall (Tsujino et al. 2017). Generally, these modeling studies demonstrated some significant discrepancies between the simulation and observation in both the location and amount of heavy rainfall. However, the Cloud Resolving Storm Simulator (CReSS) captures relatively well the location and intensity of the maximum rainfall. Lee et al. (2010) have applied CReSS in experiments of multiscale cloud and precipitation systems to simulate topographically induced localized intense rainfall over Cheju Island, South Korea. Therefore, the CReSS model was chosen for our selected case study of Typhoon Bolaven rainfall on Cheju Island.

This study was motivated by the fact that not only are there not enough studies focused on Bolaven alone, but also no research was carried out with the approach of explaining the nonlinearity of the scaling exponents of the rainfall structure by performing a multifractal analysis of Typhoon Bolaven. To jointly understand the dynamics and rainfall by multifractal space–time analysis, the Doppler S-band radars and mesoscale simulation CReSS model are compared during Typhoon Bolaven. In the following sections, a detailed description of the three-dimensional structure of Bolaven during landfall on Cheju Island measured by different observation instruments and the procedure of preparing the dataset of rainfall rate are presented. We introduce the data in section 2, and the analysis methodology is described in section 3. The results of the multifractal analysis performed on radar data and CReSS model simulation are discussed in section 4 including a comparison between the radar data and the model simulation. Last, the summary and conclusions of the results are in section 5.

## 2. Data

Typhoon Bolaven was formed as a tropical depression on 19 August 2012 and steadily intensified to a typhoon by 21 August. On 24 August, the system attained its peak intensity with a wind speed of 51.4 m s^{−1} (115 mi h^{−1}) and a barometric pressure of 910 hPa. The storm passed over Okinawa, Japan, on 26 August and it began accelerating toward the north, phasing into the weakening stage. As Typhoon Bolaven approached the Korean Peninsula, it continued to weaken and eventually made landfall on North Korea late on 28 August before transitioning into an extratropical cyclone. Even though Typhoon Bolaven was in its weakening stage, it was an extremely powerful typhoon that caused a lot of damage with severe rainfall all over the Korean Peninsula including Cheju Island. The rainfall amount was more than 250 mm in 2 days. It was regarded as the most powerful storm to strike the Korean Peninsula in nearly a decade.

Information on the typhoon location and its lowest pressure was recorded by the National Typhoon Center on Cheju Island. Figure 1a shows the track of the typhoon, and its center is indicated by red dots. Also, surface station data were collected by 22 sites on Cheju Island. Figure 1b shows the daily accumulated rainfall on 27 August. These data were collected by an automatic weather system (AWS) operated by the KMA. The maximum accumulated rainfall amount was detected in the middle of Cheju Island where the mountain Halla (1950 m above sea level) is located.

The radar dataset was obtained from 0600 to 2350 LST 27 August 2012 when the typhoon approached Cheju Island. The selected radar site was Gosan on Cheju Island operated by the KMA and covering a radius of 360 km. It recorded volume scans of reflectivity and Doppler velocity every 10 min. The data were interpolated in a Cartesian coordinate system with 1-km horizontal and 0.25-km vertical grid intervals (CAPPI), and the three Cartesian components of reflectivity were calculated. Figure 2 shows the wind field superimposed on the radar reflectivity (dB*Z*) at an altitude of 5 km on the left of the figure. The location of the Gosan radar is shown as the red star with a white circle in the middle of the figure. The vertical cross section of the red line on the left of the figure is shown on the right of the figure. Updrafts and downdrafts are clearly shown in the vertical cross section. However, there were missing radar data that were obtained as zero values at the lowest altitudes due to the minimum radar beam elevation angle of 0.5°. The multifractal analysis was performed on the same area of 256 km × 256 km at various altitudes. Further analysis was performed to cover the limitations of the missing data with a smaller domain that was 64 km × 64 km. Then, the rainfall rate was computed at separate altitudes using the reflectivity of the Gosan radar, which is one of the three radars, by using the *Z*–*R* relationship [*Z* = *aR*^{b}, with radar reflectivity factor *Z* (mm^{6} m^{−3}) and rain rate *R* (mm h^{−1})] with the values of *a* = 250 and *b* = 1.2, which are the parameters usually used for tropical convective systems.

The numerical simulation model CReSS is a three-dimensional nonhydrostatic model developed by the Hydrospheric Atmospheric Research Center (HyARC) of Nagoya University, Japan (Tsuboki and Sakakibara 2002). The initial and lateral boundary conditions were provided by the Japan Meteorological Agency global spectral model (JMA-GSM), which is a gridpoint-values reanalysis database. It has one of the highest horizontal resolutions of 0.1875° (approximately 20 km) with a time interval of 6 h. The JMA-GSM is often used for deep convection simulations, such as typhoon cases since it produces data up to 10 hPa, which contain information on the lower stratosphere and can detect the effect of significant gravity wave propagation in the upper-levels of the atmosphere.

Also, to set the surface fluxes of momentum and energy and surface radiation processes, the sea surface temperature (SST) was calculated by using a one-dimensional, vertical heat diffusion equation (Kondo 1976; Louis et al. 1981; Segami et al. 1989) that included an underground layer for ground temperature prediction. The SST at the initial time was calculated from the dataset of the North-East Asian Regional Global Ocean Observing System (NEAR-GOOS) real-time database, which was provided by the JMA. The land-use data were obtained from the dataset of the global land-cover characteristics database, which was provided by the U.S. Geological Survey.

For the simulation, the horizontal grid resolution was 1 km × 1 km with a mesh size of 936 × 1248 and a vertical grid resolution of 0.25 km containing 70 levels ranging from a near-surface level at 50 m to the top level at 17.3 km. The vertical calculation was performed with a terrain-following coordinate as the terrain effect was also considered during the calculation, and the rainfall-rate parameter was the instantaneous rainfall rate, which was obtained as one of the output parameters. It was calculated with the equation *ρV*_{t} = *dq*_{l}/*dz*. The duration of each time step was 10 min, similar to that of the radar data. The high range of the altitudes selected was due to the inclusion of the frozen precipitation at high altitudes since the scale of the typhoon was very large. A general rainfall field analysis with CreSS model data is shown in appendix. However, in the main paper, to compare the results obtained with the results of the multifractal analysis of the radar data, the same-size domains with mesh sizes of 256 × 256 and 64 × 64 covering the radar observation sites were selected.

## 3. Method

Multifractal tools are commonly used for studying rainfall fields and generalizations of fractal geometry. If the rainfall suits the fractality (Lovejoy and Mandelbrot 1985; Olsson et al. 1993), the number *N*_{λ} of nonzero rainfall rates at resolution *λ* (*λ* = *L*/$\u2113$, where *L* is the outer scale of the phenomenon and $\u2113$ is the observation scale) can be described in a scale-invariant notion as

where *D*_{F} is the fractal dimension that characterizes how much space a geometrical set occupies. Figure 3 shows the fractal dimensions of radar and CreSS simulations that were analyzed with the thresholds of 1, 5, and 10 mm h^{−1} at an altitude of 5 km to see the sparseness in the rainfall fields over the observation area with two different sizes of the domain. As is shown in all datasets in both domain sizes, the values of the fractal dimension decrease as depicted by the slope of each graph, while the threshold values increase. Each graph shows a different fractal dimension depending on the thresholds, and when we see them in detail all the lines show linear fitting (*R*^{2} > 0.9 for all black, blue, and red lines). Also, the black and white plots show the rainfall occurrences changing with the different thresholds more visibly. The existence of rain is indicated in black.

In the multifractal framework, *D*_{F} is strongly dependent on the threshold defining the occurrence of rainfall. This shows that more than one fractal dimension is needed to characterize the rainfall field (Gires et al. 2013). Therefore, given a multifractal field *ε*_{λ} at a resolution *λ*, the probability of obtaining a singularity of order greater than or equal to *γ* is (Schertzer and Lovejoy 1987)

where *c*(*γ*) = *D* − *D*_{F}(*γ*) is the codimension function and *D* is the dimension of the embedding space.

Also, multifractal fields can be described by their statistical moments of order *q*, following a power law given by

where *K*(*q*) is the scaling moment function that characterizes the scaling variability of the process studied.

Furthermore, Parisi and Frisch (1985) indicated that the two functions *c*(*γ*) and *K*(*q*) have a one-to-one relationship by applying the Legendre transform:

These show the correspondence between the orders of moments and the singularities that can also be considered as evidence of the correspondence between high values of multifractal parameters and extreme events.

In the framework of universal multifractals (UM) (Schertzer and Lovejoy 1987; Lovejoy and Schertzer 2007; Schertzer and Lovejoy 2011), the field can be described by only two “UM parameters” (*α* and *C*_{1}) when the field is conservative:

where α is the Levy index and measures how fast the intermittency evolves when considering singularities that are slightly different from the average field singularity, *C*_{1} is the codimension of the singularity of the average field and can measure its concentration, and *α*′ = *α*/(*α* − 1).

To validate UM, the theoretical *K*(*q*) function should be compared with the one obtained from the observation, which is called the empirical *K*(*q*). However, the theoretical *K*(*q*) can simulate the empirical one only up to a certain critical value of moment order. This critical value is related to what is called a multifractal phase transition (Schertzer and Lovejoy 1992). It is estimated as *q*_{c} = min(*q*_{s}, *q*_{D}), where *q*_{s} is the maximum-order moment estimated by a finite number of samples and *q*_{D} is the critical moment order of divergence.

The value of *q*_{s} is related to the maximal observable singularity *γ*_{s} using a Legendre transform, and it can be determined using the following equation:

For example, in a one-dimensional field (*D* = 1, where *D* is the dimensional field) with only one data sample available (*N*_{sample} = 1; thus *D*_{s} = 0, where *D*_{s} is the dimensional field with the number of samples), the critical value of moment order is usually *q*_{c} = *q*_{s}. This shows a linear behavior of the empirical *K*(*q*) for *q* ≥ *q*_{s}.

The moment order *q*_{D} represents the critical value of *q* for which extreme values of the field are becoming so dominant that the average statistical moment of order *q* ≥ *q*_{D} approaches infinity:

The moment order *q*_{D} can be determined from the following equation:

Determining the value of *q*_{D} can be graphically explained as the intersection between the theoretical *K*(*q*) function and the linear regression *K*(*q*) = (*q* − 1)*D* that corresponds to the theoretical *K*(*q*) with *C*_{1} = *D* and *α* = 0. Therefore, when *q*_{c} = *q*_{D}, the empirical *K*(*q*) function starts approaching infinity for *q* ≥ *q*_{D}.

First, a spectral analysis was performed to identify the occurring physical processes, such as checking the scaling behavior of the data. A scaling field shows a power-law relationship between the power spectra *E* and the wavenumber *ω* or frequency *f*:

where −*β* is the slope of the straight line that appears in the log–log plot of Eq. (10).

According to Tessier et al. (1993), if the spectrum behavior is linear with a spectral slope −*β*, the spectral exponent is linked to the degree of nonconservation *H* of the field:

where *H* is the Hurst exponent and is 0 for a conservative field, and *K*(2) is the second-moment scaling function of the conservative part of the field. For a conservative field, the estimate of *β* is lower than the dimension *D* of the embedding space. From Nykanen (2008), if *β* > *D*, the field should be fractionally differentiated before implementing a multifractal analysis.

Once the scaling behavior and conservativeness of the rainfall field have been obtained, a multifractal analysis can be performed. Two different methods, which are widely used in geophysics, can be used to assess the UM parameters indirectly and directly: the trace moment (TM; Schertzer and Lovejoy 1987) and the double trace moment (DTM; Lavallée et al. 1993) methods, respectively.

The TM method is performed on a broad range of moments *q* to determine the scaling moment function *K*(*q*), as shown in Eq. (3). It consists of taking the *q*th power of a multifractal field *ε*_{λ} at the highest resolution *λ* to repeatedly calculate the ensemble average at different scales and to represent the resulting averages $\u2329\epsilon \lambda q\u232a$ in a log–log plot as a function of *λ*. Then, following Eq. (3), the slopes of the linear regressions give the values of *K*(*q*) for each *q*. Using Eq. (6), the UM parameters can be estimated by

where *K*′(*q*) and *K*″(*q*) are the first and second derivatives of *K*(*q*), respectively.

The DTM method was specifically developed in the framework of UM and is conducted to obtain more robust estimates of *α* and *C*_{1} (Lavallée et al. 1993). It is a more reliable way of estimating UM parameters than the TM method. This technique is based on two steps. The first is to take the *η*th power of the conservative field *ε*_{λ} at the highest resolution *λ* and to normalize it:

Then, the TM method is applied to the normalized field $\epsilon \lambda \u2061(\eta )$ to obtain the scaling moment function *K*(*q*, *η*) for each *η* value:

Therefore, for a given *q* value, *K*(*q*, *η*) is plotted against *η* in a log–log plot and the slope of the curve gives an estimation of α. *C*_{1} can also be estimated from the interception of the slope and the axis log(*η*) = 0.

## 4. Results

The multifractal analysis was conducted on the rainfall rate at three different altitudes (1, 2, and 5 km) retrieved from the radar data and from the CReSS model.

### a. Multifractal analysis on radar data

The spectral and multifractal (TM and DTM) analyses were applied to the time ensemble average over the full rainfall event by considering each time step as an independent realization). For the analysis, two different sizes of domains were selected, 256 × 256 and 64 × 64, because of the missing radar data at the lowest altitudes resulting from the minimum radar beam elevation angle of 0.5°.

First of all, Fig. 4 shows the result of the spectral (ln[*E*(*k*)] − ln(*k*)) analysis of rainfall rates at different altitudes that were retrieved from radar data of the larger domain (256 × 256). Table 1 shows all the values of the spectral exponent *β* and *R*^{2}. In this case, *β* varies from 1.61 to 2.64, while *R*^{2} is always close to 1 (it varies between 0.96 and 0.98), which shows no extreme scaling break. The graph shows that the scaling behaviors of the altitudes of 1 and 2 km between ln(*k* ≥ 1) and ln(*k* ≤ 1.8) are not as linear as that of the altitude of 5 km. The TM analyses for the three altitudes with radar are presented in Fig. 5. On the left, the log–log plots of $\u2329\epsilon \lambda q\u232a\u2248\lambda K\u2061(q)$ with the resolution *λ* (in this case, it is set from 1 to 256 increasing by powers of 2) and values of *q* that were freely chosen between 0.1 and 7.0. A scaling behavior was shown up to when *q* = 4.5 at all altitudes. In Figs. 5d–f, the estimated (or empirical) scaling moment functions *K*(*q*) (in black) are compared with the semitheoretical functions that are the curves with the UM parameters α and *C*_{1} retrieved from TM (red) and DTM (green) analyses. The empirical *K*(*q*) curves are relatively fitting with the semitheoretical ones obtained by the TM method, where *q* = 0.5–2 at every altitude, and with the semitheoretical curves obtained by the DTM method, where *q* = 0.5–3. For *q* < 0.5, the altitudes of 1 and 2 km present a more similar tendency between the two functions than the altitude of 5 km. However, for *q* > 3.0, the semitheoretical function is more divergent from the empirical one at an altitude of 5 km than those at 1 and 2 km. Moreover, *q* = 3.0 can be considered as a multifractal phase transition point, as predicted by Mandelbrot (1974). This is caused mainly by one of two reasons: spatial integration or finite sample size (Schertzer and Lovejoy 1987). A DTM analysis was performed and the result is shown in Fig. 6. A clear scaling behavior as depicted in Figs. 6a–c was retrieved, along with the result of the TM analysis. For each power *η*, with a fixed value of *q* (*q* = 1.5), the slope of the linear regression gives an estimation of the scaling moment functions *K*(*q*, *η*). Figures 6d–f with S-shape curves are conditioned by the appearance of numerical limitations at smaller moments and the critical behavior of extremes at higher statistical moments. Both are characterized by the flattening of the DTM curves. The slope of the curve gives an estimation of α and *C*_{1}. The value of α increases with altitude while the value of *C*_{1} decreases, showing how concentrated the field is and how quickly the intermittency evolved. The values of α (α = 0.858; 1.072) are close to 1 and the values of *C*_{1} (*C*_{1} = 0.349; 0.255) are higher at altitudes of 1 and 2 km than at 5 km, which means they are sparse with high singularities (Biaou 2004). On the contrary, at 5 km (α = 1.421; *C*_{1} = 0.114), α is larger and *C*_{1} is smaller than at 1 and 2 km, which means smaller singularities and a more homogeneous rainfall field.

As mentioned before, because of missing radar data at the lowest altitudes, the smaller domain was selected (64 × 64). This zero field bias results in a spurious estimation of the mean intermittency (*C*_{1}) that decreases at higher altitudes when multifractality index estimates evolve in contrast. The results of spectral analysis at each altitude are shown in Fig. 7. The value of the spectral exponent *β* varies from 2.67 to 2.97, *R*^{2} is approximately 0.97 at lower altitudes but almost 1 at an altitude of 5 km (see Table 1). Regardless of the domain sizes, the graph shows that the rainfall field is the most conservative at 5 km. Figure 8 shows the result of the TM analysis in the smaller domain. The *q* values selected were the same as those in the previous TM analysis (from 0.1 to 7.0). Each of the graphs shows the empirical (black) *K*(*q*) and semitheoretical *K*(*q*) of the TM (red) and DTM (green) analyses. As $\u2329\epsilon \lambda q\u232a\u2248\lambda K\u2061(q)$ approaches infinity as *q* ≥ *q*_{D} and a finite value as *q* < *q*_{D} in Eq. (8), the moment order $\u2329\epsilon \lambda q\u232a$ converges at the region *q* < *q*_{D}. From Fig. 8e it is observed that, unlike the results from the larger domain, the empirical *K*(*q*) as a function of *q* has a smaller value than both the theoretical *K*(*q*) for the TM and DTM analyses, which means that the value of *α* does not increase dramatically after *q* = 3, as in the larger domain. However, at 5 km, the empirical *K*(*q*) is better fitted than the semitheoretical *K*(*q*) of the DTM analysis. Each value of the multifractality index from the different domains with the TM and DTM analyses are indicated in Table 2. Figure 9 shows the result of DTM analysis at the domain size of 64 × 64. A clear scaling behavior with the power of η at the fixed *q* = 1.5 is evident, as was also seen in Fig. 6.

### b. Multifractal analysis on CReSS model

The multifractal analysis of the CReSS model data was applied to the domain of the radar coverage with sizes of 256 × 256 and 64 × 64 (Figs. 10–12 ). The procedure was the same as the radar dataset, and the values of the result are also indicated in Table 2. More results of the multifractal analysis on the CReSS model simulation data in a larger domain (512 × 512) are described separately in appendix.

### c. Comparison between radar data and CReSS

From the results of the TM and DTM analyses, the difference of the estimated α values for the domain 256 × 256 was obtained as 0.4, while it was 0.003 for the estimated *C*_{1} values, whereas for the domain 64 × 64, the difference was 0.3 for α and 0.003 for the estimated *C*_{1}. This shows that both domains present different degrees of multifractality and similar intermittencies of average intensity. Comparing the results of the estimated α values from the DTM analysis, with the domain size of 256 × 256 and an altitude of 5 km for the radar (α = 1.421), it fits considerably well with the domain of multifractal analysis done with radar observations (DR) for the CReSS model (α = 1.416). Also, regarding Eq. (6), if the α values are considerably similar for two *K*(*q*) functions, these functions would be different from each other only by a ratio of *C*_{1} values, which would be, in this case, 0.105/0.114 = 0.921. For the domain size of 64 × 64, the estimated α value at an altitude of 5 km for the radar (α = 1.356) is smaller than the DR for the CReSS model (α = 1.817). However, the ratio of *C*_{1} values is 0.092/0.089 = 1.034. This supports the fact that it fits better at an altitude of 5 km (the ratio of *C*_{1} values at 1 km = 0.301 and 2 km = 0.412 with the domain size of 256 × 256, and the ratio of *C*_{1} values at 1 km = 0.472 and 2 km = 0.532 with the domain size of 64 × 64).

To obtain a more detailed comparison of linearity between the radar data and CReSS data, Fig. 13 compares each *K*(*q*) of the radar data at each altitude and *K*(*q*) of the CReSS model data (Figs. 13a–c with the domain size of 256 × 256 and Figs. 13d–f with the domain size of 64 × 64). This comparison of (*q*) functions shows that the consistency between the two different datasets is best at 5 km, because the curvature of blue dots at 5 km aligns straighter than at 1 km and 2 km, where there is a difference not only with *C*1 (departure from the bisectrix) but also with α (presence of a curvature; departure from the linear regression fit). *C*1 describes the sparseness of the mean value of the field (mean intermittency) and α describes how much it varies as it goes away from the mean value of the field (variability of intermittency). The sets of graphs on both domains demonstrate that departures from the bisectrix and linear regression fit diminish with increasing altitude.

The *q*_{s} and *γ*_{s} were calculated to see the maximal singularities, and they are illustrated in Fig. 14. It shows the tendency of each multifractality index and the consistency of the maximal singularity. It clearly shows that there is a decrease of the *α* values from the TM analysis when comparing the two domains, 256 × 256 and 64 × 64. It is due to the zero field at low altitudes. The *γ*_{s} from the DTM analysis remains close, with the same order between the two domains. The TM and DTM are independent methods but show the same tendency that there is a decrease of *γ*_{s} at the higher altitude for the radar. Also, to see the difference in detail, the absolute values of difference of singularity *γ*_{s} were taken and respectively divided by *C*_{1} of the radar and CReSS (see Fig. 15). This demonstrates that the difference in singularities from the TM analysis with the 64 × 64 domain is more than 1.5 times larger than the mean singularity. Figure 15b shows that only the 256 × 256 domain at an altitude of 5 km has a relatively similar result when the other values are much more than 1, which means the difference of singularities has become too strong.

## 5. Summary and conclusions

Mesoscale observation data of Typhoon Bolaven were collected on Cheju Island, South Korea as the strongest typhoon in a decade made landfall. An S-band radar recorded volume scans of reflectivity and Doppler velocity every 10 min and have been compared with CReSS model simulations. Detailed observational datasets from typhoons are rarely available for conducting multifractal analysis, and for the first time, a comparison of radar volume data and model simulations was performed. Despite the different analysis limitations, there is consistency in the multifractal intermittency.

For multifractal analysis, first, spectral analysis was performed to see the scaling behavior of the rainfall field. Then, TM and DTM analyses were performed in each case. From the results, regardless of the size of the domains, it is clear that there is a relatively good agreement of multifractality between radar at 5 km and CReSS in both domains (64 × 64 and 256 × 256), but not at lower altitudes. This may be because many case studies of tropical cyclones have been performed at high altitude, for example, 5 km, and used to tune the models, despite the possible presence of the bright band at these altitudes. To investigate this issue, the maximal singularities were compared between the radar and CReSS data. A comparison of the difference between singularities shows that maximal singularity tends to decrease at higher measurement altitude. These significant differences in singularities between the radar and CReSS show that the model simulates a much smoother field compared to the radar measurements at low altitudes.

This study is valuable as it was the first time a multifractal approach was applied to a large dataset for a typhoon case study. The multifractal parameters capture the vertical structure of the rainfall field in a typhoon in a simple way. Comparison of these parameters shows that the rainfall field obtained from a numerical model does not capture the detailed rainfall structure at different altitudes. With model rainfall estimates at all levels, this technique could be applied to further compare model and observation rainfall structure to help define the range of applicability of model rainfall.

## Acknowledgments

We acknowledge financial support from the Chair “Hydrology for Resilient Cities” of Ecole des Ponts, endowed by Veolia. This work is also supported by PHC STAR 2018 program in the framework of project 41606PG. This research was supported by the International Cooperation Program of the National Research Foundation of Korea (NRF) funded by the Ministry of Science and ICT (Grant 2018K1A3A1A21044390).

### APPENDIX

#### Multifractal Analysis of CReSS Model Simulation

For the multifractal analysis, the different domains were selected, as shown in Fig. A1. Similar to the radar data, to first identify the scaling behavior of the rainfall rate on the CReSS model, spectral analysis was performed, as shown in Figs. A2a and A2f. In general, the analyses among D1 and D2 show the fact that as the rainfall rate is the average of all the levels of altitude, scaling behavior is more stable than when the analysis was performed separately on each altitude of radar data. The approximate values of *β* and *R*^{2} for D1 and D2 are respectively *β* = 1.19 and *R*^{2} = 0.96 and *β* = 1.39 and *R*^{2} = 0.98, with no apparent scaling break for both of the domains.

Then TM analyses were also conducted with the CReSS model for the three domains studied. On the left side of Figs. A2b, A2c, A2g, and A2h, the scale invariance can be observed, with the resolution ranging from 1 to 512 for D1 and D2. The values of *q* are also freely set from 0.1 to 7 where the single scaling behaviors, with the imperfect but acceptable straight lines, are retrieved for all domains. On the right side of Fig. 14, the comparisons between empirical *K*(*q*) and semitheoretical *K*(*q*) from TM and DTM are well paired from *q* = 0.5 to *q* = 2.5 on both domains. The latter would be, therefore, the point of multifractal phase transition for the CReSS model data. Last, DTM analyses were performed also with the CReSS model, and the graphs are shown in Figs. A2d, A2e, A2i, and A2j. For each power *η*, the value of *q* was fixed to 1.5 as well as when analyzing radar data. Each α and *C*_{1} retrieved from DTM analysis is indicated in Table A1. The α values for both of the domains are not so different from each other, which indicates high singularities. The *C*_{1} values reflect the distribution of the rainfall. Thus, D1 is where the accumulated rainfall was distributed generally (lowest *C*_{1}), and D2 has concentratedly high accumulated rainfall (highest *C*_{1}).

Figure A3 shows the *γ*_{s} of CReSS simulation with respect to the size of the domain. The *α* from TM and DTM analyses are decreasing as the size of the domain decreases, whereas *C*_{1} is increasing. On the contrary, for the scales of 64 and 256, the estimate remains within the estimates of 512 for D2 domains. It also indicates that the bigger domain gives the better estimation, and it depends on where the domain is selected.

## REFERENCES

*Nonlinear Processes Geophys.*

*J. Korean Meteor. Soc.*

*Asia-Pac. J. Atmos. Sci.*

*Water Resour. Res.*

*J. Hydrol.*

*J. Hydrometeor.*

*Nonlinear Processes Geophys.*

*J. Appl. Meteor.*

*J. Geophys. Res.*

*Rev. Sci. Eau*

*Atmos. Res.*

*J. Meteor. Soc. Japan*

*Fractals in Geography*, L. De Cola and N. Lam, Eds., Prentice-Hall, 171–205

*Nonlinear Processes Geophys.*

*Atmos. Res.*

*IEEE Trans. Geosci. Remote Sens.*

*Workshop on Planetary Boundary Layer Parameterization*, Reading, United Kingdom, ECMWF, 59–79, https://www.ecmwf.int/file/23088/download?token=g93a5QuY

*Tellus*

*Nonlinear Dynamics in Geosciences*, Springer, 311–337

*Ann. Geophys.*

*Houille Blanche*

*J. Fluid Mech.*

*J. Geophys. Res.*

*J. Hydrometeor.*

*J. Geophys. Res.*

*Turbulence and Predictability in Geophysical Fluid Dynamics*, M. Ghill, R. Benzi, and G. Parisi, Eds., Elsevier, 111–114

*J. Korean Earth Sci. Soc.*

*Proc. Fourth Symp. on Turbulent Shear Flows*, Karlshule, Germany, University of Karlsruhe, 11.1–11.8

*Turbulent Shear Flow*, B. Launder, Ed., Vol. 4. Springer, 7–33

*J. Geophys. Res.*

*Physica A*

*Int. J. Bifurcat. Chaos*

*J. Meteor. Soc. Japan*

*J. Meteor. Soc. Japan*

*J. Meteor. Soc. Japan*

*Acta Meteor. Sin.*

*Weather Radar and Hydrology Symp.*, Exeter, United Kingdom, IAHS, 421–426

*J. Appl. Meteor.*

*High Performance Computing*, H. P. Zima et al., Eds., Springer, 243–259, https://doi.org/10.1007/3-540-47847-7_21

*J. Atmos. Sci.*

*Tellus*