Analysis of the drainage density of experimental and modelled tidal networks

Based on controlled laboratory experiments, we numerically simulate the initiation and long-term evolution of back-barrier tidal networks in micro-tidal and meso-tidal conditions. The simulated pattern formation is comparable to the morphological growth observed in the laboratory, which is characterised by relatively rapid initiation and slower adjustment towards an equilibrium state. The simulated velocity field is in agreement with natural reference systems such as the micro-tidal Venice Lagoon and the meso-tidal Wadden Sea. Special attention is given to the concept of drainage density, which is measured on the basis of the exceedance probability distribution of the unchannelled flow lengths. Model results indicate that the exceedance probability distribution is characterised by an approximately exponential trend, similar to the results of laboratory experiments and observations in natural systems. The drainage density increases greatly during the initial phase of tidal network development, while it slows down when the system approaches equilibrium. Due to the larger tidal prism, the tidal basin has a larger drainage density for the meso-tidal condition (after the same amount of time) than the micro-tidal case. In both micro-tidal and meso-tidal simulations, it is found that there is an initial rapid increase of the tidal prism which soon reaches a relatively steady value (after approximately 40 yr), while the drainage density adjusts more slowly. In agreement with the laboratory experiments, the initial bottom perturbations play an important role in determining the morphological development and hence the exceedance probability distribution of the unchannelled flow lengths. Overall, our study indicates an agreement of the geometric characteristics between the numerical and experimental tidal networks.


Introduction
Channel networks (e.g.mountainous, fluvial or tidal) exhibit a variety of morphologies and their shape and functioning has attracted considerable attention in the research sphere (Montgomery et al., 1997;Rinaldo et al., 1998Rinaldo et al., , 1999a;;D'Alpaos et al., 2005).A comprehensive understanding of channel networks' morphologies as well as their long-term evolutions is fundamental to address their response to environmental variations under increasing climate change and human pressure (French, 2006;D'Alpaos, 2011;Coco et al., 2013).Branching tidal networks are typical landscapes in back-barrier tidal basins such as the Wadden Sea (Wang et al., 2012) and the Venice Lagoon (D'Alpaos et al., 2007a).These systems consist of various geomorphological elements (Fig. 1) that are shaped by a variety of processes.Ebb deltas and barrier islands act as natural shelters for flood deltas and tidal flats (and associated salt marshes).Tidal flats are often dissected by dendritic channel networks that serve as effective nutrient pathways for flora and fauna, and provide valuable accommodating space for coastal ecosystems (Allen, 2000;Townend et al., 2011).
Different methodologies (e.g.Montgomery and Dietrich 1988;Fagherazzi et al., 1999;Passalacqua et al., 2010;Jiménez et al., 2013a) have been developed in order to ex-tract the structure of channel networks and gain insight into their geomorphic characteristics.Due to the high variability in morphologies of channel networks, it is in fact convenient to study these systems using a statistical approach (Marani et al., 2003;Passalacqua et al., 2013).Remote sensing, associated with these specifically developed techniques to analyse digital terrain maps (DTM), has enabled a more accurate measure and understanding of the geomorphic characteristics of tidal networks (Fagherazzi et al., 1999;Rinaldo et al., 1999a, b;Mason et al., 2006;Vandenbruwaene et al, 2012).
Analyses of these field data indicate that scale invariance, as widely observed in river networks, does not necessarily hold for tidal networks (Rinaldo et al., 1999a) because of the more complicated interplay between various competing landscape-forming processes operating in estuarine systems (e.g.bi-directional tidal flow and the presence of locally generated wind waves).Recently, Passalacqua et al. (2013) analysed the channel network of the Ganges-Brahmaputra-Jamuna delta using a statistical framework based on several descriptors (island area, shape factor, aspect ratio and nearest-edge distance), and their results indicated that the statistical behaviour of these descriptors could distinguish the coastal region and the area along the main rivers; hence geomorphic signatures of processes and vegetation types responsible for delta formation and evolution could be identified.
Drainage density, a measure of channelisation in watersheds, has been found to be of fundamental importance for the understanding of channel network systems since it reflects a balance between climate, geomorphology and hydrology (Moglen et al., 1998).This quantity was first defined by Horton (1932Horton ( , 1945) ) as the ratio between the total channelised stream length and the watershed area: D = L c /A.Although widely adopted, the "laws of drainage networks" developed by Horton (1945) are still under debate.Analysing unbiased samples of networks generated with a Monte Carlo method, Kirchner (1993) found that almost all possible networks follow Horton's laws, indicating the distinctiveness of an individual channel network actually cannot be revealed by Horton's theory (see also Rinaldo et al., 1998).This was also confirmed by Marani et al. (2003) through analysis of both natural networks in the Venice Lagoon and highly schematic test settings.In order to overcome the drawback embedded in Horton's laws, Marani et al. (2003) proposed another measure based on the statistical distribution of the unchannelled flow length (l), i.e. the mean distance from a point on the lagoon platform to the nearest channel.This methodology, based on a Poisson-type hydrodynamic model (Rinaldo et al., 1999a), has been used in many studies afterwards (D'Alpaos et al., 2005(D'Alpaos et al., , 2007a, b;, b;Feola et al., 2005), proving to be a useful geomorphic tool for analysing channel network systems.
Field observations, together with different geometric analyses (e.g.Rinaldo et al., 1999a;Marani et al., 2003;Novakowski et al., 2004;Vandenbruwaene et al., 2012a, b), have confirmed that tidal networks show distinctive landscapes and geometric characteristics, because of the highly complicated interactions between various processes such as waves, tides, biological activities and sea level rise.By limiting the number of active processes, controlled laboratory experiments provide an important alternative to gain insight (Stefanon et al., 2010(Stefanon et al., , 2012;;Vlaswinkel and Cantelli, 2011;Kleinhans et al., 2012;Iwasaki et al., 2013).Moreover, the data obtained from these experiments (or after proper scaling analyses) can be utilised to benchmark numerical models (Tambroni et al., 2010).Stefanon et al. (2010Stefanon et al. ( , 2012) ) successfully reproduced the major features of the evolution of tidal networks in the laboratory, which are comparable to natural systems in terms of both the morphological development and the geomorphic characteristics.Zhou et al. (2014) performed a comparative study between these laboratory experiments and numerical simulations at both laboratory and natural scales, exploring the effects of different parameterisations on the long-term morphological evolution of tidal networks.They found that eddy viscosity, sediment transport formulation and bed slope terms played an important role in the resulting morphologies of tidal networks.In this study, emphasis is given to the comparison of the general geomorphic characteristics of tidal  Specifically, we will explore whether the numerically modelled channel networks are comparable with the experimental (Stefanon et al., 2010(Stefanon et al., , 2012) ) or natural ones from the perspective of drainage density.The variation of drainage density during the morphological development will also be examined in order to shed light on tidal network ontogeny.The remaining of the paper is organised as follows: Sect. 2 describes the experimental and numerical models, Sects.3 and 4 present the results and associated discussions, and concluding remarks are made in Sect. 5.

Methods
2.1 Physical model and scaling analysis Stefanon et al. (2010Stefanon et al. ( , 2012) ) carried out experiments under a range of different settings.One of these experiments (indicated by "Run 4" in Stefanon et al., 2010) was conducted considering a constant mean water level.Run 4 is the focus of this study.The experimental domain included an initially flat lagoon, a gently sloping shelf and a deep sea (Fig. 2a).
The lagoon was connected with the shelf and the sea via an inlet which was bounded by non-erodible barriers.Initially, both the lagoon and shelf were filled with coarse, low-density non-cohesive grains (medium size D 50 = 0.8 mm and density ρ s = 1041 kg m −3 ), providing an erodible sediment layer thick enough to prevent the erosion down to the non-erodible concrete bottom.A harmonic tide, with a tidal amplitude a = 0.005 m and period T = 480 s, was simulated at the sea by a vertically oscillating weir.
The final stable network of Run 4, observed in the laboratory after approximately 8000 tidal cycles, is shown in Fig. 2b.Through scaling analysis (Stefanon et al., 2010), it was possible to "transfer" this experiment to typical natural systems where tidal networks are often observed, i.e. micro-tidal and meso-tidal basins like the Venice Lagoon and the Wadden Sea (see Table 1 for the results of the scaling analysis).The reader is referred to Stefanon et al. (2010) for a detailed description of the experiments and the scaling analysis.

Numerical model and setup
The numerical model Delft3D is used to simulate the "morphodynamic loop", which in this study consists of the depthaveraged shallow water equations, sediment transport (van Rijn, 1993) and bed level updating equations.At each hydrodynamic time step, the morphological change is calculated feeding back on the hydrodynamics (i.e."online approach").In order to reduce the computational cost, the model introduces a so-called "morphological accelerating factor" (Mor-Fac) which linearly scales up the bed level change; see Lesser et al. (2004) and van der Wegen and Roelvink (2008) for more details.(d-f) meso-tidal condition with the first bathymetry, and (g-i) meso-tidal condition with the second bathymetry.
Both the scaled-up micro-and meso-tidal cases are simulated in this study according to the scaling laws given in Table 1.A grid mesh of 48 300 rectangular cells is adopted for both cases with different cell sizes (micro-tidal case: 22.5 m × 22.5 m; meso-tidal case: 31 m × 31 m) and their input bathymetries are deduced by scaling up the bottom perturbations of the initial measured experimental bathymetry with proper multiplier factors.The values of the multiplier factors are chosen on the basis of the scaling analysis presented by Stefanon et al. (2010).In particular, the scaling ratios were derived based on the similarities of hydrodynamics (Froude similitude and the similitude of local inertia and advection) and sediment transport (similitude of Shields parameter and particle Reynolds number) between the prototype and physical model.In order to maintain those similarities, the vertical depth scaling factor for micro-tidal and meso-tidal cases is 100 and 200, respectively (Table 1).The bottom perturbation of the initial experimental bathymetry of Run 4 is approximately within ±0.004 m; hence the perturbations of the micro-tidal and meso-tidal initial bathymetries are set within ±0.4 and ±0.8 m, respectively.As for the meso-tidal case, a second initial bathymetry with different bottom perturbations is considered in order to investigate the role of different random bed perturbations on the morphological evolution of tidal networks.The magnitude of the bottom perturbations for the second bathymetry is the same as the first one, whereas the spatial distributions which were generated randomly are different.Some other key model parameters are chosen through sensitivity tests (Zhou et al., 2014), including time step (∆t = 30 s), morphological accelerating factor (MorFac = 20), transversal bed slope term (α bn = 50) and Chézy friction coefficient (C = 65 m 1/2 s −1 ).

Pattern formation and velocity field
The morphological configurations of the micro-tidal (tidal range: 1 m) and meso-tidal (tidal range: 2 m) basins after 5, 65 and 125 yr are shown in Fig. 3.Much like the laboratory observations, the initial phase of the channel development proved to be rapid.By scouring the inlet, a central channel was formed within 5 yr.The channel subsequently branched into smaller tributaries which grew wider and deeper and continued branching, meandering and elongating.After 5 yr, approximately one-third of the area of the micro-tidal and over half of the area of the meso-tidal basins (for both the bathymetries) were dissected by the expanding channel network (Fig. 3a, d and g).As the tidal prism kept increasing, the network continued to develop rapidly until a general quasi-stable tidal network took shape after around 50-60 yr (Fig. 3b, e and h).The growth of the channel network slowed down and adjusted to a stable planar configuration after approximately 120 yr (Fig. 3c, f and i).After that, only the inlet area was still subject to small erosion (leading to further inlet deepening) until the morphological equilibrium state was reached after around 200 yr.
The average depth of the inlet area of the stable tidal network observed in the laboratory is around 0.05-0.06m (Fig. 2b), and the corresponding depths of the numerically modelled micro-tidal and meso-tidal systems are 5-6 m and 10-11 m, respectively.This agrees well with the scaling analysis of Stefanon et al. (2010) (see Table 1).Overall, the average channel depths of the micro-tidal and meso-tidal simulated networks are approximately 2-3 m and 3-4 m, respectively, with the average tidal flat depths of about 0.3-0.4m and 0.5-0.6 m, respectively.These morphological features generally resemble natural examples such as microtidal Venice Lagoon and meso-tidal Dutch Wadden Sea (Stefanon et al., 2010).Although sharing a similar trend in channel network development, the detailed morphological evolution of the micro-tidal and meso-tidal simulations differs in primarily two aspects.First, these two networks display different final morphologies: the channels formed in the mesotidal case are more meandering and occupy a larger portion of the basin area than in the micro-tidal case.Second, the channel network developed more rapidly in the meso-tidal case and the channels were wider and deeper after the same amount of simulation time.This was also found by van Maanen et al. (2013), who showed that tidal networks grow more rapidly for larger tidal range.Finally, our simulations show that the initial bottom perturbations play an important role in the morphological development of tidal networks (Fig. 3d-i), which was also noticed in the laboratory experiments.
The spatial distribution of velocity magnitude at the maximum flood condition based on stable networks is shown in Fig. 4. As expected, velocities in the channels are larger than in the tidal flats (see Fig. 3c and f for morphologies).The maximum velocities are found near the inlet area, with an average value of approximately 0.6-0.7 m s −1 for the microtidal case (Fig. 4a) and 0.7-0.8m s −1 for the meso-tidal case (Fig. 4b).These values fall within a reasonable range when compared to natural lagoon systems and are also consistent with the scaling analysis (Table 1) of Stefanon et al. (2010).

Drainage density
Following the method proposed by Marani et al. (2003) in the framework of tidal systems, the drainage density of tidal networks is evaluated on the basis of the exceedance probability distribution of the unchannelled flow lengths, determined using the Poisson-type hydrodynamic model developed by Rinaldo et al. (1999a).Previous studies on natural tidal networks indicate that the probability distribution of unchannelled flow lengths is exponential (Marani et al., 2003;Feola et al., 2005;D'Alpaos et al., 2007;Vandenbruwaene et al., 2012a).The probability density function shows a linear trend when it is plotted on a semilogarithmic scale, and the inverse of the slope represents the mean unchannelled flow length.The drainage density is larger when the slope is steeper, suggesting that the channel networks are more efficient in draining tidal flats and marshes (Marani et al., 2003).
Figure 5 shows the exceedance probability distribution of tidal networks after 5, 35, 65 and 125 yr under micro-tidal conditions.After 5 yr, the slope of the semilogarithmic distribution is relatively gentle when less than one-third of the tidal flats is drained.65 and 125 yr are close since the basic structure of the tidal network has taken shape after 35 yr and only some minor morphological adjustments occur afterwards.The slope of the semilogarithmic distribution increases with time, indicating a decreasing mean unchannelled flow length and an increasing drainage density (Fig. 6).The mean unchannelled flow length decreases sharply in the first 35 yr, indicating a rapid expansion of tidal networks.The decrease of the later phase (when the system is evolving towards equilibrium) is slower (Fig. 6a).Generally, the characteristic of the probability distribution (Fig. 5), as well as the evolution of the mean unchannelled flow length and drainage density (Fig. 6), is in accordance with the trend of the morphological evolution: the initial phase of development is rapid and slows down with time, especially when the morphological equilibrium approaches.Similar results were also presented by Vandenbruwaene et al. (2012a, b), who conducted a 5 yr field study of spontaneous formation and evolution of a tidal network in the Scheldt Estuary, Belgium.
The morphological evolution of the micro-tidal and mesotidal cases differs pronouncedly (Fig. 3), which is also reflected by the exceedance probability distribution.After the same amount of time (i.e. 5 and 125 yr), the unchannelled flow length of the low exceedance probability is larger under meso-tidal conditions (Fig. 7a) than micro-tidal ones (Fig. 5), which occurs because the dimension of the mesotidal basin is larger.In order to compare the laboratoryscale experiment Run 4 with the micro-tidal and meso-tidal simulations, the unchannelled flow lengths are scaled with the corresponding basin widths (see Fig. 7b).The basin is drained more efficiently under the meso-tidal condition than micro-tidal conditions after the same time (5 and 125 yr; see also Fig. 8b) since the dimensionless exceedance probability curve of the meso-tidal case displays a larger slope.This is consistent with the smaller mean unchannelled flow length shown by Fig. 8a.The drainage densities of both the microtidal and meso-tidal simulated stable networks are larger than those found in the laboratory experiment (represented by the marker "+" in Fig. 7b).Overall, different initial conditions can result in wide differences in either the experimental (see Fig. 13 of Stefanon et al., 2010) or the numerically simulated drainage systems (Figs.7 and 9).Stefanon et al. (2010) performed two experiments using different initial bottom perturbations under the same tidal forcing and found that the final patterns obtained in the laboratory differed pronouncedly.Their finding was also reflected by the distinctive exceedance probability distributions of the unchannelled flow lengths.Our numerical model results are consistent with the laboratory observations.The evolved tidal network after 5 yr, obtained using the first bathymetry as initial condition, showed less drained area than the one obtained when considering the second bathymetry (characterised by different randomly generated bed perturbations with respect to the first bathymetry) as an initial condition, which is in agreement with the larger slope of the exceedance probability distribution (in a semilog plot) of the unchannelled flow lengths obtained starting from the second bathymetry (Fig. 9).However, the first bathymetry leads to a larger final drainage density than the second when the morphological equilibrium configuration is reached.This indicates that the lagoon of the second bathymetry can be better drained near the inlet area, while the inner lagoon is drained less efficiently when compared to the lagoon of the first bathymetry.This is simply the result of different spatially distributed random perturbations of the bottom and the morphodynamic feedbacks they generated.Analogous to the micro-tidal case shown in Fig. 6a, the mean unchannelled flow length decreases rapidly in the beginning phase and tends to be stable when morphological equilibrium is approached (blue lines in Fig. 8a).
The evolution of the mean unchannelled flow length and drainage density display asymptotic behaviour (Fig. 8), as also found by Vandenbruwaene et al. (2012a).As a primary land-forming force shaping channel networks, the tidal prism of both micro-tidal and meso-tidal simulations also shows a similar evolution (Fig. 10).The tidal prism increases sharply in the first 35 yr when the channel network expands, quickly dissecting the tidal flats, and slows down afterwards.Although sharing a similar trend, the tidal prism of the mesotidal case is larger than the micro-tidal one due to the different tidal ranges and basin sizes (Fig. 10).Initially, the tidal prism adapts to the morphological evolutions more rapidly than the drainage density does and then the two quantities vary with similar temporal scales, as also found by Stefanon et al. (2012).Both the evolution of drainage density and the tidal prism displays asymptotic behaviour (Fig. 10).Therefore, numerical simulations predict that the drainage density of a stable channel network system should correspond to a certain tidal prism.In fact, using varying water levels to simulate different sea level conditions (hence different morphological equilibrium systems), the laboratory experiments conducted by Stefanon et al. (2012) indicate an approximately linear relationship between the land-forming tidal prism and drainage density, which is also in qualitative agreement with the results of Marani et al. (2003) for natural salt-marsh basins.

Discussion
The existence of feedback mechanisms between various landscape-forming processes (e.g.tides, waves, biological effects, sea level rise and human intervention) limits our understanding of the ontogeny of tidal networks (Rinaldo et al., 1999a;Coco et al., 2013).By reducing the number of processes affecting the evolution of tidal networks, controlled laboratory experiments provide a good base for gaining insight into these systems.While tidal flows are the external driver of bathymetric changes, the slowly evolving estuarine bathymetry can also constrain the tidal flow and determine flow patterns.Therefore, estuarine landscapes are the result of the mutual adaptation between flow and morphology (i.e.morphodynamic feedback).The laboratory-observed tidal networks of Stefanon et al. (2010Stefanon et al. ( , 2012) ) are morphologically comparable to natural ones.However, their quantitative analysis almost certainly suffers from inevitable scaling effects, which can be difficult to evaluate, especially when translated to long-term numerical modelling.Based on scaling arguments (Stefanon et al., 2010), we set up numerical simulations in an attempt to understand the evolution of drainage density of tidal networks from a numerical modelling point of view and also examine to which degree numerical model results agree with the scaling laws.Quantitatively, the numerically simulated tidal networks resemble those obtained in the laboratory experiments (Stefanon et al., 2010(Stefanon et al., , 2012;;Kleinhans et al., 2012) and some natural back-barrier tide-dominated systems (e.g. the Venice Lagoon and the Dutch Wadden Sea).The flow velocities fall in a reasonable range both in the tidal channels and the flats and are in good agreement with the scaling arguments.Overall, the numerical model appears to be capable of reproducing the major features of tidal network ontogeny, i.e. from initiation to equilibrium.In this context, it is important to specify that the equilibrium (or stable) state considered in this study is not a "frozen" or "static" equilibrium state (strictly null sediment fluxes and no bed level change), which is highly unfeasible in reality (especially if one considers that external forcing is constantly changing).Instead, a different definition, sometimes indicated as "dynamical equilibrium" and characterised by null gradients in sediment fluxes (and no bed level change) over a specific timescale, is more appropriate to describe our numerical simulations.In fact our results indicate that the gradients in sediment fluxes, and for continuity bed level changes, decreased over time and approached zero after approximately 200 yr, even though sediment fluxes  were never strictly null.Our definition of dynamic equilibrium is also different from the usual definition of statistical equilibrium, which implies that small fluctuations in the system do not necessarily cease even if the overall property of the network remains unchanged.
Consistent with the field observation of the natural systems in the Scheldt Estuary, Belgium (Vandenbruwaene et al., 2012a), the drainage density of the numerically simulated tidal networks is asymptotic towards an equilibrium configuration.The exceedance probability distribution of the unchannelled flow lengths of the simulated tidal networks displays an exponential trend, in analogy with the laboratory experiments and natural tidal networks in the Venice Lagoon (Marani et al., 2003;D'Alpaos et al., 2007a).The exceedance probability distributions of the unchannelled flow lengths for the micro-tidal and meso-tidal scale are different after the same amount of simulation time, suggesting, from this perspective, the absence of scale-invariant processes.As previously indicated, uncertainties in scaling of laboratory observations could potentially cascade on the scaling adopted in the numerical simulations.Model results indicate that the drainage density of the meso-tidal basin is larger due to the larger tidal prism.The tidal prism increases rapidly during the beginning phase of tidal network development and soon reaches a relatively steady value, while the drainage density adjusts more slowly, which is consistent with Stefanon et al. (2012).
Model results also show that the random bottom perturbations play a considerable role in determining the long-term morphological development of tidal networks, and hence the overall evolution of the drainage density.Consistent with laboratory observations (Stefanon et al., 2010), the different ex-ceedance probability distributions of the unchannelled flow length (Fig. 9) reflect the pronounced difference of final configurations of tidal networks.This indicates that small differences in the initial characteristics of tidal lagoons affect the later evolution and ultimately result in pronouncedly different landscapes.Another metric to quantify the general differences between the resulting bathymetries is the hypsometric curves shown in Fig. 11a.The cyan and magenta dashed lines represent the meso-tidal starting hypsometries for bathymetry 1 and 2, respectively.The blue and red solid lines indicate the corresponding final hypsometries at equilibrium.These two numerical simulations started with input bathymetries characterised by different spatial random perturbations (the magnitude of the perturbations was the same) and nearly 95 % of the bed elevation was around the mean sea level.The existence of higher areas near the barrier islands (the red areas in Fig. 3) approximately accounts for the remaining 5 % of the lagoon (notice that in the numerical simulations these higher areas remained unchanged since most of that area was dry even at high tide).Due to the statistical similarity in the seabed perturbations, the two initial bathymetries (dashed lines) show similar hypsometries (the elevation of the first bathymetry is slightly lower).The resulting final hypsometries also share a similar morphological pattern, but some noticeable differences are present.The two hypsometric curves intersect so that the second hypsometry (red line) is higher than the first one (blue line) for areas lower than 77 %, while the opposite occurs for areas larger than 77 %.Overall, the shallower areas are less channelised, while the deep areas are more incised (second case), consistent with the final simulated morphologies (Fig. 3d-i), which show that less area is drained and the main channels are deeper in the second case.This can also be noted from Fig. 11b and c, which are the initial and final differences of bed elevation between the first and second bathymetries, respectively.The initial differences are generally within the range of −0.8 to 0.8 m, while the final differences can reach 5 m with different spatial distribution of channels and tidal flats.Nonetheless, hypsometric curves could only capture a general difference between the two final morphologies due to the aggregated nature of the method.Instead, the probability distributions of unchannelled flow lengths (Fig. 9) provided a better way to investigate the differences between the two cases.Some other metrics (e.g.evolution of lagoon bed level change, channel width-to-depth ratios and the relationship between tidal prism and cross-sectional area) have been discussed in detail elsewhere (Zhou et al., 2014), and overall, they indicate close similarity between experimental and numerical tidal networks.
Nevertheless, there are limitations in the current study that need to be addressed.Firstly, both the experiment and numerical models do not include the effects of waves, biological activities, mud (or sand-mud mixture) and possible river inflow, all of which are commonly present in estuarine environments.Therefore, more processes should be added in both the physical and numerical models so that a more realistic condition can be considered in the future.Secondly, the Poisson-type hydrodynamic model (Rinaldo et al., 1999a), as adopted to calculate the exceedance probability distribution of the unchannelled flow length (and hence the drainage density), assumes friction terms dominate over inertial and advective terms in the momentum equations.This assumption holds well for shallow basins in which the water depth is normally smaller than 5 m (e.g. the Venice Lagoon), while it may not be valid for the deep area in the model domain.
As for the current study, the effect of this limitation is minor for the micro-tidal case (water depth is generally below 5 m), whereas the effect is more pronounced for the meso-tidal case since channels are more incised with a larger tidal prism, especially near the inlet area.This limitation can be overcome by using a more sophisticated hydrodynamic model (e.g.Delft3D); however, this may be computationally demanding.The reader is referred to Jiménez et al. (2013b) for some ongoing research on this topic.

Conclusions
We adopt a state-of-the-art numerical model (Delft3D) to investigate the morphological development of tidal networks and compare it to a controlled laboratory experiment.The numerical model is able to reproduce the general behaviour of tidal network ontogeny and the modelled velocity field agrees with the scaling analysis and different natural reference systems.
The comparison primarily focuses on the analysis of the drainage density of tidal networks.Instead of using the classic Horton's method, we analyse the drainage density based on the exceedance probability distribution of the unchannelled flow lengths (the mean of this exponential probability distribution, i.e. the mean unchannelled flow length, is the drainage density).Much like laboratory experiments and natural systems, analysis of the numerically simulated tidal networks shows that the probability distribution of unchannelled lengths displays an exponential trend.In accordance with the observed morphological evolution, the simulated drainage density increases rapidly during the initial phase of tidal network development and slows down when equilibrium is approached.Because of the larger tidal prism, the meso-tidal basin is more drained than the micro-tidal case when stable channel networks are attained.The initial increase of the tidal prism is more rapid than that of the drainage density, while later the increase is slower, indicating that the tidal prism adjusts faster to a stable condition.The initial bottom perturbations are found to play an important role in determining the morphological development, which is also captured in the exceedance probability distribution of the unchannelled flow length.
Using the unchannelled flow length as a descriptor for statistical analysis, we have shown that the numerically modelled tidal networks are comparable to the experimental and natural ones.We have also shown that tidal networks display distinctive morphological features that can be influenced by various factors (e.g.tidal range, initial bottom perturbations).Finally, this study highlights the benefit (and hence need) of complementary studies between physical and numerical modelling in gaining insight into the long-term morphological evolution of estuarine systems.

Figure 1 .
Figure 1.Demonstration of a typical barrier tidal basin system, showing various geomorphological features including flood and ebb deltas, channel networks, barrier beaches, tidal flats and salt marshes.The system is Hampton Harbour inlet, located at 42 • 53 48 N and 70 • 49 06 W on the east coast of the US.The image was taken from Google Earth (dated 1 January 2010), accessed on 11 October 2013.© Image USDA, Farm Service Agency, Google Earth 2013.

Figure 2 .
Figure 2. Laboratory experimental setting and results of Stefanon et al. (2010).(a) Schematic overview of the experimental apparatus; (b) final stable channel network obtained in the lagoon of the experiment Run 4.

Figure 3 .
Figure 3. Morphological pattern formations after 5, 65 and 125 yr simulated by the numerical model.(a-c) Micro-tidal condition,(d-f) meso-tidal condition with the first bathymetry, and (g-i) meso-tidal condition with the second bathymetry.

Figure 4 .
Figure 4. Spatial distribution of velocity magnitude at the maximum flood condition based on a stable tidal network: (a) micro-tidal case and (b) meso-tidal case with the first bathymetry.

Figure 5 .
Figure 5. Semilog plot of the exceedance probability P(L > l) of the unchannelled flow length l, computed for the configurations of channel networks under micro-tidal conditions after 5, 35, 65 and 125 yr; see also Fig. 3a-c for the morphologies of corresponding networks.

Figure 6 .
Figure 6.(a) Evolution of the mean unchannelled flow length l m with time, and (b) evolution of the drainage density D of the developing tidal networks under micro-tidal conditions.

Figure 7 .
Figure 7. (a)Semilog plot of the exceedance probability P(L > l) distribution of the unchannelled flow length l, computed for the configurations of numerical tidal networks under meso-tidal conditions after 5 and 125 yr, and (b) semilog plot of the exceedance probability P(L/W > l/w) of dimensionless unchannelled flow length l/w (scaled with basin width w), computed for the configurations of stable experimental tidal networks (Run 4) and numerical tidal networks under both micro-tidal and meso-tidal conditions after 5 and 125 yr; see also Fig.3for the morphologies of corresponding networks.

Figure 8 .
Figure 8.Comparison of the evolution of micro-tidal first, meso-tidal first and meso-tidal second bathymetry with regard to (a) the mean unchannelled flow length and (b) the drainage density; see also Fig. 3.

Figure 9 .
Figure 9. Semilog plot of the exceedance probability distribution of the unchannelled flow length, computed for the configurations of channel networks with different initial bottom perturbations under meso-tidal conditions after 5 and 125 yr.

Figure 10 .
Figure 10.Evolution of tidal prism and drainage density for (a) the micro-tidal case and (b) the meso-tidal case.

Figure 11 .
Figure 11.(a) Initial and final hypsometries of the first and second meso-tidal bathymetries; (b) and (c) are the initial and final differences of bed elevation between bathymetry 1 and 2, respectively (note that the colour bar scales in panel (b) and (c) are different).For bed elevation, the upward direction is positive.

Table 1 .
Stefanon et al. (2010)nd numerical model parameters after scaling analysis ofStefanon et al. (2010).W , L and D represent width, length and depth scales of the basins; U denotes the typical velocity in tidal channels.