Spatially-smooth regionalization of flow duration curves in non-pristine basins

Abstract

The flow duration curve (FDC) is a fundamental signature of the hydrological cycle to support water management strategies. Despite many studies on this topic, its estimation in ungauged basins is still a relevant issue as the FDC is controlled by different types of processes at different time-space scales, thus resulting quite sensitive to the specific case study.

In this work, a regional spatially-smooth procedure to evaluate the annual FDC in ungauged basins is proposed, based on the estimation of the$L$-moments (mean, $L$-CV and $L$-skewness) through regression models valid for the whole case study area. In this approach, homogeneous regions are no longer required and the $L$-moments are allowed to continuously vary along the river network, thus providing a final FDC smoothly evolving for different locations on the river. Regressions are based on a set of topographic, climatic, land use and vegetation descriptors at the basin scale. Moreover, the model ensures that the mean annual runoff is preserved at the river confluences, i.e. the sum of annual flows of the upstream reaches is equal to the predicted annual downstream flow.

The proposed model is adapted to incorporate different “sub-models” to account for local information within the regional framework, where man-induced alterations are known, as common in non-pristine catchments. In particular, we propose a module to consider the impact of existing/designed water withdrawals on the $L$-moments of the FDC.

The procedure has been applied to a dataset of daily observation of about 120 gauged basins on the upper Po river basin in North-Western Italy.

Introduction

Flow duration curves (FDC) are widely used tools to represent water availability in a river basin and are thus considered for many water resources planning and management purposes, like concessions for water uses and the planning of new hydropower plants. The FDC represents the percentage of time a certain value of discharge (usually at the daily scale) is equaled or exceeded in a river section over a specified a period of time. A FDC computed for a single year is called “annual”; if the observations of multiple years are merged together, the FDC is usually referred to as “period-of-record” or “total”. The empirical FDC can be easily built up by plotting the sorted the observations versus their frequency of non-exceedance computed with the Weibull plotting position, although the FDC is more frequently represented with respect to the exceedance frequency, thus resulting as a decreasing function. The evaluation of the FDC in ungauged basins is still a major issue in hydrological modeling, despite a large body of literature available on this topic, as recently reviewed by Castellarin et al. (2013).

In this work, we develop a regional model for prediction of the FDC in ungauged basins, developed and applied in the upper Po river basin (in North-Western Italy, an area mainly characterized by alpine and piedmont environments). The statistical framework developed is able to take advantage of local information about some types of anthropic effects. The aim of the work is to provide a regional tool to estimate the mean annual FDC in a generic watershed based on morphometric and climatic descriptors.

Methods

The regional spatially-smooth (SS) statistical estimation method proposed in this paper is based on a framework developed by Laio et al. (2011) in the context of regional flood frequency analysis. Analogies between the procedures are obtained by representing the mean annual FDC by means of its $L$-moments, i.e. distribution-free statistics which describe the mean, the variability and the skewness of the empirical FDC curve. $L$-moments are linear combinations of order statistics and are widely used in statistical hydrology; here we will refer to as $\mathit{\mu }$ for the mean value, while the dimensionless coefficients $\mathit{\tau }$ ($L$-CV) and ${\mathit{\tau }}_{\mathrm{3}}$ ($L$-skewness) respectively describe the variability and the skewness of the curve (see e.g., Hosking and Wallis, 1997, and Grimaldi et al., 2011, for more details on $L$-moments).

In the SS method the $L$-moments are the variables to be regionalized, i.e. the variables that will be predicted in ungauged basins. This is not so common in the literature as most of the models use as regionalization variables directly the FDC quantiles or the parameters of an analytical function fitted to the FDC. Then, the complete FDC can be calculated by fitting a probability distribution function to the set of predicted regional $L$-moments. In this way, the errors introduced by the preliminary choice of a distribution to fit the empirical curve do not enter in the regionalization procedure. The issue of the choice of a suitable distribution to represent the estimated FDC represents the last step in the regionalization process and will be discussed in the results section.

The present approach is referred to as spatially-smooth because it does not require the delineation of homogeneous regions, i.e. groups of basins sharing the same statistical characteristics (and thus the same $L$-moments $\mathit{\tau }$ and ${\mathit{\tau }}_{\mathrm{3}}$ within the region). Rather, $L$-moments are allowed to continuously vary in the space, so the FDC can vary accordingly. Following Wagener et al. (2013), the SS approach belong to the “mapping function” paradigms in which the set of relationship between the $L$-moments and the descriptors (and not the $L$-moments themselves) remain valid for the whole area of interest.

Ultimately, this approach can be considered an extension of the index-flow approach, as the dimensionless FDC is not constant in a region but is allowed to change site by site.

To build the regional model, the sample $L$-moments of the mean annual FDC must be computed at each available gauging station; estimation in ungauged sites is then obtained by building multiple linear regression models based on some basins descriptors (i.e. morphologic, climatic, geological, etc characteristics of the catchments), also known as descriptors. For a generic variable $y$ to be regionalized, two different model structures have been considered: $$\begin{array}{}& & y=a_{\mathrm{0}}+a_{\mathrm{1}}\cdot x_{\mathrm{1}}+\mathrm{\dots }+a_{\text{d}}\cdot x_{\text{d}}+\mathit{\epsilon }& & y=a_{\mathrm{0}}\cdot x_{\mathrm{1}}^{a_{\mathrm{1}}}\cdot \mathrm{\dots }\cdot x_{\text{d}}^{a_{\text{d}}}\cdot \mathit{\epsilon }\end{array}$$ where $x_{\text{d}}$ represents the generic basin descriptor, $a_{\text{d}}$ are the coefficients to be estimated, and $\mathit{\epsilon }$ is the residual error to be minimized. Equation (1) can be solved by the Least Squares method (LS), while Eq. (2) can be still solved within the LS framework after the log-transformation of both side of the equation. After the selection of a set of suitable regression models, the prediction of $\mathit{\mu }$, $\mathit{\tau }$ and ${\mathit{\tau }}_{\mathrm{3}}$ in a generic ungauged basin can be performed.

To provide reliable results, observed data should be natural flow observations, i.e. with no alterations due to upstream water uses. However, many gauging stations are located on rivers affected by man-induced alterations: in such cases, the observations should be previously processed to provide “naturalized” values. Where the actual amount of daily derived flow is known, corrections are straightforward, but this information is seldom available. This is the case (see Fig. 1) of run-on-river power plants where the gauging station is located between the intake and the outflow of the system and no information are available about the actual flow derived by the plant.

Example of gauging station bypassed by a water flow.

[Figure omitted. See PDF]

Observed and corrected $L$-moments for gauging stations affected by relevant water intakes.

Basin	$K$	$\mathit{\mu }$ obs.	$\mathit{\mu }$ corr.	$\mathit{\tau }$ obs.	$\mathit{\tau }$ corr.	${\mathit{\tau }}_{\mathrm{3}}$ obs.	${\mathit{\tau }}_{\mathrm{3}}$ corr.
Devero at Baceno	6.069	1.73	4.51	0.385	0.218	0.374	0.161
Germanasca at Perrero	0.637	3.14	4.6	0.708	0.528	0.627	0.465
Isorno at Pontetto	4	0.75	1.93	0.485	0.285	0.677	0.315
Po at Castiglione Torinese	2.016	54.57	117.5	0.585	0.371	0.607	0.33
S.Bernardino at S.Bernardino Santino	0.414	4.35	5.82	0.728	0.576	0.67	0.55
Stura di Demonte at Vinadio	3.053	2.62	5.45	0.482	0.291	0.536	0.258
Stura di Viù at Germagnano	0.734	4.77	7.26	0.596	0.435	0.568	0.406
Toce at Domodossola	3.391	12.21	29.34	0.375	0.224	0.469	0.226
Varaita at Torrette	8.824	1.36	3.2	0.383	0.212	0.541	0.221

To overcome this problem, and exploiting also the data relative to such stations, a new methodology, developed by Ganora et al. (2013), has been used to obtain the $L$-moments of the “natural”. In essence, the $L$-moments computed on the measurements affected by changes (i.e., downstream the intake) are converted to a set of naturalized $L$-moments subsequently used in the analysis. The method only requires the maximum discharge that can be withdrawn (said $\mathrm{\Delta }Q)$ and the affected (observed) $L$-moments and does not need to assess the complete time series of the natural flow. Note that $\mathrm{\Delta }Q$ is a known characteristic of the plants, available, in our case, from a regional database of water infrastructures.

The correction method uses the parameter: $$K={\displaystyle \frac{\mathrm{\Delta }Q}{{\mathit{\mu }}^{\text{D}}}},$$ where ${\mathit{\mu }}^{\text{D}}$ is the mean observed discharge (i.e. downstream the intake; superscript D) to correct the $L$-moments by means of the equations

$$\left\{\begin{array}{l}{\mathit{\mu }}^{\text{U}}={\mathit{\mu }}^{\text{D}}\cdot \mathrm{\Psi }\left[K\right]\\ {\mathit{\tau }}^{\text{U}}={\mathit{\tau }}^{\text{D}}\cdot {\displaystyle \frac{0.5}{L_{\mathrm{2}}^{\ast }\left[K\right]}}\\ {\mathit{\tau }}_{\mathrm{3}}^{\text{U}}={\mathit{\tau }}_{\mathrm{3}}^{\text{D}}\cdot {\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{\displaystyle \frac{L_{\mathrm{2}}^{\ast }\left[K\right]}{L_{\mathrm{3}}^{\ast }\left[K\right]}}\end{array}\right..$$ In Eq. (4) the superscript $D$ indicates the downstream (observed) value, while $U$ is the upstream (naturalized) statistic. Functions $\mathrm{\Psi }\left[K\right]$, $L_{\mathrm{2}}^{\ast }\left[K\right]$ and $L_{\mathrm{3}}^{\ast }\left[K\right]$ give the correction factor for each $L$-moment and are reported in the Appendix A. Details about the correction method are discussed in Ganora et al. (2013).

Data availability

The available hydrological dataset refers to North-West of Italy (see Fig. 2) and includes 129 gauging stations with a total of 1438 station-year of daily observations, that have been selected after different quality controls. Only years with no more than 3 missing values have been included in the dataset. The available records have a mean length of 11 years, with actual durations ranging between 2 and 52 years. Many stations are quite recent and their data do not completely overlap with longer records, however they allow a larger spatial coverage of the model. More details about the dataset can be found in Ganora et al. (2013).

For each station, the annual FDCs have been computed based on the daily discharge values; the $N=365$ values have been sorted and associated to the non-exceedance frequency $F_i=\frac{i}{N+\mathrm{1}}$ , with $i=\mathrm{1}\mathrm{\dots }N$ being the index of the sorted sample (leap years curves have been resampled on the frequency domain to obtain 365 values). Finally, the mean annual FDC at each station has been computed by taking the average of each annual FDC quantile (i.e., the average frequency-by-frequency of all the annual values).

Location of gauging stations used in the analysis.

[Figure omitted. See PDF]

The database encompasses a number of basins affected by different water withdrawals for different uses and with different intake systems. In many cases, for our purposes, such anthropic effects can be neglected (e.g., negligible intake with respect to the average flow; high-elevation reservoir affecting only the upper part of the river and with negligible effect on the long-term statistics). In other cases alterations can be corrected with measurements of real intake (e.g. large irrigation canals). However, other cases require a correction of the observed streamflow statistics to provide reliable “natural” values to be used in the regional analysis (with the method described in the previous section). This is the case of 9 gauging stations located downstream the intake of run-of-the-river hydropower plants and upstream their water return. For these intakes no withdrawal data were available, as no gauging recording system was installed. Table 1 reports the main characteristics of the above mentioned stations with man-induced alterations, as well as the corrected (“naturalized”) $L$-moments.

Spatially-Smooth regional estimation

Sample $L$-moments, i.e. $\mathit{\mu }$, $\mathit{\tau }$ and ${\mathit{\tau }}_{\mathrm{3}}$, have been computed from the available average FDCs with the necessary corrections to obtain “naturalized” values in the 9 cases described above. To support the implementation of the regional model, a number of catchment descriptors (about 120), including morphologic and climatic characteristics of the basin, soil and vegetation type, land use, etc, have been computed for each basin. The most relevant descriptors have been also reported in an “Atlas of basin characteristics” (Gallo et al., 2013). The descriptors act as independent variable in the regression models and some of them are expected to provide a statistical dependence with each $L$-moments.

Concerning the mean annual flow, to implement the regression model we considered the transformed variable $$h=(\mathit{\mu }\cdot \mathrm{31\hspace{0.17em}536})/A,$$ which represents the flow in mm when $\mathit{\mu }$ in in m${}^{\mathrm{3}}$ s${}^{-\mathrm{1}}$ and $A$ the catchment area in km${}^{\mathrm{2}}$. Using $h$ in the additive model structure (Eq. 1) and considering only descriptors that represents area-averaged values (e.g. mean basin elevation, spatial average of annual rainfall, etc.), the model guarantees the congruence of the mean discharge along the stream, i.e. the mean annual flow downstream of a confluence is equal to the sum of the two mean annual flow values computed upstream. This is a notable feature of the method that allows to keep the estimates of mean annual flow congruent at the confluences.

Natural values of the $L$-moments (observed or corrected) compared with regional predictions. Colors represent record lengths; circled points are the corrected values.

[Figure omitted. See PDF]

In order to select the most appropriate regional model, we performed all the possible regressions combining 1 to 4 descriptors as independent variables to estimate $h$; each regression was tested for significance with the $t$ Student test (5 % level of significance) and for multicollinearity with the VIF test (limit value equal to 5, see e.g. Montgomery et al., 2001). Models passing the tests were finally sorted according to their adjusted coefficient of determination and the best-performing ones were further evaluated by visual check of diagnostic error plots. Regressions have been solved using the Weighted Least Square method (e.g. Montgomery et al., 2001) considering the inverse of the square root of record length as weighting coefficients.

Empirical (points) and analytical (lines) flow duration curve for the Agogna at Momo river: (a) fitting on sample $L$-moments; (b) fitting on regional $L$-moments. Probability distributions are: log-Normal (LN3), Generalized Logistic (GLO), Generalized Pareto (GPA), Gamma (GAM) and Burr XII (B12).

[Figure omitted. See PDF]

Concerning the $L$-moments $\mathit{\tau }$ and ${\mathit{\tau }}_{\mathrm{3}}$, an analogous procedure has been performed but without requiring any congruence at the confluences, and thus without restricting the set of possible descriptors. In this case, both the additive and the multiplicative model structures have been considered, thus considering a larger number of possible models (about 10${}^{\mathrm{6}}$ models have been checked for each $L$-moment and each model structure).

The final models, applicable to the whole area of interests are: $$\begin{array}{cc}& & h=-7.3605\times {10}^{\mathrm{2}}+1.2527\times \text{MAP}+3.2569\times {10}^{-\mathrm{1}}& \times z_{\text{m}}+5.2674\times {\text{fourier}}_{\text{B1}}-6.7185\times {\text{clc}}_{\mathrm{2}}\\ & \mathit{\tau }=-2.896\times {10}^{-\mathrm{1}}-2.688\times {10}^{-\mathrm{3}}\times {\text{clc}}_{\mathrm{3}}+9.643\\ & \times {10}^{-\mathrm{5}}\times a_{75}+1.688\times {10}^{-\mathrm{4}}\times \text{MAP}+2.941\\ & & \times 10\times c_{\text{int}},\end{array}$$

$${\mathit{\tau }}_{\mathrm{3}}=4.755\times z_{max}^{-0.2702}\times \text{SD}({\text{IDF}}_a{)}^{0.06869}\times {\text{CV}}_{\text{rp}}^{0.2106},$$ where MAP is the mean annual precipitation, $z_{\text{m}}$ and $z_{max}$ are the mean and maximum basin elevation, fourier${}_{\text{B1}}$ and CV${}_{\text{rp}}$ are rainfall regime parameters, clc${}_{\mathrm{2}}$ and clc${}_{\mathrm{3}}$ are land use parameters, $a_{75}$ is the 75th percentile of the hypsographic curve, $c_{\mathrm{int}}$ and $\text{SD}\left({\text{IDF}}_a\right)$ are extreme-rainfall statistics. All the parameters are computed at the basin scale; details are reported in Gallo et al. (2013).

Performances of the selected regression models can be examined in Fig. 3 where predictions in cross-validation are reported. The mean annual discharge prediction presents little scatter, while $L$-CV and $L$-skewness tend to be more dispersed. Worse performances in higher-order $L$-moments prediction were in fact expected, but it is worth nothing that the points with larger errors are in general characterized by shorter time series and thus by larger sample uncertainty. Further diagnostics on residuals (not reported) have been performed to check possible heteroskedasticity and residual correlation.

Selection of FDC analytical form

The regional statistical framework presented is related to the estimation of the $L$-moments, similarly to Laio et al. (2011). The $L$-moments can be used to fit different probability distribution functions to represent the FDC and the choice of the analytical form of the FDC can be done as the last step of the procedure. To suggest a suitable distribution to represent the FDC in ungauged basins, different distributions have been investigated, all with 3 parameters: the log-Normal (LN3) the Generalized Logistic (GLO), the Generalized Pareto (GPA) and the Gamma (GAM) (details can be found in Hosking and Wallis, 1997), and the Burr XII (also called B12) described, for example, in Ganora and Laio (2015).

A first series of test of fitting was performed with the more common distributions (LN3, GLO, GPA, GAM). The main drawback observed was that in many cases negative streamflow values in the lower tail were obtained. This occurs when the distributions have a variable lower bound that can allow (depending on the parameter set) negative predictions for high exceedance frequencies. Further constraints to the parameters of the distribution can be added to keep the result consistent, but such controls seem hardly applicable in a regional context, for two reasons: first, the adjustments in general require numerical optimization that become difficult to be implemented in large-scale models; second, it is in general preferable to avoid complex (and likely unstable) estimation methods when the model will be used by practitioners that may not have confidence with numerical techniques.

To overcome the problem of negative quantiles without optimizing the distribution's parameters, we introduce the use of the Burr probability distribution, well known in different scientific communities but rarely used in the hydrological field (Ganora and Laio, 2015). In particular, the proposed function is the Burr XII distribution, originally proposed by Burr (1942) as a two-parameter function, and later extended in the three-parameter form (cumulative frequency function): $$\begin{array}{cc}F\left(x\right)& =\mathrm{1}-{\left[\mathrm{1},-,k,{\left({\displaystyle \frac{x}{\mathit{\lambda }}}\right)}^c\right]}^{\frac{\mathrm{1}}{k}}k\ne \mathrm{0}\\ & & =\mathrm{1}-\exp \left[-,{\left({\displaystyle \frac{x}{\mathit{\lambda }}}\right)}^c\right]k=\mathrm{0},\end{array}$$ where $k$ and $c$ are the shape parameters and $\mathit{\lambda }$ is the scale parameter. The parameter $k$ can be any real number ($-\mathrm{\infty }<k<\mathrm{\infty })$, while $c>-k$ and $\mathit{\lambda }>\mathrm{0}$. The domain is always positive, being $\mathrm{0}\le x\le \mathrm{\infty }$ for $k\le \mathrm{0}$ and $\mathrm{0}\le x\le \mathit{\lambda }{k}^{-\mathrm{1}/c}$ otherwise; for this kind of application we considered only non-positive values of $k$ in order to use only unbounded (from above) curves. The condition $k=\mathrm{0}$ is the limiting case for which $F\left(x\right)$ becomes the two-parameter Weibull distribution. Another limiting case appears for $k\to -\mathrm{\infty }$ (with $c\to \mathrm{\infty }$ simultaneously) and lead to the Pareto distribution: $$\begin{array}{cc}{F}_P\left(x\right)& =\mathrm{1}-{\left({\displaystyle \frac{x_{\mathrm{0}}}{x}}\right)}^{\mathit{\gamma }}x>x_{\mathrm{0}}\\ & & =\mathrm{0}\mathrm{otherwise},\end{array}$$ where $-c/k=\mathit{\gamma }$ and $x_{\mathrm{0}}=\mathit{\lambda }(-k{)}^{-\mathrm{1}/c}$.

The Burr function results quite “flexible”, and thus suitable to fit many different shapes of empirical data, thanks to the wide range of skewness and kurtosis values it covers; nevertheless it has only three parameters, thus avoiding overfitting problems. On the other hand, despite the simple analytical form of the Burr XII, the estimation of the two shape parameters is not straightforward as it requires the joint inversion of two nonlinear equations (not shown here). For details about the estimation of the parameters with the method of $L$-moments as well for other characteristics of the distribution the reader is referred to Ganora and Laio (2015) and references therein.

To proceed with the model selection, a comparison between the different distributions is in order. This comparison has been performed by fitting the analytical curves to 365 quantiles corresponding to the non-exceedance probabilities $F_i=i/366$ with $i=\mathrm{1}\mathrm{\dots }365$. For each distribution, the difference between the empirical and the estimated quantiles has been calculated (on both untransformed and log-transformed discharge values). The RMSE provided a error statistic for each considered station.

Results show that LN3 and GPA have generally good fitting performances; GLO and B12 have slightly larger errors, but comparable with LN3 and GPA, while the GAM does not provide reliable results. However, the LN3 provides negative values in 21 % of the stations, the GLO in 51 % and the GPA in 11 %. The average duration (in days per year) affected by negative quantiles is 10, 14 and 19 days for the LN3, GLO and GPA respectively (the average is computed considering only the stations with at least one negative value). In the context of the data set analysed here the Burr distribution results the best choice, as it provides adequate fitting capabilities without the need of further optimization to avoid negative or inconsistent quantile predictions. An example of fitting is reported in Fig. 4a where the studied analytical functions have been superimposed on the empirical FDC.

Panel b of Fig. 4 shows the same record, but compared to the analytical FDCs computed from the regional $L$-moments. An overestimation of the mid-right part of the FDC is evident, while the left side of the curve seems properly recognized. In this case, although the shift could be due to an inaccurate estimation of the $L$-moments, it can also be caused by the presence of unknown (or underestimated) water intakes upstream the gauging station. Under this perspective, the regional model can be also used as a diagnostic tool to identify basins in which an improvement of the knowledge on the actual water use is recommended. This kind of analysis is currently under development.

Conclusions

A regional model for the estimation of the mean annual flow duration curve (FDC) has been presented in the paper. The model is made of three regression-based equations to estimate the $L$-moments (i.e., the mean, the $L$-CV and the $L$-skewness) of the FDC at any ungauged basin, which can be subsequently used to fit the whole curve. In particular, concerning the mean annual flow, the river network structure is explicitly accounted for: at any confluence, the mean annual flow in the downstream section is the sum of the mean values in the two upstream river reaches, thus ensuring the congruence of the prediction.

Moreover, the method allows the use in the calibration phase of records affected by water derivations, even if the actual derived flow is not known or hardly retrievable. This is possible thanks to a parsimonious corrections method that provide “naturalized” $L$-moments considering only simple “geometrical” information about the water intake and does not require the reconstruction of the whole natural time series.

Finally, among different analytical representations of the FDC, the three-parameter Burr XII distribution has been proved to be the most suitable for the case study. In fact, while it provides good fitting performances, it does not allow the estimation of negative discharge values, which frequently appear by applying other analytical forms. This is advantageous as it is not necessary to optimize the parameters of a distribution to force the predicted quantiles to positive values only, thus making the model suitable for large-scale unsupervised applications.

Due to its flexibility, the present modeling framework is easily adaptable to further extensions like, for instance, the inclusion of observations from new gauging stations (by combining local and regional $L$-moments), while its use as a diagnostic tools to investigate the actual effect of existing water derivations is already in progress.

The correction method developed by Ganora et al. (2013) is based on the hypothesis that the upstream FDC $L$-moments can be obtained by correcting the observed ones through Eq. (4) in a parsimonious way, i.e. using just simple information about the intake system. This is necessary when scarce information are available about the plant, while more sophisticated models could be adopted where the available information is more detailed.

The method assumes the intake flow is exactly corresponding to the actual available flow below the threshold $\mathrm{\Delta }Q$. In this work we also consider the plant is always in exercise, even though the original method allows to account for a number of days of stop during the year. Correction coefficients are obtained by studying the effect of water withdrawal on a virtual FDC defined by an exponential distribution; in this way the problem can be treated analytically and solutions are based o a single parameter.

Applicability of correction factors has been tested with real data (i.e. also with non-exponential FDC) showing that $\mathit{\mu }$ and $\mathit{\tau }$ can be very efficiently reconstructed while, as expected, corrections for ${\mathit{\tau }}_{\mathrm{3}}$ provide larger errors, but still useful for the scope of the analysis.

To allow the implementation of Eq. (4), from Ganora et al. (2013), we define: $$\begin{array}{cc}& d^{\ast }=\mathrm{Exp}\left[W\right[K\left]\right]\text{where}W\left[K\right]\text{is the Lambert}W\text{-function},\\ & \text{i.e. the inverse of equation}K=W{e}^W;\\ & \mathit{\psi }\left[K\right]=K/W\left[K\right];\\ & L_{\mathrm{2}}=\mathit{\psi }\left[K\right]\cdot S_{\mathrm{2}}\left[d^{\ast }\right]+k\cdot \left(d^{\ast \mathrm{2}},-,d^{\ast }\right);\\ & L_{\mathrm{3}}=\mathit{\psi }\left[K\right]\cdot S_{\mathrm{3}}\left[d^{\ast }\right]+k\cdot (-\mathrm{2}{d}^{\ast \mathrm{3}}+\mathrm{3}{d}^{\ast \mathrm{2}}-d^{\ast });\\ & S_{\mathrm{2}}\left[d^{\ast }\right]={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{d}^{\ast }\left(\mathrm{2},-,d^{\ast },+,\mathrm{2},\left(d^{\ast },-,\mathrm{1}\right),\cdot ,\mathrm{Ln},\left[,d^{\ast },\right]\right);\\ & S_{\mathrm{3}}\left[d^{\ast }\right]={\displaystyle \frac{\mathrm{1}}{\mathrm{6}}}{d}^{\ast }\left(\mathrm{6},+,d^{\ast },\left(\mathrm{4},d^{\ast },-,\mathrm{9}\right),-,\mathrm{6},(,d^{\ast },-,\mathrm{1},),(,\mathrm{2},d^{\ast },-,\mathrm{1},)\right.\\ & \left..,\cdot ,\mathrm{Ln},\left[,d^{\ast },\right]\right).\end{array}$$

Acknowledgements

Funding from the Italian Ministry of Education, Universities and Research (FIRB-RBFR12BA3Y project), and from the European Union – Horizon 2020 ERC-2014-CoG (“CWASI” project, ref. 647473) are acknowledged.