Evaluating longitudinal dispersion of scalars in rural channels of agro-urban environments

In agro-urban environments, the water resource conveyed by rural channels is susceptible to a gradual impoverishment due to the continuous combined sewer overflow release, constituting a pending and urgent issue for water management companies and the entire community. Reliable one-dimensional longitudinal dispersion coefficients D are required to model and study the hydrodynamics and water quality patterns at the scale of rural channel networks. Empirical formulas are usually adopted to estimate D but the accuracy in the prediction could be questionable. In order to identify which are the most suitable formulas to determine D in rural channels, field tracer measurements were carried out in three rural channels with typical geometry and configuration. The obtained D values were then compared with the most commonly used predicting formulas that the literature provides. The accuracy of the predictors was further checked by simulating different flow rates inside the tested channels by using a one-dimensional hydraulic model. Starting from the obtained results, indications and guidelines to choose the most suitable formulas to predict D in rural channels were provided. These indications should be followed when developing realistic quality models in the agro-urban environments, especially in those cases where direct measurements of the longitudinal dispersion coefficient D are not available. Longitudinal dispersion coefficient measurements by means of field tracer experiments in different rural channels. Field data prediction by using different formulas present in the literature. Definition of guidelines on the most reliable formulas to be used in the rural channels context. Longitudinal dispersion coefficient measurements by means of field tracer experiments in different rural channels. Field data prediction by using different formulas present in the literature. Definition of guidelines on the most reliable formulas to be used in the rural channels context.


Introduction
The correct control of agro-urban environments is becoming a leading-edge topic for water management companies. Agro-urban environments start from the edge of conurbations and extend in the surrounding countryside, including agricultural and livestock farm fields. The great complexity in the management of agro-urban environments is related to the dense interconnection between the rural channel network and the urban drainage systems of the therein towns. The interconnection occurs through spillways that, during severe rainfall events, deliver part of the sewer discharge directly into the rural network without any control or treatment, whenever the hydraulic capacity of the urban drainage system crosses a designed threshold. In this way, the rural network acts as the receiving water body of the overflows that, most of the time, comes from combined sewage systems where raw wastewater is mixed with surface stormwater runoff.
This picture is common worldwide [8,24,67,97] and, since the rural channels are normally used for farming, the frequent releases of combined sewer overflows (CSOs) could cause sanitary and environmental issues related to crops irrigation and aquifer recharge. As reported in many studies [e.g. 36,49,77], the CSO pollutant loads are comparable, in magnitude, to the emission of wastewater treatment plants (considering an annual basis). Thus, since CSO events have low occurrence frequency, this means that they release large volumes of water with a remarkable concentration of pollutants in terms of carbon, nutrients, micropollutants, heavy metals and pathogens [49,71,77].
The problem of CSOs is further exacerbated by the concomitance of two key factors: (i) the availability of space to realise appropriate treatments in high-density urban areas is constantly decreasing [45]; (ii) the mean and maximum daily precipitations are increasing because of climate change [87], activating the spillways more frequently. Moreover, the contamination risk and water quality impoverishment associated to CSO spills are more pronounced in rural channels than in rivers since the overflows represent a considerable fraction of the total flow rate.
In this challenging context, the problem of contaminants spreading in the rural channel network becomes pivotal in order to study and implement innovative solutions for the control, management and treatment of CSO water [55,72]. The correct design of solutions that exploit the intrinsic self-purification proprieties of rural channels relies on the knowledge of pollutants concentration, propagation and dispersion along the network. Besides, it is required the development of water quality models that asses the effectivness of the choosen solution.
Due to the computational effort of three-and bi-dimensional water quality models at the network scale, one-dimensional (1-D) models are commonly preferred and employed [35,48]. The application of the one-dimensional advection-dispersion equation [84,85] requires the knowledge of the longitudinal dispersion coefficient D that governs the pollutant transport. This coefficient can be directly measured or predicted by means of appropriate empirical formulas. In the first case, the measured D exactly reflects the watercourse's longitudinal dispersion conditions, although only an instantaneous picture in terms of flow rate and water stage is captured. In the second case, D can be estimated for different hydraulic conditions but it is important to be careful in choosing the formula that better reproduces the behaviour of the watercourse under examination. Empirical formulas are developed to consider the major factors influencing the longitudinal dispersion phenomenon in streamflows such as channel morphology, flow resistance and hydrodynamic regimes [34,76,78,96]. Many researchers have proposed their formula (see "Appendix"), 1 3 but the accuracy of the provided longitudinal dispersion coefficient is often questionable when it is applied in circumstances that are different from those represented by the training dataset used to determine the empirical predictor.
The difficulties associated with the determination of D represent one of the highest sources of uncertainty for hydrodynamic models simulating pollutant transport in streams. Hence, the recommended practice is to carry out field measurements to directly quantify D and then use the gained information to select the most reliable formula among the existing ones or to build a new site-specific formula [4,48,73,92].
From an experimental point of view, the longitudinal dispersion coefficient D can be obtained from tracer experiments or through velocity field measurements (together with the knowledge of the river/channel bathymetry). In the first case, D is estimated from breakthrough curves with different methods [34,40], while, in the second case, D is determined directly through its analytical definition (initially derived by Taylor [85] for pipe flows and then adapted for open-channel flows by Fischer [29]). As reported in [10,48], longitudinal dispersion coefficients obtained with these two experimental methodologies are substantially in agreement. However, velocities and bathymetry measurements required advanced instrumentation (e.g. acoustic Doppler current profilers), making this experimental procedure convenient just for large river systems, where the tracer experiments become costly and time-consuming. In common applications, tracer experiments are still the most frequently performed and different type of tracers have been successfully used (e.g. sodium chloride, fluorescent dyes or radioisotopes).
Considering the agro-urban environments, the rural channels that compose the network are usually either man-made in concrete with a trapezoidal or rectangular cross-section, or natural-shaped with an irregular cross-section and presence of abundant riparian vegetation that, interacting with the flow, affects both the flow resistance and the solute mixing [58,64,80]. The aim of this study is to perform field tracer experiments to directly measure the longitudinal dispersion coefficient in different rural channel configurations (concrete and vegetated) and, then, to understand which are the most reliable empirical formulas to predict the longitudinal dispersion coefficient within rural channels in the agro-urban environment. As water management companies in agro-urban environments urgently require trustworthy dispersion models but ad hoc field tracer measurements may not always be carried out, the present work aims to provide some guidelines to choose the most appropriate D predictor.
After this introduction, Sect. 2 describes the experimental sites chosen for the tracer experiments and the used instrumentation (Sect. 2.1 ), reports the hydraulic and geometric properties of the channels (Sect. 2.2), illustrates the tracer experiments methodology (Sect. 2.3) and shows how the longitudinal dispersion coefficients were determined (Sect. 2.4). Section 3 reports the results of the experimental campaign (Sects. 3.1 and 3.2) while an extensive analysis of the performance of empirical formulas available in the literature is presented in Sects. 3.3 and 3.4. An almost comprehensive list of the most used formulas for the determination of D is summarised in "Appendix". Finally, Sect. 4 is devoted to conclusions and summarises the principal outcomes.

Material and methods
In the following, the longitudinal, vertical and spanwise coordinate are indicated with the symbols x, y and z, respectively, while t is for the time.

Experimental sites and instrumentation
Measurement campaigns of longitudinal dispersion were performed between July and August 2020, in three different rural channels (mainly intended for irrigation and drainage of agrarian lands) in the agro-urban environment of the metropolitan city of Milan (Northern Italy): 1. the Roggia Delfinona (vegetated channel, Fig. 1a), which has an irregular cross-section with a fine-grain channel bed and highly vegetated sidewalls; 2. the Roggia Gamberina (gravel-bed channel, Fig. 1b), which has a rectangular-shaped cross-section with a coarse gravel channel bed and concrete sidewalls. The interaction between the riparian vegetation and the flow is limited; 3. the Derivatore Vittuone (concrete channel, Fig. 1c), which has a trapezoidal crosssection with both channel bed and sidewalls in concrete. There is no interaction between the riparian vegetation and the flow.
These channels were selected as they represent typical configurations that can be found at the interface between urban and rural environments in high-density anthropic areas. The field tracer experiments were performed in channel reaches 180 m long for the Roggia Delfinona and 250 m long for the Roggia Gamberina and Derivatore Vittuone. In these reaches, there are no bends, abrupt changes in the cross-sectional geometry, junctions with other rural channels or hydraulic structures. The field experiments were carried out in summer when only a steady-state irrigation flow is present in the channels unless heavy rain events occur.
For each study channel, the mean longitudinal velocity was measured in a representative cross-section ( Fig. 2a-c) by means of a portable magnetic-inductive flow meter (Flow Sensor NAUTILUS C 2000, OTT, USA), whose accuracy is ± 1% of the measured values and it is particularly suitable in shallow water situations (minimum operating flow depth equal to 3 cm.). The geometry of the cross-sections was evaluated with a graduated rod, whose precision is ± 1 mm, while the channel slope was measured with a portable Global Positioning System (GPS) (GRS-1, TOPCON, Japan). The used GPS provides an accuracy equal to ± 2.5 mm along the vertical position. In order to reduce the measurement error, we took five measurements for each point and we considered the vertical position equal to the average of the five values. We computed the bottom slope according to the following procedure: the point with the lowest elevation (the thalweg) was identified in both the upstream (the injection point indicated in Fig. 3) and downstream (the third test section indicated in Fig. 3) cross-sections; then, the bottom slope was computed as S = y∕L 3 , where y is the height difference between the two cross-sections and L 3 is their longitudinal distance (see Fig. 3); the obtained slope S was validated by comparing the values obtained through this procedure for five different sub-reaches within the considered channel. For all the cases, the discrepancy among the obtained slope values was negligible ( ± 3%).
Concerning measurements of dispersion coefficient, the tracer involved in each trial was sodium chloride (NaCl) and its concentration in the flow was measured with four electrical conductivity (EC) meters (PCE-PHD 1, PCE, Italy). The instrument precision is ± 2% of the measurement range, while the frequency rate is 0.5 Hz. We decided to use NaCl mainly because it is easily available, inexpensive, colorless and it is generally not toxic at low concentrations. The idea underpinning the employment of EC meters in tracer experiments is that the ion concentration, naturally present in the water, increases after the injection of NaCl tracer. Thus, the background electrical conductivity (i.e. the baseline concentration C 0 ) experiences a sudden growth that is read by the EC meters. It is known that the measured electrical conductivity is directly proportional to the NaCl concentration over a wide range of concentrations [15,41]. For this reason, the EC meters were calibrated against laboratory samples of known salinities between 0 and 1 g/l of NaCl (i.e. the typical concentrations encountered during the field experiments). The employed EC meters also Fig. 1 View of the three study channels: a Roggia Delfinona (vegetated channel); b Roggia Gamberina (gravel-bed channel); c Derivatore Vittuone (concrete channel). The first two channels are located in the Gaggiano municipality while the third channel is inside the Sedriano municipality. The cyan arrows indicate the flow direction measure the water temperature with a precision ± 0.8 C • , allowing us to quantify the water kinematic viscosity and density .

Hydraulic characterisation of the rural channels
Thanks to the regularity of the channels' geometry and the vegetation distribution (see Fig. 1a-c), as well as the long distance from potential sources of disturbance, the study c reaches can be plausibly supposed to be in quasi-uniform flow conditions [14], a situation commonly encountered in rural open-channel flows [13,25,26]. In this condition, the energy slope, the free-surface slope and the bottom slope of the channel are assumed to be equal.
The flow rate in each rural channel was estimated by using the so-called conventional current meter methodology [99]. More in details, the velocities were measured in more than ten vertical transects for each representative cross-section (Fig. 2a-c) by placing the current meter in three selected positions along the vertical direction, i.e., starting from the bottom, equals to 20%, 40% and 80% of the local water height y(z). This method is called the three-point method and it is suitable in situations where the velocities in the vertical direction are unusually distributed, situation that can be encountered in vegetated or partially-vegetated channels. Once acquired the current meter velocity data, the value of the flow rate Q was estimated by means of the mid-section method: in the i-th vertical transect, the mean velocity U i was computed as: where U 20 i , U 40 i and U 80 i are the velocity measured at the three aforementioned elevations. The velocity U i was then assumed to be equal to the mean velocity in a polygon area that extends from halfway the distance between the adjacent transects (on the lateral direction), and from the bottom to the free surface (on the vertical direction). In this way, the flow rate was given by: where n is the total number of verticals and i is the area of the i-th polygon. Figure 2 reports the vertical transects and the elevations used in the current meter measurements for each representative cross-section of the rural channels. The y and z coordinates were normalised by means of the channel width W and the maximum flow depth h, respectively, while the longitudinal mean velocity by using the bulk velocity U b . Fig. 2 also shows the distribution of the normalised longitudinal mean velocity in the cross-sections. It is worth noting how the presence of the riparian vegetation in Roggia Delfinona ( Fig. 2a) creates preferential flow pathways where the velocity reaches its maxima value. This is in agreement with other field studies pursued in vegetated or partly-vegetated open-channels [e.g. 25,26]. The isovels are influenced by the bed morphology [60] in Roggia Gamberina (Fig. 2b). In Derivatore Vittuone, the isovels instead resemble those ones reported in Tominaga et al. [88] for flume experiments in trapezoidal smooth-bed open-channel flows (Fig. 2c). Table 1 lists the hydraulic parameters that were determined by using the aforementioned procedure, along with the geometric features of the rural channels. As reported by Nezu and Nakagawa in [59], the aspect ratio = W∕H dictates the dimension and the intensity of the cellular secondary currents that are generated in the proximity of the sidewalls and corners because of turbulence anisotropy. Looking at Table 1, the Roggia Delfinona and the Derivatore Vittuone present a similar value of , which is comparable to the threshold ≃ 5 − 6 identified by Nezu and Nakagawa [59] and discriminates two different categories of openchannel flows: the narrow open-channels ( ≲ 5 − 6 ) where the presence of secondary currents strongly affects the flow across its entire width, and the wide open-channels ( where the effects of secondary current are bounded near the sidewalls and the flow in the mid-cross-section zone is little altered. However, field and laboratory experiments in openchannel flows [1,59,60] reveal that spanwise heterogeneities in the riverbed trigger the formation of secondary currents also within the flow cross-section. The occurrence of secondary Table 1 Hydraulic conditions of the rural channels According to [14], the bulk velocity is U b = Q∕A ; the hydraulic radius is R h = A∕P , where P is the wetted perimeter; the mean cross-sectional flow depth is H = A∕W ; the shear velocity is  Sketch on the EC meters disposition in a channel for the NaCl tracer measurements. It is also reported the distance from the injection point for all EC meters within each channel currents also modifies the size of the large-scale turbulent coherent structures [9,66,101] that are responsible for the transport of a large portion of the turbulent kinetic energy and momentum [54]. In light of these considerations, we can assume that the Roggia Delfinona and the Derivatore Vittuone probably behave as narrow open-channels, while the Roggia Gamberina as a wide open-channel. The non-dimensional Chézy coefficient = U b ∕u * (the dimensional Chézy coefficient is defined as C = √ g ) represents the other fundamental non-dimensional quantity that is usually taken into account for the hydraulic characterisation of rural channels. The values in Table 1 indicate that, while the Roggia Gamberina and the Derivatore Vittuone have a very similar Chézy coefficient, the Roggia Delfinona has a value indicating greater resistance to the flow motion due to the presence of the riparian vegetation (it is important to recall that ≈ 20 indicates hydraulically smooth-bed configurations while lower values are referred to hydraulically rough-bed conditions). The values obtained are in accordance with other field studies [e.g. 81].

Description of the tracer measurements
The tracer measurements involved four EC meters displaced in three different test sections ( Fig. 3). Three of them were placed at the centreline of the channel, while the fourth was placed 0.5 m from the right sidewall in the third test section. All the EC meters were mounted on wood beams that allowed for their correct positioning. The sensors of the EC meters were located around 0.1 m below the free-surface for each experiment. The distances between the injection point ( x = 0 ) and the three test sections are reported in Fig. 3 and indicated with L 1 , L 2 and L 3 , respectively.
The tracer (NaCl) was injected inside the flow in two different manners: (i) by introducing 400 mg/l of a saline solution made of 2 kg of NaCl dissolved in 5 l of water (relative salt dilution method); (ii) by introducing 1 kg of very fine dry NaCl (dry injection method). Great attention was paid during the injection phase, i.e. it was poured quickly and perpendicular to the flow direction, in order to simulate a slug injection [19]. All the methodologies were repeated three times for each channel so to guarantee six observations for all the channels.
To check if the tracer was uniformly mixed within the monitored channel reach, the criterion proposed by Fischer et al. [34] was used. This criterion identifies the distance L m (the so-called mixing length) required to obtain a complete mixing of the tracer across the entire cross-section. If the data are taken beyond L m , all the distributions in the spanwise direction are likely similar among them. The mixing length can be estimated as [34]: where a is a coefficient that depends on the position of the injection within the crosssection and z = bHu * is the transverse mixing coefficient. Usually, a = 0.1 for injections that take place in the mid-channel [11,34], although some authors [47,79] report that L m results may be overestimated by using that a value. Considering the transverse mixing coefficient z , a theoretical framework for its prediction is still missing [39]. Therefore, to evaluate the coefficient b, we had to rely on experimental works conducted in laboratory flumes or in natural rivers/channels. By taking into account literature results for uniform straight channels, b varies between 0.08−0.26 while its experimental mean value is approximately equal to 0.15 [12,20,34,74]. By using the values reported in Table 1 gives L m equal of about 40 m, 1275 m and 105 m for Roggia Delfinona, Roggia Gamberina and Derivatore Vittuone, respectively. This implies that the third test section in Roggia Delfinona and Derivatore Vittuone was placed beyond the mixing length L m , whereas the third test section in Roggia Gamberina was within L m . Due to the network configuration of the rural channels in the agro-urban environment (i.e. characterised along their path by a dense interconnection with other rural channels and spillways from the urban drainage system), it was not possible to carry out measurements with a test section placed more downstream, especially for the Roggia Gamberina, otherwise, the steady-state conditions would have been no longer applicable. The quality of the tracer measurements was checked by means of the tracer mass recovery ratio R r , expressed as [19,21]: is the mass recovered and M is the mass injected into the channel. The times t a and t b are the elapsed time of the initial and ending edge of the concentration distribution, respectively. They were chosen as the time where C is greater than 3% with respect to the baseline concentration C 0 of the water in the channel. This threshold was adopted because it delimits the natural oscillation (and noise) of the electrical conductivity from the increase due to the passage of the tracer cloud. It was found that, at the third test section, the mean recovery ratio are 0.80, 0.93 and 0.86 for the Roggia Delfinona, Roggia Gamberina and Derivatore Vittuone, respectively. These values are within the acceptable range of 0.8−1.2 indicated by Kadlec and Wallace [43]. The tests carried out with the dry injection method reveal the worst R r values, which are lower by about 4−9% than the values obtained with the relative salt dilution method.

Determination of the longitudinal dispersion coefficient
The shear dispersion process, formulated in the milestone works of Taylor [83][84][85], is the combination of longitudinal advection and lateral diffusion. After a sufficient time (or space), where the vertical and transversal mixing has homogenised the concentration in the entire cross-section, the shear dispersion can be modelled by means of the so-called 1-D advection-dispersion equation for steady flows that is usually presented in the form of a Fickian-type diffusion equation [34,74]: where C and U are the cross-sectional average concentration and velocity, while D is the longitudinal 1-D dispersion coefficient. The specific conditions under which Eq. 5 could be used are [83][84][85]: (i) the solute is non-reactive; (ii) the solute has travelled a distance approximately equal to the mixing length; (iii) the flow is in steady-state conditions; (iv) the turbulence is statistically stationary; (v) the geometry of the cross-section is constant in the flow direction. Under these conditions, the solution of Eq. 5 for a slug injection yields: The quantification of a trustworthy longitudinal dispersion coefficient usually relies on field-work measurements. The literature provides different approaches to obtain D from local tracer measurements, and the most frequently used among those are: the reduction peak method, the direct application of the analytical solution (Eq. 6), the method of moments [28], and the routing procedure [30]. All these approaches have pros and cons, but their estimates of the dispersion coefficients are generally comparable [40]. The dispersion coefficient for Roggia Delfinona and Derivatore Vittuone were derived by comparing the measured concentration distribution in the third test section with the analytical solution (Eq. 6). The calibration involved both D and U, as often done in the literature [see e.g. 3,92,40], and the best fit between the experimental data and the analytical solution was obtained through a least squares minimization.
For what concerns the estimation of the dispersion coefficient of the Roggia Gamberina, whose measurements were carried out within the mixing length, we tried to extrapolate a plausible and realistic dispersion coefficient from the available data. To achieve this goal, we adopted the procedure described in the following. Considering the Eq. 6, the decay of the concentration maximum evolves along the x coordinate as: Consequently, D can be estimated as: Thus, we first computed the mixing length for the Roggia Gamberina by using Eq. 3, founding L m ≈ 1275 m. Successively, we set U = U b and x = L m , so that the only unknown parameter in Eq. 8 was C max (L m ) . We extrapolated this value from the measured C max (L i ) , where i = 1, 2, 3 depend on the test section, by fitting the measured values with an appropriate curve. Since different fitting curves return different pairs of the values C max and D, we selected the curves that provide the upper and lower limits for C max and D. These curves are a power law of the form Y = a 1 X b 1 and a hyperbole of the form of Y = 1∕ � a 2 + b 2 √ X � , see Fig. 4a. In this way, the real values should be accurately bounded by the obtained upper and lower values.
To assess the goodness of this estimation, we compared the concentration maxima of our experiments with those of Sukhodolov et al. [81] in a normalised fashion (Fig. 4b). This work was selected because of the similarity between the hydraulic and geometric characteristics of its studied rivers and the rural channels herein reported (in Sukhodolov et al. [81] = 4.66 − 100.67 and = 2.85 − 11.78 ). The results highlight consistency among our data, Sukhodolov's ones, and the trend for the decay of concentration maxima predicted by the theory (Eq. 7). The extrapolated mean value for the Roggia Gamberina deviates from the linear trend of around 25%. It should be noted that our data collapse well with Eq. 7 starting from x∕W ≈ 30 , somewhat earlier than those of Sukhodolov et al. [81] (that linearly collapse in the log-log plot starting from x∕W ≈ 80 ). This could be attributed to the greater regularity in terms of geometries and flow conditions presents in rural channels when compared to rivers, indicating that C max starts to follow Eq. 7 somewhat earlier.

Monitoring campaign results
An example of the longitudinal concentration distributions C obtained with the relative salt dilution method is presented in Fig. 5 for each channel. Comparing the concentration distributions in the last test section (i.e. the last blue bell and the red bell present in each panel of Fig. 5), we can understand the evolution of the tracer cloud in the spanwise direction. The Roggia Delfinona (Fig. 5a) and Derivatore Vittuone (Fig. 5c) have almost identical concentration distributions, i.e. the difference between the concentration maxima C max is less than 8%, the root mean squared error RMSE is less than 0.0044, and the coefficient of determination R 2 is greater than 0.91 (see [17] for the definition of RMSE and R 2 ). Furthermore, the areas under the time-concentration curves observed in the centreline and near the right sidewall differ by 5% for Roggia Delfinona and 1.5% for Derivatore Vittuone. This indicates that the lateral gradient of concentration has efficiently diffused the tracer in the z-direction. The Roggia Gamberina (Fig. 5b), instead, presents somewhat different central and lateral concentration distributions, revealing a not uniform mixing in the whole cross-section, as also predicted by means of Eq. 3. Figure 6 shows an example of the match between the analytical solution (Eq. 6) and the experimental results for a representative trial (Derivatore Vittuone, relative salt dilution Normalised concentration maxima as function of the non-dimensional distance x/W. Our C max data (the coloured ones) derive from the average of all the six trials, except for the Roggia Gamberina extrapolated value (coloured-hollow symbol) that is the mean of the upper and lower limit, estimated with the fitting procedures, as described in the text. The dashed grey line is Eq. 7 fitting the data of Sukhodolov et al. [81]  Vittuone, respectively. The slightly worse RMSE and R 2 values for the Roggia Delfinona may be attributed to the fact that, when the concentration C is measured at a fixed point during the time span t, the concentration distribution has enough time to evolve, exhibiting a skewed distribution with respect to the theoretical Gaussian one [53]. Thus, since the  (Table 1), the concentration distribution of the former presents a more pronounced skewed shape than the latter's one (this is also noticeable in Fig. 5a, b).
The difference between the measured bulk velocity U b and the velocity U coming from the application of Eq. 6 to the tracer data differ by approximately 7% at most. This discrepancy has the same order of magnitude as the one typically associated with the salt dilution gauging measurements of bulk quantities (i.e. velocity and discharge) in streamflows [18,68]. This further confirms that the last EC meter is placed beyond the mixing length L m since, otherwise, the discrepancy between the bulk velocity of the flow and velocity of the NaCl cloud should have been much higher (e.g. up to a 30% [19]).

Longitudinal dispersion coefficients
The longitudinal dispersion coefficients D obtained from the field measurement campaign by means of the procedures described in Sect. 2.4 are listed in Table 2. The reported values are the average of the three trials for each methodology (i.e. relative salt dilution method and dry injection method). The longitudinal dispersion coefficients  obtained with the dry injection method are systematically higher than those obtained with the salt dilution method. The percentage variation between the first and second row in Table 2 amounts to 11-12% for the Roggia Delfinona and Derivatore Vittuone, while it amounts to 30% for the Roggia Gamberina. Since it is not possible to assess which methodology provides the most accurate D estimation, in the following, we will consider the averaged D values (Table 2), while considering the individual values obtained with the two methodologies as uncertainty bounds.
The obtained D values are in agreement with those obtained for similar rural and irrigational channels, as reported in [34,74].
In a natural open-channel flow, the outer scales used for normalisation purpose are the mean flow depth H (related to the size of the large-scale wedge-type vortices that populate the flow [59]) and the shear velocity u * (notoriously related to the turbulent velocity fluctuations imposed by the flow [69]). These outer scales are commonly used to scale the longitudinal dispersion coefficient D [34], yielding the non-dimensional longitudinal dispersion coefficient D∕Hu * . By looking at the scaled values of D in Table 2, D∕Hu * increases as the aspect ratio increases. This is presumably due to different factors such as: the increase of the variation of the spanwise mean velocity profiles, and the decrease of the intensity of the secondary currents (that causes the reduction of the transverse dispersion coefficient). The increase of D∕Hu * with is in agreement with past studies [48,78,102], although the literature reports no obvious linear or non-linear relationships between these two quantities [89].

Comparison with the existing formulas
The values of D presented in Table 2 are associated to particular hydraulic and geometric conditions (Table 1). For this reason, it is common practice to use equations in order to determine D in situations other than those encountered during the experimental campaigns. In particular, the building-up of reliable water quality models needs a suitable quantification of D. The literature provides a plethora of formulas to calculate D (or D∕Hu * ), but the applicability to rural channels in the agro-urban environment is not so straightforward since these formulas have been derived using different methodologies (i.e. mathematical, theoretical, empirical) and fitted to field data coming from somewhat different situations (i.e. ranging from very wide rivers [e.g. 62, 100] to small-steep streams [e.g. 21,81]).
To identify the most accurate formulas, 31 expressions from 24 studies were selected and compared to the field measurements ("Appendix" reports an exhaustive, albeit not complete, description for each formulation in Table 4). The comparison between the empirical predictions ( D predicted ) and the experimental observations ( D observed ) are displayed in Fig. 7. As well explained by Fischer et al. [34], the empirical formulations for the estimation of the longitudinal dispersion coefficient do not claim to be highly accurate, but, rather, they try to identify the order of magnitude of the phenomenon. Nevertheless, and quite surprisingly, many of the existing formulas overestimate by even two orders of magnitude the field results (Fig. 7).
The enlargement in Fig. 7 allows us to better identify which formulas give the best estimates while Table 3 quantitatively reports the values of the estimates and the percentage relative errors R.E., which is defined as:  Fig. 7 Comparison between the normalised measured ( D observed ) and predicted ( D predicted ) longitudinal dispersion coefficient. Green, orange and light-blue symbols stand for the Roggia Delfinona (vegetated channel), Roggia Gamberina (gravel-bed channel) and Derivatore Vittune (concrete channel), respectively. In the enlargement, the dashed lines represent a relative error (R.E.) equal to ± 50% . For details about the formulas, see Table 4 in "Appendix" From Fig. 7 and Table 3, it emerges that the empirical formulas able to predict the experimental results for the three studied rural channels with a R.E. lower than (or equal to) ± 50% are:  [52]) experiments. The presented work aims to test these formulas over a different dataset, whose hydraulic conditions (see Table 1) are not expected to completely fall within the range of the formulas' training datasets. The majority of the formulas fail to properly predict the longitudinal dispersion coefficient D in rural channels and this is likely due to the fact that these formulas try to include and integrate the effects of different natural factors-such as the presence of dead-zones [16,90], river bedforms [6], hydraulic structures [98], meanders [30,32], compound Table 3 Formulas that predict at least one of the studied rural channels with a relative error (R.E.) of ≈ 50% For the sake of clarity, only the name of the first author of the study is reported in the first column. For more details about the formulas, see Table 4 in "Appendix"

Roggia Delfinona
Roggia Gamberina Derivatore Vittuone cross-sections [5], river confluences [70] and riparian vegetation [58,64,80]-in a sort of universal formula for natural streams. Furthermore, Table 4 in "Appendix" shows that almost all the formulas have been developed by using field data taken from the same few works [10,21,31,37,38,56,62,74,86,100], i.e. a similar dataset is used for the development of the empirical predictors. When it is required to predict D in different watercourses with respect to the dataset (as the rural channels herein analysed), the empirical formulations give a considerable error and only those with a more robust theoretical background offer a reliable estimate. Some valuable and more recent D dataset coming from field studies [e.g. 22,73,40,48] should be included in the future, especially when soft computing techniques are used to generate new predictors [96]. Further focusing on the formulas that give the best results, we can find the following indications: (i) for rural channels with a low aspect ratio ( ≲ 10 ), Parker 1961 [63] and Koussis & Rodríguez-Mirasol 1998 [47] give the best estimates, and the Chézy coefficient seems not to be very relevant; (ii) for rural channels with no or minimal riparian vegetation ( ≳ 10 ), Deng  Furthermore, Parker 1961 [63] provides a reasonable estimation in all situation but the relative error increases up to 70%. In particular, [63] underestimates the value of D as the aspect ratio increases, as it can be noted by manipulating the original formula as follows: As increases, R h → H and D H u * → 20.195 . Consequently, this drawback of Eq. 10 should be taken into account in real applications.
The formulation of Magazine et al. 1988 [52] seems to be excessively influenced by the type of roughness presents in the channel, which is expressed both by and n w (i.e. the Manning's coefficient of the sidewalls of the channel, see Sect. 3.4). Thus, it gives roughly the same value of D for both Roggia Gamberina and Derivatore Vittuone, whereas they have similar channel's roughness but rather different cross-sectional characteristics. Therefore, its practical applicability is discouraged in those contexts where the geometry of the watercourses changes considerably (such as in rural channel networks).

Prediction of dispersion with different flow rates
Since the flow rate inside rural channels sensibly changes during flood events, it is important that the chosen predictor give reasonable estimates of D also when the flow rate Q changes. In order to study how the selected formulas (Table 3) work against a variation in the flow rate, we performed a series of hydraulic simulations using the 1-D numerical model HEC-RAS (Hydrologic Engineering Center-River Analysis System) [7].
The simulations for the three study channels were carried out considering a reach 250 m long with constant cross-sections that is equal to the representative ones reported in Fig. 2 but extended in the vertical direction according to the geometric data gathered during the field campaigns.
A first set of simulations was done in order to calibrate the Manning's coefficients, which usually represents the most significant element in model calibration [57,65], both for the bed and sidewalls with respect to the steady-state hydraulic conditions encountered in the field campaigns (Table 1). For the steady-state simulations, a downstream boundary condition together with a steady flow rate Q were required since the flow regime is always subcritical ( Table 1). The used downstream boundary condition was the energy line slope assumed equal to the bottom slope S, i.e. assuming uniform flow conditions. The resulting Manning's coefficients for the sidewalls ( n w ) and the bed ( n b ) were equal to: n w = 0.098 s/m ( [7] and in [14] for similar conditions. After this preliminary calibration, a series of steady-state simulations with different Q values was performed by assuming constant Manning's coefficients as the flow rate change. The results are displayed in Fig. 8 and show that the evolution of the mean flow depth H (Fig. 8a) and the bulk velocity U b (Fig. 8b) is in agreement with the past literature [e.g. 14] as the flow rate Q increases. Besides, the curves of the aspect ratio (Fig. 8c) show that, starting from Q∕Q max ≈ 0.1 (where Q max is the bankfull flow rate), Roggia Delfinona and Derivatore Vittuone have < 10 , while Roggia Gamberina has always > 10 . Finally, for what concerns the flow resistance expressed by the Chézy coefficients (Fig. 8d), Roggia Gamberina and Derivatore Vittuone show an increasing value of as the flow rate increases. In contrast, starting from Q∕Q max ≈ 0.1 , Roggia Delfinona has an almost constant value of .  Figure 9 reports the evolution of the longitudinal dispersion coefficient as the flow rate increases for each study channel. As reported in the literature [40,74,94], it is expected that D increases with Q. This behaviour is caught by almost all the selected formulas listed in Table 3. Only the formula proposed by Koussis & Rodríguez-Mirasol 1998 [47] does not  Table 3 follow this trend, giving an almost constant trend starting from Q∕Q max > 0.3 for Roggia Delfinona (Fig. 9a) or a decreasing trend for Derivatore Vittuone (Fig. 9c). For this reason, the use of this formula should be avoided when different flow rates are considered. It is worth noting that some formulas are more sensitive to changes in the cross-sectional shape than the others. Considering Roggia Delfinona, which presents an irregular cross-section (Fig. 2a), the formulas of Parker 1961 [63] and Sukhodolov et al. 1997 [81] remain smooth as Q increases, while Liu 1977 [51], Deng et al. 2001 [20] and also Noori et al. 2017 [61] are more winding, following the cross-sectional change. Furthermore, all the values of D given by the selected predictors for the vegetated channel are bounded by the formulas of Noori et al. 2017 [61] and Sukhodolov et al. 1997 [81], which can be seen as the upper and lower boundaries, respectively. Finally, Fig. 9b shows the different behaviour of the formulas of Deng et al. 2001 I [20] and Wang & Huai 2016 I [95] as the flow rate increase, highlighting comparable results up to Q∕Q max < 0.4 and significant divergence beyond this value. This occur also for the Derivatore Vittuone (Fig. 9c), although the two formulas start to deviate earlier ( Q∕Q max ≈ 0.1).

Conclusions
The longitudinal dispersion coefficient in typical rural channels (Roggia Delfinona, Roggia Gamberina and Derivatore Vittuone) located in an agro-urban environment was studied in order to understand which empirical formulas are the most suitable when carrying on onedimensional models at the network scale. The tracer measurements were done using two different methodologies, namely the relative salt dilution and the dry injection methods. The longitudinal dispersion coefficients D in Roggia Delfinona and in Derivatore Vittuone were quantified by directly applying the analytical solution for slug injection, while in Roggia Gamberina an extrapolated value was estimated since the tracer measurements were done within the mixing length. The two tracer measurement methodologies give comparable results, although tracer mass recovery ratio R r is lower in the case of the dry injection method, likely because a small fraction of the salt settles down during the process (around a 5-10%). On the contrary, the longitudinal dispersion coefficient D results to be higher (around 10%) when the salt dilution method is adopted. By averaging the values of D given by the two methodologies the following longitudinal dispersion coefficient values were estimated: 0.437 m 2 /s for Roggia Delfinona (completely vegetated channel), 1.005 m 2 /s for Roggia Gamberina (rectangular channel with a gravel bed) and 0.545 for Derivatore Vittuone (completely concrete channel).
The comparison between the field data and 31 formulas taken from the literature reveals that only 9 of those formulas are able to predict D, with a relative error R.E. = ± 50% , at least in one rural channel. Furthermore, the simulations with growing flow rates carried out with the hydraulic model HEC-RAS reveal that not all the selected formulations yield increasing values of D as Q increases, as it is instead ascertained by the literature [40,74,94].
After this examination, we indicate the adoption of the following predictive formulas when rural channels are present in the study domain: (i) if the aspect ratio is ≲ 10 , Parker 1961 [63] gives the best estimate with expected R.E. < 20% ; (ii) for man-made rural channels with Chézy coefficient ≳ 10 , both Deng  Most of the formulas that reasonably predict our field results are theory-based, whereas the empirical formulas, which are instead based on statistical approaches or soft computing techniques, seem to be less adequate to envisage the longitudinal dispersion coefficient in agro-urban environments. This could be attributed to the fact that, in most cases, the training dataset of these latter kinds of formulas is the same and based on tracer experiments accomplished in wide rivers in the U.S.A. (see Table 4). As a consequence, other waterway configurations, such as rural channels, are not properly taken into account.
As discussed by Wallis and Manson [94], a possible way to obtain better predictors could be to build different formulas for different hydraulic/geometric situations rather than try to produce a single and universal formula. An example of this approach is given by Wallis and Heron [93], who propose an empirical formula for small rivers with satisfactory results. The realisation of focused and precise formulas for the determination of the dispersion coefficients, each refined for different situations, could be a possible way to improve our knowledge and the reliability of water quality models.
Finally, we mention another issue that future works should address and that regards vegetation maintenance. As reported by Västilä et al. [91], the extensive vegetation cutting practised to increase the flow conveyance inside rural channels also impacts the transport of pollutants, especially in the longitudinal direction, leading to the passage of faster and higher concentration peaks. Hence, not only the spatial variability of rural channel characteristics (e.g. channel morphology and hydrodynamic regime) but also the seasonal variability due to vegetation maintenance practices have to be considered in evaluating the longitudinal dispersion coefficient in rural channels.

Appendix: Formulas for the estimation of the longitudinal dispersion coefficient
Many authors have proposed formulas to estimate the one-dimensional longitudinal dispersion coefficient D so far, and Table 4 reports a list of the most famous and used ones. The 24 selected studies offer 31 formulas that are suitable for different hydraulic and geometric configurations of channels and natural rivers.
In the column 'Derivation', the method underpinning the formula proposed by the authors is reported as: Mathematical when referring to analytical derivations that are only based on working hypotheses (e.g. the Reynolds analogy, logarithmic distribution of the longitudinal mean velocity profile, etc.); Semi-Theoretical when referring to those formulas that are based on physically-sound assumptions together with coefficients coming from experimental experiences; Empirical when referring to those formulas that are mainly driven by the experimental data. This last case is further subdivided in Statistical, which indicates more conventional methods to obtain the searched equations (e.g. dimensional analysis, regression methods, analysis of variance, etc.), and Soft Computing, which includes different techniques based on artificial intelligence and machine learning (e.g artificial neural network, genetic algorithms, model tree, particle swarm optimization, etc.).