A new description of probability density distributions of polar mesospheric clouds

In this paper we present a new description of statistical probability density functions (pdfs) of polar mesospheric clouds (PMCs). The analysis is based on observations of maximum backscatter, ice mass density, ice particle radius, and number density of ice particles measured by the ALOMAR Rayleigh–Mie–Raman lidar for all PMC seasons from 2002 to 2016. From this data set we derive a new class of pdfs that describe the statistics of PMC events that is different from previous statistical methods using the approach of an exponential distribution commonly named the g distribution. The new analysis describes successfully the probability distributions of ALOMAR lidar data. It turns out that the former g-function description is a special case of our new approach. In general the new statistical function can be applied to many kinds of different PMC parameters, e.g., maximum backscatter, integrated backscatter, ice mass density, ice water content, ice particle radius, ice particle number density, or albedo measured by satellites. As a main advantage the new method allows us to connect different observational PMC distributions of lidar and satellite data, and also to compare with distributions from ice model studies. In particular, the statistical distributions of different ice parameters can be compared with each other on the basis of a common assessment that facilitates, for example, trend analysis of PMC.

1 Introduction

First studies of probability distributions of polar mesospheric clouds (PMCs) were reported by Thomas (1995) using data from the UVS (ultraviolet spectrograph) instrument on board the Solar Mesosphere Explorer (SME) satellite and from the Solar Backscatter Ultraviolet (SBUV) instrument on the Nimbus-7 satellite over the period 1978–1986, measuring scattered limb albedo at 265 nm and nadir albedo at 273.5 nm, respectively. Thomas (1995) introduced empirical measures in the statistical analysis of PMC brightness distributions. He showed that the frequency distribution of PMC albedo derived from both SME and SBUV satellite data can be approximated by an (normalized) exponential probability function, see Fig. 3 in Thomas (1995). Secondly, the author also proposed to use cumulative frequency numbers (the so-called g function) of clouds, g(A), exceeding a certain albedo A in order to better represent the exponential populations. Examples of g distributions are plotted on a semilogarithmic scale in Fig. 4 in Thomas (1995), clearly indicating an approximately linear behavior of cumulative frequencies in a logarithmic format.

In the following years many observational PMC analyses of seasonal statistics have been published frequently using the g function, e.g., reports from Wind Imaging Interferometer (WINDII) and Polar Ozone and Aerosol Measurement II data (Shettle et al., 2002), SBUV data (Deland et al., 2003), Student Nitric Oxide Explorer (SNOE) data (Bailey et al., 2007), ice water content data derived from SBUV (DeLand and Thomas, 2015), or ALOMAR lidar data (Fiedler et al., 2017). Also model analyses have used the g function investigating trends and long-term changes in PMC parameters (Lübken et al., 2013; Berger and Lübken, 2015).

The g-function approach has been relatively successfully applied to many kinds of different PMC parameters as brightness, albedo, maximum backscatter ratio, integrated backscatter, ice water content, ice mass densities, ice particle size, or ice particle number density since frequency histograms of all these parameters have sometimes a nearly, at least piecewise, exponential shape. Furthermore, sometimes PMC data seem to fit almost perfectly to exponential distributions, particularly when using cumulative standardizations of data (Thomas, 1995). An example of a good exponential fit is the frequency distribution of ALOMAR backscatter data that are discussed in Sect. 3.1.1. On the other hand, in some statistical applications it is obvious that the exponential approach describes the data rather insufficiently, see examples of ice mass density, ice radius, and ice number density in Sect. 3.1.2. Therefore it is a desirable task to provide some more aspects on the theory of PMC statistics.

This paper makes an attempt to investigate in more detail the statistics of probability density functions (pdfs) of PMC climatology for various ice parameters. In the following we analyze a PMC data record of maximum backscatter, ice mass density, ice particle radius, and number density from the period 2002–2016 measured by the ALOMAR Rayleigh–Mie–Raman (RMR) lidar. From the analysis of these ALOMAR data, we derive a new class of pdfs of PMC distributions that, as we will show, modifies and improves the exponential (g-function) approach as introduced by Thomas (1995).

2 Description of ALOMAR lidar data

The data set obtained by the ground-based RMR lidar, located at the Arctic station ALOMAR (69^∘ N, 16^∘ E), consists of occurrence frequency, brightness, and altitude of PMC. The RMR lidar is in operation on a routine basis during the summer seasons (PMC season: 20 May to 20 August) since 1997. Since summer 2002 the lidar system has the general capability to run in a multiple wavelength (3-color) mode. We briefly summarize the 3-color lidar technique: laser pulses at three separated wavelengths (355, 532, 1064 nm) are emitted, scattered back by air molecules and ice particles in the atmosphere, and collected by telescopes. The received light is recorded by single photon counting detectors with an integration time of 15 min. After separation of the ice particle and molecular backscatter signal, we extract three vertical profiles of so-called backscatter coefficients, which are a measure of height-dependent brightness of the ice cloud. At the height of maximum backscatter (MBS) at 532 nm we calculate three MBS values. We assume that at the altitude of MBS, typically located near 83 km, the actual shape of the ice particle distribution can be described by a normal distribution. Then we derive from the three measured MBS values the characteristics of the normal distribution with mean ice radius, ice number density, and variance (Baumgarten et al., 2007). Finally, we also estimate from these ice parameters the actual ice mass density (IMD) at the MBS height. Such a Gaussian assumption has been widely used in PMC data processing of lidar and satellite data, e.g., ALOMAR lidar (Baumgarten et al., 2010) and AIM satellite with SOFIE/CIPS instruments (Hervig and Stevens, 2014; Bailey et al., 2015). Also microphysical model studies show strong evidence of Gaussian-distributed ice particles at the height of maximum brightness of PMC, e.g., Berger and von Zahn (2002) and Rapp and Thomas (2006).

In this paper we will analyze the climatology of all ice seasons from 2002 until 2016 merging all 15 seasons to one data record. Within this combined data set we then get a total number N of 8597 observations, which is sufficiently numerous in order to avoid excessive statistical irregularities in a frequency histogram of the data.

3 The exponential probability distribution (g function)

In general, the seasonal climatology of PMC events with measured ice parameters, such as integrated backscatter, maximum backscatter, column ice mass, albedo, or ice mass density, has been supposed to follow an exponential distribution that we name E(x) with ice parameter variable x. In the following we summarize the general characteristics of the exponential distribution that allows us to compute a numerical test for exponentially distributed data. The properties of the exponential probability distribution will be also compared with the characteristics of our new probability distribution approach introduced in Sect. 4.

The general form of the exponential distribution E(x) with scale parameter α>0 is defined as a pdf given by $E\left(x\right)=\mathit{\alpha }\exp (-\mathit{\alpha }x)$ that fulfills the normalization condition of a pdf with ${\int }_{\mathrm{0}}^{\mathrm{\infty }}E\left(x\right)\mathrm{d}x=\mathrm{1}$. Thomas (1995) defined the g function g(x) as the cumulative probability E_cum with

$$\begin{array}{}\text{(1)}& g\left(x\right)=E_{\mathrm{cum}}\left(x\right)=\underset{x}{\overset{\mathrm{\infty }}{\int }}\mathit{\alpha }{e}^{-\mathit{\alpha }x}\mathrm{d}{x}^{\prime }=e^{-\mathit{\alpha }x}.\end{array}$$

Taking the logarithm of E yields a straight line $\ln \left(E\right)=\ln \mathit{\alpha }-\mathit{\alpha }x$. For a given class of values [x₁;x₂] the likeliness of this class is proportional to the area enclosed by the continuous probability distribution and is obtained by integrating E on the segment length (bin size) $\mathrm{\Delta }x=x_{\mathrm{2}}-x_{\mathrm{1}}$ as ${\int }_{x_{\mathrm{1}}}^{x_{\mathrm{2}}}E\mathrm{d}x=-e^{-\mathit{\alpha }{x}_{\mathrm{2}}}+e^{-\mathit{\alpha }{x}_{\mathrm{1}}}$.

A statistical analysis of ice parameters has to take into account the aspect of specific sensitivities of different instruments. For example the ALOMAR lidar is generally sensitive to a backscatter signal larger than a threshold about 2–$\mathrm{3}\times {\mathrm{10}}^{-\mathrm{10}}$ m⁻¹ sr⁻¹ (Fiedler et al., 2017). When considering a threshold (x_th) the exponential pdf E(x) is normalized according to $A{\int }_{x_{\mathrm{th}}}^{\mathrm{\infty }}E\left(x\right)\mathrm{d}x=\mathrm{1}$ with a scaling factor A=exp (αx_th). We summarize the properties of the exponential distribution taking into account a threshold in Appendix A.

For a threshold of zero (x_th=0) we get the regular exponential distribution E(x) that has the mean $\mathit{\mu }=\mathrm{1}/\mathit{\alpha }$; median $\mathit{\nu }=\ln \left(\mathrm{2}\right)/\mathit{\alpha }$; mode η=0, which is the value that occurs most frequently in the data sample; variance ${\mathit{\sigma }}^{\mathrm{2}}=\mathrm{1}/{\mathit{\alpha }}^{\mathrm{2}}$; and standard deviation $\mathit{\sigma }=\mathrm{1}/\mathit{\alpha }$. Note that the exponential distribution has the unique property that the mean μ and standard deviation σ are identical, see also Eqs. (A1) and (A4). In combination with the median (Eq. A2), these equations form a simple statistical constraint, namely

$$\begin{array}{}\text{(2)}& \mathit{\mu }-x_{\mathrm{th}}=\mathit{\sigma }=(\mathit{\nu }-x_{\mathrm{th}})/\ln \left(\mathrm{2}\right).\end{array}$$

This allows us to test whether a given observational data sample shows good conformity with an exponential (g-function) distribution.

For a given data sample x_i ($i=\mathrm{1},\mathrm{\dots },N$) assuming a threshold x_i>x_th we use the common estimates of mean $\stackrel{\mathrm{‾}}{m}$ and variance s² (standard deviation s) with

$$\begin{array}{}\text{(3)}& \stackrel{\mathrm{‾}}{m}={\displaystyle \frac{\mathrm{1}}{N}}\sum _i^N{x}_i,s^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{N-\mathrm{1}}}\sum _i^N{\left(x_i-\stackrel{\mathrm{‾}}{m}\right)}^{\mathrm{2}},x>x_{\mathrm{th}}.\end{array}$$

In addition we also calculate the median $\stackrel{\mathrm{̃}}{m}$ and mode $\widehat{m}$. Hence testing a data sample to be exponentially distributed means that mean, median, and standard deviation of the sample have to fulfill the following identity:

$$\begin{array}{}\text{(4)}& \mathit{\mu }-x_{\mathrm{th}}=\mathit{\sigma }={\displaystyle \frac{\mathit{\nu }-x_{\mathrm{th}}}{\ln \left(\mathrm{2}\right)}}⟶\stackrel{\mathrm{‾}}{m}-x_{\mathrm{th}}=s={\displaystyle \frac{\stackrel{\mathrm{̃}}{m}-x_{\mathrm{th}}}{\ln \left(\mathrm{2}\right)}}.\end{array}$$

We will use this condition to analyze the ALOMAR data with respect to possible exponential (g-function) distributions.

3.1 Analysis of ALOMAR data on exponential distributions (g function)

3.1.1 Analysis of maximum backscatter data

We investigate the frequency distribution of MBS at 532 nm in units of 10⁻¹⁰ m⁻¹ sr⁻¹. We assume a threshold of 3 that corresponds to the instrumental sensitivity of the ALOMAR lidar. Then we sort the x = MBS data to a bin size of one per class starting from the threshold value and calculate a frequency histogram. Finally, we normalize the histogram so that the sum of all frequency classes equals one.

Figure 1 (a) Logarithm of frequency distribution of maximum backscatter (x = MBS) in units of 10⁻¹⁰ m⁻¹ sr⁻¹ (gray points x>3; black circles $\mathrm{0}<x<\mathrm{3}$). The bin size is Δx=1. The straight line (solid red) has been derived from a least-squares fit to MBS data with x>3. (b) Same as (a), but original nonlogarithmic frequency distribution (gray bars x>3; black bars $\mathrm{0}<x<\mathrm{3}$). The exponential fit derived from panel (a) is shown as a red curve. Values of mean, median, and standard deviation are given to compare fit and original data taking into account a threshold of x_th=3. The relative error given in percent describes the quality of exponential fitting, see text for details.

Figure 1a shows the frequency distribution of x = MBS data in a semilogarithm diagram. The first impression is that the data points are almost perfectly approximated by a linear regression besides some statistical noise. This indicates that an exponential function describes the distribution of data with a high accuracy. Figure 1b shows the distribution histogram in an original nonlogarithmic representation. We see that the exponential fit matches the data histogram with a high precision. The good quality of the fit is characterized by a small relative error of 6.5 % that is calculated as a sum of $\mathrm{100}\mathit{\%}\cdot {\sum }_j^M|E_j-{\mathcal{X}}_j|$ for x>3 with theoretical exponential frequencies E_j and normalized frequencies 𝒳_j of data x per class j with a total of M classes. The high quality of the fit is also supported by the fact that theoretical mean, median, mode, and standard deviation (μ, ν, η, σ) using Eqs. (A1)–(A4) and estimates of mean, median, mode, and standard deviation ($\stackrel{\mathrm{‾}}{m}$, $\stackrel{\mathrm{̃}}{m}$, $\widehat{m}$, s) from the data sample derived from Eq. (3) all coincide within their error bars. Now we perform the proposed exponential (g-function) test with $\stackrel{\mathrm{‾}}{m}-x_{\mathrm{th}}=s=(\stackrel{\mathrm{̃}}{m}-x_{\mathrm{th}})/\ln \left(\mathrm{2}\right)$, see Eq. (4), and insert the values from the data sample of mean ($\stackrel{\mathrm{‾}}{m}=\mathrm{12.0}\pm \mathrm{0.3}$), median ($\stackrel{\mathrm{̃}}{m}=\mathrm{9.0}\pm \mathrm{0.4}$), and standard deviation ($s=\mathrm{9.2}\pm \mathrm{0.5}$). The error uncertainties have been estimated with bootstrap methods. We find that $\stackrel{\mathrm{‾}}{m}-x_{\mathrm{th}}=\mathrm{12.0}-\mathrm{3}=\mathrm{9}\pm \mathrm{0.3}$, $s=\mathrm{9.2}\pm \mathrm{0.5}$, and $(\stackrel{\mathrm{̃}}{m}-x_{\mathrm{th}})/\ln \left(\mathrm{2}\right)=(\mathrm{9.0}-\mathrm{3})/\mathrm{0.69315}=\mathrm{8.7}\pm \mathrm{0.6}$. Hence the identity is fulfilled when allowing for uncertainties introduced by statistical errors. We conclude that lidar MBS data are very likely exponentially distributed and follow a g function.

3.1.2 Analysis of ice mass density, ice radius, and ice number density data

Figure 2 (a) Logarithm of frequency distribution of ice mass density (y = IMD) in units of mg m⁻³ (gray points y>20; black circles $\mathrm{0}<y<\mathrm{20}$). The bin size is Δy=2. The straight line (solid red) has been derived from a least-squares fit to IMD data with y>20. (b) Same as (a), but original nonlogarithmic frequency distribution (gray bars y>20; black bars $\mathrm{0}<y<\mathrm{20}$). The exponential fit derived from (a) is shown as a red curve. Values of mean, median, and standard deviation are given to compare fit and original data taking into account a threshold of y_th=20. The relative error given in percent describes the quality of exponential fitting. (c, d) Same, but for ice radius r in units of nm with bin size Δr=1 and threshold r_th=20. (e, f) Same, but for ice number density n in units of 1 cm⁻³ with bin size Δn=10 and threshold n_th=30.

Now we investigate other ice parameters from the ALOMAR data set with respect to exponential distributions, namely the frequency distributions of IMD in units of mg m⁻³ (threshold 20, bin size of 2), ice radius r in units of nm (threshold 20, bin size of 1), and ice number density n in units of cm⁻³ (threshold 30, bin size of 10). We will show that these parameters do not follow an exponential distribution (g function). In Fig. 2a we plot the frequency distribution for y = IMD data in a semilogarithmic diagram. We show in the following that the data points have no dominant linear shape. There exist systematic deviations from data and the theoretical exponential fit. In comparison to the fit curve, data points are systematically smaller at y=20–40. Vice versa, data points substantially exceed fit values in the range y=40–90. Also, frequencies in all classes below the threshold are significantly smaller than a proposed exponential fit. Indeed, the frequency histogram in the nonlogarithmic frame shows these systematic deviations between data and exponential fit even more pronounced, see Fig. 2b. With a relative error of about 19 % the exponential curve fails to satisfactorily fit the data. Also, significant differences exist between fit and data parameters of mean, median, mode, and standard deviation (Fig. 2b). Finally, we apply the exponential (g-function) test for IMD data and get the following results: finding the mean ($\stackrel{\mathrm{‾}}{m}=\mathrm{62.5}\pm \mathrm{1.3}$), median ($\stackrel{\mathrm{̃}}{m}=\mathrm{53.5}\pm \mathrm{1.4}$), and standard deviation ($s=\mathrm{35.2}\pm \mathrm{1.2}$) directly calculated from the data sample, we get $\stackrel{\mathrm{‾}}{m}-y_{\mathrm{th}}=\mathrm{62.5}-\mathrm{20}=\mathrm{42.5}$ not equal to s=35.2 not equal to $(\stackrel{\mathrm{̃}}{m}-y_{\mathrm{th}})/\ln \left(\mathrm{2}\right)=(\mathrm{53.5}-\mathrm{20})/\mathrm{0.69315}=\mathrm{48.3}$. Hence the condition of identity is not satisfied even allowing for uncertainties introduced by statistical errors again calculated from bootstrap methods. That is why we have to conclude that the lidar IMD data are very likely not exponentially distributed. When we investigate a possible exponential distribution for ice radius r and ice number density n data, see Fig. 2c–f, we even see larger discrepancies between data and exponential fits with, for example, relative errors of about 29 %, also indicating that both r and n are very likely not exponentially distributed. This is supported by the fact that the test of mean, median, and variance fails again and shows large inequalities.

We summarize that ice mass density, ice radius, and ice number density do not follow an exponential distribution in contrast to maximum backscatter. In the following section we will show that this is reasonable and is based on the fact that a functional link between MBS and the other data sets of IMD, r and n does miss a linear relationship.

3.2 Test of linearity between maximum backscatter and ice mass density, ice radius, ice number density data

Figure 3 (a) Maximum backscatter (x = MBS) versus ice mass density (y = IMD) in a logarithmic frame for all data with correlation coefficient R and regression parameters c and d, see text for more details. Regression points $p_{\mathrm{1}/\mathrm{2}}=[m\pm \mathrm{\Delta }m;n\pm \mathrm{\Delta }n]$ are calculated with $m=\mathrm{1}/N{\sum }_i^N\ln x_i$, $n=\mathrm{1}/N{\sum }_i^N\ln y_i$, $\mathrm{\Delta }m=\sqrt{\mathrm{1}/N{\sum }_i^N(\ln x_i-m{)}^{\mathrm{2}}}$, and $\mathrm{\Delta }n=\sqrt{\mathrm{1}/N{\sum }_i^N(\ln y_i-n{)}^{\mathrm{2}}}$. Mean and median are calculated from original nonlogarithmic data. The solid line shows the mean regression defined by regression points and corresponding c and d values. Dashed lines result from simple regression analysis of $y\left(x\right):x\Rightarrow y$ and $x\left(y\right):x\Rightarrow y$. (b) Same for x = MBS versus ice radius r.

Linearity between MBS and IMD, ice radius r and ice number density n data, is a necessary and sufficient condition that also shows that IMD, r and n data, samples are exponentially distributed, see next section. In the following we will test this constraint. Figure 3a shows a scatter plot in a logarithmic frame for simultaneously measured MBS and IMD data. In order to test a linear relationship between x = MBS and y = IMD we introduce a general fit function described by a power law condition as

$$\begin{array}{}\text{(5)}& y\left(x\right)=c{x}^d\iff x\left(y\right)=(y/c{)}^{\mathrm{1}/d}\end{array}$$

with the two constants c (linear constant) and d (power constant). Only for d=1 we expect a perfect linear dependence between x and y. First of all, the logarithmic values of x and y do not yet have a high linear correlation (R=0.54), see Fig. 3a. Since the correlation coefficient R is unequal one, the two regression lines resulting from $y\left(x\right):\ln x\mapsto \ln y$ and $x\left(y\right):\ln y\mapsto \ln x$ differ from each other. This means that the best choice of a regression fit is determined by a straight line through the two regression points, which are defined by the means plus/minus standard deviations of logarithmic x and y data. Note that the positions of regression points also relate to the half-width of the angle that is spanned by the two regression lines y(x) and x(y). For our mean regression line we estimate d=0.873. The statistical error for d is $\mathrm{\Delta }d=\pm \mathrm{0.012}$ with a confidence level of 95 %, which indicates a significant nonlinearity. Hence we conclude that the pdf describing the distribution of IMD data is very likely not an exact exponential function and its cumulative distribution does not follow precisely a g-function description because the criterion of “linearity” is violated. In Fig. 3b we show a second example for the correlation between MBS and ice radius n. Again the correlation is about R=0.55, but now the power constant is even much smaller with d=0.497, which is far away from unity. Finally, we investigated the linearity between MBS and ice number density n where we find a weak negative correlation of $R=-\mathrm{0.15}$ (not shown here). A best fit analysis yields a power value of $d=-\mathrm{0.534}$, which again fails significantly the constraint of unity. Hence we conclude that ice radius and ice number density distributions should also not follow exponential (g-function) distributions.

4 A new probability density function for PMC parameters

In this section we will present the major part of the new statistical approach in order to describe frequency distributions of different PMC parameters.

There exists a general mathematical method (“integration by substitution”) that provides the opportunity to transform between pdfs with different statistical variables. This is done by the following procedure: assuming a given pdf P(x) with variable x, then the transformation from x to a new variable y(x) with a new pdf Q(y) is specified by

$$\begin{array}{}\text{(6)}& Q\left(y\right)=∥\partial x/\partial y∥\cdot P\left(x\right(y\left)\right)\end{array}$$

with x(y) being the inverse function of y(x). Here the absolute value of the derivative $\partial x/\partial y$ has to be calculated so that the new pdf Q is defined positively everywhere. In order to apply this approach one needs generally two requirements (1) Any transformation between the two pdfs, P and Q, needs an initial guess in one of the two pdfs, either P or Q. (2) An analytic formula of a forward and backward model must be available that describes the functional dependence between the two statistical ice variables x and y. In the following we discuss how we satisfy these two requirements.

We apply this method for two ice parameters, namely MBS with variable x and an unknown ice parameter named u (e.g., this unknown ice parameter might be ice particle radius). For condition (1), we use the hypothesis that the distribution of maximum backscatter data (MBS) is perfectly represented by an exponential pdf and its cumulative distribution is described by a g function according to Eq. (1). For condition (2), we assume a power form of a fit function used in Eq. (5) that also allows us to analytically calculate the inverse function. We discuss a suitable justification of this assumption in Sect. 6.2. Hence the forward model is u(x)=cx^d and the backward model is $x\left(u\right)=(u/c{)}^{\mathrm{1}/d}$. Then the new distribution 𝒰 for the arbitrary ice parameter u using Eq. (6) is given by

$$\begin{array}{ll}\text{(7)}& {\displaystyle }\mathcal{U}\left(u\right)=& {\displaystyle }∥\partial x/\partial u/|\cdot E(x\left(u\right)){\displaystyle }=& {\displaystyle }\left.∥{\displaystyle \frac{\mathrm{1}}{\mathrm{d}u}}\cdot {\left({\displaystyle \frac{u}{c}}\right)}^{\mathrm{1}/d}\right|\cdot \mathit{\alpha }{e}^{-\mathit{\alpha }(u/c{)}^{\mathrm{1}/d}}.\end{array}$$

Equation (7) can be simplified to a more general form with

$$\begin{array}{}\text{(8)}& {\displaystyle }\mathcal{U}\left(u\right)=\stackrel{\mathrm{̃}}{a}\left|\stackrel{\mathrm{̃}}{b}\right|u^{\stackrel{\mathrm{̃}}{b}-\mathrm{1}}{e}^{-\stackrel{\mathrm{̃}}{a}{u}^{\stackrel{\mathrm{̃}}{b}}},\stackrel{\mathrm{̃}}{a}=\mathit{\alpha }{\left({\displaystyle \frac{\mathrm{1}}{c}}\right)}^{\mathrm{1}/d},\stackrel{\mathrm{̃}}{b}={\displaystyle \frac{\mathrm{1}}{d}}.\end{array}$$

In a next step we introduce, in an arbitrary manner, a third ice parameter named z for which we assume again the same power law (Eq. 5) now valid between z and u as

$${\displaystyle }z\left(u\right)=\stackrel{\mathrm{̃}}{c}{u}^{\stackrel{\mathrm{̃}}{d}}\iff u\left(z\right)=(z/\stackrel{\mathrm{̃}}{c}{)}^{\mathrm{1}/\stackrel{\mathrm{̃}}{d}}.$$

Again we apply Eq. (6) and calculate the unknown pdf 𝒵(z):

$$\begin{array}{ll}{\displaystyle }\mathcal{Z}\left(z\right)=& {\displaystyle }∥\partial u/\partial z∥\cdot \mathcal{U}\left(u\right(z\left)\right)\\ {\displaystyle }=& {\displaystyle }∥{\displaystyle \frac{\mathrm{1}}{\stackrel{\mathrm{̃}}{d}z}}\cdot {\left({\displaystyle \frac{z}{\stackrel{\mathrm{̃}}{c}}}\right)}^{\mathrm{1}/\stackrel{\mathrm{̃}}{d}}∥\cdot \stackrel{\mathrm{̃}}{a}\stackrel{\mathrm{̃}}{b}{\left((z/\stackrel{\mathrm{̃}}{c}{)}^{\mathrm{1}/\stackrel{\mathrm{̃}}{d}}\right)}^{\stackrel{\mathrm{̃}}{b}-\mathrm{1}}{e}^{-\stackrel{\mathrm{̃}}{a}(z/\stackrel{\mathrm{̃}}{c}{)}^{\stackrel{\mathrm{̃}}{b}/\stackrel{\mathrm{̃}}{d}}}\\ {\displaystyle }=& {\displaystyle }∥{\displaystyle \frac{\stackrel{\mathrm{̃}}{a}\stackrel{\mathrm{̃}}{b}}{\stackrel{\mathrm{̃}}{d}}}∥\cdot z^{-\mathrm{1}}{\left((z/\stackrel{\mathrm{̃}}{c}{)}^{\mathrm{1}/\stackrel{\mathrm{̃}}{d}}\right)}^{\stackrel{\mathrm{̃}}{b}}{e}^{-\stackrel{\mathrm{̃}}{a}(z/\stackrel{\mathrm{̃}}{c}{)}^{\stackrel{\mathrm{̃}}{b}/\stackrel{\mathrm{̃}}{d}}}.\end{array}$$

At first glance the algebraic expression for 𝒵 looks particularly complex, but 𝒵 can be transformed to a general form with $a=\stackrel{\mathrm{̃}}{a}(\mathrm{1}/\stackrel{\mathrm{̃}}{c}{)}^{\stackrel{\mathrm{̃}}{b}/\stackrel{\mathrm{̃}}{d}}$ and $b=\stackrel{\mathrm{̃}}{b}/\stackrel{\mathrm{̃}}{d}$ as

$$\begin{array}{}\text{(9)}& \mathcal{Z}\left(z\right)=a\left|b\right|z^{b-\mathrm{1}}{e}^{-a{z}^b}(a>\mathrm{0},b\ne \mathrm{0}).\end{array}$$

Equation (9) represents our final result. The pdf 𝒵(z) describes the general form of the new statistical distribution. Note that the algebraic expressions of Eqs. (8) and (9) formally coincide. This means that any probability distribution of a new ice parameter that is connected to other ice parameters through our functional power law (Eq. 5) can be described by the general pdf given by Eq. (9). The constants a and b represent two free parameters in the 𝒵 distribution, which we name the scale parameter a and the shape parameter b. Obviously, the 𝒵 pdf is identical with an exponential pdf (or g function) in the limit b=1. This shows the close interconnection of the new 𝒵 pdf to the commonly used exponential (g-function) approach. We will show in the following that any distribution from so different ice parameters, such as maximum backscatter, ice mass density, ice radius, and number density of ice particles, can be described on a uniform basis with a high accuracy by 𝒵. Vice versa this indicates that these ice parameters are connected depending on each other by the uniform power law relation (Eq. 5), more details are discussed in Sect. 6.2.

5 Application of the 𝒵 distribution to real data

5.1 General properties of the 𝒵 distribution

In this section we first show some general characteristics of the new 𝒵 distribution. From these properties we derive conditions and constraints that will allow to estimate the specific values of the two free constants in 𝒵, the scale parameter a and shape parameter b, for a given data sample.

First we show that 𝒵 is a correct pdf satisfying the normalization condition $\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\mathcal{Z}\mathrm{d}z=\mathrm{1}$:

$${\displaystyle }{\displaystyle }\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\mathcal{Z}\mathrm{d}z=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}a\left|b\right|z^{b-\mathrm{1}}{e}^{-a{z}^b}\mathrm{d}z={\left.-{\displaystyle \frac{\left|b\right|}{b}}{e}^{-a{z}^b}\right|}_{\mathrm{0}}^{\mathrm{\infty }}=\mathrm{1}.$$

The definition range of $\mathcal{Z}(z;a,b)$ is z≥0, a>0, and b≠0 with $\mathcal{Z}(z<\mathrm{0})=\mathrm{0}$ and

$$\begin{array}{ll}\text{(10a)}& {\displaystyle }& {\displaystyle }\mathcal{Z}=a\left|b\right|z^{b-\mathrm{1}}{e}^{-a{z}^b},{\displaystyle }& {\displaystyle }\ln \left(\mathcal{Z}\right)=\ln \left(a\right|b\left|\right)+(b-\mathrm{1})\cdot \ln \left(z\right)-a{z}^b,b>\mathrm{0}.\end{array}$$

For a negative b the distribution 𝒵 is described by

$$\begin{array}{ll}\text{(10b)}& {\displaystyle }& {\displaystyle }\mathcal{Z}=a\left|b\right|{\left({\displaystyle \frac{\mathrm{1}}{z}}\right)}^{|b-\mathrm{1}|}{e}^{-a(\mathrm{1}/z{)}^{\left|b\right|}},{\displaystyle }& {\displaystyle }\ln \left(\mathcal{Z}\right)=\ln \left(a\right|b\left|\right)+|b-\mathrm{1}|\cdot \ln (\mathrm{1}/z)-a(\mathrm{1}/z{)}^{\left|b\right|},b<\mathrm{0}.\end{array}$$

The cumulative form of 𝒵 for b>0 is given by

$$\begin{array}{ll}\text{(11a)}& {\displaystyle }& {\displaystyle }{\mathcal{Z}}_{\mathrm{cum}}\left(z\right)=\underset{z}{\overset{\mathrm{\infty }}{\int }}\mathcal{Z}\mathrm{d}{z}^{\prime }=e^{-a{z}^b},{\displaystyle }& {\displaystyle }\ln \left({\mathcal{Z}}_{\mathrm{cum}}\right)=-a{z}^b,\ln \left(\right|\ln \left({\mathcal{Z}}_{\mathrm{cum}}\right)\left|\right)=\ln \left(a\right)+b\ln \left(z\right).\end{array}$$

For b<0 we have to choose the cumulative calculation in reverse order starting the integration at zero. Naming the reverse cumulative with index zero as ${\mathcal{Z}}_{\mathrm{cum}}^{\mathrm{0}}$ we get

$$\begin{array}{ll}\text{(11b)}& {\displaystyle }& {\displaystyle }{\mathcal{Z}}_{\mathrm{cum}}^{\mathrm{0}}\left(z\right)=\underset{\mathrm{0}}{\overset{z}{\int }}\mathcal{Z}\mathrm{d}{z}^{\prime }=e^{-a(\mathrm{1}/z{)}^{\left|b\right|}}=e^{-a{z}^b},{\displaystyle }& {\displaystyle }\ln \left({\mathcal{Z}}_{\mathrm{cum}}^{\mathrm{0}}\right)=-a{z}^b,\ln \left(\right|\ln \left({\mathcal{Z}}_{\mathrm{cum}}^{\mathrm{0}}\right)\left|\right)=\ln \left(a\right)+b\ln \left(z\right).\end{array}$$

Only the cumulative descriptions from Eqs. (11a) and (11b) allow us, in principle, to roughly estimate the constants a and b for a given data sample using the double-logarithmic functional dependence, whereas the direct logarithm of 𝒵 (Eqs. 10a, 10b) offers no possibility to solve for a and b. However, the method calculating the double-logarithmic cumulative is not recommended. Several numerical tests showed that a stable estimation of a and b from noisy data applying this double-logarithmic approach is an almost impossible task. Instead, we propose two different methods that rely on much more powerful principles (see next Sect. 5.2). Additionally, we have to take care of a possible negative value of b that can be only identified using Eq. (11b). In fact such a case occurs in the analysis of ALOMAR data. In Sect. 5.3 we will give an example where only a negative slope parameter describes the distribution of number density of ice particles.

Generally, the 𝒵 distribution has the ability to characterize many different types of distributions, see Fig. 4. Especially the shape parameter b determines the shape of the 𝒵 distribution describing nonlinear exponential, exponential, right-skewed, left-skewed, or symmetric curves. For $\mathrm{0}<b<\mathrm{1}$ the pdf increase is nonlinear and exponentially accelerated to infinity as z approaches zero. For b=1 the pdf is exactly an exponential distribution having a positive finite value for z equal zero. For b>1, the function tends to zero as z approaches zero. When b is between 1 and 2, the function is right-skewed and rises to a peak quickly, then decreases for large z. When b has an approximate value between 3 and 4, the function becomes symmetric and bell-shaped like a normal distribution. Note that exact symmetry is given for a skewness equal to zero, which is true at z=3.60232. For b values larger than approximately 5, the function becomes again asymmetric changing the skewness to the left. For b<0, the function is skewed to the right and decreases steeply towards zero as z approaches zero. Note that 𝒵 is never negative and owns a local maximum described by the mode whenever b is negative or larger than 1. Finally we see that a double-logarithmic presentation of cumulative functions describes linear shapes with slope b, see Fig. 4j–l.

It is interesting to note that our new 𝒵 distribution is closely related to a more general Weibull distribution (Wilks, 1995). Nevertheless there is a difference concerning the shape parameter b, which in our case is not only defined for positive values but also for negative values. Such a case is disregarded by a classical 2-D Weibull distribution.

Figure 4 (a–c) Examples of 𝒵(z) function with different parameter values a and b, see Eq. (10). (d–f) Same but for ln (𝒵) from Eq. (10); (g–i) same but for 𝒵_cum from Eq. (11); (j–l) same but for $\ln \left(\right|\ln \left({\mathcal{Z}}_{\mathrm{cum}}\right)\left|\right)$ from Eq. (11).

Now we shortly summarize the mathematical descriptions of median, mode, mean, variance, and standard deviation parameters of 𝒵 for the case of a zero threshold. The calculations are described in detail in Appendix B for the general case of a nonzero threshold.
Median:

$$\begin{array}{}\text{(12a)}& \mathit{\nu }={\left({\displaystyle \frac{\ln \mathrm{2}}{a}}\right)}^{\mathrm{1}/b}.\end{array}$$

Mode:

$$\begin{array}{ll}\text{(12b)}& {\displaystyle }& {\displaystyle }\mathit{\eta }={\left({\displaystyle \frac{b-\mathrm{1}}{ab}}\right)}^{\mathrm{1}/b}\text{for}b>\mathrm{1},b<\mathrm{0};{\displaystyle }& {\displaystyle }\mathit{\eta }=\mathrm{0}\text{for}\mathrm{0}<b\le \mathrm{1}.\end{array}$$

Mean:

$$\begin{array}{}\text{(12c)}& \mathit{\mu }={\displaystyle \frac{\mathrm{\Gamma }\left(\frac{b+\mathrm{1}}{b}\right)}{a^{\frac{\mathrm{1}}{b}}}}.\end{array}$$

Variance and standard deviation:

$$\begin{array}{}\text{(12d)}& {\mathit{\sigma }}^{\mathrm{2}}={\displaystyle \frac{\mathrm{\Gamma }\left(\frac{b+\mathrm{2}}{b}\right)}{a^{\frac{\mathrm{2}}{b}}}}-{\mathit{\mu }}^{\mathrm{2}},\mathit{\sigma }=\sqrt{{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{b+\mathrm{2}}{b}\right)}{a^{\frac{\mathrm{2}}{b}}}}-{\mathit{\mu }}^{\mathrm{2}}}.\end{array}$$

The expressions of mean and variance use the gamma function $\mathrm{\Gamma }\left(t\right)=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{x}^{t-\mathrm{1}}{e}^{-x}\mathrm{d}x$. Notice that the gamma function is defined for all real values of t except t=0 and all negative integer values of t. Note also that median, mode, mean, variance, and standard deviation parameters of 𝒵 coincide with those of an exponential distribution in the limit as b equals 1.

5.2 Two computational methods to estimate the free parameters a and b of 𝒵 from a given data sample

In this section we present two numerical methods to calculate the scale parameter a and shape parameter b describing the new 𝒵 distribution. First of all, since any measurement depends on a specific instrumental sensitivity, we have to introduce a threshold that we name z_th. The remaining data sample consists of N observations z_i with z_i>z_th. Then we calculate the mean $\stackrel{\mathrm{‾}}{m}$ and standard deviation s of data z_i using Eq. (3), and also the median value $\stackrel{\mathrm{̃}}{m}$ from data z_i.

Method (1): we investigate the corresponding theoretical moments from 𝒵. In Appendix B we derive the theoretical mean μ (Eq. B5) and median ν (Eq. B3) for the 𝒵 distribution with a threshold constraint. Taking the estimates of mean $\stackrel{\mathrm{‾}}{m}$ and median $\stackrel{\mathrm{̃}}{m}$ from the sample as best proxies for the theoretical mean μ and median ν values of 𝒵, we get the following equations:

$$\begin{array}{}\text{(13a)}& \mathit{\nu }={\left({\displaystyle \frac{\ln \mathrm{2}}{a}}+z_{\mathrm{th}}^b\right)}^{\mathrm{1}/b}⟶a={\displaystyle \frac{\ln \mathrm{2}}{{\stackrel{\mathrm{̃}}{m}}^b-z_{{\mathrm{th}}^b}}},\end{array}$$

$$\begin{array}{ll}\text{(13b)}& {\displaystyle }\mathit{\mu }=& {\displaystyle }-{\displaystyle \frac{A\left|b\right|\mathrm{\Gamma }\left(\frac{b+\mathrm{1}}{b},a{z}_{\mathrm{th}}^b\right)}{b{a}^{\mathrm{1}/b}}}⟶\mathrm{0}{\displaystyle }=& {\displaystyle }\stackrel{\mathrm{‾}}{m}\cdot b{a}^{\mathrm{1}/b}+A\left|b\right|\mathrm{\Gamma }\left({\displaystyle \frac{b+\mathrm{1}}{b}},a{z}_{\mathrm{th}}^b\right).\end{array}$$

Note that the use of a threshold constraint involves the introduction of a scaling factor $A=\exp \left(a{z}_{\mathrm{th}}^b\right)$, which is present in Eq. (13b). Inserting the algebraic term of a (right side of Eq. 13a into the right side zero-equation Eq. 13b) and using the threshold value of z_th yields an equation only for b, which has to be computed iteratively. Once a numerical value of b has been estimated with a sufficient accuracy, we insert this b value into the upper right equation to get the numerical value for a.

We note that in classical statistics the method of moments determines a and b from the mean and variance equations. In principle this approach should be possible here too, but in practise the algebraic structure of the variance equation is too complicated, see Eq. (B6) in Appendix B. This means that the variance equation, if at all, is only iteratively solvable, whereas the use of the median equation offers an analytical transformation to a. Generally, we recommend to apply the proposed method using the mean and median equations. This straight-forward method is easy to program and produces reliable estimates of parameters a and b.

Method (2): we also present a second method using a maximum likelihood approach, see Appendix C. The parameters are again calculated from two equations (Eq. C5) with

$$\begin{array}{ll}\text{(14)}& {\displaystyle }& {\displaystyle \frac{\mathrm{1}}{a}}=-z_{\mathrm{th}}^b+{\displaystyle \frac{\sum z_i^b}{N}},{\displaystyle }& {\displaystyle }\mathrm{0}={\displaystyle \frac{\mathrm{1}}{b}}+a\cdot \ln z_{\mathrm{th}}\cdot z_{\mathrm{th}}^b+{\displaystyle \frac{\sum \ln z_i}{N}}-a\cdot {\displaystyle \frac{\sum \ln z_i\cdot z_i^b}{N}}.\end{array}$$

Interestingly, the left equation includes a term $\mathrm{1}/N\sum z_i^b$, which is the mean of the sample values weighted by power b, whereas the right equation includes the mean $\mathrm{1}/N\sum \ln z_i$ of logarithmic data and $\mathrm{1}/N\sum z_i^b\ln z_i$. This shows a similarity to the computation of regression points used in Fig. 3. We insert a into the right equation that yields a unique equation for b, which again can be solved iteratively. Once b is fixed, the left equation allows us to determine a. In the following we will test our lidar data samples with these two procedures and we will show that both methods produce almost identical results.

Figure 5 Frequency distributions and Z-function analysis of ALOMAR data. Parameters a and b have been estimated with the mean and median methods. The relative error given in percent describes the quality of the 𝒵-function fit; (a) maximum backscatter data x. (b) ice mass density y; (c) ice particle radius r; (d) ice number density n, see text for more details.

5.3 𝒵 distributions applied to ALOMAR data

Applications of the 𝒵 distribution to ALOMAR data of maximum backscatter (x), ice mass density (y), ice particle radius (r), and ice number density (n) are shown in Fig. 5. Note that thresholds have been computed from the regression functions (Eq. 5) described in Sect. 3.2 on the basis of $x_{\mathrm{th}}=\mathrm{3}\times {\mathrm{10}}^{-\mathrm{10}}$ m⁻¹ sr⁻¹, resulting in y_th=22 mg m⁻³, r_th=22.3 nm, and n_th=662 cm⁻³. The values of scale parameter a and shape parameter b have been calculated with the method of mean and median equations (method 1). Then the theoretical curves of 𝒵 and theoretical values of mean, median, mode, and standard deviation have been calculated by inserting the values of a, b, and threshold z_th into Eqs. (B1)–(B6). Obviously the pdf 𝒵 sometimes has no simple exponential shape, which is the case for ice mass density, ice radius, and ice number density. As we see in Fig. 5 all 𝒵-pdf curves (in blue) match the original data histograms with a high accuracy. The relative error is in a range of about 6 %–10 % except that ice number density has a relative error of 15 %. When we compare the mean, median, mode, and standard deviation derived from the theoretical distribution and corresponding estimates from data samples, we see a precise coincidence of mean and median values. Not surprising this is due to the fact that the parameters a and b have been computed by the mean and median method, which guarantees the preservation of mean and median values. Nevertheless standard deviation and mode also always show a good agreement within the error range. A closer look to the maximum backscatter distribution shows that MBS data are almost perfectly exponentially distributed with b=0.931, which is not too far away from b=1 for an exact exponential pdf. As we had already shown, see Sect. 3.1.1, MBS data are very likely exponentially distributed, now the 𝒵-distribution analysis confirms this result. Hence we conclude that the commonly used exponential (g-function) analysis might only be a reasonable statistical method in the case of analyzing MBS lidar data.

In contrast to MBS, the 𝒵 distribution of IMD shows a function that converges rapidly to zero for small IMD values. The distribution is described with b=1.355, which significantly deviates from b=1 for a precise exponential function. Note that the mode of the data sample at 40 mg m⁻³ differs from the theoretical mode of 23 mg m⁻³ because of a relatively high statistical noise in the data. But mean, median, and standard deviation values agree almost perfectly. Similar to IMD, the ice radius distribution indicates a significant nonexponential behavior with b=1.833. The distribution converges to zero as the radius approaches zero. The curve is skewed to the right and has a maximum at r=25.8 nm, which differs only slightly from the mode of the data sample at r=27.8 nm. Again mean, median, and standard deviation values agree almost perfectly.

The sample of ice number density shows a completely different behavior with a slope parameter that is negative with $b=-\mathrm{0.819}$. The physical meaning is that the parameter ice number density is negatively correlated with all other ice parameters. For example, large ice numbers n correspond to small ice radii, IMD and MBS values. As a consequence this leads to a threshold of n in the reverse direction, that is from large values to small values defined by $n<n_{\mathrm{th}}=\mathrm{662}$ cm⁻³. One can see this feature in the right tail of 𝒵(n) plotted as a dashed curve, see Fig. 5d. The reverse behavior is also present for small values of n. Small values of n are measured for very bright PMC events with large MBS that have small occurrence rates. Therefore, the number of small ice particles has a relatively high uncertainty due to their low occurrence frequency, and it is this statistical error that produces some deviations from the fit curve to the data in the range of n=0–80 cm⁻³. We note that the numerical procedure computing the pair (a, b) from the method of mean/median (Eqs. 13a, 10b) has automatically detected the existence of a negative slope parameter b without any a priori information.

Figure 6 Same as Fig. 5, but parameters a and b have been estimated from the maximum likelihood method. (a–c) Maximum backscatter, ice mass density, and ice radius: gray bars indicate values larger than the threshold and black bars indicate values smaller than the threshold. (d) Ice number density: gray bars indicate values smaller than the threshold and black bars indicate values larger than the threshold.

Now we repeat the analysis using method 2. Figure 6 summarizes the (a,b) values and statistical moments calculated from the method of maximum likelihood estimators. As can be seen the maximum likelihood approach computes almost identical results for all ice parameters. We have also added in Fig. 6a–c the histogram bars (in black) for all data being smaller than the threshold. Please keep in mind that the calculation of theoretical distribution curves is based exclusively on data larger than the threshold. Hence, decreasing or increasing a threshold will change the specific values of a and b. Figure 6d shows the ice number density distribution where we have added in the histogram (in black) all data being larger than the threshold. Again, also the maximum likelihood method has automatically detected the existence of a negative slope parameter b for the ice number density distribution.

6 Discussion

6.1 Construction of artificial data

In the derivation of the 𝒵 distribution we used the assumption that all ice parameters of maximum backscatter x = MBS, ice mass density y = IMD, ice particle radius r, and ice number density n are connected with one another by the power law given in Eq. (5). In Sect. 5.3 we showed that the 𝒵 pdf describes with a high accuracy each distribution of these ice parameters, which, in turn, means that indeed there exists at least an approximative power law between ice parameters. We discuss a suitable justification of this power law relation in more detail in Sect. 6.2. In the following we will show that the use of a 𝒵 distribution allows us to construct artificial unknown data samples of various ice parameters that approximate true data to a high degree. We think that such an application is one of the most beneficial outcomes from the new 𝒵-distribution approach. We explain the numerical procedure by the help of a practical example.

We already showed a linear dependance in the logarithmic frame using linear regression (LR) for maximum backscatter and ice particle radius, see Fig. 3b. Hence we can compute artificial ice radius proxies ${\stackrel{\mathrm{̃}}{r}}_i$, named as LR proxy of true data r_i, as a function of MBS-data x_i from the regression power law function (Eq. 5) with ${\stackrel{\mathrm{̃}}{r}}_i=c{x}_i^d$ and with power law coefficients c=13.509 and d=0.497. Figure 7a shows a comparison between LR proxy and original ice radius data where we test the identity of the two data samples. The correlation coefficient is the same as shown in Fig. 3b with R=0.55. Mean and median values of proxy and original data are almost identical, and a regression analysis shows a perfect identity (c=1.000 and d=1.000). Now we calculate the frequency histogram of LR proxy ice radii, see Fig. 7b, and compare the histogram with the original 𝒵 distribution of ice radii already shown in Fig. 5c. We find that the LR proxy approximates the mean, median, mode, and standard deviation values of the original 𝒵 distribution with an relative error of 9.5 % comparable to the original error of 9.1 %. We conclude that a linear regression analysis of logarithmic data offers a good opportunity to approximate data, provided that a pair of data samples exists that allows for the calculation of power law coefficients c and d from regression methods.

Figure 7 (a) Proxy p of ice radius versus original ice radius data without any threshold. The proxy has been derived from maximum backscatter data using the fit function that has been estimated by linear regression (LR proxy) between original logarithmic MBS and ice radius data, see Fig. 3b. (b) Frequency distribution of LR proxy (gray and black histogram) with a threshold r_th=23.3 nm. For comparison we also plot the original 𝒵-pdf curve (blue) from the analysis of original ice radius data, see Fig. 5c. The relative error describes the accuracy between LR-proxy data and original 𝒵-function fit. (c) Same as (a), but for 𝒵-proxy data resulting from the 𝒵-pdf analysis of MBS data, see text for more details; (d) same as b with 𝒵-proxy data.

Now the 𝒵-distribution approach offers a more general possibility to derive artificial data samples without any knowledge of correlation and regression coefficients. Indeed we will show that results from the 𝒵 approach are very close to results from a regression analysis. Again, our goal is to approximate ice radius data from a given maximum backscatter data sample. But now we suppose that no data of ice particle radii r exist, hence any correlation and regression analysis is not possible. First, we assume that a data sample of x = MBS of number N exists and also its 𝒵 distribution $\mathcal{Z}(x,a_x,b_x)$ with a_x=0.140 and b_x=0.931 is well known, see Fig. 5a. Secondly, we assume that we know a priori the form of the 𝒵 distribution $\mathcal{Z}(r,a_r,b_r)$ of ice radius r, e.g., with values of parameters a_r and b_r from Fig. 5c ($a_r=\mathrm{1.269}\times {\mathrm{10}}^{-\mathrm{3}}$, b_r=1.833). Please keep in mind that such information about scale and shape parameters of the ice radius distribution could be also provided from independent satellite measurements that are capable of measuring ice particle radii, e.g., AIM-SOFIE.

Our new proxy method (𝒵 proxy) requires the following transformations. We first transform the x_i values ($i=\mathrm{1},\mathrm{\dots },N$) into the z domain with $x_i\mapsto z_i:z=a_x{x}_i^{b_x}$ followed by a second transformation with $z_i\mapsto r_i:r_i=(z_i/a_r{)}^{\mathrm{1}/b_r}$ resulting in

$$\begin{array}{}\text{(15)}& {\displaystyle }{x}_i\mapsto z_i\mapsto r_i:r_i={\left[{\displaystyle \frac{a_x\cdot x_i^{b_x}}{a_r}}\right]}^{\mathrm{1}/b_r}=c\cdot x_i^d\end{array}$$

with

$${\displaystyle }c={\left({\displaystyle \frac{a_x}{a_r}}\right)}^{\mathrm{1}/b_r},d={\displaystyle \frac{b_x}{b_r}}.$$

Note that the derivation of c and d in Eq. (15) is based on the same mathematical steps when we developed the 𝒵 distribution from Eqs. (8) to (9). Inserting the a_x, b_x, a_r, and b_r values into Eq. (15) determines the power law coefficients for 𝒵 proxy $\stackrel{\mathrm{̃}}{r}$ with c=12.994 and d=0.508. These values do not exactly coincide with c and d values obtained from the regression method, see above, but the identity test between 𝒵 proxies and true ice radii shows a very good coincidence, see Fig. 7c. Again, mean and median values of proxy and original data are practically identical, and a regression analysis shows an almost perfect identity (c=1.095 and d=0.980). Finally we calculate a frequency histogram of 𝒵 proxies, see Fig. 7d, and find a good agreement between proxies and true pdf. Mean, median, and standard deviations of 𝒵-proxy data correspond perfectly to original ice radius data, and the relative error has now even decreased to 9.2 %.

We summarize that we present a new method in order to construct artificial data samples provided 𝒵-descriptions of these data sets exist. By means of a consecutive arranging of ice parameters starting at a given data sample, this method allows us to construct any artificial data sample within $(x,y,r,n)$. This method can also be applied to other data sets, e.g., ice parameter measurements from satellite observations. For example, a data sample of ice water content (IWC) obtained from satellite measurements might be analyzed in terms of a 𝒵 distribution estimating the scale and shape parameters a and b of the IWC distribution. This would allow us to establish a connection between satellite IWC data to lidar data samples $(x,y,r,n)$ through Eq. (15), hence the satellite IWC data could be transferred to lidar maximum backscatter, ice mass density, ice particle radius, and ice number density. Vice versa the knowledge of a satellite IWC 𝒵 distribution would allow us to transform lidar observations into IWC proxies and compare these with the original IWC observed by the satellite. We think that our proposed transformation method could be very helpful to connect different ice parameter data from different instruments, either from satellite observations or ground-based measurements. We also think that this new approach might be important in trend analysis of PMC.

In the next section we will discuss the power law assumption (Eq. 5) and the physical meaning of the shape parameter b, which might be introduced as a new trend variable in the analysis of PMC long-term changes.

6.2 Discussion of the power law assumption between PMC parameters

In this section we discuss some theoretical aspects of the power law dependence on ice parameters in order to validate the justification of Eq. (5). We use again the assumption as already discussed in Sect. 2 that at the altitude of maximum brightness (MBS) and ice mass density (IMD) there exists in the real atmospheric background an ice particle distribution that is perfectly Gaussian-distributed (𝒩_i) as

$${\displaystyle }{\mathcal{N}}_i\left(\mathit{\tau }\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_i\sqrt{\mathrm{2}\mathit{\pi }}}}\exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{\mathit{\tau }-r_i}{{\mathit{\sigma }}_i}}\right)}^{\mathrm{2}}\right).$$

We also assume that the geometric shapes of these ice particles are spheres with ice radii τ with mean radius r_i and variance ${\mathit{\sigma }}_i^{\mathrm{2}}$. Again index $i=\mathrm{1},\mathrm{\dots },N$ relates to the ith measurement in a given data sample of number N. 𝒩_i is normalized to ${\int }_{\mathrm{0}}^{\mathrm{\infty }}{\mathcal{N}}_i\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau }=\mathrm{1}$. When we assume an ice number density of n_i particles per cubic centimeter we get the expression ${\int }_{\mathrm{0}}^{\mathrm{\infty }}{n}_i\cdot {\mathcal{N}}_i\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau }=n_i$. Note that mean ice radii r_i and ice number densities n_i are elements of our lidar data climatology, which we have introduced in Sect. 2.

Furthermore, we assume from the analysis of lidar observations (three-color measurements) that the relation of mean radius and variance is according to σ_i=0.4 r_i for r_i<37.5 nm and σ_i=15 nm for r_i≥37.5 nm (Baumgarten et al., 2010). This assumption has also been applied in the analysis of AIM/SOFIE-CIPS PMC satellite data (Lumpe et al., 2013; Hervig and Stevens, 2014). In order to simplify calculations we apply this relation σ_i=0.4 r_i also for r_i larger than 37.5 nm. We now investigate the question of which backscatter lidar and ice mass signals result from such an ice distribution.

We compute the mass of a spherical ice particle with radius τ as 4∕3πρ_iceτ³ with density of ice ρ_ice=932 kg m⁻3. The backscatter signal from a single ice particle is calculated as a_Lτ^5.8 with the lidar constant $a_{\mathrm{L}}=\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{11}}$ m². Then the maximum backscatter x_i and ice mass density y_i are estimated by an integration of the radius distribution from zero to infinity as

$$\begin{array}{ll}{\displaystyle }{x}_i=& {\displaystyle }{a}_{\mathrm{L}}n\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\mathit{\tau }}^{\mathrm{5.8}}{\mathcal{N}}_i\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau }\\ {\displaystyle }=& {\displaystyle }{a}_{\mathrm{L}}n\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\mathit{\tau }}^{\mathrm{5.8}}{\displaystyle \frac{\mathrm{1}}{\mathrm{0.4}{r}_i\sqrt{\mathrm{2}\mathit{\pi }}}}\exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{\mathit{\tau }-r_i}{\mathrm{0.4}{r}_i}}\right)}^{\mathrm{2}}\right)\mathrm{d}\mathit{\tau }\\ {\displaystyle }{y}_i=& {\displaystyle }{\displaystyle \frac{\mathrm{4}}{\mathrm{3}}}\mathit{\pi }{\mathit{\rho }}_{\mathrm{ice}}n\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\mathit{\tau }}^{\mathrm{3}}{\mathcal{N}}_i\left(\mathit{\tau }\right)\mathrm{d}\mathit{\tau }\\ {\displaystyle }=& {\displaystyle }{\displaystyle \frac{\mathrm{4}}{\mathrm{3}}}\mathit{\pi }{\mathit{\rho }}_{\mathrm{ice}}n\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\mathit{\tau }}^{\mathrm{3}}{\displaystyle \frac{\mathrm{1}}{\mathrm{0.4}{r}_i\sqrt{\mathrm{2}\mathit{\pi }}}}\exp \left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{\mathit{\tau }-r_i}{\mathrm{0.4}{r}_i}}\right)}^{\mathrm{2}}\right)\mathrm{d}\mathit{\tau }\end{array}$$

assuming a constant number density n of ice particles. Only the integral of y_i is analytically computable with a solution in which the error function defined by the integral erf$\left(x\right)=\mathrm{2}/\sqrt{\mathit{\pi }}{\int }_{\mathrm{0}}^x\exp (-t^{\mathrm{2}})\mathrm{d}t$ is part of the solution:

$$\begin{array}{ll}{\displaystyle }{x}_i=& {\displaystyle \frac{\mathrm{4}}{\mathrm{3}}}\mathit{\pi }{\mathit{\rho }}_{\mathrm{ice}}n\left[{\displaystyle \frac{\mathrm{37}}{\mathrm{50}}}{r}_i^{\mathrm{3}}\mathrm{erf}\left({\displaystyle \frac{\mathrm{5}\left(\mathit{\tau }-r_i\right)}{\sqrt{\mathrm{8}}{r}_i}}\right)\right.\\ {\displaystyle }-& {\displaystyle \frac{\sqrt{\mathrm{2}}}{\mathrm{125}\sqrt{\mathit{\pi }}}}{r}_i\left(\mathrm{25}{\mathit{\tau }}^{\mathrm{2}}+\mathrm{25}{r}_i\mathit{\tau }+\mathrm{33}{r}_i^{\mathrm{2}}\right)\\ {\displaystyle }& {\displaystyle }{\left.\exp \left({\displaystyle \frac{\mathrm{25}\mathit{\tau }}{\mathrm{4}{r}_i}}-{\displaystyle \frac{\mathrm{25}\left({\mathit{\tau }}^{\mathrm{2}}+r_i^{\mathrm{2}}\right)}{\mathrm{8}{r}_i^{\mathrm{2}}}}\right)\right]}_{\mathrm{0}}^{\mathrm{\infty }}.\end{array}$$

The integral for x_i includes the term τ^5.8 that arises from Mie-scatter theory for light scattering of a wavelength of 532 nm (ALOMAR RMR lidar) at spheres in a range of radii with 1–100 nm. The exponential value of 5.8 approximates exact Mie-scatter calculations with a relative error less than 0.5 % in this radii range. Unfortunately, the integral can only be solved analytically if the exponent is an integer number as 5 or 6. Nevertheless, we are able to solve this integral by means of numerical methods with the specific exponent of 5.8. In a next step, we construct analytical approximations f_x(r_i) and f_y(r_i) for both integral solutions using a typical value of n=200 cm⁻³ with

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{f}_x\left(r_i\right)=x_i=a_{\mathrm{1}}{a}_{\mathrm{L}}n{r}_i^{\mathrm{5.8}},\\ {\displaystyle }& {\displaystyle }{f}_y\left(r_i\right)=y_i=a_{\mathrm{2}}{\displaystyle \frac{\mathrm{4}}{\mathrm{3}}}\mathit{\pi }{\mathit{\rho }}_{\mathrm{ice}}n{r}_i^{\mathrm{3}}.\end{array}$$

The linear constants a₁ and a₂ with values a₁=4.20 and a₂=1.47 are optimal dimensionless parameters. f_y approximates the analytical solution of y_i with a relative error less than 0.7 % in the range $r_i=[\mathrm{0},\mathrm{45}$ nm], and less than 1.2 % in the range $r_i=[\mathrm{45},\mathrm{70}$ nm]. A precise solution of x_i resulting from numerical methods of integration is approximated by f_x with a relative error less than 0.3 % in the range $r_i=[\mathrm{0},\mathrm{45}$ nm], then the relative error increases linearly to a maximum error of 5 % at r_i=70 nm. We find that the solutions f_x and f_y approximate the general power law condition p=cq^d (Eq. 5) inside a small error range. Hence these analytical examples show that MBS is a function of ice radius proportional to ∼r^5.8 (d=5.8), the same is also true for IMD (∼r³, d=3). It also follows that MBS and IMD are consequently connected through a power law condition with MBS ∼ IMD^5.8∕3 ($d=\mathrm{5.8}/\mathrm{3}=\mathrm{1.93}$).

But the new form of the 𝒵-distribution technique opens up whole new perspectives for the validation of the analytical examples based on the ALOMAR lidar data samples. We transform the z distribution of IMD into the MBS domain using Eq. (15) with $y_i\mapsto z_i:z_i=a_y{y}_i^{b_y}$ and $z_i\mapsto x_i:x_i=(z_i/a_x{)}^{\mathrm{1}/b_x}$ that gives

$$\begin{array}{}\text{(16)}& {\displaystyle }x={\left[{\displaystyle \frac{a_y\cdot y^{b_y}}{a_x}}\right]}^{\mathrm{1}/b_x}=c\cdot y^d\end{array}$$

with

$${\displaystyle }c={\left({\displaystyle \frac{a_y}{a_x}}\right)}^{\mathrm{1}/b_x},d={\displaystyle \frac{b_y}{b_x}}.$$

We insert into Eq. (16) the values of a_x=0.140, b_x=0.931, $a_y=\mathrm{4.321}\times {\mathrm{10}}^{-\mathrm{3}}$, and b_y=1.355 from Fig. 5a and b and get c=0.023 and d=1.46. The power constant (d=1.46) derived from the shape parameters b_x and b_y of the 𝒵-distribution analysis of real ALOMAR IBS and IMD data is significantly different from the power estimate (d=1.93) belonging to the analytical example that necessitates various assumptions, e.g., Gaussian-distributed ice particles at the height of maximum backscatter, constant ice particle number, or spherical shape of ice particles. Hence, we conclude that the determination of shape parameters b from a 𝒵-distribution analysis of observational data therefore provides a qualitative indication of the actual microphysical state that controls real ice formation processes. This leads to the idea that as a future task long-term changes in PMC formation might be characterized by potential long-term changes in b that indicate long-term changes in atmospheric background conditions and microphysical ice constraints of ice formation.

7 Summary and conclusions

In this study we present a new method to describe statistical probability density functions (pdfs) for different ice parameters of PMC. We analyze a climatology of ice seasons from 2002 until 2016 as measured by the ALOMAR lidar. From this data set we derive ice cloud parameters of maximum backscatter, ice mass density, ice radius, and ice number density whose occurrence frequencies are investigated with respect to exponential distributions. We show that only maximum backscatter follows an exponential distribution, whereas ice mass density, ice radius, and ice number density frequencies fail to fit satisfactorily to an exponential distribution. The reason for these deviations from exponential behavior is based on the fact that these ice parameters are not linearly dependent on each other.

We introduce a new probability density distribution (𝒵 function, see Eq. 9) that instead assumes a general power law relation among ice parameters, see Eq. (5). The new 𝒵 distribution is described by two free constants with scale parameter a and the shape parameter b. We point out that the new distribution is closely related to a more general Weibull distribution. The new distribution has been applied to maximum backscatter, ice mass density, ice radius, and ice number density data from the ALOMAR data set. As a result all data distributions are described with a high accuracy by 𝒵. We discuss that the exponential distribution (g function) is a special case of the more general 𝒵 function with shape parameter b=1. We present two numerically stable methods (method of mean and median, method of maximum likeliness) that allow to derive the values of free constants a and b describing the actual 𝒵-function shape for a given data sample.

Perhaps the most important application of the new method is the possibility to construct unknown data sets for different ice parameters that approximate true data to a high degree. We show in Sect. 6.1 that a linear regression analysis in a logarithmic data frame offers a good opportunity to approximate data provided that a pair of data samples exists that allows for the calculation of power law coefficients c and d from regression methods. The 𝒵-distribution approach offers a more general possibility to derive artificial data samples without any knowledge of correlation and regression coefficients. This allows for the connection of different observational PMC distributions of lidar and satellite data, and also with distributions resulting from ice model studies. In particular, the statistical distributions of different measured ice parameters can be compared with each other on the basis of a common assessment that again should be helpful in combining trend analysis of PMC long-term time series from different observational data sets.

Appendix A: Properties of the exponential distribution (g function)

When considering a threshold (x_th) the exponential pdf E(x) is normalized according to $A{\int }_{x_{\mathrm{th}}}^{\mathrm{\infty }}E\left(x\right)\mathrm{d}x=\mathrm{1}$ with a scaling factor A=exp (αx_th). It follows that the mean μ is then given by

$$\begin{array}{ll}{\displaystyle }\mathit{\mu }=& {\displaystyle }\underset{x_{\mathrm{th}}}{\overset{\mathrm{\infty }}{\int }}x\cdot AE\left(x\right)\mathrm{d}{x}^{\prime }\\ {\displaystyle }=& {\displaystyle }A\underset{x_{\mathrm{th}}}{\overset{\mathrm{\infty }}{\int }}x\cdot \mathit{\alpha }{e}^{-\mathit{\alpha }x}\mathrm{d}{x}^{\prime }\\ {\displaystyle }=& {\displaystyle }{\left.-A{\displaystyle \frac{\left(\mathit{\alpha }x+\mathrm{1}\right)e^{-\mathit{\alpha }x}}{\mathit{\alpha }}}\right|}_{x_{\mathrm{th}}}^{\mathrm{\infty }}\\ {\displaystyle }=& {\displaystyle }\mathrm{0}-(-)e^{\mathit{\alpha }{x}_{\mathrm{th}}}{\displaystyle \frac{\left(\mathit{\alpha }{x}_{\mathrm{th}}+\mathrm{1}\right)e^{-\mathit{\alpha }{x}_{\mathrm{th}}}}{\mathit{\alpha }}}.\end{array}$$

This yields for the mean

$$\begin{array}{}\text{(A1)}& \mathit{\mu }=x_{\mathrm{th}}+\mathrm{1}/\mathit{\alpha }.\end{array}$$

The median ν denotes the boundary of separating the higher half from the lower half of the distribution with

$$\begin{array}{ll}{\displaystyle }\mathrm{0.5}=& {\displaystyle }\underset{\mathit{\nu }}{\overset{\mathrm{\infty }}{\int }}A\cdot E\left(x\right)\mathrm{d}{x}^{\prime }=A\underset{\mathit{\nu }}{\overset{\mathrm{\infty }}{\int }}\mathit{\alpha }\exp (-\mathit{\alpha }x)\mathrm{d}{x}^{\prime }={\left.-A{e}^{-\mathit{\alpha }x}\right|}_{\mathit{\nu }}^{\mathrm{\infty }}\\ {\displaystyle }=& {\displaystyle }\mathrm{0}-(-)e^{\mathit{\alpha }{x}_{\mathrm{th}}}\cdot e^{-\mathit{\alpha }\mathit{\nu }}.\end{array}$$

The equation is solved for the median with

$$\begin{array}{}\text{(A2)}& \mathit{\nu }=x_{\mathrm{th}}+\ln \left(\mathrm{2}\right)/\mathit{\alpha }.\end{array}$$

The mode is the value η at which AE(x) takes its maximum value

$$\begin{array}{}\text{(A3)}& \mathit{\eta }=x_{\mathrm{th}}.\end{array}$$

The variance σ² in an exponential distribution considering a threshold is calculated with

$$\begin{array}{ll}{\displaystyle }{\mathit{\sigma }}^{\mathrm{2}}=& {\displaystyle }\underset{x_{\mathrm{th}}}{\overset{\mathrm{\infty }}{\int }}(x-\mathit{\mu }{)}^{\mathrm{2}}\cdot AE(x)\mathrm{d}{x}^{\prime }=A\underset{x_{\mathrm{th}}}{\overset{\mathrm{\infty }}{\int }}(x-\mathit{\mu }{)}^{\mathrm{2}}\cdot \mathit{\alpha }{e}^{-\mathit{\alpha }x}\mathrm{d}{x}^{\prime }\\ {\displaystyle }=& {\displaystyle }{\left.-{\displaystyle \frac{A\left(\mathit{\alpha }\left(x-\mathit{\mu }\right)\left(\mathit{\alpha }\left(x-\mathit{\mu }\right)+\mathrm{2}\right)+\mathrm{2}\right)e^{-\mathit{\alpha }x}}{{\mathit{\alpha }}^{\mathrm{2}}}}\right|}_{x_{\mathrm{th}}}^{\mathrm{\infty }}\\ {\displaystyle }=& {\displaystyle }\mathrm{0}-(-)e^{\mathit{\alpha }{x}_{\mathrm{th}}}{\displaystyle \frac{\left(\mathit{\alpha }\left(x_{\mathrm{th}}-\mathit{\mu }\right)\left(\mathit{\alpha }\left(x_{\mathrm{th}}-\mathit{\mu }\right)+\mathrm{2}\right)+\mathrm{2}\right)e^{-\mathit{\alpha }{x}_{\mathrm{th}}}}{{\mathit{\alpha }}^{\mathrm{2}}}}.\end{array}$$

Inserting $\mathit{\mu }=x_{\mathrm{th}}+\mathrm{1}/\mathit{\alpha }$ simplifies the algebraic expression and shows that the variance σ² (standard deviation σ) is independently from a given threshold:

$$\begin{array}{}\text{(A4)}& {\mathit{\sigma }}^{\mathrm{2}}=\mathrm{1}/{\mathit{\alpha }}^{\mathrm{2}}.\end{array}$$

Appendix B: Properties of 𝒵 distribution

In the following all quantities take into account a threshold z_th. We introduce a scaling factor $A=\exp \left(a{z}_{\mathrm{th}}^b\right)$. Setting the threshold to zero means a scaling factor A=1 and gives the regular expressions for cumulative pdf, median, mode, mean, and variance, see Eqs. (12a)–(12d).

Probability density function $\mathcal{Z}(z\ge z_{\mathrm{th}},a>\mathrm{0},b\ne \mathrm{0})$:

$$\begin{array}{ll}\text{(B1)}& {\displaystyle }& {\displaystyle }\mathcal{Z}\left(z\right)=A\cdot a\left|b\right|z^{b-\mathrm{1}}{e}^{-a{z}^b},{\displaystyle }& {\displaystyle }\mathrm{1}=\underset{z_{\mathrm{th}}}{\overset{\mathrm{\infty }}{\int }}{e}^{a{z}_{\mathrm{th}}^b}\cdot a\left|b\right|z^{b-\mathrm{1}}{e}^{-a{z}^b}\mathrm{d}z.\end{array}$$

Cumulative form of 𝒵_cum:

$$\begin{array}{ll}\text{(B2)}& {\displaystyle }{\mathcal{Z}}_{\mathrm{cum}}\left(z\right)=& {\displaystyle }A\underset{z}{\overset{\mathrm{\infty }}{\int }}a\left|b\right|z^{b-\mathrm{1}}{e}^{-a{z}^b}\mathrm{d}z{\displaystyle }=& {\displaystyle }{\left.-A{e}^{-a{z}^b}\right|}_z^{\mathrm{\infty }}=A\cdot e^{-a{z}^b},z\ge z_{\mathrm{th}}.\end{array}$$

Median ν:

$$\begin{array}{ll}{\displaystyle }\mathrm{0.5}=& {\displaystyle }\underset{\mathit{\nu }}{\overset{\mathrm{\infty }}{\int }}\mathcal{Z}\mathrm{d}z=\underset{\mathit{\nu }}{\overset{\mathrm{\infty }}{\int }}Aa\left|b\right|z^{b-\mathrm{1}}{e}^{-a{z}^b}\mathrm{d}z\\ {\displaystyle }=& {\displaystyle }{\left.-A{e}^{-a{\mathit{\nu }}^b}\right|}_{\mathit{\nu }}^{\mathrm{\infty }}=A{e}^{-a{\mathit{\nu }}^b}\end{array}$$

$$\begin{array}{ll}\text{(B3)}& {\displaystyle }⟶\mathit{\nu }=& {\displaystyle }{\left({\displaystyle \frac{\ln \left(\mathrm{2}A\right)}{a}}\right)}^{\mathrm{1}/b}{\displaystyle }=& {\displaystyle }{\left({\displaystyle \frac{\ln \mathrm{2}+\ln \left(e^{a{z}_{\mathrm{th}}^b}\right)}{a}}\right)}^{\mathrm{1}/b}\\ {\displaystyle }=& {\displaystyle }{\left({\displaystyle \frac{\ln \mathrm{2}}{a}}+z_{\mathrm{th}}^b\right)}^{\mathrm{1}/b}.\end{array}$$

Mode η:

$$\begin{array}{ll}\text{(B4)}& {\displaystyle \frac{\partial \mathcal{Z}}{\partial z}}=& {\displaystyle }\mathrm{0}=-A\cdot a\left|b\right|z^{b-\mathrm{2}}\left(ab{z}^b-b+\mathrm{1}\right)e^{-a{z}^b}{\displaystyle }⟶\mathit{\eta }=& {\displaystyle }{\left({\displaystyle \frac{b-\mathrm{1}}{ab}}\right)}^{\mathrm{1}/b}(b>\mathrm{1},b<\mathrm{0}).\end{array}$$

Mean μ:

$$\begin{array}{ll}\text{(B5)}& {\displaystyle }\mathit{\mu }=& {\displaystyle }\underset{z_{\mathrm{th}}}{\overset{\mathrm{\infty }}{\int }}z\mathcal{Z}\mathrm{d}z=A\underset{z_{\mathrm{th}}}{\overset{\mathrm{\infty }}{\int }}a\left|b\right|z^b{e}^{-a{z}^b}\mathrm{d}z{\displaystyle }=& {\displaystyle }-{\displaystyle \frac{A\left|b\right|\mathrm{\Gamma }\left(\frac{b+\mathrm{1}}{b},a{z}_{\mathrm{th}}^b\right)}{b{a}^{\frac{\mathrm{1}}{b}}}}.\end{array}$$

Details of calculation:
Substitute

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }u=a^{\frac{b+\mathrm{1}}{b}}{z}^{b+\mathrm{1}}⟶\mathrm{d}z={\displaystyle \frac{\mathrm{1}}{(b+\mathrm{1})a^{\frac{b+\mathrm{1}}{b}}{z}^b}}\mathrm{d}u\\ {\displaystyle }& {\displaystyle }\int Aa\left|b\right|z^b{e}^{-a{z}^b}\mathrm{d}z={\displaystyle \frac{Aa\left|b\right|}{b{a}^{\frac{\mathrm{1}}{b}+\mathrm{1}}+a^{\frac{\mathrm{1}}{b}+\mathrm{1}}}}\cdot \int e^{-u^{\frac{b}{b+\mathrm{1}}}}\mathrm{d}u\end{array}$$

We solve

$${\displaystyle }\int e^{-u^{\frac{b}{b+\mathrm{1}}}}\mathrm{d}u=-{\displaystyle \frac{(b+\mathrm{1})\mathrm{\Gamma }\left(\frac{b+\mathrm{1}}{b},u^{\frac{b+\mathrm{1}}{b}}\right)}{b}}$$

Inserting

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{Aa\left|b\right|}{b{a}^{\frac{\mathrm{1}}{b}+\mathrm{1}}+a^{\frac{\mathrm{1}}{b}+\mathrm{1}}}}\cdot \int e^{-u^{\frac{b}{b+\mathrm{1}}}}\mathrm{d}u\\ {\displaystyle }& {\displaystyle }=-{\displaystyle \frac{Aa\left|b\right|(b+\mathrm{1})\mathrm{\Gamma }\left(\frac{b+\mathrm{1}}{b},u^{\frac{b+\mathrm{1}}{b}}\right)}{b\left(b{a}^{\frac{\mathrm{1}}{b}+\mathrm{1}}+a^{\frac{\mathrm{1}}{b}+\mathrm{1}}\right)}}\end{array}$$

Re-substitute

$$\begin{array}{ll}{\displaystyle }=& {\displaystyle }-{\displaystyle \frac{Aa\left|b\right|(b+\mathrm{1})\mathrm{\Gamma }(\frac{b+\mathrm{1}}{b},a{z}^b)}{b\left(b{a}^{\frac{\mathrm{1}}{b}+\mathrm{1}}+a^{\frac{\mathrm{1}}{b}+\mathrm{1}}\right)}}\\ {\displaystyle }=& {\displaystyle }-{\displaystyle \frac{A\left|b\right|\mathrm{\Gamma }\left(\frac{b+\mathrm{1}}{b},a{z}^b\right)}{b{a}^{\frac{\mathrm{1}}{b}}}}.\end{array}$$

Here we use the gamma function $\mathrm{\Gamma }\left(a\right)=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{t}^{a-\mathrm{1}}{e}^{-t}\mathrm{d}t$ and the incomplete gamma function $\mathrm{\Gamma }(a,x)=\underset{x}{\overset{\mathrm{\infty }}{\int }}{t}^{a-\mathrm{1}}{e}^{-t}\mathrm{d}t$. Notice that the gamma function is defined for all real values of a except a=0 and all negative integer values of a. The same applies to Γ(a,x) with x≥0.
Variance σ²:

$$\begin{array}{ll}{\displaystyle }{\mathit{\sigma }}^{\mathrm{2}}& {\displaystyle }=\underset{z_{\mathrm{th}}}{\overset{\mathrm{\infty }}{\int }}(z-\mathit{\mu }{)}^{\mathrm{2}}\mathcal{Z}\mathrm{d}z\\ {\displaystyle }& {\displaystyle }=A\underset{z_{\mathrm{th}}}{\overset{\mathrm{\infty }}{\int }}ab(z-\mathit{\mu }{)}^{\mathrm{2}}{z}^{b-\mathrm{1}}{e}^{-a{z}^b}\mathrm{d}z\\ {\displaystyle }& {\displaystyle }={\left.-A{\mathit{\mu }}^{\mathrm{2}}{e}^{-a{z}^b}-{\displaystyle \frac{A\left[\mathrm{\Gamma }\left(\frac{b+\mathrm{2}}{b},a{z}^b\right)-\mathrm{2}{a}^{\frac{\mathrm{1}}{b}}\mathit{\mu }\mathrm{\Gamma }\left(\frac{b+\mathrm{1}}{b},a{z}^b\right)\right]}{a^{\frac{\mathrm{2}}{b}}}}\right|}_{z_{\mathrm{th}}}^{\mathrm{\infty }}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{A\left[\mathrm{\Gamma }\left(\frac{b+\mathrm{2}}{b},a{z}_{\mathrm{th}}^b\right)-\mathrm{2}{a}^{\frac{\mathrm{1}}{b}}\mathit{\mu }\mathrm{\Gamma }\left(\frac{b+\mathrm{1}}{b},a{z}_{\mathrm{th}}^b\right)\right]}{a^{\frac{\mathrm{2}}{b}}}}-{\mathit{\mu }}^{\mathrm{2}},\end{array}$$

for

$$\begin{array}{ll}\text{(B6)}& {\displaystyle }A=\mathrm{1},z_{\mathrm{th}}=& {\displaystyle }\mathrm{0}:{\mathit{\sigma }}^{\mathrm{2}}={\displaystyle \frac{\mathrm{\Gamma }\left(\frac{b+\mathrm{2}}{b}\right)}{a^{\frac{\mathrm{2}}{b}}}}-{\mathit{\mu }}^{\mathrm{2}}{\displaystyle }=& {\displaystyle \frac{\mathrm{\Gamma }\left(\frac{b+\mathrm{2}}{b}\right)-{\mathrm{\Gamma }}^{\mathrm{2}}\left(\frac{b+\mathrm{1}}{b}\right)}{a^{\frac{\mathrm{2}}{b}}}}\\ {\displaystyle }=& {\displaystyle \frac{\mathrm{2}\mathrm{\Gamma }\left(\frac{\mathrm{2}}{b}\right)-\frac{\mathrm{1}}{b}{\mathrm{\Gamma }}^{\mathrm{2}}\left(\frac{\mathrm{1}}{b}\right)}{b{a}^{\frac{\mathrm{2}}{b}}}}.\end{array}$$

Appendix C: Estimation of parameters (a,b) using the maximum log-likelihood method

For a single observation, the likelihood function l of Z is calculated from Eq. (B1). Given a sample of N observations with threshold z_th, the likelihood function $l(a,b)={\prod }_{i=\mathrm{1}}^N\mathcal{Z}\left(z_i\right)$ is

$$\begin{array}{}\text{(C1)}& l(a,b)=e^{Na{z}_{\mathrm{th}}^b}{\left(a\left|b\right|\right)}^N\prod _{i=\mathrm{1}}^N{z}_i^{b-\mathrm{1}}{e}^{-a{z}_i^b}.\end{array}$$

Taking the logarithm of l yields the log-likelihood function $L(a,b)=\ln \left(l\right)={\sum }_{i=\mathrm{1}}^N\ln \left(\mathcal{Z}\right(z_i\left)\right)$

$$\begin{array}{ll}\text{(C2)}& {\displaystyle }& {\displaystyle }L(a,b)={\displaystyle }& {\displaystyle }Na{z}_{\mathrm{th}}^b+N(\ln a+\ln |b\left|\right)+(b-\mathrm{1})\sum _{i=\mathrm{1}}^N\ln \left(z_i\right)-a\sum _{i=\mathrm{1}}^N{z}_i^b.\end{array}$$

The derivative with respect to parameter a is

$$\begin{array}{}\text{(C3)}& {\displaystyle \frac{\partial L(a,b)}{\partial a}}=N\cdot z_{\mathrm{th}}^b+{\displaystyle \frac{N}{a}}-\sum _{i=\mathrm{1}}^N{z}_i^b,\end{array}$$

and for parameter b

$$\begin{array}{ll}\text{(C4)}& {\displaystyle }& {\displaystyle \frac{\partial L(a,b)}{\partial b}}={\displaystyle }& {\displaystyle }N\cdot a\cdot \ln z_{\mathrm{th}}\cdot z_{\mathrm{th}}^b+{\displaystyle \frac{N}{b}}+\sum _{i=\mathrm{1}}^N\ln z_i-a\sum _{i=\mathrm{1}}^N{z}_i^b\ln z_i.\end{array}$$

Setting each of the derivatives equal to zero yields for a and b.

$$\begin{array}{ll}\text{(C5)}& {\displaystyle }& {\displaystyle \frac{\mathrm{1}}{a}}=-z_{\mathrm{th}}^b+{\displaystyle \frac{\sum z_i^b}{N}},{\displaystyle }& {\displaystyle }\mathrm{0}={\displaystyle \frac{\mathrm{1}}{b}}+a\cdot \ln z_{\mathrm{th}}\cdot z_{\mathrm{th}}^b+{\displaystyle \frac{\sum \ln z_i}{N}}-a\cdot {\displaystyle \frac{\sum \ln z_i\cdot z_i^b}{N}}.\end{array}$$

These are the maximum-likelihood estimators for scale parameter a and shape parameter b.