A GRASS GIS application for vertical sorting of sediments analysis in River Dynamics

The extreme versatility in different research fields of GRASS GIS is well known. A tool for the vertical sorting of sediments in river dynamics analysis is illustrated in this work. In particular, a GRASS GIS python module has been written which implements a forecasting sorting model by Blom&Parker (2006) to analyze river bed composition’s evolution in depth in terms of grainsize. The module takes a DEM and information relative to the bed load transport composition as input. It works in two different and consecutive phases: the first one uses the GRASS capabilities in analyzing geometrical features of the river bed along a chosen river reach, the second phase is the "numerical" one and implements the forecasting model itself, then executes statistical analyses and draws graphs, by the means of matplotlib library. Moreover, a specific procedure for the import of a laser scanner cloud of points is implemented, in case the raster DEM map is not available. At the moment, the module has been applied using flumes data from Saint Anthony Falls Laboratory (Minneapolis, MN) and some first results have been obtained, but the "testing" phase on other flume’s data is still in progress. Moreover the module has been written for Geoinformatics FCE CTU 2011 39 Minelli A. et al: A GRASS GIS application for vertical sorting of sediments analysis in River Dynamics GRASS 65 on a Ubuntu Linux machine, even if the debugging of a GRASS 64, Windows version, is also in progress. The final aim of this work is the application of the model on natural rivers, but there are still some drawbacks. First of all the need of a high resolution DEM in input, secondly the number and type of data in input (for example the bed load composition in volume fraction per each size considered) which is not easily obtainable, so the best solution is represented by testing the model on a well instrumented river reach to export in future the forecasting method to un-instrumented reaches.


Introduction
The vertical sorting of sediments in gravel and sandy-gravel bed rivers is one of the factors directly involved in the evolution processes of river morphology [1].In order to obtain a better understanding and modeling of the main physical processes occurring in vertical sorting and due to the large amount of field work behind grainsize analyses on a real river, some prediction methods for the study of vertical sorting of sediments have been formulated in time [2]; [3]; [4]; [5]; [6]; [7].
In this work, the Equilibrium Sorting Model [8] has been implemented on the Open Source Geographical Informative System GRASS [9], and a specific tool, which can be integrated in the main GIS structure, has been created.The specific choice to implement the theoretical model in GIS has been realised in order to achieve the aim of computerizing the model prediction, simplifying the field work of sampling not only but also sieving and classifying a terrain.In addition to this an other main scope is represented by obtaining the connection between the predictive capabilities of the model to the GIS capabilities in analyzing terrain morphology.
In fact, the model application requires a priori an accurate knowledge and extraction of bedforms characteristics on the river bed surface, that GRASS GIS can automatically provide.Moreover, it is interesting testing the application of GIS at microscopical scale (i.e.grainsize specific scale).
First of all, starting from the hypothesis that, there is equilibrium in entrainment and deposition fluxes of sediments in steady conditions, a resolution method for the model has been formulated.Starting from few geographical information (such as for example a Digital Elevation Model) a geographical and geometrical analysis concerning bed form dimensions, is performed.Secondly some statistical analyses are executed on the main geometrical features of the bedforms.Once the geometry characteristics of the river bed and statistical distribution of the model's parameters are analyzed, the prediction can be performed: the resolution method uses an iterative procedure to find the vertical grainsize composition values for each point and for each different value of depth (over the bed surface) examined.
The code written (v.eqsm.py) is a python script for GRASS GIS and uses python scientific, graphical and mathematical libraries.Currently the module is in testing phase.The predictive processes have been verified step-by-step during the implementation and, along the following calibration phase, experimental data from the main channel of the Saint Antony Falls Laboratory (Minneapolis, MN, USA) have been used [10].This choice has been lead due to the necessity of calibrating the tool in a controlled environment with less external in-Minelli A. et al: A GRASS GIS application for vertical sorting of sediments analysis in River Dynamics fluences than possible.At present time, the model predictions are being verified using flume data although the module has to be applied on real river data as soon as possible, due to the necessity of having well developed bed forms to formulate a reasonable prediction on.
Currently the tool implementing the Equilibrium Sorting Model takes in input, as required by the model, both data relative to grainsizes to formulate predictions on it and bedload composition for the grainsize examined.The tool generated as output the prediction in depth of the sediment sorting for one or more specific points investigated and for a specified number of layers in depth.

The Equilibrium Sorting Model in GIS
The equilibrium sorting model is a vertical sorting model of sediments formulated by A. Blom and G. Parker (2006) [11] as a steady simplification of the main framework for the vertical sorting phenomena modeling [12].The model supposes the continuity of terrain in depth [13] in the way that not only a terrain "active layer" is involved in the entrainment and deposition processes, below the bed surface, but each grain, of a specific size, at a specific depth has a certain probability to be held by the current (Fig. 1).Moreover, the model is grainsize specific, in the way that different grainsizes can be examined: for this first implementation a bimodal mixture has been considered.In addition to this, the GRASS GIS use, allows to analyze the river bed, considering the bedform presence along the river reach examined, as accurately as the elevation map resolution is and allows the user to apply the model starting by any geographical-referred data.Deepening into the model, the evaluation of vertical sorting is done point-to-point along a river reach; it is possible to simplify the river reach profile as reported in Fig. 2. Suppose to separate each bedform in an ascending (stoss) and a descending (lee) face, it is possible to hypothesize that on lee face only deposition occurs and on stoss face both entrainment and deposition occur (Fig. 2).The equation depth specific, simulating the mass conservation through entrainment and deposition fluxes over an infinitesimal long reach segment in time changes as follows: • C i is the mass concentration of sediments of grainsize i;  • c b is the total mass concentration of sediments in the bed surface;

Geoinformatics FCE CTU 2011
• P s is the probability density function that a sediment can be held by the current; • F i is the volume fraction of sediments of grainsize i; • D ei is the deposition rate of a specific grainsize; • E ei is the entrainment rate of a specific grainsize.
All terms are averaged over a series of bedforms (reach length).
Moreover the probability density function (P s ) depends on the mean elevation of bed river, averaged over a series of bedforms ( ηa ).
If the framework is reduced at steady conditions and the equilibrium between entrainment and deposition rates is hypothesized, the prediction on lee faces is formulated and the only unknown variable is represented by F i : where: • F leelociE is the value of volume fraction of grainsize i sediments in the lee deposit (downstream the bedform); • F topi is the volume fraction of sediments of grainsize i overtaking the bedform; • J is the Heaviside step function; • η bmin and η bmax are the minimum and maximum value recorded for the bottom of each bedform in the river reach; • λ and ∆, as shown in Fig. 2, are the length and drop of a single bedform; • ω i is the lee sorting function, able to predict the depth where we can find, in the lee deposit, the i-th grainsize; • pbE is the probability density function for the bed surface not to be greater than a certain value, depends on the mean bed level.

Minelli A. et al: A GRASS GIS application for vertical sorting of sediments analysis in River Dynamics
These parameters depend both on the bedload composition (F topi ) and on geometrical features characterizing the river bed (e.g.λ, ∆ and η), where the first (bedload composition) is given in input to the module, while the second (geometry) is calculated by GRASS GIS examining the DEM.The solution method and consequently the procedure implemented (considering two different grainsizes) allows to solve iteratively a set of three equations involving both sediment transport and bedform geometrical parameters to find the "right values" for F i at each fixed depth value in each fixed point.The equations to be solved in sequence for each point examined are: 1. φ mleeE = N i=1 φ i * F leeiE -which gives as result the mean composition of the lee deposit of the bedform where the point, chosen to perform prediction, falls (φ mleeE ); F leeiE is the volume fraction of sediments of i-th grainsize in the lee deposit; τ bE -which gives as result the values of lee sorting parameter (δ iE ) for each grainsize examined, this parameter is involved in the definition of lee sorting function (ω iE ); ρ, ρ s are the density of water and material (respectively), τ bE is the bed shear stress (averaged on the river reach); λ∆ pbE * dη b which gives as result the volume fraction of elements of i-th grainsize in the lee deposit (F leeiE ) this value is converted subsequently in the corresponding F i value using the relation: F aiE is the volume fraction of sediments of i-th grainsize in bedload transport material, ν is the angle of repose of material, z is the nondimensional height (relative to the mean height of the specific bedform) of the point where prediction has to be performed.
The analysis starts hypothesizing an equal volume concentration for each grainsize examined (F 1 = F 2 = 50%) then, the three equations are solved and the values obtained for Fi are compared to the hypothesized values.If values are not matching, for each volume fraction, the mean value between hypothesized and obtained is taken, another iteration is performed, a second set of values for F i are obtained and so on, until the values match.

v.eqsm.py -workflow and GUI:
v.eqsm.py is the new GRASS module which allows to predict the vertical sorting of sediments for two different grainsizes (bimodal mixture), starting by DEM and some bedload transport information.The module is written in python language and it is available for download at the following link: http://svn.osgeo.org/grass/grass-addons/vector/v.eqsm/.It uses scientific, graphical and mathematical python libraries i.e. scipy, math and pyplot.The module works as shown in Fig. 3 diagram and it is described step by step as follows: Minelli A. et al: A GRASS GIS application for vertical sorting of sediments analysis in River Dynamics • input information are relative to: grainsizes to perform the prediction on, the digital elevation map of the river reach and about the current and the bedload transport (sediment concentration of studied grainsizes, and density of the bedload); • the first phase consists on calculation by the means of GIS of geometrical features of the river reach for the entire reach, starting from DEM as input data; • statistical analyses on geometrical parameters values are executed and graphical representation concerning how these values are distributed are produced; • GRASS GIS individuates the lee faces (where the model can give a prediction) and there's a interactive procedure to choose the point (or points) where to perform the prediction; • in the last phase, the iterative procedure is used to perform the prediction on selected point(s), the model; since a good convergency is found, results are consequently written in a text report, in addition to this graphs illustrate the vertical sorting in depth for the two grainsizes, and vector file of points is also provided where prediction is stored.
The module Graphical User Interface is reported in Figure 4. To provide the prediction in depth, it is possible to specify a maximum prediction depth ("dmax" parameter) and the number of steps to perform the analysis on ("nstep" parameter).By using this approach the user can control the accuracy of the analysis by defining a more or less continuous pattern for the prediction.Moreover, if a raster map of elevation information is not available, the module can take in input a laser scanner cloud of points ("DEM importation" section on the GUI).

v.eqsm.py -an experiment:
The v.eqsm.pytesting phase is still running for the following reasons: first of all it was necessary to test the module on trusted data, secondarily it was necessary to have data from a controlled environment, as it is the flume one, to have higher accuracy than the possible.The model works efficiently only if bedforms are well defined so, it was quite difficult to find laboratory data (laboratory with a flume sufficiently extended, where water conditions are similar to natural ones) which had these characteristics.At last, to export the model on real river, it is necessary to have high accuracy DEM For these specific grainsizes, we have to obtain the following volume fraction values for the bed surface: F 1 , F 2 = 0.799, 0.211.Then, the DEM is imported and the concentration of chosen grainsizes is read in the recirculation sample.The aim of the experiment is to predict the grainsize fraction on the channel bed surface, starting from bedload composition informations and channel bed morphology.Resuming, the input parameters are: • grainsize dimensions to investigate: φ 1 , φ 2 = −3.41,−0.79; • bedload material density: 1744 N/mc; • average bed shear stress (averaged on the river reach length): 10.47 Pa; • volume fraction (%) of chosen grainsizes (φ 1 , φ 2 ) in the bedload: 0.14556, 0.70888; • maximum depth of prediction: 10 mm; • nr. of steps in depth to investigate: 5.
The command syntax given from terminal is: eqsm.py dem=dem vector=l1 dmax=10 nstep=5 res=10 lee_out=lee_out \ point_out=point_out fione=-0.79428fitwo=-3.41009ros=1744 to=12.752526971\ faone=0.659383856fatwo=0.340616144n_iter=20 txtout=texout graphout=graph Other input data are DEM, the vector file of the river reach axis (for this test a reach of 10 meters length has been considered) and the name for the vector file of points where prediction are stored, the text file of report where predictions are written and the name for graphics in output.

Minelli A. et al: A GRASS GIS application for vertical sorting of sediments analysis in River Dynamics
The module takes some minutes to work and then gives in output some graphs, reports and vector files (Fig. 6): • the graphic of river reach profile (Fig. 6a), reporting the concavity and convexity rates, respecting the mean bed level; • a graphic of delta parameter (total height of each bedform) raw and logarithmic distribution (Fig. 6b and 6c); • the graphics of the probability density distribution (using the Weibull distribution, Fig. 6d and 6e) and the Weibull plots (Fig. 6f and 6g) for maximum and minimum height of each dune, respecting the mean bed level; • the graphic reporting the volume fraction for each grainsize in depth (Fig. 6h); • an output text file reporting both an input data summary and the output of the prediction (Fig. 6i); • the vector file of the lee faces along the entire river reach (Fig. 6j): in each feature some information about the morphology of the lee face and the entire bedform are stored, e.g. the lee face length, the dune length and total height, the position of the dune respecting the mean bed level etc.; • the vector file of point where informations about the prediction are stored (Fig. 6k).
Figure 6a: the river reach profile.
From the graphic in Figure 6a, reporting the reach profile extracted (using r.profile) from the DEM in Figure 5, it is clear that the profile is mostly convex, so the bed surface is often higher than the mean bed level, and downstream, near the end of the reach, there's a ca.8cm peak which is probably an instrumental recording error (laser scanner bad reflection), it is reasonable not to consider into the prediction that specific bedform.That peak is highlighted also in Figure 6b, where there's only one occurrence of a 8cm delta value, at anyway the frequency delta distribution remains logarithmic (a lot of low delta values) and the graph in  Fig. 6c fits well that distribution.The presence of low delta values is due to the not-so-well developed bedforms on the channel, and this is due to, on turn, the fact that the experiments collected in Streamlab 2006 dataset were not finalized to this specific kind of tests.So more runs with higher water discharge values should be performed to obtain a better developed bed surface (with dunes instead of ripples - [14]; [15]; [16]).From graphs in Figures 6d and 6e it is clear an almost perfect mirroring of etat and etab values, this means that the river profile is, for a high percentage, "average": where the top of the bedform (etat) lies over the mean bed level and the bottom (etab) lies below the mean bed level.In Figures 6f  and 6g the Weibull plot is reported and the envelope fits correctly the extreme values of the parameters [17].From the graphic reporting the prediction output (Fig. 6h) we can evidence  an almost linear behavior of grainsize composition with depth.This result has been also found by Blom (2008) [18].Figures 6j and 6k, show how results are stored in GRASS GIS, the advantage to have these vector files in output is that it is always possible to export the results and to deal with the complete output of the grainsize prediction and morphological analysis.

Conclusions and future developements:
Comparing the values obtained from the forecasting to the ones expected (recorded in the dataset), we can see by the graphics, that there are very low differences.It is evident that  more accurate studies on similar cases have to be done, but these first results are encouraging in any case.It has to be said that it's not easy to find data from appropriate structures as the main channel of SAFL, which is sufficiently large in the way that it has a behavior that can be compared, in some aspect, to a real river one and which allows to test efficiently the module, since there are few structures like that in the world and the river bed morphology has to be particular (bedforms have to be well evidenced and defined, no ripples).At anyway, using available data, a statistical analysis on river reach morphology and geometrical parameters distribution is performed.Moreover the prediction results are given in three different formats (text file, graphics and vector files).In the end, the Equilibrium Sorting Model and the module written (eqsm.py)can help researcher and worker in sampling and classification work  because, if we can apply this model to real rivers, a prediction of river bed composition is formulated in the way that the sampling work can be limited to the minimum and necessary.This is a particular field of application for GIS, most of all because the scale of the analysis is very large and units are millimeters instead of meters (normal unit of measure for GIS data).So the application to grainsize distribution in depth is quite original, but the results obtained are encouraging.Some limitation is still present both in the module and in the model.For example, this model can predict the bed composition with maximum depth equal to the total bedform height because of the peculiar definition of some model parameters, consequently the model gives "right" prediction at higher depth, only if the bedforms are well developed; currently, this represents an open task for the model.By the other hand, the module examines only two grainsizes and it is written for GRASS 6.5 and Linux systems only, even if a Windows version is in debug phase.Moreover, to use the module on real river we need some high resolution data e.g. the DEM, which must have a resolution at least  comparable to that one of the elements investigated (bedform dimensions) and a so accurate kind of data is often difficult to find, unless specific and/or deep studies.

Figure 1 :
Figure 1: the real behavior of terrain erosion in depth (left) and simplified schema supposing erosion only on a strata with L a thickness (right) -image courtesy of G. Parker, C. Paola and S. Leclair.
Minelli A. et al: A GRASS GIS application for vertical sorting of sediments analysis in River Dynamics

Figure 2 :
Figure 2: (on the left) the river reach simplified profile and some geometrical parameters describing dune profile: lengh of the dune and of stoss and lee face (λ), relative height of the dune (∆), total height of bottom of each dune and mean bed level (η); (on the right) entrainment and deposition rates on stoss and lee faces with solid discharge crossing the top of the dune (q top ) -image courtesy of A. Blom and G. Parker.
(or bathymetries) to identify and geometrically define bedforms characteristics.Data from Saint Anthony Falls Laboratory (Minneapolis, MN, USA) has been used for this work.Particularly some experiments have been done to test the reliability of the module using the large amount of data in Streamlab Dataset (2006), which is freely accessible by the following URL: https://repository.nced.umn.edu/browse.php?dataset_id=8.This dataset has been collected from March 2006 to November 2007 and it derives from a series of experiments executed in the main channel of the laboratory (dimensions: 84 x 2.75 x 1.8 m).Additional runs have been executed with different purposes (e.g.testing sieving methods or biological parameters under specific conditions) and both bed sample composition in depth and recirculation sample composition are available.Moreover laser scanner outputMinelli A. et al: A GRASS GIS application for vertical sorting of sediments analysis in River Dynamicsare available and each sample or cloud of points have a temporal reference, consequently, some bed samples have been selected in the way that, for each bed sample (which represents the expected result -it can be used to verify the prediction) it is possible to employ the nearest (in time) recirculation sample and cloud of points as input data.Data recorded in August 2006 have been used for this work.In particular, three samples from the main channel of Streamlab 2006 dataset are available for testing v.eqsm.py(Fig.5).

Figure 5 :
Figure 5: Digital Elevation Model of the main channel (unit is 1mm); in red, samples available for testing eqsm.py.Sample #1 has been chosen.
Figure 6b: total height frequency values.
Figure 6h: graphic of volume fraction distribution in depth for examined grainsizes.

Figure 6i :
Figure 6i: text file with resume of input data and prediction results.

Figure 6j :
Figure 6j: in the GRASS Monitor -the extracted lee faces vector file (in white) plotted over the entire studied reach (in black) and the DEM (in green); in the output Form of the query -geometrical informations stored in each lee face; in the bottom view -a 3D visualization by NVIZ of the DEM and the studied river reach (in black).

Figure 6k :
Figure6k: in the GRASS Monitor -the vector points file (in yellow) where prediction is executed and results are stored, plotted over the lee faces (white) and reach (black) vector file; in the output Form of the query -prediction results stored in the points vector file; in the bottom view -the river reach profile (in blue) and the vector points file (in yellow) where prediction is executed.