Application of Asymmetrical Statistical Distributions for 1D Simulation of Solute Transport in Streams

: Analytical solutions of the one-dimensional (1D) advection–dispersion equations, describing the substance transport in streams, are often used because of their simplicity and computational speed. Practical computations, however, clearly show the limits and the inaccuracies of this approach. These are especially visible in cases where the streams deform concentration distribution of the transported substance due to hydraulic and morphological conditions, e.g., by transient storage zones (dead zones), vegetation, and irregularities in the stream hydromorphology. In this paper, a new approach to the simulation of 1D substance transport is presented, adapted, and tested on tracer experiments available in the published research, and carried out in three small streams in Slovakia with dead zones. Evaluation of the proposed methods, based on di ﬀ erent probability distributions, conﬁrmed that they approximate the measured concentrations signiﬁcantly better than those based upon the commonly used Gaussian distribution. Finally, an example of the application of the proposed methods to an iterative (inverse) task is presented.


Introduction
Solute transport in streams and rivers is strongly related to river characteristics, such as mean flow velocity, velocity distribution, secondary currents, and turbulence features. These parameters are mainly determined by the river morphology and the discharge conditions. Most natural channels are characterized by the relevant diversity of morphological conditions. In natural channels, changing river width, curvature, bed form, bed material, and vegetation are the reason for this diversity. In rivers, which are regulated by man-made constructions such as spur dikes, groins, stabilized bed, and so on, the morphological diversity is often less pronounced, and thus, flow velocities are more homogeneous.
In natural channels, some of these morphological irregularities, such as small cavities existing in sand or gravel beds, side arms and embayments, bigger obstacles, bank vegetation, and uprooted trees, can produce recirculating flows, which occur on different scales on the riverbanks and the riverbed. These irregularities act as dead zones for the current flowing in the main stream direction. In regulated rivers, groyne fields are the most important sort of dead zones. Groyne fields can cover large parts of the river, significantly affecting its flow field. Dead-water zones or dead zones can be defined as geometrical irregularities existing at the river periphery, within which the mean flow velocity in the main stream direction is approximately equal to zero (or even negative) [1,2].

Theory
The general form of the one-dimensional advection-dispersion equation, which is valid once cross-sectional mixing has been achieved, can be found in the literature, e.g., [19] in the form where t is the time (s), C is the concentration of a substance (kg.m −3 ), D x is the dispersion coefficient in the longitudinal direction (m 2 s −1 ), v x is water velocity in the longitudinal direction (m.s −1 ), M s is a function representing the sources of solute (kg.m −3 s −1 ), x is the spatial coordinate distance (m).
It is well known that downstream of the release of a tracer into a turbulent flow, two zones are established: the advective zone and the Gaussian (or equilibrium) zone [20]. The former is defined as the distance it takes a tracer particle to sample the entire flow field, and it is characterized by a skewed concentration profile due to the velocity profile in the channel. After a sufficient time, defined by a Lagrangian timescale [21], the differential longitudinal advection and the cross sectional mixing attain a balance and the skewness decays until the concentration profile is Gaussian and the variance of the tracer could increase linearly with the time [20]. However, in some cases advection dominates dispersion and the concept of frozen cloud is used by assuming that no longitudinal dispersion occurs during the time taken for the tracer to pass a sample site [20,22].
For contaminant transport in streams with transient storage, a term expressing the exchange between the transient storage (dead zones) and the main stream is usually added in (1). Furthermore, an equation for the transport inside the dead zone is needed. If we can assume the water flow and transport parameters as constant, the resulting equations can be written in the form [9,23,24]: where C s is the concentration of the substance in the storage zone (kg.m −3 ), m is the mass exchange coefficient between the main stream and the storage zone (s −1 ), and n is the ratio between the storage zone and the main stream cross-sectional area (-). If m or n becomes zero, the (2) reduces to (1). A general form of solution for diffusion Equation (1) by [19], and eventually by [25,26] could be written as where A is the cross-sectional area of the stream (m 2 ) and f is the unknown function ("similarity solution ") [19]. The unknown function f can be determined in two different ways: 1.
Based on experiment derive a curve fit to real data.

Analytical Solution
Using the second method, we get the one-dimensional analytical solution of the advection-dispersion equation for simplified conditions and immediate solute input in the form [26] C( where v is the average flow velocity in the stream. The simplified conditions comprise a prismatic streambed, steady and uniform flow, and conservative substance, as well as assumption of lateral and vertical homogeneous concentration of the flow (1D), so theoretically, this equation can be applied after achieving cross-sectional mixing in the river. Initial and boundary conditions are typically zero concentration for the solute in stream and no lateral inflows, but using the superposition principle and dividing the stream to sections, different (nonzero) initial and boundary conditions could also be imposed. By comparing (4) with (5), the analytical solution of the function f could be interpreted as having a form of Gaussian normal distribution with parameters of normal distribution (e.g., standard deviation σ, etc.). In fact, this curve does not have the exact shape of a Gaussian normal distribution, rather it has an asymmetric form, caused by substitution of the standard deviation σ (which is constant in Gaussian distribution) by the time-dependent term √ 2D x t. Analytical solution of advection-dispersion equation (5) is currently considered as the standard analytical solution, but in our opinion the applicability of this solution is limited to flow without barriers, with the assumption of a symmetrical spreading of particles in a longitudinal direction. In real streams, however, this assumption is not always right; movement of particles is slowed down (delayed) by the accumulation of particles in dead zones of the streams.
For this reason, in the case of flow in real conditions with occurrence of dead zones in the stream, it is not appropriate to use the Gaussian normal distribution to approximate the function f from (4). Instead a different statistical distribution with asymmetric shape should be used.
After reviewing the literature, we selected four statistical distributions for the present study. Our selection was based mainly on the criterion of simplicity, easy mathematical formulation, and the number of distribution parameters. The first distribution is the "classical" Gaussian distribution (GaussD), according to (5). As mentioned earlier, this is considered the "standard solution". Inspired by the work of [18] we selected two other distribution types-the Gumbel distribution (GumbelD) and the log-normal distribution (LogNormD)-modified for three parameters. The log-normal distribution was also used in previous studies and it was stated that instream advection and dispersion processes can produce only log-normal concentration distributions [12]. Finally, the generalized extreme values distribution (GEVD) was considered. Notably, the first two distributions (GaussD, GumbelD) are two-parametric, whereas the last two (LogNormD and GEVD) are three-parametric distributions. It should also be mentioned that the Gumbel distribution is a special case of the GEV distribution.
Other potential distributions (models) were the TOSS model [27] and, finally, the generalized inverse Gaussian distribution [28]. However, these are too complicated with many embedded functions (e.g., modified Bessel functions of the second kind, etc.) and complex expressions for the distribution of specific parameters and mean, median, or mode determination. Another possible distribution was the inverse Gaussian distribution [29], but this has two parameters only, so the accuracy would presumably be similar to the GumbelD approximation method. Finally, another potential approach found in the literature is that based on signal deconvolution to obtain the residence time distribution in a stream [30].

The GumbelD Approximation Method
As an appropriate approximation of the function f from the Equation (4), the formula of the Gumbel distribution is proposed. The general equation of the Gumbel distribution is [31] where p is the distribution probability (density), µ is the location parameter, and σ is the scale parameter. The parameters from Equation (6), (considering Equation (4)) can be defined as follows: where D x,G is the dispersion coefficient in the longitudinal direction (m 2 s −1 ) used in the proposed model. To be dimensionally consistent, z is a dimensionless parameter and σ has the dimension of length (m). To get the results in concentration units (kg m −3 ) the term M/A must be included into (6). Substituting parameters from (8) and (9) into (6), the one-dimensional analytical solution has the form:

Three-Parametric LogNormD Approximation Method
We consider, as an appropriate approximation of the function f from the (4), the expression of the three-parametric log-normal distribution. The general equation of this distribution is [32] , f or x > θ where p is the distribution probability (density), µ is the location parameter (peak location), σ is the variance (scale) parameter, and θ is the shift (lower boundary) parameter. The time t max , in which the maximum (peak) concentration occurs is defined as As the peak time t max is known from field measurements, Equation (12) yields It could be assumed that the parameter t 0 (start time of nonzero concentrations) can be expressed as a proportion of the peak (maximal concentration) time, i.e., t 0 = k t max , where 0 < k < 1.

Generalized Extreme Value Distribution (GEVD) Method
Another approximation of the function f from (4) is the generalized extreme value (GEV) distribution. The general equation of the GEV distribution is [33] where p is the distribution probability (density), µ is the location parameter (peak location), σ is the variance (scale) parameter, and ξ is the shape parameter. The peak time t max is defined as In practical computations, we can almost always assume the case ξ 0, so then The variance (scale) parameter σ can be defined similar to previous cases (see (8)) Then

Literature Data Use
The literature dataset from [34] was used to evaluate all approximation methods described above. They presented measurements from 51 field tracer experiments performed between 1959 and 1973. In the present study, we excluded the tracer experiments performed on short lengths i.e., with lengths less than 10 kilometers. It can be assumed, that in such cases the results are distorted by the influence of the initial transversal and vertical mixing. We also omitted all the experiments with noncentral tracer inputs (i.e., tracer input at the bank, input at multiple points, etc.) because of the lack of data for these particular tracer experiments.
As a result, 29 tracer experiments with overall measurements in 141 profiles (there were up to 8 measuring profiles in each tracer experiment) were selected. Statistical evaluation was performed for each profile using each of the applied approximation methods (GaussD, GumbelD, LogNormD, and GEVD). This focused on finding the best fit between the measured values and the approximation with the considered method, i.e., to find the optimal set of method parameters providing the best goodness of fit based on the minimal root mean square error (RMSE) between the measured and approximated values, i.e., where RMSE is the root mean square error (kg.m −3 ), c m,t is the measured value (concentration) at the time t (kg.m −3 ), c a,t is the approximated value in the time t (kg.m −3 ), t 1 is the measurement start time (sec) and t 2 is the measurement end time (sec), NRMSE (Eq. 24) is the normalized root mean square error (-), and c max and c min are the maximal and minimal concentration in each particular experiment (kg.m −3 ). The optimal set of parameters were analyzed and set up using the built-in Excel Solver function. The target value was the RMSE and the variable parameters were the peak time (mode), the scale parameter (to replace the term M/A in corresponding equations), and the corresponding method parameters (dispersion coefficient variance, shape parameters). The target value was the main indicator for the optimization process, looking for a minimum value for this indicator to achieve the maximum conformity of the measured and modelled data for the evaluated parameters. Limits for the peak time for the optimization process were typically set up from 0.8 up to 1.2 of the real (experimental) peak time to give some variability for the optimization process. Theoretically, there should be no limits, but by setting such limits, the convergence of the optimization procedure was secured. The scale parameter limits were set up to have a maximal difference of 0.001 (1% ) difference between the sum of the tracer in modelled and experimental values. There were just formal limits for other parameters, e.g., the dispersion coefficient should a nonzero and positive value, the lower boundary coefficient k (see (17) should be positive and less than 1 (0 < k < 1). The considered rivers exhibited a large variability of hydraulic conditions. The length of the experiments was in the range from 2.5 up to 294 km (1.6-183 miles), while the range of discharges and velocity was from 0.48 up to 6824 m 3 .s −1 (17-241000 cfs-cubic feet per second) and from 0.07 up to 1.86 m.s −-1 , respectively. All the data are listed in Table 1. * means the numbering according to [34] Examples of the graphical evaluation of the considered distributions are shown in the Figures 1-3.

Field Experiments Data Use
The ability of the methods to approximate the 1D dispersion in extreme conditions of the channel with submerged vegetation in a large portion of the cross-sectional area was also tested.
For these tests, experimental data from three tracer studies in Malá Nitra, Šúrsky, and Malina streams conducted between 2012 and 2016 were used. These tracer experiments focused on the determination of longitudinal and transversal dispersion coefficients and on pollutant spreading within the primary mixing zone (prior to the transversal and vertical homogenization). In the first two streams, kitchen salt was used as tracer and the water conductivity was measured. For the experiment on the Malina stream we used a coloring agent (E133, brilliant blue) as a tracer, with subsequent spectrophotometric detection of agent concentration. In all the cases evaluated in this paper, an instantaneous input of the tracer was applied.
The discharge in all experiments was determined using the cross-section velocity method. For velocity determination the device SonTeK Flow 3D tracker was used. Water level slopes were determined using the GPS precise levelling method (with ±10 mm accuracy), and based on the water level slope, the roughness (Manning) coefficient was calculated using the Manning's formula [35].
The Malá Nitra stream is located within the village of Veľký Kýr (N48.181799°, E18.155373). The experiments were performed in two reaches with lengths of 785 and 1340 meters. The first reach of

Field Experiments Data Use
The ability of the methods to approximate the 1D dispersion in extreme conditions of the channel with submerged vegetation in a large portion of the cross-sectional area was also tested.
For these tests, experimental data from three tracer studies in Malá Nitra, Šúrsky, and Malina streams conducted between 2012 and 2016 were used. These tracer experiments focused on the determination of longitudinal and transversal dispersion coefficients and on pollutant spreading within the primary mixing zone (prior to the transversal and vertical homogenization). In the first two streams, kitchen salt was used as tracer and the water conductivity was measured. For the experiment on the Malina stream we used a coloring agent (E133, brilliant blue) as a tracer, with subsequent spectrophotometric detection of agent concentration. In all the cases evaluated in this paper, an instantaneous input of the tracer was applied. The discharge in all experiments was determined using the cross-section velocity method. For velocity determination the device SonTeK Flow 3D tracker was used. Water level slopes were determined using the GPS precise levelling method (with ±10 mm accuracy), and based on the water level slope, the roughness (Manning) coefficient was calculated using the Manning's formula [35].
The  Figure 4.
The main characteristics of all three streams used for tracer experiments in Slovakia are listed in Table 2. In all the three rivers used for field experiment data, the length was smaller than the length of rivers from the literature data. Therefore, the question of the 1D well-mixed assumption should be questioned.
Vertical mixing in a river needs, typically, a length of L mix-vert = 100*h [20], where h is the water depth. In the three Slovakian rivers, water depth h was always < 1 m (on average, 0.5, 0.6, and 0.88 m). So, L mix-vert < 100 (from 50 to 88 m). This means in the field measurement the assumption of vertical mixing holds.
The estimation of the length for transversal mixing L mix-transv is more difficult because no established theory exists to predict transverse mixing rate [36]. This rate is related to several boundary conditions of the mixing process (riverbank/riverbed irregularities, vegetation, meandering, bends, islands, confluences). The vegetation is expected to affect Manning coefficient [37], which is related to the friction factor, the shear velocity, and, ultimately, to the transverse mixing rate.  For scaling analysis, again, the length for transverse mixing is L mix-transv = c 1 *(W 2 /*h), where W is channel width and c 1 is a numerical constant related to the rate of the process of transversal mixing [38]. In the three Slovakian rivers, c 1 could be in the order of 20-40, so the L mix-transv should be in the order of 150-200 m, which is quite a lot shorter than the minimum distance of the field experiments (Table 2). Hence, the assumption of complete mixing in the cross-section should be valid.

Literature Data
The evaluation of the results from the approximation methods is based on comparison of the normalized root mean square error (NRMSE), which is the standard root mean square error (RMSE) divided by the maximal difference of measured values (c max − c min ) [39]. Since the minimal concentration in all experiments was set up to be zero, the NRMSE = RMSE/c max . For all experiments, the NRMSE was determined to evaluate the fit between the measured and approximated values. Results are listed in Table 3. If the NRMSE for the GaussD method will be regarded as 100%, values less than 100% demonstrate that the approximation method was more accurate than the GaussD method, whereas values over 100% demonstrate that the GaussD method was more accurate. Such relative results are shown in Table 4. The number of values exceeding the assigned accuracy limit are listed in Table 5. An accuracy limit, e.g., 40%, means that the NRMSE of the particular method was less than 40% of the NRMSE of the GaussD method (=100%). So, for example, for the 40% accuracy limit, the GEVD method was more precise than the GaussD method in 76 of total 141 experiments (53.9%). Accuracy limit of 100% means the same precision as the GaussD method. Further evaluation focused on the relationship (correlation) between the Gaussian (GaussD method) dispersion coefficient and the dispersion coefficient for examined approximation methods. The results are listed in Table 6.  Table 7. Similarly, the correlation between the stream conditions and the values of parameters used in each approximation method were determined (Table 8).  The ranges of dispersion parameters for each approximation method are presented in Table 9.

Field Experiments
Tables 10-13 list the results of the proposed analysis for the field experiments. The comparison was that which was already used in the previous section. The results of the approximation methods for the Malina stream, where a large portion of submerged and emerged vegetation was observed, are shown in Figures 5 and 6. It is also noted that the considered methods were able to simulate the 1D dispersion process in such conditions.  Table 13. Correlation between the dispersion coefficients. The results of the approximation methods for the Malina stream, where a large portion of submerged and emerged vegetation was observed, are shown in Figures 5 and 6. It is also noted that the considered methods were able to simulate the 1D dispersion process in such conditions.

Literature Data
As it can be seen from Table 3, the approximation accuracy of the examined approximation methods was generally better than that of the GaussD method. An average value of NRMSE for each approximation method from all experiments was determined. If we assume the average NRMSE of the GaussD method to be 100%, the relative NRMSE average of the GumbelD method was 53.1%, the LogNormD and GEVD methods were even better-40.6 and 38.2%, respectively (see Table 4). The extreme (min, max.) values of the NRMSE were caused mainly by the extreme accuracy (or inaccuracy) of the GaussD approximation method with the measured values, i.e., in cases where the NRMSE was lower than 2% or larger than 10%.
The standard deviation (the 68.3 percentile) was also significantly better for the LogNormD and GEVD methods-51.3 and, , 46.2%, respectively vs. 67.6% for the GumbelD approximation method. Table 3 shows that all the three approximation methods were better than the GaussD method, but, in turn, the LogNormD and the GEVD were clearly better than the GumbelD approximation method. The GEVD approximation method was better almost in the whole limit range (0%-100%) than the LogNormD method. In particular, the LogNormD and GEVD methods were less accurate than the GaussD method in only two cases (1.4%) and three cases (2.1%), respectively. This happened in cases, when the time-concentration curve was perfectly symmetric (presumably a prismatic streambed with no obstacles), so the GaussD method could be fully applied.
Further results were obtained in the correlation examination between the specific approximation methods. Table 6 shows a very high correlation between the dispersion coefficients of GaussD, GumbelD, and GEVD methods (r 2 > 0.99), whereas all the correlations of the LogNormD method were quite weak (r 2 was about 0.2). This means that the dispersion coefficients for the GumbelD and GEVD methods were very close to our understanding of the role of the dispersion coefficient, as identified from the standard Gauss-based analytical solution of the ADE (Equation (5)). Even the values of the dispersion coefficient between the GaussD and GEVD methods were quite similar, probably due to a linear dependency between them.
The weak correlation between the LogNormD and the GaussD method should be considered a major disadvantage of the LogNormD method, as it will create practical difficulties when using this method. This is due to the dependence of the form of the approximation curve on the parameter t0 (see (15), whereas the dispersion coefficient is very difficult to estimate due its relatively small range

Literature Data
As it can be seen from Table 3, the approximation accuracy of the examined approximation methods was generally better than that of the GaussD method. An average value of NRMSE for each approximation method from all experiments was determined. If we assume the average NRMSE of the GaussD method to be 100%, the relative NRMSE average of the GumbelD method was 53.1%, the LogNormD and GEVD methods were even better-40.6 and 38.2%, respectively (see Table 4). The extreme (min, max.) values of the NRMSE were caused mainly by the extreme accuracy (or inaccuracy) of the GaussD approximation method with the measured values, i.e., in cases where the NRMSE was lower than 2% or larger than 10%.
The standard deviation (the 68.3 percentile) was also significantly better for the LogNormD and GEVD methods-51.3 and, 46.2%, respectively vs. 67.6% for the GumbelD approximation method. Table 3 shows that all the three approximation methods were better than the GaussD method, but, in turn, the LogNormD and the GEVD were clearly better than the GumbelD approximation method. The GEVD approximation method was better almost in the whole limit range (0%-100%) than the LogNormD method. In particular, the LogNormD and GEVD methods were less accurate than the GaussD method in only two cases (1.4%) and three cases (2.1%), respectively. This happened in cases, when the time-concentration curve was perfectly symmetric (presumably a prismatic streambed with no obstacles), so the GaussD method could be fully applied.
Further results were obtained in the correlation examination between the specific approximation methods. Table 6 shows a very high correlation between the dispersion coefficients of GaussD, GumbelD, and GEVD methods (r 2 > 0.99), whereas all the correlations of the LogNormD method were quite weak (r 2 was about 0.2). This means that the dispersion coefficients for the GumbelD and GEVD methods were very close to our understanding of the role of the dispersion coefficient, as identified from the standard Gauss-based analytical solution of the ADE (Equation (5)). Even the values of the dispersion coefficient between the GaussD and GEVD methods were quite similar, probably due to a linear dependency between them.
The weak correlation between the LogNormD and the GaussD method should be considered a major disadvantage of the LogNormD method, as it will create practical difficulties when using this method. This is due to the dependence of the form of the approximation curve on the parameter t 0 (see (15), whereas the dispersion coefficient is very difficult to estimate due its relatively small range (average 0.413, std. dev. 0.1). On the other hand, the parameter t 0 can be very easily measured in tracer experiments as it is the time of the first nonzero concentration measured in the stream.
It should be noted that the values of the dispersion coefficient for the GumbelD and the GEVD methods were very similar. In fact, the Gumbel distribution, having the shape parameter ξ equal to zero (see (18)), is a special case of the GEV distribution. The results showed that the nonzero value of the shape parameter allows better fit of the approximation curve to the measured data, both for standard and distorted hydraulic conditions in the stream (flood, dead zones, etc.).
As it can be seen in Table 7, there are no noticeable or significant dependencies between the hydraulic parameters of examined streams and the particular method accuracy expressed by the NRMSE parameter. On the other hand, some correlations can be observed between the longitudinal dispersion coefficients D x , D x,G , and D x,GEV (Gauss, Gumbel, and GEV methods) on stream distance from injection point, stream velocity, and stream width (see Table 8). The parameter ξ, used in the GEV method (see (22)) has only a weak correlation to the stream hydraulic parameters, but graph analysis revealed that higher values of this parameter correspond with the long "tail" of the time-concentration curve. Hence, it is possible to assume that this parameter characterizes the extent of dead zones in the stream. Very weak (almost none) correlations were found between the dispersion parameters obtained using the LogNorm approximation method, which can be regarded as a serious disadvantage of this method (difficulty with parameters prediction).

Field Experiments
The results of the statistical test performed on field data are presented in Tables 10-13 The results were generally very similar to those of the literature data and confirmed the above conclusions.
Higher accuracy of the LogNormD and GEVD approximation methods was observed, i.e., a better fit of the modelled data to the experimental data. In addition, the average precision of all examined methods increased almost by 100%. This statement also applies to the minimal and maximal values, as well as for the standard deviation (see Tables 10 and 11).
The exceeding values test also documents the higher precision of proposed approximation methods using the field experiment data. Table 12 shows that no experiments were less precise than 60% and 50% of the GaussD method using the GumbelD, the LogNormD, and GEVD methods, respectively. The correlation results between the dispersion coefficients were almost identical to the cases using literature data (Table 13). The enhanced statistical analysis (similar as in Tables 7-9) was not possible for the field experiments carried in Slovakia, due to the low variety of hydraulic parameters.

Application of the Proposed Approximation Method to Inverse Tasks
The initial idea to explore different approximation methods for use in solute transport tasks originated within an ongoing research grant focused on localization of pollution sources on rivers (VEGA 1/0805/16, see Funding), i.e., solving the inverse task. One of the methods examined for solving such a task was the so-called "brute force method", which is based on simulation of all possible scenarios (contaminant distribution in time, distances, etc.) and comparison of the results with measured data. The best fit between the simulated and measured data set up the position of the pollution source. Such an approach requires a large number of simulations to be performed and the results are highly dependent on the simulation precision and adaptation of local conditions. Use of the GaussD method (5) led to unsatisfactory results (overall precision about 30-50%), so the idea was to use much more precise simulation methods, which are fast and able to describe the pollutant breakthrough curve more precisely, in conditions of real streams with transient (dead) zones. The solute (tracer) source localization experiment was performed on real data from the Malina stream experiment. Both numerical experiments were performed using our own simulation software LOPS (LOcalisation of Pollution Sources) with the time step of 10 minutes and space (longitudinal) step of 50 m using two different simulation methods: GaussD and the GumbelD methods. The longitudinal dispersion coefficients were determined by the results of the tracer experiment and the measured breakthrough curve was intentionally shifted in time by 3600 seconds. In the first case (GaussD method, (5)), the best fit between the measured and modelled data was achieved at a distance of 1750 m, whereas the exact distance was 1,435 m. The inaccuracy was quite large (22%) and mainly due to the inaccuracy of the model results (using GaussD method) compared with the field data.
The same conditions, but different approximation method-GumbelD (Equation (11)) were used in the second case. The fit between the modelled and measured data was much better than in the previous case-the difference sum was 0.229 in the first case, in the second case it was 0.024, i.e., approximately ten times less. Also, the solute source distance was determined almost exactly at 1450 m.

Conclusions
The application of alternative forms of the 1D analytical solution for the hydrodynamic dispersion in rivers with dead zones was explored in this study. Besides the "classical" Gauss distribution, Gumbel, log-normal, and generalized extreme value distributions were considered to evaluate their ability to approximate the dispersion process from an instantaneous injection. For this evaluation, both literature data and experimental data collected in three small streams in Slovakia were applied.
For both datasets, the results of this study show that all three approximation methods provided a higher accuracy than the Gaussian-based method. The highest accuracy was obtained using the LogNormD and the GEVD methods, having very similar precision. In particular, the LogNormD and GEVD methods were less accurate than the GaussD approximation method in only 0.7% and 2.1% of the considered cases, respectively. The GumbelD method was the least accurate when compared with the LogNormD and GEVD methods, but still more accurate than the GaussD method by around 53%.
A partial comparison of the presented methods with the dead zone model was presented in [40], where the GaussD and the GumbelD were compared with the aggregated dead zone (ADZ) model, represented by the OTIS-P model [8]. As referred in this paper, the performance of both methods (GumbelD and ADZ) was comparable, so it can be assumed that it will be similar also in the case of the LogNormD and GEVD methods. In our view, the relative ease of use of the methods used is clearly in favor of the proposed methods, for their use is sufficiently simple in a commonly available spreadsheet editor, no specialized SW is needed.
All the three examined approximation methods are suitable, especially for streams with significant influence of dead zones, which cause a strong asymmetry in the shape of the concentration distribution curves over time. These approximation methods can be shortly characterized as follows: • The Gumbel method is very easy to compute, but is less accurate because this method does not allow justification of the time concentration curve shape expressing the extent of dead zones along the examining stream reach; • The LogNormD approximation method is more precise than Gumbel and Gauss methods, but its parameters are not directly determinable and there is considerable uncertainty in their determination (there is a simple estimation of the t 0 parameter, but it has difficult predictable dispersion and location parameters (µ, D x,LN ), which have no direct physical meaning); • The GEVD approximation method is the most precise method, the dispersion coefficient (D x,GEV ) has similar meaning and way of determination as the dispersion coefficients D x and D x,G . It also has relatively predictable shape coefficient ξ (higher values correspond with longer "tails" of the time concentration curves, i.e., it expresses the extent of the dead zones). The correlation between the dispersion coefficient used in this method (D x,GEV ) and the dispersion coefficient D x , used in Equation (5) can be a subject of further research, with the aim of determination of the D x,GEV value, as well as the shape coefficient (ξ) value.
The application of the Gaussian distribution (GaussD method) in cases of streams with large a extent of dead zones is problematic and its use leads to significant differences between measured and simulated values. More detailed and precise methods for the dispersion simulation in streams with dead zones can be found in the literature [10], but these methods are quite difficult to apply as they are based on a solution of partial differential equations and require a high number of parameters.
An added advantage of the novel approach herein proposed is its simplicity and low computational effort. This allows its application for fast and simple simulations of contaminant dispersion in rivers and streams using the discrete advection-diffusion equation (ADE) technique, as described by [41], and applying the superposition principle, which is a common feature of all analytical solutions. These features have been used successfully to solve inverse tasks, where a large number of simulations are required. As documented in the discussion, such approach has significantly increased the accuracy of the inverse task solution.