Bedforms-ATM, an open source software to analyze the scale-based hierarchies and dimensionality of natural bed forms

Bedforms-ATM (Bed forms analysis toolkit for multiscale modeling) is a software designed to hierarchize and quantify the dimensionality of natural bed forms fields. It comprises four modular applications, namely: (1) wavelet analysis, (2) Hovmöller analysis, (3) multiscale discrimination, and (4) threedimensionality analysis. Bedforms-ATM also provides insights on bed form systems dynamics and their interrelationship with the surrounding hydrodynamic characteristics. The software structure encourages its expandability via the collaboration from the community of users. Both fluvial and synthetic bed form data accompany Bedforms-ATM. © 2018 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Although the need of discriminating bed forms scales (e.g., bars, dunes and ripples) that are commonly present in unidirectional, oscillatory and combined flows has been highlighted by many researchers [e.g. [1][2][3], no standard nomenclature and procedure exist to systematically identify the scale and magnitude of such sedimentary features [3]. As a result, the discrimination between ripples and dunes is taken as obvious and many authors provide * Corresponding author at: Universidad del Norte, Barranquilla, Colombia.
E-mail address: rgutierrezll@uninorte.edu.co (R.R. Gutierrez). limited or no explanation regarding the criteria they followed for their discrimination in field and laboratory studies [1]. Reference [3] introduced a spectral-structural technique that allowed for a quantitative characterization of three bed form hierarchies ( Fig. 1-a and 1-b) even for curved bed form fields (Fig. 1-c). However, it is not able to differentiate between two dimensional (2D) and three dimensional (3D) bed form features. The structural component of the technique presented by [3] obeys the scalebased classification scheme proposed by [4]. It defines ripples as features having wavelengths less than 0.6 m; dunes as features characterized by wavelength thresholds of 5, 10, and 100 m determining small, medium and large dunes; and bars as features having wavelengths higher than 100 m. The spectral component of the technique proposed by [3] is based on the coupled application of robust spline filters and the one dimensional continuous wavelet transforms frame developed by [5].
Bed forms in sedimentary environments are predominantly 3D in planform [3]. This characteristic induces marked differences in local turbulence, drag coefficient, and dispersal patterns of sediments [6][7][8]. Therefore, a complete bed forms description requires distinguishing between 2D and 3D bed forms. To this end, [9] proposed a quantitative continuum metric, named 3-dimensionality index (3DI), which varies between 0 (strictly 2D) and 1 (strictly 3D). 3DI (cf. Fig. 1) quantifies the normalized autocorrelation between the fluctuations of a bed form field η(x i , y i ) respect to the local cross-sectional meanη x (y i ), and the fluctuations of the local stream-wise meanη y (x i ) respect to the bed form field meanη xy (Eqs. (1)-(4)). Thus, a strictly 2D bed form has a flat crest perpendicular to the flow, and the more irregular the cross-sectional profiles are, the greater the value of 3DI.
Most computer programs designed to analyze bed forms fields (e.g., Bed form tracking tool by [2]) are not currently publicly available or are discontinued. Herein we present Bedforms-ATM (Bed forms analysis toolkit for multiscale modeling), a MATLAB open source software that performs the variability analysis of bed form fields by quantitatively describing both their hierarchies and dimensionality.
This paper presents the technical aspects of the development and capabilities of Bedforms-ATM and is organized as follows: Section 2 describes the data processing methods used in this study, Section 3 provides a functional description of Bedforms-ATM, Section 4 focuses on its implementation, Section 5 presents case histories of its application on bed forms and synthetic bed forms, and Section 6 discusses on its potential impact. Finally, Section 7 outlines the main conclusions of this study. Fig. 1 shows an approximately equally spaced (i.e., ∆x ≈ constant) elevation bed form profile η(x i , y J ) ∈ η(x i , y j ), also represented herein by the vector η y,J (x i ). The software performs the wavelet analysis (Application 1) of each bed form profile in η(x i , y j ) by adapting the one dimensional continuous wavelet analysis package developed by [5]. The continuous wavelet transform of η y,J (x i ) is mathematically defined as the convolution of it with a continuous wavelet function. By performing a Monte Carlo analysis, [5] estimated empirical factors that allow for determining both the local wavelet power spectrum (i.e., the modulus of the wavelet transform) and its cone of influence (i.e., the divide that discriminates reliable from spurious local power spectrum results). This is done for continuous wavelet functions such as Morlet, Paul, Marr and the derivative of the Gaussian, at a given confidence level.

Background
Thus, the spatial distribution of the wavelengths along η y,J (x i ) is obtained. Subsequently, the wavelet spectrum provides information about the dominant constituent wavelengths of η y,J (x i ). The reader is kindly referred to [5] and [3] for technical details of the application of wavelet transforms on bed forms variability analysis.
Application 2 (power Hovmöller analysis) is an extension of Application 1 [5]. Basically, the Power Hovmöller analysis allows for locating bed form features in η y,J (x i ) that lie in a specific wavelength interval by passing the bed form profile with a continuouswavelet-based high-pass filter and low-pass filter, which are defined by such interval. Thus, by applying that procedure to each bed form profile in η(x i , y i ), the locations of such features are mapped on the bed form field domain [3].
Application 3 (scale-based discrimination) allows for obtaining three hierarchies from η y,J (x i ) by setting two wavelength thresholds, θ 1,3 and θ 3,3 . The smallest scales are comprised in Hierarchy 1, which are determined following the structural criteria presented by [3]. For instance, if Hierarchy 1 corresponds to ripples, θ 1,3 = 0.6 m and if Hierarchy 1 corresponds to small dunes, θ 1,3 = 5 m.
The largest scales are comprised in Hierarchy 3, and are retrieved from the wavelet spectrum according to a threshold θ 3,3 typically set >100 m. Hierarchy 2 is determined by subtracting scales in Hierarchy 1 and Hierarchy 3 from η y,J (x i ). A given situation could present more than three scales of interest, i.e ripples, small dunes, intermediate dunes, large dunes, and bars. In such cases, the described procedure would first discriminate ripples, dunes and bars as hierarchies 1, 2 and 3, respectively. Further, the discrimination procedure can be applied again on the intermediate scales in Hierarchy 2, distinguishing small, intermediate and large dunes.
It is important to note that the methods used for Applications 3 are optimized versions of those proposed by [3] in the following aspects: learning curves to allow for an efficient bias-variance trade-off at each hierarchy are used, and the bias-variance error is distributed among them, and thus the discriminated hierarchies are denoised at a better extent. Likewise, when the bed form data is provided as a water depth field (e.g. h(x i , y j ) in Fig. 1), the application verifies the expected weak correlation between h 1,3 (x i , y i ), namely a Hierarchy 1 field, and the mean water depth (e.g.h x (y J ) in Fig. 1) by using the proxy parameter named scale-variance-ratio (SVR 1,2 ) proposed by [3] representing the ratio between the standard deviation of Hierarchy 1 and that of Hierarchy 2. The reader is kindly referred to [3] and the flow diagrams in the User Manual of Bedforms-ATM for details on the methodology of Applications 1, 2, and 3.
Subsequently, the average 3DI of a resolution cell, which is limited by nodes η(x i , y j ) and η(x i+1 , y j+1 ) (Figs. 1-a and 1-b), is obtained by estimating T b over a spatial window of arbitrary fixed length Lx by Ly being centered in such resolution cell. Finally, contours of T b over the whole bed form field are obtained. It is important to note that the aforementioned process excludes the resolution cells of the bed form field periphery; although, the algorithm minimizes the number of such cells.

Software architecture
Bedforms-ATM applications are built to run in successive sequence. This stems from the fact that some of Application 2 input data is Application 1's output data, both Application 3 and Application 4 use Application 1's output, and Application 4 inherits input  [3] (when j = J, bed form profile η(x i , y J ) ∈ η(x i , y j ) is also represented by vector η y,J (x i ) in this paper). By following the nomenclature proposed by [3], vector-wise Hierarchies 1, 2, and 3 from such bed form profile are η y,J 1,3 (x i ), η y,J 2,3 (x i ), η y,J 3,3 (x i ), respectively. This representation also applies to bed form fields h(x i , y j ) defined by water depth measurements. (c) Curved bed form field which is analyzed by transforming rectangular coordinates into streamwise s and normal n local coordinates, such that all bed form profiles have a constant sinuosity. Thus, resolution cells are trapezoidal, where ∆s j+1 > ∆s j and ∆n = n i,j+1 − n i,j is constant, provided that ∆s j = s i,j − s i+1,j , and thereby, all bed form profiles have the same sample size m. Synthetic bed form data for both rectangular and curved fields are provided as accompanying information of this paper (see Appendix A, Supplementary File 1, Section 2). data from Application 1 and 3 (see Appendix A, Supplementary File 1, Section 1). However, all of these applications are presented in a single MATLAB platform. The input data is collected from the Application 1's interface, and in order to avoid data inconsistencies, it is collected as a set of streamwise equally spaced vectors representing bed form profiles from either linear or curved elevation or water depth bed form fields, provided that all vectors are sampled at the same sampling period, approximately (Bedforms-ATM verifies that not a single sample period is higher than four times the average sample period, and that the number of sample period outliers is less than 5% of the vector length. If these conditions are not met, the program is prevented to run.) These vectors can be described by either relative (i.e., local streamwise coordinate and elevation or water depth) or absolute coordinates (i.e., Easting, Northing, and elevation or water depth) being saved in TXT or MAT file extensions. Additionally, the input data can be provided as a cell array MAT file containing a bed form profile vector in each cell.

Software functionalities
Application 1's interface (Appendix A, Supplementary File 1, Section 1) lists all the steps to create, select input data, preview the input data, perform a run time test (for large data), and run the application. It also highlights the steps that have been successfully completed. We foresee that many potential users of Bedforms-ATM may not be familiar with the wavelet analysis theoretical background; therefore, Application 1's interface lists all the necessary parameters to be set for successfully performing the wavelet analysis of bed forms profiles (advanced users can however edit them to obtain parameter specific estimations.) Thus, it guides the user to the appropriate selection of both the significance level and the frequency increment for the wavelet analysis; likewise, it allows for the use of two sets of wavelet functions, namely: (1) the derivatives of the Gaussian, among which the second derivative, i.e. the Ricker wavelet function or so called Mexican hat wavelet function, has been widely used in the analysis of geophysical signals [3]; and (2) the Morlet function which is most suited for analyzing bed forms variability [3]. If the latter wavelet function is selected, the interface guides the user to the proper selection of the Morlet's nondimensional frequency parameter, k 0 . Additionally, it allows the user to set the output settings (e.g., figures title, signal name, signal unit, sample name, and sample unit). Finally, the user can select the file format (e.g., MATLAB fig, jpg, tiff, or PDF) in which the application's output will be saved. When the run time test indicates that it will take a long period of time, the application suggests the user to print the output to a MAT file and provides a function to print it later, which results in a reduced computational time.
Application 2 performs a scale-averaging of the wavelet power spectra at multiple locations (Fig. 5 in Appendix A, Supplementary File 1), and thus, the spatial variability of a field of bed form transects can be assessed. This application estimates the appropriate wavelength intervals to perform the power Hovmöller analysis, which is presented in its interface, and thereby the analysis of low density scale populations or non-existent scales are avoided. Lastly, Application 2's interface provides the user with two alternatives, namely: (1) to export the output in TecPlot R ⃝ format, or (2) to print the output to file in jpg, tiff or PDF formats. Since the code is open, advanced users can modify it to include other formats such as eps format.
Application 3 presents a plot of the structural classification scheme suggested by [4], and once the input files are selected (i.e., Application 1's output), the application depicts the scales over the aforementioned plot and suggests the user the three scales of analysis, namely: Hierarchy 1, Hierarchy 2, and Hierarchy 3. Additionally, when the whole analysis is performed over a water depth field, this application allows the user for selecting the cross correlation analysis between SVR 1,2 and SVR 2,3 , and the mean water depth at each bed form profile as presented in Fig. 2. This operation represents a physically based assessment of the discrimination output. Finally, the application prints the output to file in both PDF and MAT file extensions. Additionally, it can be exported to a Tecplot R ⃝ input format file. Application 4 allows the user to select: (1) the analysis field, namely: the bed form field η(x i , y j ), or one of the hierarchies being obtained from Application 3; and (2) a sub-domain where the analysis will be performed. The application suggests the window size over which T b is integrated for the user-selected hierarchy as well as bed form field. Additionally, the user can choose printing the output in PDF, MAT or Tecplot R ⃝ file formats.

Implementation
The program along with the corresponding user manual, and test data from the Paraná River (150 bed form profiles) are available free of charge at https://sourceforge.net/projects/bedforms-atm/ and can be easily downloaded and installed. Additionally, the same information can be obtained from https://github.com/rgutierrezl/ BedformsATM. Synthetic data from both curved and rectangular plots to test the program are provided by [10]. Please refer to Appendix A, Supplementary File 1, Section 2 for technical details on the aforementioned data. An explicative video is available at www.pucp.edu.pe/tnbVV8 and can also be found at Appendix A, Supplementary File 2. Bedforms-ATM is divided into four interfaces. To call them, the user should execute the following scripts from the MATLAB console:

Illustrative examples
Bedforms-ATM was tested by using field data from the Paraná River from [11], and synthetic data from [10]. The application of the software on the former data effectively provides information of the spatial distribution of the bed forms wavelength spectrum (Application 1 and Application 2) as shown in Appendix A, Supplementary File 1, Section 3. Fig. 2 presents the output of Application 3 for the Paraná River. Fig. 2-a shows the centered first hierarchy (i.e., small dunes, h 1,3 ), and Fig. 2-b presents the centered second hierarchy (i.e., medium and large dunes, h 2,3 ) which shows less noise when compared to the output presented by [3] as the algorithm distributes the resulting bias-variance error of the method. Similarly, Fig. 2-c shows the third hierarchy (i.e., bars, h 3,3 ) exhibiting a manifestly less noisy output than that presented by [3].
The scale-variance-ratio between the first and second hierarchy (SVR 1,2 , Fig. 2-d) confirms the expected weak correlation with water depth, which is explained by the fact that h 1,3 is mainly induced by local turbulence. On the contrary, SVR 2,3 ( Fig. 2-d) shows a strong correlation with water depth, due to the fact that larger bed form scales are controlled by hydraulic conditions in fluvial environments [3]. Fig. 3-a shows Application 4 output for the Paraná River when the analysis is performed over h(x i , y j ), in which the effective area of the output is smaller than that of h(x i , y j ). Different integration window sizes are used in the analysis, and consequently different degrees of smooth outputs are obtained. This stems on the fact that, similarly to the application of Parzen windows, the bigger the window the smoother the output. Thus, no single window width is ideal overall because some sizes allow to gain information in some areas, whereas losing in others [12]. The user can visually select the output he/she considers the best based on his/her interest as shown in Fig. 3-b, and thus a localized analysis can be performed ( Fig. 3-c).
Three-dimensional analysis on the Paraná River bed forms indicates that high values of T b , namely markedly 3D bed forms, are mostly concentrated at the deeper areas, possibly because three dimensional bed forms control at some extent the Reynolds stresses, the drag coefficients, and the dispersal patterns of sediments [7].

Impact
Even though meteorological models have been usually tested with observation using objective scientifically rigorous methods, many morphological models have been evaluated by means of the modeler's subjective judgment, or by mean squared score (i.e. the Brier skill score) in which the initial bed is chosen as reference [13,14]. In this contribution, the efficiency of Application 3 is analyzed and presented on the basis of objective statistics that are commonly used to evaluate meteorological models. It was successfully applied on both synthetic and fluvial bed forms, [15] reported the effectiveness of Bedforms-ATM to discriminate seabed form features from the Grådyb tidal channel inlet in Denmark, which were previously studied by [16]. Thus, in the light of the results (see Appendix A, Supplementary File 1, Section 4), by coupling the capabilities of wavelet and Hovmöller analysis and 3DI, Bedforms-ATM provides an effective frame to analyze the variability of bed forms features from different sedimentary environments. It was coded following the best practices of object-oriented and userfriendly software, which allows for analyzing bed form data sets in a single platform without the need of passing from one software to another. In the instances that users are not familiar with the wavelet analysis, all the applications guide them to select the appropriate parameters to provide reasonable results. Likewise, the graphical outputs of the software can be exported as vector graphics in publication-quality resolution.
We designed Bedforms-ATM as a modular and straightforward to modify, and thus allowing for the expandability of the software. We expect and encourage new applications to be developed by the users community. In this context, any feedback (e.g., bug reports, suggestions for further development, and pieces of contributed code or even ready plug-ins) benefits the users community [17].
The usability of the software was tested by releasing a beta version in .p MATLAB files, which unlike .m format files, cannot be  edited. An error report form was also provided. Few reports were received, and the present version, corrects the reported errors. Finally, the results presented herein provide valuable information for both scientific and practical applications. For the former application, for instance, when coupled to river morphodynamic models, it could enhance our understanding of the interrelationship between the three hierarchies discriminated herein. For the latter application, it could potentially be used to estimate the input of each hierarchy in the global roughness of a reach or perform a bathymetry reduction which are relevant for flood and navigability modeling.

Conclusions
Most of the computer programs aimed to analyze bed forms in both fluvial and marine environments are not public available or are discontinued. Bedforms Analysis Toolkit for Multidimensional Modeling (Bedforms-ATM) potentially fills this niche, and, in the light of our results from applying the software on fluvial and synthetic bed form fields, we believe it has the potential to be used as an initial standard frame towards bed forms dimensionality analysis and scale-based discrimination, and provide valuable insights on the understanding of bed form systems dynamics and their interrelationship with the surrounding hydrodynamic characteristics. Bedforms-ATM is presented herein as an open source code software whose structure has been designed to encourage its expandability via the collaboration from the users community. The software will greatly benefit from future MATLAB code for parallel computing and GPU which unfortunately are not provided in this version.