Multifractal evaluation of simulated precipitation intensities from the COSMO NWP model

. The framework of universal multifractals (UM) characterizes the spatio-temporal variability in geophysical data over a wide range of scales with only a limited number of scale-invariant parameters. This work aims to clarify the link between multifractals (MFs) and more conventional weather descriptors and to show how they can be used to perform a multi-scale evaluation of model data. The ﬁrst part of this work focuses on a MF analysis of the climatology of precipitation intensities simulated by the COSMO numerical weather prediction model. Analysis of the spatial structure of the MF parameters, and their correlations with external meteorological and topographical descriptors, reveals that simulated precipitation tends to be smoother at higher altitudes, and that the mean intermittency is mostly inﬂuenced by the latitude. A hierarchical clustering was performed on the external descriptors, yielding three different clusters, which correspond roughly to Alpine/continental, Mediterranean and temperate regions. Distributions of MF parameters within these three clusters are shown to be statistically signiﬁcantly different, indicating that the MF signature of rain is indeed geographically dependent.Thesecond part of this work is event-based and focuses on the smaller scales. The MF parameters of precipitation intensities at the ground are compared with those obtained from the Swiss radar composite during three events corresponding to typical synoptic conditions over Switzerland. The results of this analysis show that the COSMO simulations exhibit spatial scaling breaks that are not present in the radar data, indicating that the model is not able to simulate the observed variability at all scales. A comparison of the operational one-moment microphysical parameterization scheme of COSMO with a more advanced two-moment scheme reveals that, while no scheme systematically outperforms the other, the two-moment scheme tends to produce larger extreme values and more discontinuous precipitation ﬁelds, which agree better with the radar composite.

Abstract. The framework of universal multifractals (UM) characterizes the spatio-temporal variability in geophysical data over a wide range of scales with only a limited number of scale-invariant parameters. This work aims to clarify the link between multifractals (MFs) and more conventional weather descriptors and to show how they can be used to perform a multi-scale evaluation of model data.
The first part of this work focuses on a MF analysis of the climatology of precipitation intensities simulated by the COSMO numerical weather prediction model. Analysis of the spatial structure of the MF parameters, and their correlations with external meteorological and topographical descriptors, reveals that simulated precipitation tends to be smoother at higher altitudes, and that the mean intermittency is mostly influenced by the latitude. A hierarchical clustering was performed on the external descriptors, yielding three different clusters, which correspond roughly to Alpine/continental, Mediterranean and temperate regions. Distributions of MF parameters within these three clusters are shown to be statistically significantly different, indicating that the MF signature of rain is indeed geographically dependent.
The second part of this work is event-based and focuses on the smaller scales. The MF parameters of precipitation intensities at the ground are compared with those obtained from the Swiss radar composite during three events corresponding to typical synoptic conditions over Switzerland. The results of this analysis show that the COSMO simulations exhibit spatial scaling breaks that are not present in the radar data, indicating that the model is not able to simulate the observed variability at all scales. A comparison of the operational one-moment microphysical parameteriza-tion scheme of COSMO with a more advanced two-moment scheme reveals that, while no scheme systematically outperforms the other, the two-moment scheme tends to produce larger extreme values and more discontinuous precipitation fields, which agree better with the radar composite.

Introduction
Validation of precipitation fields simulated by a numerical weather prediction model is a delicate task as reference data (rain gauges, radar scans) are typically available at a different spatial and temporal resolution than the model. Traditional point-based verification scores are generally unable to provide sufficient information about the forecast quality, as they do not take into consideration the spatial structure and are affected by the so-called "double penalty" (Gilleland et al., 2009). Indeed, small displacements in the simulated forecast features will be penalized twice, once for missing the observation and again for giving a false alarm. The impact of this double penalty is related to the variability in the simulated fields, which tends to increase with the resolution of the model. Numerous methods have been proposed in recent years to address this issue. Some methods rely on the use of traditional scores but applied on filtered fields, estimating the forecast performance as a function of scale and precipitation intensity (e.g., Ebert, 2008;Mittermaier et al., 2013), while others detect specific features of forecast and verification fields and compare these features based on their attributes (e.g., Davis et al., 2006;Wernli et al., 2008). Yet another type of validation technique relies on the separation Published by Copernicus Publications on behalf of the European Geosciences Union.
of scales with the use of space-frequency tools such as the 2-D wavelet transform (Vasić et al., 2007).
Multifractals (MFs) offer a convenient way to analyze the variability in complex geophysical systems globally over a wide range of scales. In the context of MFs, the statistical properties of a field are related to the resolution by a power law (Schertzer and Lovejoy, 1987). Universal multifractals (UMs) are a framework of MFs based on the concept of multiplicative cascades, which allows for analysis and simulation of a high variability across scales with only a small number of parameters with physical meaning (e.g., Schertzer and Lovejoy, 1987;Lovejoy and Schertzer, 2007). In meteorology, UMs have been used to study a large variety of complex natural phenomena such as the distribution of rainfall intensities at the ground (e.g., Marsan et al., 1996;Gires et al., 2015a, b), atmospheric turbulence (e.g., Parisi and Frisch, 1985a;Schertzer and Lovejoy, 2011) or climate change (e.g., Schmitt et al., 1995;Royer et al., 2008). Gires et al. (2011) used the UM framework to compare simulations of Meso-NH, a non-hydrostatic numerical weather prediction (NWP) model developed by Météo-France, with composite radar images during a heavy convective rainfall event. This comparison showed that the radar quantitative precipitation estimation (QPE) and the model simulations were generally characterized by similar ranges of scaling and agreed quite well with a simple space-time scaling model. Gires et al. (2011) focused on a rather flat area and a single event. It is therefore relevant to study the MF behavior of model simulations and radar observations over a more complex terrain and in a broader meteorological context. This work is thus divided into two parts. The first part aims to illustrate the use of MFs for characterizing regional patterns of precipitation and to relate MFs to meteorological and topographical features. To this end, a large-scale analysis of 5 years of simulated precipitation intensities from the COSMO (Consortium for Small-scale Modeling) numerical weather prediction model is conducted. The MF properties of the corresponding climatology of precipitation intensities are then studied and related to several regional and meteorological descriptors. The second part of this paper extends the work of Gires et al. (2011) over Switzerland for three different meteorological situations (snowstorm, convective summer precipitation and stratiform rainfall) by using simulated precipitation intensities from the COSMO model run with different microphysical parameterizations.
This article is structured as follows: in Sect. 2, the COSMO model and the Swiss radar composite are briefly described. The studied events, radar datasets and model variables are then described in detail. Section 3 provides a short summary of the UM framework. In Sect. 4, a climatological analysis of precipitation intensities simulated by COSMO is performed with the UM framework in relation to external geographical and meteorological descriptors. In Sect. 5, a spatial and temporal analysis of precipitation intensities on the ground simulated by COSMO is performed during three characteristic events. The results are then compared with the UM analysis of the radar composite. Finally, Sect. 6 gives a summary of the main results and concludes this work.
2 Description of the data

The COSMO model
The COSMO model is a mesoscale limited-area numerical weather prediction model initially developed as the Lokal Modell (LM) at the Deutscher Wetterdienst (DWD). It is now operated and developed by various weather services in Europe, including Switzerland. Besides its operational applications, it is also used for scientific purposes in weather prediction and for regional climate simulations. The COSMO model is a non-hydrostatic model based on the fully compressible primitive equations integrated using a split-explicit third-order Runge-Kutta scheme (Wicker and Skamarock, 2002). The spatial discretization is based on a fifth-order upstream advection scheme on an Arakawa C-grid with Lorenz vertical staggering. Height-based Gal-Chen coordinates are used in the vertical (Gal-Chen and Somerville, 1975). The model uses a rotated coordinate system where the pole is displaced to ensure approximatively horizontal resolution over the model domain. Subgrid-scale processes are taken into account with parameterizations. In particular, grid-scale clouds and precipitation are parameterized operationally with a onemoment scheme similar to Rutledge and Hobbs (1983) and Lin et al. (1983), with five hydrometeor categories: rain, snow, graupel, ice crystals and cloud droplets. Snow is assumed to be in the form of rimed aggregates of ice crystals that have become large enough to have an appreciable fall velocity. Cloud ice is assumed to be in the form of small hexagonal plates that are suspended in the air and have no appreciable fall velocity. The particle size distributions (PSDs) are assumed to be exponential for all hydrometeors except for rain where a gamma PSD is used: where N is expressed in units of m −3 mm −1 , D is the equivolume diameter in mm, N 0 is the intercept parameter (m −3 mm −1−µ ), the slope parameter (mm −1 ) and µ the dimensionless shape parameter.
In the operational one-moment scheme, the only free parameter of the PSDs is the slope parameter which can be obtained from the prognostic moment of order three (mass concentrations). The intercept parameter N 0 is either assumed to be constant, or in the case of snow, to be temperature dependent. The scale parameter µ is set to zero (exponential PSDs) for all hydrometeors except for rain where it can be chosen a priori and is set to 0.5 by default.
A more advanced two-moment scheme, which adds hail as a sixth hydrometeor category, was developed for COSMO Atmos. Chem. Phys., 17, 14253-14273, 2017 www.atmos-chem-phys.net/17/14253/2017/ by Seifert and Beheng (2006). In this scheme all PSDs are assumed to be gamma distributions where the intercept and slope parameters are free parameters that can be obtained from the prognostic moments of order zero (number concentrations) and order three (mass concentrations). As this scheme significantly increases the overall computation time, it is currently not used operationally.
In COSMO, the interaction of various microphysical processes and their feedback on the simulated flow fields are represented by a system of budget equations for q x , the specific mass fraction in kg x per kg air of hydrometeor x.
where S x represents the microphysical sources and sink per unit mass of moist air, F x the turbulent fluxes and P x the precipitation or sedimentation fluxes defined by P T is the terminal fall velocity of hydrometeor j . The precipitation intensity at the ground is then simply the sum of the sedimentation fluxes of all hydrometeors at the lowest model level. In terms of terminal velocities COSMO generally assumes power laws: v T = aD b , where D is the particle equivolume diameter. Note that in the two-moment scheme a more refined empirical relation by Rogers et al. (1993) is used for raindrops. Numerically, this system of differential equations is treated with a time splitting method, in which the advection terms v · ∇q x are first integrated over a COSMO time step (20 s) and the budget equations are then solved for the microphysical source terms and sedimentation only. In the operational microphysical scheme the source terms include (1) nucleation and depositional growth of cloud ice, (2) autoconversion of cloud water to rain, (3) collection mechanisms, (4) diffusional growth of rain and snow and (5) melting and freezing mechanisms. Details about the parameterization of all these source terms can be found in Doms et al. (2011).
In the operational setup, the COSMO model uses a prognostic turbulent kinetic energy (TKE) closure at level 2.5 for the parameterization of atmospheric turbulence. This scheme is similar to Mellor and Yamada (1982), the main difference being the use of variables that are conserved under moist adiabatic processes: total cloud water and liquid water potential temperature. Additionally, a so-called "circulation term" is included, which describes the transfer of nonturbulent subgrid kinetic energy from larger-scale circulation towards TKE. The reader is referred to Baldauf et al. (2011) and the model documentation (Doms et al., 2011) for a more in-depth description of the various COSMO subgrid parameterizations.

Climatological study
In the first part of this work, a MF characterization of all precipitation intensities simulated by COSMO during 5 years is performed. The computed MF parameters are then com-pared with various descriptors. A total of 43 115 hourly time steps of COSMO simulations in analysis mode covering a period of 5 years (2011 to 2016) were retrieved from the Me-teoSwiss archives. These data are available at a 2 km resolution (0.02 • angular resolution) using a setup known operationally as "COSMO-2" (COSMO, 2015), which was used for forecast and analysis until the beginning of 2016. The first and last 100 km of the COSMO-2 domain along both longitudinal and latitudinal directions were discarded from the analysis in order to avoid border effects. The remaining area is shown in Fig. 2 (Domain 1). This area was then divided into 209 subsquares of size 64 × 64 (128 km× 128 km), with an overlap of two-thirds between consecutive squares. This ratio has been chosen as a trade-off between representativity (total number of squares) and computation time.
Besides the MF parameters, which will be described later on, 11 local descriptors of the geography and meteorology were computed from the COSMO data within every square (Table 1). Spatial maps of all descriptors are given in Fig. A1 in Appendix A.

QPE comparison 2.3.1 Simulation of events
In the second part of this work, three precipitation events were simulated with COSMO and compared to the radar QPE in terms of the MF properties of their rainfall intensity fields. These events correspond to three typical meteorological situations observed over Switzerland. A brief description of the events is given in Table 2, and 500 hPa geopotential and temperature charts are shown in Fig. 1. To simulate these events, COSMO was used in its version 5.01 with the standard "COSMO-2" setup (COSMO, 2015). As was done in similar studies (Bohme et al., 2009), a spin-up time of 12 h was used to account for the cold start of the model. For the initial and boundary conditions, analysis forcings of MeteoSwiss, obtained with the COSMO-7 model run at a 7 km spatial resolution, were used. This allows for forcing the model with the most accurate information available at the time of simulation. In addition, the events were also simulated using the non-operational two-moment scheme, while keeping all other namelist parameters unchanged. For all simulations, model outputs were written every 5 min of simulated time, which corresponds to the temporal resolution of the Swiss radar composite.

QPE data
In Sect. 5, precipitation intensities at the ground simulated by the COSMO are compared with the QPE product from the Swiss operational radar composite. The Swiss radar composite consists of plane position indicator measurements of  Strong summer convection triggered by the presence of very warm and wet subtropical air over Switzerland the four 1 operational polarimetric C-band radars. The QPE product of MeteoSwiss is computed in the following way: the linear equivalent radar reflectivity measurements at up to six 1 • × 1 • × 83 m clutter-free radar bins are corrected for partial beam-blocking and averaged to derive polar 1 • × 1 • × 500 m radar bins. Reflectivity measurements are then converted to equivalent precipitation intensity with a Marshall and Palmer (1948) Z-R relationship. The precipitation estimation at the ground is extrapolated from multi-radar observations aloft using a weighting function that depends on the altitude above the ground and the radar visibility. The correction for the vertical profile of reflectivity (VPR) is done with an average profile based on aggregation over a few hours and over the visible part of the area located less than 70 km from the radar. More information on the MeteoSwiss QPE estimation can be found in Germann et al. (2006). Note that the Plaine-Morte radar was only installed in 2014 and was thus not available during the first event (26 March 2010). The Swiss radar composite extends radially up to 250 km from every single radar (Fig. 2). However, the quality of the product improves with the proximity to the radar and in the areas where the radar scanning domains overlap. To perform a comparison of rain 1 The Weissfluhgipfel radar was not yet installed at the time of the considered events intensities, a smaller field of 128 × 128 km 2 was chosen in the center of the domain where the quality of the product is optimal (Domain 2 in Fig. 2). For the second event (8 April 2014) the domain was moved slightly to the left in order to better follow the evolution of the precipitation event (Domain 3 in Fig. 2).

Multifractality
Let be a normalized (divided by its mean) conservative field, which can be one-or two-dimensional (time series or spatial map). The field also needs to have the same length in all dimensions, and this length must be a power of two.
The fractal dimension D f of indicates how a corresponding binary field (where all values larger than a given threshold are set to 1) scales with the resolution. The resolution is defined by the ratio between outer scale L and observation scale l (λ = L/ l). where N λ is the number of positive samples (rainy pixels for example) at a given resolution, which can be obtained with the help of box counting. It is possible to interpret this result in a probabilistic way. Indeed, consider a line or cube of size l, then p is the probability that it intersects the field. This probability scales with the following resolution: where D is the dimension of the field (1 for a time series, 2 for a spatial field) and c f = D − D f is called the fractal codimension of the field. It is clear that the values of D f (and c f ) depend on the threshold that is used. Several thresholds and corresponding values of D f are thus required to characterize the field. In the MF framework, this characterization is performed with scale-independent thresholds (Schertzer and Lovejoy, 1987), which yields a MF version of Eq. (4).
where c(γ ) is the codimension function, which is convex and increasing and γ is a so-called singularity, which is independent of scale. λ γ can thus be seen as a scale-invariant threshold. Since c(γ ) takes several values, MFs can be understood as a hierarchy of embedded fractal sets. It can be shown that Eq. (5) implies that the statistical moments q of the field scale with the resolution (Schertzer and Lovejoy, 1987): where K(q) is the moment scaling function which is related to c(γ ) by a Legendre transform (Parisi and Frisch, 1985b). For a conservative field, λ = 1. The quality of the scaling can be studied with the trace moment (TM) method which consists of a log-log plot of the up-scaled fields as a function of the resolution λ for each moment q, the slope being the moment scaling function.
In the UM framework (Schertzer and Lovejoy, 1987), the moment scaling function for a conservative field K c (q) can be fully characterized with only two parameters, α and C 1 .
C 1 is the mean intermittency codimension and measures the clustering of the average intensity at increasing scales. C 1 is equal to zero when the field is homogeneous. α is the multifractality index and measures the clustering variability with respect to the intensity level; α ∈ [0, 2]. The size of the sample limits the insight one can get of a statistical process. For a MF processes, if N s samples are available this will result in a maximum singularity γ s (and corresponding moment order q s ) beyond which the values of the statistical estimates of the codimension and moment scaling function are not considered as reliable (e.g., Schertzer and Lovejoy, 1987;Lovejoy and Schertzer, 2007). It can be shown that in the MF framework, q s and γ s are equal to the following: where 1 α + 1 α = 1 and D s is the sampling dimension defined by N s = λ D s .
An example of the use of γ s can be found in Royer et al. (2008) who investigated the impact of climate change on rainfall extremes with a climate model. They observed an increase in γ s over time, which could result in a possible increase in the intensity of rainfall extremes over the next hundred years. Douglas and Barros (2003) and Hubert et al. (1993) also used the maximum singularity γ s in the estimation of probable maximum precipitation.
In this work, the UM parameters are estimated with the double TM method (Lavallée et al., 1993). This method relies on the fact that, in the context of UM, the moment scaling function K(q, η) of the field (η) λ , obtained by raising the field at a power η and up-scaled at resolution λ, can easily be expressed as a function of α (Lavallée et al., 1993).
α is thus the slope of the linear part of K(q, η) as a function of η in a log-log plot.

Non-conservative fields
In the case of a non-conservative field φ, we have φ λ = 1. One way to consider non-conservative fields within the UM framework it to assume that they can be expressed as the following: where H is the non-conservation parameter (H = 0 for conservative fields) and is a conservative field characterized by a moment scaling function K c (q) with parameters C 1 and α.
The moment scaling function of the non-conservative field φ λ is then given by H can be related to the spectral slope β by where β is the exponent of the power law that characterizes the relation between power spectrum and wave numbers: Hence the larger the value of the slope β, the shorter the decorrelation range. If β is larger than the dimension of the field, the field is non-conservative. λ can be estimated from φ λ with a fractional integration (for H < 0) or differentiation (for H > 0) of order H , which is equivalent to a multiplication by k H in Fourier space. In practice however, for H > 0, particularly when H > 0.5, λ max (the field at the maximum resolution) is often approx-  When α = 0, the field is mono-fractal which means that a single fractal dimension is sufficient to fully characterize the field. The higher α, the larger the variability within areas with precipitation. C 1 mean intermittency Measures how concentrated the average field is, C 1 = 0 for a uniform field. The larger C 1 , the larger the intermittency of the field. High C 1 associated with high α implies strong extremes. β Negative of the spectral slope The larger β the shorter the decorrelation range of the data. If β = 0, the power spectral density is the one of white noise, meaning there is no decorrelation. When β is large, large-scale phenomena have a large contribution to the variability in the data, which means that the fields have a larger correlation range (smoother fields). H non-conservation parameter Scale-independent proportionality factor relating a conservative field to the nonconservative field (a field for which the average is not scale independent). Can be seen as a smoothing factor. If H is positive, the field is too smooth and one needs to differentiate it to retrieve a conservative field. If H is negative, the field is too discontinuous and one needs to integrate it to retrieve a conservative field. γ s maximum singularity Maximum observable singularity (scale independent threshold) from the data. Large γ s implies that stronger extremes are present in the data. R 2 TM coefficient of determination Coefficient of determination of the relation between a given moment of order q and the scale (Eq. 6) on a log-log plot. If the field is MF, the plot is a straight line with slope K(q). In practice, q = 1.5 is often used as a reference to determine R 2 . R 2 can be seen as an estimation of the quality of scaling and of the viability of the MF approach for the data. D f Fractal dimension Factor relating the number of rainy observations at a given scale to the corresponding scale. The larger D f , the more uniform the binary precipitation field (obtained with a threshold of zero) is.
imated by the renormalized absolute fluctuations of the field. Table 3 provides an overview of all MF parameters and their interpretation, while Appendix B provides some visual examples of how varying the MF parameters affects the spatial structure of randomly generated MF fields.

Spatio-temporal analysis
The MF analysis of time series of two-dimensional fields, such as the ones considered in this study, can be performed both in space, by considering an ensemble of twodimensional fields (one sample for every time step), or in time, by considering an ensemble of one-dimensional time series (one sample for every coordinate in the twodimensional field). A simple spatio-temporal scaling model (e.g., Marsan et al., 1996;Deidda, 2000;Macor et al., 2007;Radkevich et al., 2008) is based on the hypothesis of an anisotropy coefficient between space and time: where H t is the anisotropy coefficient between space and time, which in the theory of Kolmogorov (e.g., Kolmogorov, 1962;Marsan et al., 1996) is equal to one-third. This result implies identical α and proportional C 1 and H parameters.

Climatological analysis of MF parameters
Within all 209 selected squares (Sect. 2.2), the MF parameters α, C 1 , β, H , γ s , R 2 and D f were computed both in time and in space by performing an ensemble average, i.e., by taking the mean parameters over all the available realizations of the process. In space, every time step is considered as a realization of a two-dimensional geophysical process, while in time every COSMO grid point is considered as a realization of a one-dimensional geophysical process (time series). For the analysis in space, this implies 43 115 realizations of a 64 × 64 2-D field and for the analysis in time, 4096 realizations (64 × 64) of a time series of length 32 768 (2 15 , which is the closest power of 2 to 43 115). Analysis of the overall temporal power spectrum (not displayed) revealed the presence of a strong peak centered around a resolution of 3 h, which is very likely to be caused by the nudging scheme of COSMO (assimilation within the model is performed every 3 h). To remove the effect of the nudging in the estimation of MF parameters, only the larger timescales (from 8 h to 5 years) were considered. Even though some areas are characterized by values of H > 0.5 (non-conservative fields), it was decided to treat all areas in a consistent way, by working on the original fields instead of the fluctuations. Indeed, it was observed that using fluctuations for areas where H > 0.5 was causing important discontinuities in the spatial structure of MF parameters. Indeed, using fluctuations is only a crude way to address the non-conservativity of the fields. A proper correction using fractional integration would allow for a much smoother correction, since it is proportional to H , but is computationally intractable because of the very large dataset that is used. In order to test the effect of zeros on the overall analysis, the MF parameters were also estimated in space 2 by using only the fields where there is precipitation over at least 50 % of the surface. This did, however, not impact the main conclusions in terms of correlations and spatial structure of the MF parameters. Hence the subsequent study was performed on the raw precipitation fields without any kind of filtering.

Correlation study
In the first step, the relationships between MF parameters evaluated in time and space and the descriptors detailed in Table 1 were analyzed by looking at the non-parametric (Spearman) correlations. Figure 3 displays the correlation plots for the MF parameters in space and time.
The following conclusions can be drawn from the correlation plots. D f is strongly correlated with the latitude, indicating that the fractal dimension of rain is larger at higher latitudes (i.e., cooler climates). D f is also quite strongly correlated to the wet fraction: the more often it rains, the higher the fractal dimension. Similarly, C 1 also has a strong latitudinal trend and tends to decrease in regions with a high wet fraction (less intermittency), which is typically the case in cooler climates. α also seems to show a latitudinal trend, though not as strong as for C 1 and D f . The link between α and the standard deviation is not obvious since the standard deviation is a second-order statistic, while the UM parameters are based on all moments. C 1 and γ s tend to decrease with the altitude, while β and H tend to increase. This could be due to orographic effects, precipitation over mountains being both dominated by large-scale circulations and generally abundant.
The correlation values are roughly consistent in time and in space, although the ones in time are generally higher. In time the correlation between β and H is lower than in space. This can be explained by the lower values of β in time. Be-2 In time, it would not be possible to filter out non rainy time steps, as it would break the continuity of the times series.   Note that the part on the right of the correlation plot in time has been truncated since it is the same as on the correlation plot in space. Correlations that are not statistically significant (for a significance level of 5 %) and correlations which are not consistent in sign between time and space are left in white.
cause of this, and as a consequence of Eq. (12), H becomes more correlated with α and C 1 .

Hierarchical clustering
A hierarchical clustering of all 209 areas was performed based on the value of their descriptors (Table 1), using the Ward linkage method (Ward, 1963). Investigation of the dendrogram gives an optimal number of clusters of either three Atmos. Chem. Phys., 17, 14253-14273, 2017 www.atmos-chem-phys.net/17/14253/2017/ or six. Figure 4 shows the resulting classification for three different clusters. As could be expected, the clustering of the meteorological and topographical descriptors results in a meaningful spatial distribution. Indeed, all clusters are spatially very coherent, with Cluster 1 corresponding mostly to the Alpine regions, from the Mediterranean Sea to Austria, Cluster 2 corresponding mostly to the cooler temperate regions in the east of France and south of Germany, and Cluster 3 corresponding to the warmer Mediterranean regions in the south of France, in Italy and the Balkans. Note that over land areas, a somewhat similar classification can be obtained by aggregating clusters of the famous Köppen (1936) climate classification. A comparison between the meteorological classification and the Köppen classification is performed in Appendix C.
The distributions of MF parameters within these clusters, as illustrated in Fig. 5, highlight some obvious discrepancies in MF parameters between clusters. In space, α, H , β and D f seem to be lower in cluster 3 than in the others, while C 1 and γ s are higher. Clusters 1 and 2 do not exhibit such marked differences but differ nonetheless by stronger values of D f and slightly larger values of H , β and α. In time, cluster 1 (Alpine regions) is characterized by high values of α and β and low values of C 1 , which indicates frequent rainfall and a large variability in intensities. Cluster 2 differs by its low values of β and H , indicating that small temporal frequencies play a bigger role in the overall variability in precipitation.
The statistical significance of these discrepancies was confirmed both with the MANOVA (multivariate analysis of variance) and the non-parametric Kruskal-Wallis statistical tests. All tests were performed with a significance level of 2.5 %. The MANOVA reveals that the multivariate means of MF parameters (both in time and space) are significantly different between the three clusters, as well as between all three pairs of clusters taken separately (1 vs. 2, 2 vs. 3 and 1 vs. 3). The non-parametric Kruskal-Wallis test, performed separately for all MF parameters, reveals that distributions of all MF parameters are significantly different between the three clusters. A pairwise comparison (1 vs. 2, 2 vs. 3 and 1 vs. 3) reveals that, in time, all MF parameters are significantly different between all pairs of clusters, except for H which is not statistically different between clusters 1 and 2. In space, the situation is similar, but this time it is C 1 and γ s which are not significantly different between clusters 1 and 2.
To summarize, the statistical analysis shows that the MF parameters of precipitation intensities are significantly different within the three climatological clusters. Figure 6 shows the spatial structure of MF parameters in space and reveals that α is particularly large over Bavaria, the Piemont region of Italy and the Champagne region of France. However, these "clusters" of large values are quite difficult   Table 1. to relate to other trends, especially in terms of descriptors. Generally, the Swiss Alps are characterized by relatively low values of α. In terms of C 1 , there is a clearer trend, with a strong latitudinal and elevational negative gradient, which can be related to an increase in overall precipitation amounts. In terms of β, the regions over Italy show smaller values, which tend to indicate that precipitation events exhibit strong high-frequency components (such as for example short convective storms). Over the Alps and the more continental regions, where the precipitation systems are mostly frontal, β increases. The trend is similar for H , which is related to β. For γ s , we observe a similar trend to C 1 with a decrease with altitude and elevation. This indicates that precipitation extremes are stronger in southern regions. For R 2 it is difficult to identify a proper spatial structure, except that the foothill regions of the Alps generally have large R 2 values. Finally for D f , there is mostly a latitudinal trend: the value of D f increases with the latitude, which can be related to a decrease in the number of zeros (smaller C 1 and generally larger α).

Spatial structure of MF parameters
In time, similar conclusions can be drawn for C 1 , D f , γ s and β. For the other parameters, the latitudinal trend seems to be much more visible than in space, especially for H and R 2 . R 2 seems to be generally larger in time than in space. Unlike in space, for α in time there is a clearer dependency on the altitude. α in time tends to be larger in mountainous regions, indicating a larger temporal variability in precipitation intensities in these regions.
To summarize, the MF signature of precipitation is related to the topography, the climate and the typical meteorological conditions. As such, MFs can be used as a way to characterize precipitation fields and assess the realism of simulated atmospheric variables.

Comparison of simulated precipitation with radar QPE
The previous section presented the MF approach in a climatological context, which helps to link meteorology/geography and MF parameters. In the present section, we evaluate the quality of the precipitation simulated by COSMO with two different microphysical schemes, by comparing it with quantitative precipitation estimations from the Swiss radar network using the UM framework. This comparison requires the COSMO model to be run at the radar temporal resolution (5 min) in a very expensive setup (twomoment scheme). As such, a climatological comparison over several years is not feasible from a computational point of view. Hence, the comparison is now conducted on the event scale, which also makes the intercomparison of the microphysical parameterizations easier. The power law Z-R relationship by Marshall and Palmer (1948), which is used by MeteoSwiss to derive the QPE, does not correspond with the Z-R relationships derived from the COSMO microphysical parameterizations. This explains part of the discrepancy between radar QPE and simulated precipitation intensities, but in general should not impact the validity of the MF comparison. Indeed, if the Z-R relationship derived from the COSMO parameterizations can be approximated by a power law, then the correction needed to account for discrepancies in Z-R relationships is itself a power law: R corr = aR b , where R corr is the precipitation intensity one would obtain by first converting COSMO precipitation intensities to reflectivities and then back to precipitation intensities using the radar QPE Z-R relationship. It can be shown  Figure 6. Spatial representation of the MF parameters estimated in space for all areas. The special colormap for H has been chosen to separate positive and negative values and highlight non-conservative areas where |H | > 0.5. For R 2 , circles represent all zones where R 2 < 0.9, whereas squares represent all zones with good scaling (R 2 ≥ 0.9). As previously mentioned, the size of the represented squares is not to scale. The colors in the background correspond to the hierarchical classification.  that in the context of UM, the corrected field will have the same value of α and the same scaling properties as the original field, while C 1 will be multiplied by b α . Moreover, for the one-moment scheme, it was observed that while the intercept parameter a changes significantly, the exponent parameter b is almost the same: R corr = 0.68R 0.98 . As the exponent 0.98 is close to unity, this implies an almost direct proportionality, and as such even C 1 should barely be affected. Note that this power law was derived by using the T-matrix method (Mishchenko et al., 1996) to compute radar cross sections at the C-band.
For the two-moment scheme, things are more complicated as no one-to-one relationship exists between rain rate and reflectivity. However, a rough estimation of the error in C 1 was derived by considering a representative set of rainy time steps from all events. The estimated values of b for the two-moment scheme varied between 0.81 and 1.23, which would imply a maximum relative error in C 1 of 51 % on spatial fields and 23 % on time series.
Overall, correcting precipitation fields for discrepancies in the Z-R relation is challenging, especially in the solid phase and for the two-moment scheme, where deriving a Z-R relation from the model parameterization is difficult. However, in the MF context, only C 1 should be affected and only with significant solid precipitation or when using the two-moment scheme. This should be kept in mind when interpreting C 1 values.

Scaling analysis
A MF comparison of the precipitation fields simulated by COSMO in its one-moment and two-moment schemes with the QPE product from the Swiss radar composite was performed. As a first step, a spectral analysis was performed both in time (ensemble of one-dimensional time series of precipitation intensities) and in space (ensemble of twodimensional maps of precipitation intensities). Figure 7 shows the spectral analysis in space for all events and data. A best-fit line is shown for the radar QPE from which the value of β is computed. β is equal to −m, where m is the slope of the best-fit line.
For the 26 March 2010, we observe a single scaling regime for the radar QPE, with a good scaling both on large (16-64 km) and small scales (2-16 km), as the spread around the line is relatively small. For the model intensities, we observe strong discrepancies with the radar QPE in terms of spectral slope at smaller scales (2-8 km), which are not well represented. A possible explanation for this break in scaling properties of the model is the fact that large scales are dominated by the dynamics of the model (primitive equations of the atmosphere), whereas smaller scales are dominated by the parameterizations of subgrid phenomena (turbulence, convection). However, even at larger scales (8-64 km), the agreement between radar QPE and model simulations is still quite poor in terms of spectral slope. Obviously, for this rainfall event, COSMO is not able to recreate the spatial structure of precipitation observed by the radar.
For the 8 April 2014, the scaling is similar between radar and model precipitation intensities, possibly indicating that for this stratiform rain event, parameterizations and dynamics match better. Both radar and simulations show a weak scaling break at around 8 km.
For the last event, we observe again a good scaling for the radar QPE and a much worse scaling for the model precipitation intensities. However, in contrast with the first event, the larger scales (> 16 km) are not well represented. Indeed, inspection of the time series of precipitation shows that COSMO is not able to locate accurately the convective cells of precipitation and generally overestimates their extent. In terms of microphysical parameterizations, we observe that the spectral slopes of the one-moment scheme are gener-128 Figure 7. Spectral analysis in space of the QPE products during the three events. The displayed lines are best-fit lines that take into account a possible scaling break. The associated value of β is given in the legend. Note that the maximal represented observation scale, which corresponds to the Nyquist frequency, is twice the maximal resolution.
ally closer to the ones obtained from the radar QPE. This is especially visible for the last (convective) event, where the COSMO simulations show weak scaling (β close to zero). This implies that the simulated rainfall intensities are dominated by small-scale features, while large-scale features are underestimated. Note also that for large-scale features, the power density function of COSMO simulations corresponds to white noise, indicating that the COSMO model has a shorter decorrelation range than the radar data.
The spectral analysis in time (not displayed) generally shows similar results, but with larger values of β and overall better scaling (less spread). Table 4 displays the non-conservation parameter H and the fractal dimension D f evaluated for time series of precipitation intensities (analysis in time) and for spatial fields of precipitation intensities (analysis in space), for both the radar QPE (in regular font), the COSMO one-moment scheme (in Atmos. Chem. Phys., 17, 14253-14273, 2017 www.atmos-chem-phys.net/17/14253/2017/ bold), and the COSMO two-moment scheme (in italic) and for all events. A value of H larger than 0 indicates that the field is smoother than the observed field from a direct MF cascade process, and a value of H smaller than 0 indicates that the field is too discontinuous. The larger D f , the more uniform the binary precipitation field is. In space D f = 2 implies that it rains everywhere, while in time D f = 1, implies that it rains all the time. Taking the radar as a reference, one sees that the convective event is characterized by the largest values of H , followed by the snowfall event and the stratiform event. When comparing H between COSMO and the radar QPE, one observes that H in time is always larger with the one-moment scheme than on the radar QPE, indicating that the temporal structure of the simulated fields is likely to be too smooth. In space, the trend is not as obvious and the match between the radar QPE and the precipitation intensities simulated with the one-moment scheme seems better. In terms of fractal dimension, it can be seen that for all events, both in space and time, the radar QPE has the most discontinuous binary precipitation field due to its smaller values of D f . COSMO simulations are characterized by larger values of D f , indicating a wider coverage of precipitation and less convoluted precipitation fields and time series. It is interesting to notice that the two-moment scheme gives values of D f that are closer to those of QPE, which indicates that it is better at simulating small-scale variations in the temporal and spatial occurrence of precipitation. Overall, it is worth noticing that the two-moment scheme almost always has smaller H values than the one-moment scheme, which indicates that it is more discontinuous both in time and space.
In order to account for the fact that the fields are mostly non-conservative (|H | > 0.5), and to treat all fields in a consistent way, all further analysis were performed on fluctuations of the original fields (Eq. 14). Note that while this does not result in perfectly conservative fields, it still makes them more conservative since all values of |H | are smaller than 0.5 after taking the fluctuations. Figure 8 shows the TM analysis in space and time for the three events and for a value of q = 1.5. The value of K(q = 1.5) is the slope of the best-fit lines. Repeating the TM analysis for various values of q allows us to characterize K(q) and to estimate α and C 1 . For the first two events, a scaling break can be observed at large scales for the COSMO precipitation intensities (64-128 km). These scales were excluded from the analysis due to the limited number of points in this scale range. For consistency with the observations of the spectral analysis, the scale range of 2-8 km was excluded from the analysis for the first event. This scale range does not scale well on COSMO simulations when compared with radar observation in this event. For the last event, two scaling regimes are observed for the COSMO intensities (2-16 and 16-128 km), which were studied separately. On the contrary, for the radar QPE no scaling break is observed in space. In time, a weak scaling break can be observed both for radar and COSMO intensities at a resolution of around 160 min. Hence, results are discussed only for the timescales between 5 and 160 min (smaller scales).

Spatio-temporal analysis
Values of α, C 1 and γ s obtained with an analysis in time and in space of the three events are given in Fig. 9. For the first two events, all parameters are computed only on the smaller scales (up to 64 km in space and up to 160 min in time), in order to account for the observed scaling break. For the last event both scale ranges are considered.
For the first event, both COSMO microphysical schemes give very similar MF parameters and the discrepancy with the radar QPE is quite important. In space, it can be observed that α is slightly smaller in the COSMO simulations than on the radar QPE. It is clear as well that the simulated C 1 is too small compared with the radar observations. This tends to indicate that COSMO is underestimating the spatial intermittency. Generally, the observed discrepancies in α and C 1 tend to indicate that the spatial structure of the simulated fields is too smooth and lacks the variability observed by the radars. In time, the agreement is better for C 1 , but COSMO has clearly higher values of α indicating a larger temporal variability than the radar QPE. For this event, there is a noticeable discrepancy between the maximum singularity γ s in space obtained from the radar QPE (0.721) and the γ s obtained from the model (around 0.6 for both schemes). This indicates that during this event COSMO had a tendency to underestimate extreme values, which might be caused by its difficulty to accurately simulate snowfall events, since COSMO does not consider partially melted snow (Frick and Wernli, 2012). Note that QPE of snow is very difficult and it is likely that the radar QPE itself is already underestimating precipi- tation intensities (Speirs et al., 2016), which would make this difference in γ s even more noteworthy.
For the stratiform rain event, the MF parameters of the COSMO simulations are in better agreement with the radar QPE. In time, the two-moment COSMO scheme gives values that are in relatively close agreement with the radar QPE and, in this regard, outperforms the one-moment scheme. COSMO simulations show generally smaller values of α and C 1 than the radar QPE, which is a trend that is observed for all events.
For the last convective event, two scaling regimes are considered in space: larger scales (16-128 km) and smaller scales (2-16 km). As already observed in the spectral analysis, there is a better agreement between the radar observations and the simulations with the one-moment scheme at smaller spatial scales. In time, however, the temporal intermittency of COSMO is smaller than for the radar QPE, which can be explained by the fact that COSMO generally overestimates the extent of the convective systems. Compared with the one-moment scheme and the radar QPE, the two-moment scheme has a smaller α in space but a larger α in time, as well as a smaller intermittency in time and space.
In summary, the observations of the spatio-temporal analysis are consistent with the spectral and scaling analysis where (1) a strong discrepancy in scaling behavior was observed between COSMO and the radar QPE at small scales for the first event, (2) a better scaling of the model precipitation intensities was observed for the second event, (3) a discrepancy in scaling at large scales was observed between COSMO (especially for the two-moment scheme) and the radar QPE for the third event.
Overall, it can be observed that except for the first event where both schemes give similar values, the two-moment scheme is usually characterized by a larger C 1 than the onemoment scheme, both in time and space, whereas in terms of α there is no recurring trend. However, one must keep in mind that this difference in C 1 is within the expected range of uncertainty for this parameter (see second paragraph of Sect. 5). For the MF parameters α and C 1 , there is generally a good agreement between radar observations and simulations in the range of scales where the model exhibits a good scaling behavior. In terms of α and C 1 , none of the two microphysical schemes seem to perform significantly better than the other. The two-moment scheme is, however, generally characterized by a slightly larger maximum singularity γ s , indicating a better capacity to simulate extreme values. This is especially visible in the last convective event. In terms of space/time ratios, the observed ratios differ significantly from the theoretical model: the α space/time ratio is always larger and the C 1 space/time ratio always smaller than the theoretical values (1 and 1.44 respectively). Gires et al. (2011) found different breaks for a Cevenol event (strong precipitation events occurring in fall in the south of France): roughly 16 km in space and 1 h in time, as well as a better agreement with a simple space-time model, but only for large scales, which are not the primary focus  Figure 9. α, C 1 and γ s parameter values obtained with an analysis in time and space on the fluctuations of the precipitation intensities for the three events. For the last event, both the parameters at large and small spatial scales are displayed. The numbers in blue are the space/time ratios for α and C 1 . of this study. These differences could be associated with the fact that the topography of the area analyzed in this paper is more pronounced than in Gires et al. (2011). It should also be noted that, in our case, the values of UM parameters α and C 1 exhibit a better agreement between observations and model simulations in the relevant range of scales.

Time series of UM parameters
To compare time series of UM parameters α and C 1 , the focus is put on the third (convective) event which shows the largest temporal variability. It was observed that the conclusions drawn for the third event in terms of discrepancies between radar and model MF parameters can be generalized to all events. Figure 10 shows the time series of α and C 1 throughout the third (convective) event for the COSMO and the radar QPE precipitation intensities, as well as some illustrative precipitation fields that will be discussed.

August 2015
During the convective event, four different phases can be identified. In the first short phase (12:00-14:00 UTC), observations and simulations agree relatively well in α and C 1 . This period corresponds to the initial stages of the event when only a few isolated cells are present (Fig. 10a). In the second phase (14:00-17:00 UTC), a large convective system is crossing the domain on the radar observations, which causes a strong increase in α and a decrease in C 1 . This convective system is, however, located more in the south in the simulation and enters the domain only at around 15:30 UTC (Fig. 10b). During the third phase (17:00-21:00 UTC), the large convective system is visible on the simulated field, whereas on the observed radar fields, the most intense convective cells are already out of the domain. This causes a larger α in the simulations than in the observations (Fig. 10c). Finally, during the last phase (21:00-24:00 UTC), a new convective system is visible on the observed field but is more or less absent on the simulated fields. This causes a discrepancy, the simulated fields having a smaller α and a larger C 1 than the observations (Fig. 10d). For this event, the spa- tial and temporal shifts between the convective system simulated by COSMO and the real convective system observed by the radar network are the main causes of the bad scaling observed at larger scales. Note that the succession of phases detailed before is also clearly visible in the time series of wet area fraction displayed in Fig. 11.
For the two other events, the conclusions are similar: discrepancies in MF parameters between simulated and observed precipitation intensities are caused primarily by tem-  Figure 11. Fraction of wet area during the event of the 13 August 2015.
poral and spatial shifts in the simulated precipitation patterns. The effect of such shifts on the MF analysis hints at the possibility of a further analysis based not on a fixed study domain but on a study domain which follows the precipitation system, in a way similar to Nykanen and Harris (2003). Using a mobile window would make the discrepancies in MF parameters depend much more on the small-scale structures of simulated precipitation intensities, since it would strongly reduce the effect of misplacement of the simulated precipitation systems.

Conclusions
In this work a spatial and temporal analysis of precipitation intensities simulated by the COSMO NWP model was performed using the UM framework which allows for the representation of the variability across scales with a limited number of parameters.
The first part of this work focused on a MF analysis of the climatology of precipitation intensities simulated by COSMO in its operational analysis mode. Analysis of the correlations between MF parameters and external meteorological and topographical descriptors revealed that the fractal dimension (D f ) and the mean intermittency (C 1 ) are strongly correlated to the fraction of rainy simulations. Additionally, the fractal dimension tends to increase and the mean intermittency tends to decrease with latitude, which indicates that rainfall fields are more homogeneous at higher latitudes. The effect of topography is visible in the values of C 1 and the maximum singularity γ s (related to extreme values) which tend to decrease with altitude, as well as in H and β which tend to increase with altitude. This indicates a smaller intermittency and less rainfall extremes in mountainous regions, as well as smoother rainfall intensity fields, which can be linked to the dominance of large-scale orographic effects. A hierarchical clustering was performed based on the meteorological and topographical descriptors. The resulting classification into three clusters was shown to correspond well with the famous Köppen (1936) climate classification. Distributions of MF parameters within these three clusters were found to be statistically significantly different, indicating that the MF signature of rain is indeed climate dependent. Finally, investigation of the spatial structure of MF fields confirmed the conclusions of the correlation analysis, namely that the values of β and H are mostly influenced by the altitude (simulated precipitation tends to be smoother at higher altitudes) and D f and C 1 are mostly influenced by the latitude (the intermittency decreases with latitude).
The second part of this work focused on three different events, one cold front associated with heavy snowfall, one stationary front associated with stratiform rain and a stable atmosphere, and one summer convection event with heavy rain. All events were simulated at a 2 km resolution with both the standard operational one-moment microphysical parameterization of COSMO and a more advanced twomoment microphysical scheme. A comparison of the precipitation intensities at the ground simulated by COSMO and the Swiss radar composite was performed in terms of their MF signature. Although the radar data show a single scaling regime over the studied spatial-scale ranges (1-128 km), the COSMO simulations display scaling breaks for the first and the last event. It can be observed that during the snowstorm event, COSMO is unable to properly reproduce radar observations at small scales, which might be caused by the intrinsic difficulty of simulating solid precipitation. During the last convective event, the opposite can be observed and COSMO is struggling to reproduce the larger scales, due to its difficulty in locating accurately the convective system in time and space during this event. In the temporal scales, a scaling break is observed both for the radar data and the COSMO simulations at around 3 h. Comparisons of the one-moment and two-moment COSMO microphysical parameterizations show that the fields simulated by the two-moment scheme tend to display a larger intermittency and variability than the one-moment scheme. This does not generally translate into a better agreement of MF parameters α and C 1 with the radar composite, except during the stratiform event where the twomoment scheme performs slightly better. However, the twomoment scheme gives a consistently better agreement with the radar QPE in terms of spatial and temporal fractal dimensions, which measure how convoluted the precipitation occurrence signal is.
Ultimately, the MF framework can be used to identify the scale ranges in which the model is able to simulate realistic fields of water contents, and as such this technique can be used as a diagnostic tool for model evaluation.
Data availability. The data used in this work are property of Me-teoSwiss and cannot be shared without their explicit authorization. The codes used in this work can be shared on request to the first author.  Atmos. Chem. Phys., 17, 14253-14273, 2017 www.atmos-chem-phys.net/17/14253/2017/ Appendix B: Visual example of the effect of MFs on the structure of a field Figure B1 illustrates the effect of varying α and C 1 on randomly generated isotropic conservative MF fields. One can see how increasing α increases the variability within nonzero intensity regions, whereas increasing C 1 decreases the intermittency and makes the field look more spatially homogeneous. Figure B2 illustrates the effect of H on isotropic MF fields, with constant α and C 1 . H can be considered as a kind of smoothness parameter that denotes the order of integration (H < 0) or differentiation (H > 0) needed to obtain the observed field from a direct MF cascade process. Figure B1. Illustration of the effect of α and C 1 on randomly generated fields. Blue pixels correspond to zero intensity, whereas pixels with non-zero intensity are shown with a grayscale colormap. Taken from Lovejoy (2017). Figure B2. Illustration of the effect of H on randomly generated fields for α = 1.2 and C 1 = 0.05. Blue pixels correspond to zero intensity, whereas pixels with non-zero intensity are shown with a grayscale colormap. Taken from Lovejoy (2017). Figure C1 shows the Köppen classification of the 209 subsquares used in Sect. 4. Indeed, the classification is quite close to the one obtained in Sect. 4.2 with topographical and meteorological descriptors. The main differences are the fact that Cluster 3 gets smaller and Cluster 2 larger, especially in the southwest, that some areas in the north of Italy are in Cluster 2, and that sea regions are absent, since the Köppen classification does consider only land areas.

Appendix C: Comparison with the Köppen classification
Distribution of MF parameters within Köppen areas (not displayed) show only minor deviations from those obtained with the meteorological classification (Fig. 5). However, some differences are visible in space for the distributions of α and C 1 , which tend to be less distinguishable between Clusters 1 and 2 than on the original meteorological classification. As a consequence, in the Köppen classification, differences in the distributions of α in space are not statistically significant according to the Kruskal-Wallis test. All other statistical tests, however, give similar outcomes than with the original classification, but often with a larger p value.