The iFlow modelling framework v2.4 A modular idealized process-based model for flow and transport in estuaries

. The iFlow modelling framework is a width-averaged model for the systematic analysis of the water motion and sediment transport processes in estuaries and tidal rivers. The distinctive solution method, a mathematical perturbation method, used in the model allows for identiﬁca-tion of the effect of individual physical processes on the water motion and sediment transport and study of the sensitivity of these processes to model parameters. This distinc-tion between processes provides a unique tool for interpret-ing and explaining hydrodynamic interactions and sediment trapping. iFlow also includes a large number of options to conﬁgure the model geometry and multiple choices of turbulence and salinity models. Additionally, the model contains auxiliary components, including one that facilitates easy and fast sensitivity studies. called simple way, the model complexity

Abstract. The iFlow modelling framework is a widthaveraged model for the systematic analysis of the water motion and sediment transport processes in estuaries and tidal rivers. The distinctive solution method, a mathematical perturbation method, used in the model allows for identification of the effect of individual physical processes on the water motion and sediment transport and study of the sensitivity of these processes to model parameters. This distinction between processes provides a unique tool for interpreting and explaining hydrodynamic interactions and sediment trapping. iFlow also includes a large number of options to configure the model geometry and multiple choices of turbulence and salinity models. Additionally, the model contains auxiliary components, including one that facilitates easy and fast sensitivity studies.
iFlow has a modular structure, which makes it easy to include, exclude or change individual model components, called modules. Depending on the required functionality for the application at hand, modules can be selected to construct anything from very simple quasi-linear models to rather complex models involving multiple non-linear interactions. This way, the model complexity can be adjusted to the application. Once the modules containing the required functionality are selected, the underlying model structure automatically ensures modules are called in the correct order. The model inserts iteration loops over groups of modules that are mutually dependent. iFlow also ensures a smooth coupling of modules using analytical and numerical solution methods. This way the model combines the speed and accuracy of analytical solutions with the versatility of numerical solution methods.
In this paper we present the modular structure, solution method and two examples of the use of iFlow. In the exam-ples we present two case studies, of the Yangtze and Scheldt rivers, demonstrating how iFlow facilitates the analysis of model results, the understanding of the underlying physics and the testing of parameter sensitivity. A comparison of the model results to measurements shows a good qualitative agreement.
iFlow is written in Python and is available as open source code under the LGPL license.

Introduction
The dynamics of estuaries and tidal rivers is characterized by the complex interplay of mutually interacting processes related to the water motion (i.e. tidal propagation, river runoff), salinity and sediment dynamics, transport of nutrients and bathymetric changes. In many estuaries and tidal rivers these processes are subject to constant change due to human interventions, such as dredging and canalization, or to natural changes, such as sea level rise or changing river discharge. These changes may lead to practical problems. Focussing on the hydrodynamics and sediment dynamics, examples are increasing risks of flooding related to tidal amplification or reflection (e.g. Friedrichs and Aubrey, 1994;Winterwerp et al., 2013;Schuttelaars et al., 2013) and deteriorating ecosystems due to a decreased light penetration caused by increasing suspended sediment concentrations (e.g. Colijn, 1982;Cloern, 1996;De Jonge et al., 2014). Many systems face several simultaneous natural and anthropogenic changes, which each affect multiple processes. Therefore the understanding of these processes and their interrelations through models, in combination with observational evidence, is of paramount importance in anticipating the effect of future natural and anthropogenic change.
A wide range of process-based models has contributed to the present-day understanding of flow and transport processes. These models range from linear one-dimensional along-channel models to non-linear three-dimensional numerical models. One way of classifying models is to describe their position in the spectrum ranging from exploratory to complex models (Murray, 2003). On one end of this spectrum, exploratory, or idealized, models typically include a limited number of processes that are thought to be important for the particular phenomenon that is studied. These models come in many forms, ranging from one-dimensional to three-dimensional and from analytic to numeric. The common property of these models is their excellent ability to quickly investigate the sensitivity to parameter variations and to systematically study individual physical processes. Since they are often custom-built, the applied solution techniques do not allow for an easy extension to more processes or complex model domains. Therefore the comparison between these models and real-life systems has to be qualitative and one needs to consider carefully the effect of the underlying assumptions. On the other side of the spectrum, complex models aim at a quantitative comparison of the model results with observations in a wide range of real systems. This requires the implementation of most known processes and their mutual interactions through state-of-the-art parametrizations. As a result, such models are typically numerical and nonlinear, and computation times are relatively long. This makes complex models less suitable for identifying the essential processes and conducting extensive sensitivity studies.
The aim of the iFlow modelling framework is to combine the strengths of both approaches identified above, that is, to represent some of the complex processes and interactions contained in complex models while retaining the ability to analyse these processes and study their sensitivity. iFlow is a width-averaged model for hydrodynamics and sediment transport processes in single-branch estuaries and tidal rivers, focussing on global estuarine processes. Within this context, the model is able to cover a wide range of complexity, reaching out to both the idealized and complex model types. This requires a structured and systematic approach. This approach starts from the exploratory model of Chernetsky et al. (2010), which solves for a specific subset of hydro-and sediment dynamical processes using a combination of analytical and semi-analytical solution methods. The power of iFlow lies in its ability to extend this basic model by adding more complex and realistic interactions, which can either be included in or excluded from the model depending on the application. These extensions can often only be resolved numerically and sometimes require iterative methods. The model thus naturally consists of a set of coupled and mutually interacting components that solve for different processes using different solution methods. These model components are called modules in iFlow. Modules form code-independent entities that can be developed independently and can be easily added to the model without requiring changes to other modules. The iFlow core takes care of the coupling of modules through a simple standardized input/output protocol, thus facilitating interactions between processes in different modules. This allows for a natural development of the model by implementing new processes or different implementations of already existing ones, motivated by the needs for the application at hand.
iFlow currently includes several modules that allow for the computation of the flow and suspended sediment transport. Most of these modules focus on identifying the effect of individual processes and to this end use a perturbation approach. This approach has been successfully applied before in the context of estuarine research by e.g. Ianniello (1977Ianniello ( , 1979, Chernetsky et al. (2010), Cheng et al. (2010) and Wei et al. (2016). The perturbation approach is used to identify processes that balance at different orders of magnitude. Under suitable assumptions of weakly non-linear flow, the leading-order flow and sediment balances reduce to linear equations describing the propagation of the tide and tidal re-suspension of sediment. These balances match classical exploratory model results (e.g. Prandle, 1982;Hansen and Rattray, 1965;Friedrichs and Aubrey, 1994, and references therein). However, non-linear processes and other processes that are not of leading order are not neglected. Rather, linear estimates of the non-linear processes are taken into account at the first and higher orders. Because of the linearity, the effects of each process on the flow and sediment concentration can be evaluated separately. In this way, the fully non-linear solution can theoretically be approximated to any degree of accuracy, while the effects of individual processes and interactions can still be analysed. Practically, it turns out that the qualitative properties of the solution are often well described by only a limited set of orders and processes.
Summarizing, the iFlow philosophy revolves around three central ideas: 1. the model is easily extendible by new processes; 2. the model allows for the combination of different solution methods for different processes, including analytical and numerical solution methods; 3. it is possible to identify the effects of individual physical forcing mechanisms and interactions.
This paper is structured into three main parts. Firstly, Sect. 2 discusses the modular model structure in detail using a basic example involving four modules. This section ends with a list of all modules currently included in the model in Sect. 2.3. This forms the introduction to the second part of this paper, which discusses the specific modules that form iFlow's current functionality in Sects. 3-5. Section 3 presents the model domains and numerical grids currently allowed. Section 4 then provides a discussion of the modules for hy-drodynamics and sediment dynamics, focussing on the assumptions and options in these modules. A short outline of the other main modules, including the various turbulence closures, salinity models and sensitivity modules, is provided in Sect. 5. Section 6 presents two examples of model applications, to the Yangtze and Scheldt rivers. The paper ends with conclusions and a guide to the code availability in Sects. 7 and Code availability. While this paper provides an overview of the model features and methods, an in-depth user manual and a full technical description of the model are provided in the Supplement.

Modular structure
In order to satisfy the three criteria set in the introduction (extendibility, interchangeability and ease of analysis), the structure of iFlow has to be modular. Modules are separate model entities that implement certain physical processes or perform auxiliary tasks, such as plotting or initiating a sensitivity study. A module may use any approach to obtain the required variables, for example solving a set of equations, loading measured or modelled data from a file or even linking to another modelling suite. Modules are code-independent, meaning that the interaction between different modules is only on input and output level, not on code level. This allows an independent development of modules by different developers while ensuring seamless interaction between different modules. It also allows easy interchangeability of modules that compute the same variables but that differ in the physical processes taken into account or the type of implementation used.
Depending on the problem at hand, users can select which variables to save, which physical processes to include and which auxiliary tasks to perform by selecting a set of modules. These modules are listed in an input file, together with the input parameters required by these modules. Upon the start of a simulation, iFlow will read the input file and start an automated two-step process: ordering the modules into a call stack and then calling the modules in this order. Below, these steps are explained and illustrated using the example displayed in Fig. 1, which gives a simplified demonstration of the computation of the leading-order flow velocity (i.e. linear propagation of the tide) through a set of four interacting modules.

Building the call stack
As a first step, iFlow reads the input file (Fig. 1a) and compiles a list of the modules. In order to determine the order in which to call these modules, iFlow needs information on the input required and output returned by each module. This information is documented in a registry file (Fig. 1b), which is provided with the modules and does not need to be given on input. The call stack is made by matching the output pro-vided by each module to the input required by the other modules, such that the required input is available at the moment a module is called.
The input file lists four modules with a specific task each: RegularGrid for making a grid, Geometry2DV for setting the model geometry, Hydrolead for computing the leading-order hydrodynamics and KEFitted as turbulence closure. At the end, the input file lists the variables that are required by the user, e.g. for saving or plotting; here these variables are the leading-order velocity u 0 and eddy viscosity A v (more information on these variables and the underlying equations will be provided in Sects. 4 and 5). The registry file (Fig. 1b) contains the same modules with their input and output variables. Using the registry file, iFlow assesses that the outputs of the Hydrolead module, u 0 , and of KEFitted, A v are needed to obtain the required variables. iFlow then constructs the call stack by determining the modules needed in order to run Hydrolead and KEFitted. Focussing on Hydrolead, it follows from the registry that this module requires nine input variables. These variables may be provided in the input file, by the output of other modules or in a configuration file (not shown here; see the manual for details). Three of these input variables, A 0 , phase 0 and Q 0 , are provided in the input file, while the other six follow from the output of other modules. By matching all the input for and output of the four modules, iFlow constructs the call stack depicted in Fig. 1c.
The call stack shows a loop between Hydrolead and KE-Fitted, which is necessary as both require each other's output as input. This interdependency is resolved by defining KEFitted as an iterative module. Behind the keyword inputInit in the registry of KEFitted it can be seen that this module does not require the flow velocity u 0 , computed by Hydrolead, for its first run. In subsequent runs of the iteration, u 0 is required. iFlow recognizes the interdependency and constructs the smallest possible iteration loop, here involving the two interdependent modules only. The number of iterations follows dynamically from a convergence criterion that is implemented in the KEFitted iterative module.
As a consequence of the way that iFlow constructs the call stack, the model will not use modules that are not needed to compute the required variables. A notification of this is given when running a simulation. Similarly, a notification is given if the call stack cannot be completed, because certain input variables are missing.
The example discussed here can easily be extended, e.g. by adding modules for computing additional variables or by adding auxiliary modules for saving the output or plotting it. To allow for more flexibility, the input and output files allow for a number of additional options that are beyond the scope of this paper, such as submodules and input-dependent output requirements. Details on this are provided in the iFlow manual in the Supplement. files for a model with four modules. The core uses the input and registry files to make a call stack (c) with the correct order of the modules. The output of each module is stored in the data container to be used as input to other modules.

Running and data management
After construction of the call stack, the modules are called sequentially in the determined order. As modules are required to be code-independent, they are not allowed to communicate directly with each other. Instead, the iFlow core regulates the distribution of the required input data and collection of the resulting output. The management of these data is facilitated by the DataContainer in the iFlow core. It collects the module's output upon completion and handles the input data requests by each module; see Fig. 1c. To simplify the interchangeability of modules and the analysis of data, the DataContainer supports various data types and data decompositions, as is discussed more elaborately below.
Different modules used within one simulation can have widely different degrees of complexity and are allowed to use different solution methods. Therefore the requested input and resulting output data can be of different types, including scalars, multi-dimensional arrays and analytical function descriptions. In our example, Geometry2DV sets a constant depth H , which is saved as a scalar value (see also Fig. 1c). Other implementations of the depth allow for depths varying over the horizontal x coordinate according to prescribed analytical functions or data on a grid. This difference in the way the depth is prescribed should not influence the functioning of other modules. The DataContainer allows this by providing a uniform interface to all data types. This means that there is one command for a module to retrieve H (or any other variable) regardless of the underlying data type. The DataContainer handles this command based on the data type. For example, the RegularGrid module requests H on grid points. If H is stored as a scalar, the DataContainer automatically extends this scalar value to all requested points. If H is stored as an analytical function description, this function is evaluated at the grid points. Data stored on numerical grids may as well be used as input to analytical functions. If the numerical data are requested at other coordinates than the grid points, the DataContainer automatically interpolates these data to the requested coordinates. Similarly, a module can access the derivative of a variable. iFlow sees whether an analytical function or numerical data for this derivative is provided and, if not, will automatically perform numerical differentiation.
Since iFlow is designed to improve the understanding of physical processes, modules may offer decompositions of data into contributions resulting from different physical components. The method of decomposition is the responsibility of individual modules. An example of this using the perturbation method will be discussed in Sect. 4 for the hydrodynamics and sediment dynamics modules. Within iFlow's philosophy, it should be possible to interchange these modules with others that do not make decompositions or make decompositions in different components, without affecting other modules. The DataContainer supports this using subvariables. This is illustrated in Fig. 1c for the flow velocity variable u 0 . This has contributions induced by the tide and by the river discharge, such that the sum of both yields the total flow velocity u 0 . The KEFitted turbulence model does not require this decomposition and does not necessarily need to be aware that such a decomposition exists. It can therefore simply request u 0 and iFlow will automatically sum the tide and river contributions. Alternatively a module may request a list of all the sub-variables of u 0 and request each of these contributions separately.
The DataContainer as an interface for different data types and decompositions of data thus ensures that modules with different (e.g. analytical and numerical) solution methods can be used together. Additionally, a module can easily be replaced by a different module that results in the same out-put variables through other processes, without requiring any changes to other modules.

iFlow standard modules
The iFlow modelling framework includes a number of standard modules that may be used to simulate and analyse the water motion and sediment dynamics in estuaries and tidal rivers. Together, the standard modules provide a full model for hydrodynamics and sediment dynamics that may be used in different combinations to model various levels of complexity. The modules are organized into four packages, general, analytical2DV, numerical2DV and semi_analytical2DV, containing auxiliary modules and modules using analytical, numerical or semi-analytical (i.e. largely analytical, with numerical components) solution methods respectively. All included standard modules and the location where they can be found are listed in Table 1.
A short introduction to many of these modules is provided in Sects. 3-5. Section 3 introduces the standard module for geometry and grid. The standard modules for hydrodynamics and sediment dynamics are introduced in Sect. 4. A short explanation of other modules related to salinity, turbulence, reference level and sensitivity analyses is given in Sect. 5.

Model domain and grid
The iFlow core has a flexible definition of the model dimensions that allows for anything from one-dimensional to three-dimensional models. In this paper we will discuss the standard modules in iFlow version 2.4, which are only for a two-dimensional width-averaged (2DV) model. The alongchannel axis is defined as the x coordinate and the vertical axis is defined as the z coordinate. The length of the estuary is thus measured by following the channel between the seaward boundary x = 0 and the landward boundary x = L and can be freely chosen. The width, B, and bed level, H , of the estuary can be provided as arbitrary smooth functions of x; see Fig. 2. The bed level H is relative to the mean sea level at the mouth (MSL) defined at z = 0. iFlow contains several built-in functions describing the depth and width, including polynomial and exponential functions. These functions and their derivatives are computed analytically to obtain maximum accuracy. Alternatively, the depth and width may be provided as a list of numerical data on a grid.
The surface level relative to z = 0 is denoted by R + ζ , where R denotes the reference level and ζ denotes the surface elevation. The reference level R is a quick estimate of the local mean surface level, such that H + R is always positive and is a good approximation of the mean water depth. By default, R = 0, but the use of a non-zero reference level is required if the river bed is above MSL over parts of the domain. A non-zero reference level is also useful when the mean surface elevation above MSL becomes of the same or- Figure 2. Model domain. The model is two-dimensional in the along-channel (x) and vertical (z) directions and is width-averaged. The depth and width are allowed to vary smoothly with x. Evaluate the result of a sensitivity analysis using a cost function that compares model results to data and plot the result numerical2DV RegularGrid Create a 2DV standard grid and output grid.

HydroLead
Leading-order hydrodynamics using fully numerical methods HydroFirst First-order hydrodynamics using fully numerical methods HydroHigher Higher-order hydrodynamics up to any order using fully numerical methods HigherOrderIterator Auxiliary module for higher-order computations (i.e. above first order) ReferenceLevel Computation of a sub-tidal reference level based on the river-induced set-up SedDynamicLead Leading-order sediment dynamics using fully numerical methods SedDynamicFirst First-order sediment dynamics using fully numerical methods SedDynamicSecond Second-order sediment dynamics restricted to river-induced resuspension of sediment, using fully numerical methods StaticAvailability Sediment transport and trapping. Closure module for SedDynamicLead, SedDy-namicFirst and SedDynamicSecond. SalinityLead Dynamic leading-order salinity computation using fully numerical methods SalinityFirst Dynamic first-order salinity computation using fully numerical methods semi_analytical2DV HydroLead Leading-order hydrodynamics. Fully analytical in the vertical direction and numerical in the horizontal direction HydroFirst First-order hydrodynamics. Fully analytical in the vertical direction and numerical in the horizontal direction SedDynamic Leading-, first-and second-order sediment dynamics and transport/trapping using analytical solutions, but with numerical integration. The second-order sediment dynamics is restricted to river-induced resuspension.
analytical2DV Geometry2DV Create a two-dimensional geometry with arbitrary depth and width SaltHyperbolicTangent Diagnostic (i.e. prescribed) well-mixed salinity field according to a tanh function SaltExponential Diagnostic (i.e. prescribed) well-mixed salinity field according to an exponential function TurbulenceUniform Prescribed vertically uniform eddy viscosity and roughness TurbulenceParabolic Prescribed eddy viscosity with a parabolic vertical profile and constant roughness KEFittedLead KEFittedFirst KEFittedHigher KEFittedTruncated Set of modules for a vertically uniform eddy viscosity depending on the local velocity and depth, and for the roughness depending on the local velocity. The dependency between the eddy viscosity and roughness is drawn from relations obtained from a k − ε model.  der of magnitude as the depth. In such cases, the bed level alone is not a good estimate of the mean water depth. More details on the computation of R are provided in Sect. 5.3.
As discussed in Sect. 2, iFlow modules can use a combination of analytical and numerical solution methods. Each of these modules and solution methods may or may not require a numerical grid and grids may serve different purposes. Apart from using grids for (partly) numerical computations, a grid may be used to save or plot variables as numerical data. iFlow allows for using different grids in different modules or omitting a grid altogether. As a result, computations in different modules may use grids with different resolutions and the output may be stored on yet a different grid. Automatic linear interpolation of data between different grids ensures a smooth coupling of modules using different grids. Here, the standard grid module of iFlow, called Reg-ularGrid, is discussed. RegularGrid defines two grids: one computational grid used in all numerical modules and one potentially different output grid. In many cases it is useful to have an output grid with a low resolution to limit the size of the output data, while using a higher resolution computational grid for the benefits of the model accuracy. iFlow grids are curvi-linear and may be non-equidistant in both the x and z direction. More details can be found in the iFlow manuals, attached as a Supplement.

Equations and solution methods for hydrodynamics and sediment dynamics
The standard modules for computing the hydrodynamics and sediment dynamics fit particularly well in the iFlow philosophy as they allow for a separate analysis of the physical contributions to the total result. These analysis properties re-sult from the perturbation approach that is used to solve the continuity, momentum and sediment balances. The steps taken in the perturbation analysis are listed in Fig. 3, which also forms the outline of this section. After presenting the basic width-averaged equations (Sect. 4.1), these are reduced in complexity via a scaling analysis (Sect. 4.2), perturbation approach (Sect. 4.3) and harmonic decomposition (Sect. 4.4). The perturbation approach and harmonic decomposition allow for a particularly good analysis under a set of standard forcing assumptions, which will be discussed in Sect. 4.5. Finally we will discuss the two solution methods (semi-analytical and fully numerical) implemented in the standard modules (Sect. 4.6). Throughout the whole section we will focus on the assumptions made in this procedure and the way in which this approach helps to analyse the model results.

Equations
The water motion is described by the Reynolds-averaged width-averaged shallow water equations that solve for the water level elevation ζ (x, t), horizontal velocity u(x, z, t) and vertical velocity w(x, z, t). Here, t denotes time. We neglect the effects of Coriolis and assume that density variations are small compared to the average density, allowing for the Boussinesq approximation. The resulting momentum equation reads (e.g. Chernetsky et al., 2010) Here, g is the acceleration of gravity, ρ is the density with reference density ρ 0 and the vertical eddy viscosity is de-noted by A v . The subscripts x, z and t in the equations denote derivatives with respect to these dimensions. The background horizontal eddy viscosity A h has been neglected. The momentum equation has a no-stress boundary condition at the free surface and a partial slip condition at the bed The parameter s f denotes the partial slip roughness coefficient. For s f → ∞, the partial slip condition reduces to a noslip condition u = 0. The partial slip law becomes a quadratic bottom friction law if s f is made dependent on the local velocity (see also Sect. 5.1).
In addition we use the width-averaged, depth-integrated continuity equation, which reads Here A = A(t) is the time-dependent tidal forcing at the seaward boundary and Q is the river discharge imposed on the landward boundary. Finally the width-averaged continuity equation reads with a non-permeability condition at the bed The sediment dynamics is described by the widthaveraged sediment mass balance equation, which solves for the sediment concentration c(x, z, t) in the model domain. The sediment is assumed to consist of non-cohesive, fine particles that have a uniform grain size (i.e. constant settling velocity) and are transported primarily as suspended load. At the surface we do not allow for transport of sediment through the water surface and at the bottom we assume that the diffusive flux equals the erosion flux E. The resulting equation is (e.g. Chernetsky et al., 2010) with vertical boundary conditions In Eq. (9), w s is the settling velocity and K h and K v are the horizontal and vertical eddy diffusivity. Boundary condition 11 is valid under the assumption that H x is much smaller than one. We assume that K v is related to the vertical eddy viscosity coefficient A v as K v = A v /σ ρ , where σ ρ is the Prandtl-Schmidt number that converts viscosity to diffusivity. The erosion flux E is related to the so-called reference concentration c through E = w s c . In turn, the reference concentration is defined as where ρ s is the density of sediment, τ b (x, t) = ρ 0 A v u z is the bed shear stress (again assuming H x 1), g = g(ρ s − ρ 0 )/ρ 0 is the reduced gravity, d s is the mean grain size, and a(x) is the availability of easily erodible fine sediment.
The dimensionless sediment availability function a(x) describes how the sediment is distributed over the system. This function is unknown and can be determined by imposing the so-called morphodynamic equilibrium condition (Friedrichs et al., 1998;Huijts et al., 2006;Chernetsky et al., 2010). This condition implies that the total amount of sediment in the estuary varies on a timescale that is much longer than the timescale at which the easily erodible sediment is redistributed over the system. In other words, it is assumed that the amount of sediment in the system is a constant, and we will look for the equilibrium distribution of this sediment in the estuary. Equilibrium in this context means that there is a balance between the tidally averaged erosion and deposition at the bottom or, equivalently, that the tidally averaged transport of sediment is divergence free. The latter is described by the morphodynamic equilibrium condition where uc is the advective sediment transport and K h c x is the diffusive sediment transport. As a boundary condition to this expression it is prescribed that the upstream river carries negligibly little sediment, so that there is no sediment influx from the landward boundary at x = L. The resulting morphodynamic equilibrium condition can be written as (Chernetsky et al., 2010) This implies no net (i.e. tidally averaged) transport of sediment in the entire domain. As the concentration c in Eq. (14) depends on the availability a(x), the above condition is an equation for a. This determines the availability up to a constant factor a * , which should be prescribed on input. This factor determines the total amount of sediment in the system. As the amount of sediment in the system directly affects the concentration in the water column, the absolute magnitude of the concentration may be calibrated directly by changing a * .
The relevant result of the sediment model therefore consists of the relative differences between concentrations at different locations along the estuary instead of the absolute magnitude of the concentration.

Scaling and assumptions
The first step in the perturbation approach is the scaling of the equations. This approach uses a systematic mathematical procedure to determine the relative importance of the different terms in the equations for water motion and sediment dynamics. The most dominant terms will be called leadingorder terms. Terms that are significantly smaller than these leading-order terms will be further categorized according to their relative importance. The most dominant terms, after separating leading-order terms, are called first-order terms. This categorization continues, with all terms of second or higher order generally referred to here as higher-order terms.
The scaling requires four crucial assumptions. Firstly we assume i.e. the ratio of the typical water level amplitude to the depth is much smaller than unity. The small parameter ε is used to define of which order a term is. A term is defined to be of first order if its typical relative magnitude is of order ε compared to the leading-order terms. Similarly, an nth-order term is of order ε n with respect to the leading-order terms. Secondly, it is assumed that the typical tidal wave length and the typical length-scale of bathymetric variations are of the same order of magnitude as the length of tidal influence into the estuary. This implies that sudden local bathymetric variations are not allowed. Rather, bathymetric changes should be smooth over the length of the estuary. Likewise, the method is restricted to long waves, such as tides. Short waves, such as wind waves, are not accounted for. As a consequence of this assumption, the non-linear advection term uu x + wu z in Eq. (1) and uc x + wc z in Eq. (9) scale with ε. It is found that, by these two assumptions, the leading-order equations are all linear, while all non-linearities in the velocity, concentration and water level elevation only appear as first-order or higher-order effects.
Thirdly, it is assumed that the horizontal density gradient is small. More precisely, the internal Froude number should be of order ε or, equivalently, ρ x L tide /ρ 0 should be of order ε 2 , where L tide is the length of tidal influence. As a consequence, the baroclinic pressure term g Finally, the horizontal diffusion term (K h c x ) x is assumed to be of order ε 2 .

Perturbation approach and decomposition
Instead of neglecting first-and higher-order non-linear effects, as is done in conventional linearization techniques, the perturbation approach expands these non-linearities into a series of linear estimates. To this end, the solution variables u, w, ζ an c are written as an asymptotic series ordered in the small parameter ε, i.e.
where [·] 0 denotes a quantity at leading order, [·] 1 denotes a quantity of order ε, [·] 2 order ε 2 , etc. In addition, the eddy viscosity and diffusivity, density, tidal forcing, river discharge and fall velocity are written as similar series. These series are substituted into the equations. The resulting equations are still equivalent to the original system of equations.
The analysis up to this point has merely identified what terms in the equations are of leading and higher orders. The perturbation approach is illustrated here for the momentum and depth-averaged continuity equations for the hydrodynamics, which may be used to compute u and ζ . A first approximation of the equations for the hydrodynamics can be made by neglecting all terms of first and higher orders. The leading-order momentum equation is formulated as with boundary conditions The leading-order depth-averaged continuity equation reads with boundary conditions Compared to the original equations, these leading-order equations omit the non-linear advection, density forcing and all occurrences of ζ in the integration boundaries. These terms feature in the first-and higher-order equations. As a result the leading-order equations have become linear and contain two forcing terms, which are named in the equation: the tidal forcing and river discharge. The linearity is a powerful property, as it allows for applying the principle of superposition. This means that the effect of each of the tidal forcing and river discharge may be evaluated separately and independently and may be summed to obtain the total solution. This is the principle that allows iFlow to make a decomposition of the physics into the responsible forcing mechanisms. An improved approximation of the solution results from constructing the balance of first-order terms. Again focussing on the momentum and depth-averaged continuity equations, these consist of a linear set of equations for u 1 and ζ 1 . The first-order momentum equation is given by with boundary conditions The first-order depth-averaged continuity equation reads with boundary conditions The forcing terms to these first-order equations are defined as the known terms that do not depend on u 1 or ζ 1 and are again marked by a name in the equation. The forcing mechanisms are the first-order tidal forcing A 1 at the entrance, first-order river discharge Q 1 , the density forcing and linear estimates of the non-linearities acting on the flow. These non-linearities include the effects of momentum advection, the tidal return flow and velocity-depth asymmetry. The tidal return flow is the flow that compensates for the mass transport due to correlations between the tidal velocity and surface variation.
The velocity-depth asymmetry accounts for the effect that the velocity profile differs between ebb and flood due to different water levels. Finally, temporal or spatial variations of the leading-order eddy viscosity may be included at first order, so that the interactions between these variations and the leading-order flow appear as a forcing at the first order. Note that some of these mechanisms appear in multiple places in the equations. As the equations are again linear, the principle of superposition allows iFlow to compute the effect of each of these forcing mechanisms separately and independently and sum them to obtain the total result. All forcing mechanisms to the leading-and first-order equations are summarized in Table 2. A similar approach for the sediment balance also results in linear equations at the leading and first order, forced by different physical mechanisms. The leading-order sediment balance describes a local balance between vertical turbulent mixing and the settling of sediment. It is forced at the bed, where sediment is locally resuspended by the leading-order erosion flux E 0 . This erosion rate involves the leading-order bed shear stress τ 0 b , which is derived from the leading-order velocity. The leading-order concentration thus is the concentration locally resuspended by the leading-order tide. The first-order equation describes a similar balance between vertical diffusion and the settling of sediment, but is forced by different components. Firstly, it is forced at the bed by the first-order erosion rate E 1 , which represents the erosion due to the first-order bed shear stress. This involves the firstorder velocity and therefore the flow caused by all mechanisms that act on the first-order hydrodynamics. Secondly, the first-order balance is forced by horizontal sediment advection uc x + wc z , which results in what is known as spatial settling lag effects (Postma, 1954;Van Straaten and Kuenen, 1957;De Swart and Zimmerman, 2009). Thirdly, the firstorder balance involves a forcing from the covariance between the sediment concentration and the surface elevation. Finally, if the eddy diffusivity and fall velocity have first-order contributions, their covariances with the leading-order concentration appear as first-order effects as well. All forcing mechanisms on the leading-and first-order sediment balances are summarized in Table 2. Similar to the hydrodynamics, all the contributions to the sediment concentration by different forcing terms can be evaluated separately and independently due to the principle of superposition. The ordered sediment equations are described in the manuals.
Similar to the approach outlined above for the first-order terms, higher-order approximations of both the hydrodynamics and sediment dynamics can be made by composing a balance of the terms on second, third and higher orders. It is assumed that all external forcing terms (i.e. external tidal forcing, river discharge) act on the leading and first orders. The second and higher orders therefore only contain estimates of non-linear interactions of lower order contributions. The sum of all estimates of the non-linear terms at all orders should return the total solution to the original non-linear system of equations. If the scaling assumptions are satisfied, it follows that the contributions at higher order rapidly become smaller. The solutions at leading and first order then provide a fairly accurate estimate of the total solution. The higher-order systems are nevertheless useful in cases where Table 2. Separate forcing mechanisms to the water and sediment motion and the order at which these mechanisms appear.

Short name Explanation Order
Hydrodynamics Tide Tidal amplitude forced at the seaward boundary 0 and 1 River Constant river discharge at the landward boundary 0 (numerical) or 1 Baroclinic Forcing by the along-channel baroclinic pressure gradient 1 Advection Effect of momentum advection uu x + wu z 1 Tidal return flow The return flow required to compensate for the mass flux induced by tidal correlations between the velocity and water level elevation 1 Eddy viscosity Effect of higher-order eddy viscosity contributions. 1 Velocity-depth asymmetry Correction for the alteration of the velocity profile due to the application of the no-stress boundary condition at z = R instead of the real surface z = R + ζ 1 Sediment dynamics Erosion Local resuspension at the bed 0 and 1 Spatial settling lag Effect of sediment advection uc x + wc z 1 Surface correction Correction because the transport across the time-dependent water surface is specified at z = R instead of the real surface z = R + ζ 1 Fall velocity correction Effect of higher-order variations of the fall velocity 1 Mixing correction Effect of higher-order variations of the eddy diffusivity 1 the scaling assumptions are only marginally satisfied or when studying a particular process that involves a non-linear interaction that appears at higher order.

Harmonic decomposition
The external forcing of the hydrodynamics in iFlow consists of a sub-tidal flow and a limited number of tidal constituents.
In the remainder of this paper we will assume that these tidal constituents are the M 2 tide and its overtides, as these are the most common. In general, one can choose any single tidal base mode and its overtides in the model. The solution to the non-linear system of equations also consists of a sub-tidal component, the M 2 tide and possibly infinitely many overtidal components. As the sediment dynamics is forced by the hydrodynamics, the sediment concentration is described by the same components. This means that the solution can be written as a sum of the sub-tidal component and these tidal constituents. However, instead of accounting for infinitely many components, the signal is truncated after p components, where p can be chosen arbitrarily. As an example, for the velocity u 0 we then write Re û 0 n e niωt , whereû 0 n is the complex amplitude of the nth component of u 0 , where n = 0 denotes the sub-tidal component, n = 1 the M 2 component, n = 2 the M 4 component, etc. A similar decomposition is made for all quantities that vary on the tidal timescale.
As a consequence of this harmonic decomposition, the equations are solved for each frequency component. This eliminates the need to solve the equations by time-stepping. This is a major advantage when computing (dynamic) equilibrium states of the hydrodynamics and sediment concentration, as iFlow can compute these states immediately. This is in contrast to time-stepping models, which often need many time steps and a large computational time to go from an initial state to the equilibrium state.
Details of the equations per frequency component can be found in the manuals. For the case where the leading-order eddy viscosity, eddy diffusivity, partial slip parameter and fall velocity are constant in time, this procedure is the same as in Chernetsky et al. (2010), also see the manual on the semi-analytical model implementation. If these assumptions do not hold, the matrix-solution procedure suggested by Dijkstra (2014) is followed, also see the manual on the numerical model implementation.

Standard forcing
Under certain assumptions about the external forcing, the resulting frequency components of the solutions form an especially well-analysable set. We will call these assumptions the standard forcing assumptions. These are the same as in e.g. Chernetsky et al. (2010) and are the following: 1. the leading-order hydrodynamics is only forced by an M 2 constituent; 2. the first-order is forced only by an M 4 constituent; 3. the river discharge only appears at first order; 4. the eddy viscosity and partial slip parameter do not vary on the tidal timescale; 5. the fall velocity and eddy diffusivity do not vary on the tidal timescale; and 6. the leading-order density variation only contains a subtidal and M 4 component.
Under these assumptions the leading-order hydrodynamics describes the linear propagation of the M 2 tide and only consists of an M 2 frequency. The first-order hydrodynamics consists of a sub-tidal component forced by the river discharge and an M 4 component forced by the external tidal forcing. The density-induced flow and non-linear components appearing at first order are also described by sub-tidal and M 4 components. The first-order flow therefore describes the sources of tidal asymmetry, both caused by external forcing and internal generation.
Assuming the standard forcing assumptions hold, the leading-order sediment dynamics contains the sub-tidal, M 4 , M 8 , etc. components. The first-order sediment dynamics conversely contains the M 2 , M 6 , M 10 , etc. components. In many examples, the leading-order and first-order concentrations are truncated after the M 4 tidal component. This is because the higher harmonics beyond the M 4 component are unimportant for the net transport of sediment and are therefore of less interest.
The main advantage of the standard forcing assumptions is their effect on the morphodynamic equilibrium condition, Eq. (14). This forms a sub-tidal balance of sediment transport terms at second order, which reads We can distinguish between three types of transport terms. The first describes the covariance between the velocity and concentration, i.e. 0 −H uc dz . The dominant covariance terms that result in a sub-tidal transport are R −H u 0 c 1 dz and R −H u 1 c 0 dz . The term u 0 c 1 only generates a sub-tidal transport due to the covariance between the leading-order M 2 flow and M 2 variation of the first-order concentration. The term u 1 c 0 generates transport due to M 4 -M 4 covariance and the product of both sub-tidal contributions. As the model computes the effect of different physical mechanisms contributing to u 1 and c 1 (see Table 2), the transport terms can be subdivided further into the transport caused by particular physical mechanisms. This way, we obtain a subdivision of 0 −H u 1 c 0 dz , with components named after the different contributions to u 1 . Likewise, the components in the subdivision of 0 −H u 0 c 1 dz are named after the contributions to c 1 . One exception to this is the "erosion" contribution to c 1 , which is again subdivided further into the u 1 velocity contributions that cause the erosion.
In addition to these terms, the model includes the sub-tidal transport by R −H u 1 river c 2 river-river dz , i.e. the covariance between the river-induced velocity and the river-induced sediment resuspension. This transport is a fourth-order term according to the scaling and therefore formally does not belong in this balance. However, it typically becomes the dominant term near the end of the tidal influence where all tidally induced transport mechanisms vanish. It is therefore an important mechanism to avoid an unrealistically high degree of sediment trapping at the upstream boundary.
The second type of transport term is the covariance between the velocity, concentration and the varying water surface elevation, with dominant contribution u 0 c 0 ζ 0 . No further subdivision of this term can be made. This term represents the drift of sediment with the moving surface and is largely compensated for by the tidal return flow, which is part of R −H uc dz . Therefore we will consider the transport due to this drift and the tidal return flow together as one term under the name "tidal return flow".
The final type of transport terms are the terms involving the horizontal eddy diffusivity, K h c 0 x and K h c 2 river-river,x . It is assumed that the horizontal diffusivity is constant in time, so that the term K h c 1 x is zero-averaged over the tide. The diffusive transport thus describes horizontal background diffusion of the tide-and river-induced resuspended sediment. Physically, this background diffusion is caused by unresolved flows such as lateral circulation.
Under the standard assumptions, the morphodynamic equilibrium condition thus yields an extensive set of sediment transport terms, which together should sum to zero. By investigating the separate transport terms, it can be inferred which of these mechanisms promote sediment export and which promote sediment import. An example of this is provided in Sect. 6.2.

Semi-analytical vs. numerical solution method
The iFlow hydrodynamics and sediment dynamics modules offer two ways of solving the equations: semi-analytical and numerical. The semi-analytical method follows Chernetsky et al. (2010) and uses fully analytical formulations for the vertical velocity and sediment profiles, but uses a numerical method to solve for the water level elevation. This solution method is fast and accurate, but may only be applied if the forcing satisfies certain conditions. The required conditions are the standard forcing assumptions above, together with the requirement that the eddy viscosity, eddy diffusivity and fall velocity are uniform over the water column.
The numerical method was introduced, because the assumptions on the forcing in the semi-analytical method can be too restrictive for specific applications. The numerical method allows for arbitrary vertical profiles of the eddy viscosity, eddy diffusivity and fall velocity. The numerical Table 3. Allowed forcing and turbulence options in the semi-analytical and numerical solution methods.

Semi-analytical Numerical
Orders hydrodynamics Leading and first Any Orders sediment dynamics Leading and first Leading and first Eddy viscosity/diffusivity Vertically uniform, sub-tidal in leading order and M 2 frequency in first order Vertical variations and leading-order and first-order time variations allowed Bottom boundary condition Partial slip with constant roughness Partial slip with time-varying roughness or no-slip Leading-order forcing tidal components M 2 any First-order forcing tidal components M 4 any River discharge first order leading or first order Fall velocity Vertically uniform, sub-tidal in leading order and none in first order Vertical variations and leading-order and first-order time variations allowed method also allows for releasing the standard forcing assumptions. It allows any number of tidal constituents as long as they are overtides of a base component, often the M 2 tide. These tidal constituents may be imposed at either the leading or the first order depending on the situation. The river flow may additionally be imposed at the leading order, if appropriate. The eddy viscosity, eddy diffusivity, partial slip parameter and fall velocity are also allowed to vary in time at leading or first order. This means that the numerical model may be used with the same restrictions as the semi-analytical method, but these restrictions may be relaxed for further functionality. This is at the cost of potentially larger computational times and lower accuracy, depending on the numerical grid resolution. An overview of the differences between the restrictions in the semi-analytical and numerical methods is provided in Table 3. Some of the additional functionality of the numerical method affects the sediment transport balance. The possible addition of more harmonic components leads to additional transport terms, such as a transport contribution due to the M 6 -M 6 covariance between the velocity and concentration. When a sub-tidal or M 4 velocity is entered at the leadingorder velocity, e.g. through the river discharge or externally prescribed M 4 tide, the covariance between the leading-order velocity and concentration, 0 −H u 0 c 0 dz , yields a sub-tidal contribution. According to the scaling, this contribution dominates over all transport contributions in Eq. (29), so that those contributions should no longer be considered. The term 0 −H u 0 c 0 dz in the new balance can again be subdivided according to the physical mechanisms that contribute to the velocity and concentration. However, the balance now only concerns the leading-order velocity and concentration, for which the model computes only one or two contributions (see Table 2). The subdivision of the transport therefore leads to much fewer terms and typically provides less insight into the underlying physics.
The choice to keep a simulation within the restrictions of the semi-analytical method or to extend it to the full possibil-ities of the numerical method thus has a direct effect on the ability to analyse the results. This is an example of the classical trade-off between model complexity and ability to analyse the results as was mentioned in the introduction. A major strength of iFlow is that it offers one software environment where one can experiment with the degree of complexity required for a simulation for a specific application.
5 Introduction to the modules for turbulence and salinity

Turbulence models
iFlow provides a number of modules to parametrize the eddy viscosity and roughness parameter (see also Table 1), referred to as the turbulence model. The simplest turbulence model available is implemented in the TurbulenceUniform module and assumes a vertically uniform eddy viscosity and constant partial slip roughness parameter, which may only vary with the depth (Friedrichs and Hamrick, 1996;Schramkowski and De Swart, 2002), according to with A v0 , s f,0 , m and n provided as input to the model. The input parameters A v0 and s f,0 may include time-variations (in combination with the numerical hydrodynamics only). This turbulence model is the usually chosen highly simplified turbulence model and was applied in several studies, including Chernetsky et al. (2010), De Jonge et al. (2014 and Wei et al. (2016). However, Schramkowski et al. (2016) showed that multiple values of the calibration parameters A v0 and s f,0 result in equivalent results. Therefore there is a degree of arbitrariness to the calibration parameters in this turbulence model.
In order to resolve this arbitrariness, iFlow includes a set of modules named KEFitted. These models depend only on one calibration parameter and include more physical dependencies of the eddy viscosity. These KEFitted turbulence modules define parametrizations for A v and s f derived by fitting the results of a one-dimensional numerical model with k − ε closure for a large number of barotropic tidal model configurations. The turbulence closures provide a number of options. The most important option is the choice of roughness parameter to provide on input. If the roughness parameter s f,0 is provided, the turbulence model uses the relation This model only has the calibration parameter s f,0 and requires a choice for n. It thus eliminates the need to calibrate A v0 and m. To leading order, because it is assumed that ζ H + R, this model is the same as Eqs. (30)-(31) with m = 1+n and A v0 = 0.5s f (H (x = 0)) m . This model is recommended over Eqs. (30)-(31), as it only has a single calibration parameter and thus leads to a definite best calibration parameter setting. However, note that this relation is derived for a unidirectional flow and its is assumed that any flow in another direction does not affect this relation.
Alternatively, the KEFitted turbulence models may be provided with a roughness parameter z * 00 . The formulations for the eddy viscosity and partial slip roughness then read as where u * is the bed friction velocity, which may be related to the depth-averaged velocity (see Burchard et al., 2011). The parameters z * 00 and n should be provided as input, and α, β, f 1 (z * 0 ) and f 2 (z * 0 ) are known hard-coded parameters and functions obtained by fitting results of the k − ε model (see the manual for details). This model therefore also contains only one calibration parameter z * 00 and requires a choice for n. These formulations relate the vertically uniform eddy viscosity and partial slip parameter to the local bed shear stress velocity and water depth. As a result, the bottom boundary condition for the hydrodynamic Eq. (3) has become a quadratic friction law. This model introduces non-linearity, as there now is a mutual relation between the flow velocity and the water surface elevation on the one hand and eddy viscosity and the partial slip parameter on the other hand. This non-linearity is resolved by an iteration loop over the turbulence and hydrodynamic modules, which is automatically constructed by the iFlow core as exemplified in Sect. 2. Due to the non-linearity, this model introduces more complexity compared to the previous models and is therefore only recommended when the case at hand requires this complexity, for example because of large variations in u * in space or time. An example of this is given in Sect. 6.1.
iFlow implements four modules that implement the above KEFitted relations. The KEFittedLead, KEFittedFirst and KEFittedHigher modules make an ordering of the above equations to determine the leading-order, first-order and higher-order eddy viscosity and partial slip parameter. The KEFittedTruncated module uses the sum of all computed orders of the velocity and water surface elevation to compute a total eddy viscosity and roughness parameter without ordering (i.e. a truncation method).
Finally, the TurbulenceParabolic turbulence model is similar to Eqs. (30)-(31), but assumes the eddy viscosity to have a parabolic profile in the vertical direction. This turbulence model assumes s f → ∞, so that the bottom boundary condition for the hydrodynamics reduces to a no-slip law. The roughness is instead described by a roughness height z 0 . The formulations for A v and z 0 read as The parameters A v0 and z * 0 = z 0 (x = 0)/H (x = 0), m and n are provided as input, z * = z/(H + R) and the dimensionless surface roughness z * s is determined by the model such that A v equals 10 −6 m 2 s −1 at the surface, i.e. approximately the molecular viscosity. The parabolic eddy viscosity profile represents a more realistic shape in barotropic flows and therefore results in more realistically shaped velocity profiles. However, this model faces a similar degree of arbitrariness in the choice of A v0 and z * 0 as in Eqs. (30)-(31) and may only be used in combination with the numerical solution method.

Salinity
The iFlow standard modules include two types of salinity models: diagnostic (i.e. prescribed) and prognostic (i.e. resolved). The diagnostic modules prescribe a sub-tidal vertically uniform (well-mixed) salinity that varies in the alongchannel direction. The SalinityHyperbolicTangent module formulates this as (see also Warner et al., 2005;Talke et al., 2009) and SalinityExponential formulates this as The prognostic salinity model (modules SalinityLead, SalinityFirst) follows work done by McCarthy (1993) and Wei et al. (2016). The model is based on the perturbation approach, where it is assumed that the leading-order salinity consists of a sub-tidal vertically uniform (well-mixed) salinity. Vertical and temporal variations of the salinity appear at higher orders. For more information we refer to Wei et al. (2016).

Reference level
The hydrodynamic module relies on the water depth being positive and much larger than the time varying surface elevation (see assumption 1 in Sect. 4.2). The model fails or becomes inaccurate if the bottom lies above or close to MSL. In many cases this problem can be resolved by the iFlow ReferenceLevel module. This module computes a quick estimate of the sub-tidal water level elevation based on the riverinduced set-up. This is often sufficient, because the river is often the dominant flow term in the most upstream reach, where the bottom level is highest.
The river-induced set-up is estimated numerically using the leading-order momentum and depth-averaged continuity equations, assuming it is purely forced by a constant discharge Q and the resulting water level elevation is given by R. These equations read as This system is non-linear in R as the integral in the second equation contains R in the integration boundary and u, which depends on R according to the first equation. Nevertheless, the system can be solved without iterating by starting at the mouth and working upstream. At the mouth (x = 0), R = 0 by definition. Therefore R x can be computed from the above system of equations. The value of R at the next grid point x = x follows from a simple first-order routine: R( x) = R(0) + R x (0) x. The total reference level follows by repeating this procedure for all horizontal grid cells. More accurate computations of the river-induced set-up follow from the hydrodynamic modules, so that the relatively low numerical accuracy of the reference level computation will not reduce the precision of the overall result. The reference level still depends on the eddy viscosity. If a KEFitted turbulence model is used, the eddy viscosity in turn depends on the reference level. To resolve this interdependency efficiently, without needing to iterate between the turbulence model and reference level module, the KEFitted turbulence models have a built-in routine to compute the reference level. Therefore, the ReferenceLevel module can be omitted when the KEFitted module is used.

Sensitivity analysis module
iFlow's standard sensitivity analysis module Sensitivity provides a powerful analysis tool, by easily allowing a user perform a full model simulation for various values of one or more input variables. On input, the user provides the names of the variables to loop over, as well as a list with the values for these variables. A final input parameter indicates whether all combinations of parameter values should be tested or whether the values of all variables should be changed simultaneously. The iFlow core then automatically decides which modules should be included in the loop and runs these modules for all prescribed parameter settings, saving the results to a file after each loop. The sensitivity analysis is therefore a general tool that may be combined with any set of modules to loop over any set of variables and values. An example of the use of the sensitivity module is given in the model evaluation in Sect. 6.2.

Model evaluation
The use of a 2DV perturbation approach for hydrodynamics or sediment dynamics similar to iFlow's semi-analytical method has been demonstrated before by e.g. Chernetsky et al. (2010) and Wei et al. (2016). An application of iFlow itself has been presented before by Dijkstra et al. (2017) to identify the exchange flow caused by eddy viscosity-shear covariance (ESCO).
Here, iFlow is applied to two case studies. The aim of these cases is to show the application of iFlow, demonstrate ways it can be used to analyse the results and qualitatively compare the model results to measurements. While this aim requires discussing some of the physical mechanisms observed in the model, these physical mechanisms are not the focus of this section. The first case study is an assessment of hydrodynamic effects of the river discharge on the tidal propagation in the Yangtze River, China. This case demonstrates some of the advanced hydrodynamic settings in iFlow, including the use of the reference level, leading-order river discharge and velocity-dependent eddy viscosity. Due to the inclusion of the leading-order river discharge, this can only be done using the numerical implementation of the modules in iFlow. The calibration of this model is also demonstrated. The second case study presents an assessment of the estuarine turbidity maximum (ETM) in the Scheldt River estuary. Here, we will use the standard assumptions on the forcing (see Sect. 4.5) and vertically uniform profiles for the eddy viscosity, eddy diffusivity and fall velocity, so that we can apply the semianalytical solution method. This method is preferred over the numerical method as it is faster and more accurate. This case also features a demonstration of the sensitivity analysis module of iFlow.
As a result of iFlow's flexible modular structure, the modules used are different from application to application. The modules used in the two applications presented below are shown in Fig. 4. Both cases use the modules for generating the model domain and grid. The first case then makes use of a module for the reference level and velocity-dependent eddy viscosity and partial slip parameter. As the latter module requires the leading-order velocity, it automatically iterates over the leading-order hydrodynamics module until the result has converged. The calibration routine calibrates on both the leading-order and first-order hydrodynamics using the roughness parameter in the turbulence model. It therefore automatically constructs an iteration loop over these modules. The second case is a linear sequence of modules without any need for iteration loops. Only the sensitivity analysis initiates a loop. This loop is kept as small as possible, so that a loop over the discharge and the externally prescribed M 4 tidal phase only requires a loop over the first-order hydrodynamics and sediment dynamics. long from Wusong to Datong. In order to ensure that the tidal wave damps out, the model domain is extended to 1500 km, of which only the first 560 km are analysed. Measurements of the width-averaged bed level and nearsurface width are provided by Guo et al. (2014). The bed level is characterized by large variations caused by local width variations and river bends. Smoothing this profile, the bed level is well characterized by a horizontal bed with a depth of 10 m. The width is strongly converging, from 25 km at the mouth to a fairly constant 3 km between 200 and 500 km. This width profile is approximated by the exponent of a rational function given by B = 1000 · exp 3.8 × 10 −5 x + 10 8.8 × 10 −11 x 2 + 2.5 × 10 −5 x + 3.2 . We will distinguish between two forcing conditions: wet and dry season conditions. For both we will assume average tidal conditions, for which the primary forcing components are a leading-order M 2 tide with amplitude 1.09 m and a firstorder M 4 tide with amplitude 0.22 m and a phase difference of −44 • (Guo et al., 2016). We assume a representative discharge of 50 000 m 3 s −1 for the wet season and 15 000 m 3 s −1 for the dry season. In both conditions, the river is assumed to force the water motion at leading order. The effects of salinity or sediment on the flow are not considered for simplicity. We will only consider the leading-order eddy viscosity computed by the turbulence module KEFittedLead with roughness parameter z * 00 (see Eqs. 34-36). The leading-order eddy viscosity is assumed uniform in the vertical and the leadingorder eddy viscosity and partial slip parameter are assumed constant in time and dependent on the leading-order velocity. The eddy viscosity and partial slip roughness parameter are therefore a function of the leading-order M 2 tide and the river discharge.
The model is calibrated by adjusting the roughness parameter z * 00 , such that the computed water levels match the observed water levels for the wet season. The model is calibrated through the sensitivity analysis module. This module constructs a loop over the hydrodynamic modules for a range of different z * 00 values (see also Fig. 4a). The results of each computation are compared to the measurements by using the cost function introduced by Jones and Davies (1996). The result is plotted in Fig. 6, which shows the value of the cost function for the M 2 tide and M 4 tide as a function of z * 00 . The actual value of the cost function is not displayed, since there is no interpretation to this value. The best fit to the measurement is found for the smallest cost, which is for z * 00 = 9.6 × 10 −5 for the M 2 tide and z * 00 = 1.3 × 10 −4 for the M 4 tide. Only one value for z * 00 can be chosen. Since these values are close, we proceed with a rounded value of z * 00 = 1 × 10 −4 . The same roughness value is used for the dry season case.

Results
The resulting water level amplitude and phase are plotted in Fig. 7. The lines show the model results for the sub-tidal flow, M 2 and M 4 tides in the wet season (dashed line) and dry season (solid line). The dots and crosses indicate measurements presented by Guo et al. (2016) for the dry and wet seasons respectively. We find a good correspondence between the measured tidal water level amplitude and phase. This is even true for the dry case, for which the model has not been recalibrated. We additionally find a good correspondence between the measured and modelled sub-tidal water levels, even though the model has also not been calibrated for these. Most importantly, the model captures the correct trends between the wet and dry season, such as increased tidal damping of the M 2 and M 4 tide in the wet season and a slower propagation (i.e. larger phase differences) of the M 2 tide in the wet season.
We will demonstrate how the model can be used to uncover the main processes that cause the differences in tidal propagation between the dry and wet seasons. The non-linear first-order processes that involve tide-river interactions are advection, tidal return flow and velocity-depth asymmetry. Through the perturbation method, these processes can be inferred directly from the model results. The correction to the M 2 tidal amplitude from each of these terms is only of the order of some centimetres in both the dry and wet seasons. The higher-order non-linear terms (not shown here) pose even smaller corrections to the M 2 tide.
The model is rerun without the reference level module (i.e. R = 0) to see the effect of the difference in reference depth between the dry and wet seasons. The resulting water level amplitude is plotted in Fig. 8. The first striking observation is that the sub-tidal water levels are now unrealistically high in the wet season (> 100 m). Since the river-induced water level set-up is large, the assumption that ε 1 (i.e. the ratio of the water level and depth is small) is violated, leading to an unreliable computation. Nevertheless, the M 2 and M 4 tides are hardly affected by this and show the same characteristics of the tide-river interactions as before. We can thus conclude that the reference level is an essential model feature in model cases with a large river-induced set-up, but does not seem to be essential in tide-river interactions.
The effect of the river flow on the eddy viscosity and partial slip parameter can be assessed by switching off the dependency of these quantities on the river flow. The KEFitted turbulence module provides an option to switch off any physical mechanism that can be separated explicitly from the solutions. The river flow can therefore be disregarded in the computation of the eddy viscosity and partial slip parameter, while the dependence on the M 2 tidal velocity is still accounted for. The resulting water level amplitude is plotted in Fig. 9. The M 2 and M 4 tides now do not experience a sufficient degree of damping to vanish at the 560 km point. Also, the differences in the M 2 tidal amplitude between the dry and wet seasons have nearly vanished. The main effect of the river run-off on the tidal amplitude is thus through the effect the river flow has on the eddy viscosity and partial slip parameter.
6.2 ETM location in the Scheldt River

Model settings
The tidal Scheldt River, situated in the south-west of the Netherlands and north-west of Belgium, runs from its mouth in the North Sea to a set of locks and sluices near Ghent, 160 km upstream; see Fig. 10. The river serves as the main shipping channel to the port of Antwerp, which is located approximately 75 km upstream from the mouth. Dredging activities for maintaining and deepening the shipping channel are one of the main reasons to study the development of the fine sediment dynamics in the Scheldt River. The estuary is over 6 km wide and averages a depth of 15 m at its mouth. The estuary converges to a width of about 50 m and an av-S T Figure 9. Water level amplitude for the Yangtze case when the effect of the river discharge on the eddy viscosity and partial slip parameter is omitted. The figure shows results for a dry season situation (solid line and dots) and a wet season situation (dashed line and crosses); see Fig. 7 for more explanation. The M 2 tide is damped less and there is a smaller difference between the M 2 tidal amplitude in the wet and dry seasons. erage depth of about 3 m at the end of the tidal influence. To obtain a schematized depth and width profile, a polynomial is fitted through the geometry data of 2013 (Coen et al., 2015). The depth profile is approximated by a fifth-order polynomial given by H (x) = −2.9 × 10 −24 x 5 + 1.4 × 10 −18 x 4 − 2.4 × 10 −13 x 3 + 1.7 × 10 −8 x 2 − 5.2 × 10 −4 x + 15.3, and the width profile is approximated by an exponent of a rational function given by The model is forced by a leading-order M 2 tidal amplitude of 1.77 m and a first-order M 4 tidal amplitude of 0.14 m, which is −1.3 • out of phase with the M 2 tide. The eddy viscosity A v is computed using the KEFittedLead module using s f,0 as the input parameter and using n = 0 (see Eqs. 32-33). Therefore the partial slip parameter is constant in space and time and the eddy viscosity is assumed to be uniform over the vertical, linearly dependent on the depth and constant in time. The salt water influence typically reaches up to the port of Antwerp (i.e. 75 km). It is assumed that the salinity is well mixed and can be described by a prescribed horizontal salinity profile, which is obtained by fitting a tangent hyperbolic function to summer and winter measurements and taking the mean as the representative profile (see Warner et al., 2005;Talke et al., 2009;Schramkowski et al., 2015). This results in the following expression for the salinity profile: The river discharge varies over the seasons, with an average summer discharge of 30 m 3 s −1 and an average winter discharge of 80 m 3 s −1 . Sediment concentrations found in the Scheldt are moderate, with near-surface concentrations only occasionally and locally exceeding 200 mg L −1 . Based on yearly averaged data of the suspended matter concentration , the ETM is typically found between 100 and 140 km upstream from the mouth. However, monthly averaged data indicate that for moderately high discharges, the ETM can be found around 60-70 km upstream from the mouth or even disappear entirely.
The model for the Scheldt Estuary is calibrated by tuning the partial slip roughness parameter s f 0 to measured water level data. The calibration procedure is similar to that for the Yangtze River, but is not shown here (for details, see Dijkstra et al., 2016). A best fit was found for s f 0 = 0.0048 m s −1 , which results in a good agreement with the M 2 water level and phase and the M 4 phase, but leads to an overestimation of the M 4 water level.

Results
Using the above parameter values and settings we compute the tidally averaged sediment concentration for average summer and winter conditions; see Fig. 11. Since we are interested only in the location and relative magnitude of the ETM, the concentration is scaled by the maximum concentration in the domain. It follows that, for average summer conditions, two ETMs are present: a strong ETM near the weir at approximately 150 km and a weaker one at approximately 120 km. For average winter conditions the ETM is pushed in the downstream direction to approximately 100 km from the mouth. These results are in qualitative agreement with observations and thus suggest that the model captures the most Figure 11. Tidally averaged sediment concentration for summer and winter conditions. Parameter values corresponding to the Scheldt Estuary. important physical mechanisms underlying ETM dynamics in the Scheldt Estuary.
In order to further investigate the underlying physical mechanisms of the ETM dynamics in the Scheldt Estuary, we look closer at the individual processes that contribute to the sediment transport. As explained in Sect. 4.6, iFlow allows investigating the transport contribution due to the individual contributions to the sediment concentration and the flow velocity. Five of the, in this case, most important transport contributions are shown in Fig. 12 for average summer and winter conditions. The individual transport terms are scaled with the maximum value. Note that the total transport, obtained by adding all contributions, equals zero by definition of the morphodynamic equilibrium. For both summer and winter conditions the main up-estuary transport is due to the internally generated asymmetries of the velocity and depth and tidal return flow. During winter conditions, the spatial settling lag (i.e. sediment advection) is important as well. The down-estuary transport is mainly due to the river flow. The external M 4 tide additionally promotes export of sediment in winter conditions.
To illustrate iFlow's capacity to easily perform an extensive sensitivity analysis, we further analyse the influence of the external M 4 tidal component on the ETM dynamics in winter conditions. It is not likely that this component changes on a short time-scale and we select it here purely for illustrative purposes. The sensitivity study comprises 361 different values of the external M 4 tidal phase ranging between −180 and 180 • . The results of all individual simulations are postprocessed and the results are shown in Fig. 13. The result shows that the ETM can shift between approximately 70 and 110 km from the mouth depending on the phase of the external M 4 tide. The M 4 tidal transport induces maximum sediment export at a phase of approximately 50 • , whereas it induces minimum export for a phase of approximately −100 • . For phases between −60 and 140 • , the model also indicates the existence of a concentration minimum. Such a minimum is characterized by a divergence of the sediment transport.

Conclusions
We have demonstrated that iFlow provides a flexible and versatile modular environment for modelling flows and sedi- ment transport in estuaries and tidal rivers. The model focusses on idealized approaches that allow the systematic analysis of physical processes and the sensitivity of these processes to model parameters. Due to the modular nature, iFlow offers a software environment where users can easily adjust the processes included in a simulation, thereby allowing users to adjust the degree of complexity, computational time and ability to analyse the results to a specific application. The iFlow core supports these adjustments by automatically taking care of the communication between modules, order of modules and smooth coupling of modules that use different solution methods. iFlow version 2.4 additionally includes a number of standard modules especially designed to analyse individual processes affecting the flow and sediment transport.
This has been illustrated in two applications of iFlow to the study of non-linear hydrodynamic interactions in the Yangtze River and sediment trapping mechanisms in the Scheldt River. By a systematic approach of switching particular processes on and off and by the decomposition of the results according the underlying physical processes, the model allows for a unique insight into the physics of these systems. As iFlow allows for a quick set-up and calibration of a model and quick sensitivity study, the model is especially well suited to gain a first insight into the essential processes of a system and response of the system to changes. The comparison of the model results with observations in these systems should be mainly interpreted qualitatively, focussing on the relative importance of processes, behaviour and sensitivity. Nevertheless, in the shown applications, there is a good quantitative correspondence between the model result and observations considering the degree of schematization in the model.
Both case studies used different modules and interactions, so that the model complexity suits the analysis relevant to the application. This extendibility, interchangeability and ease of analysis form the main ideas of iFlow. These ideas are reflected in the architecture of the modular set-up, data man-agement and the modules offered within this version of the model. Wrapped around this is a set of tools and modules to support input, output and visualization of the results to make the model user-friendly.
As the structure of iFlow can be adapted and modules can be added easily by new users, there is no such thing as a single iFlow model. Also, the provided default modules for hydrodynamics, turbulence and sediment dynamics may be replaced if this is useful for a particular application. For example, these modules may be replaced by a coupling to a complex model (e.g. as demonstrated for turbulence by Dijkstra et al., 2017) or observations. By coupling such module replacements to other modules one can construct unique model set-ups for studying a certain process or for comparing different model implementations within one modelling framework.
The future ambitions for the model involve further developments of modules for turbulence and morphology and for the transport of sediment, salinity and nutrients. Users are encouraged to contribute to this development by developing and sharing modules or sharing model applications.
Code availability. When using iFlow in any scientific publication, technical report or otherwise formal writing, authors are strongly requested to cite this paper and mention the name iFlow.
The iFlow code is property of the Flemish Dutch Scheldt Committee (VNSC) and is licensed under LGPL (GNU Lesser General Public License). In summary, this means that the code is open source and may be used freely for non-commercial and commercial purposes. Any alterations to the existing iFlow source code (core and modules) must be licensed under LGPL as well and are therefore open source. However, new modules or a coupling between iFlow and other software may be published under a different licence. Nevertheless, users of iFlow are encouraged to make their own developed modules and model applications open source as well.