Introduction

Tidal networks, characterized by repeating bifurcations and a tree-like structure (Fig. 1), serve as the arteries for the dynamic circulation of matter and energy within tidal flats and wetlands worldwide1,2. These networks are integral to the functioning of ecosystems that provide invaluable ecosystem services, including carbon storage, habitat provision, and coastal defense3,4. Understanding the structure and allometry of tidal networks is essential to unraveling the ecology of these environments5, and hence to facilitating effective restoration and conservation in the face of environmental change6,7.

Fig. 1: Branching tidal networks worldwide.
figure 1

a Overview of the collected tidal networks covering 21 sites worldwide, encompassing a diverse range of biomorphodynamic conditions, i.e., from micro- to macrotidal ranges and from bare to vegetated flats. b–v Remote sensing data of branching tidal networks. The randomly colored patches with white borderlines represent the extracted tidal networks, while the thin gray lines indicate the channel skeletons. Panels (j), (n), and (r) are derived from digital elevation model data, and the remaining panels are extracted from satellite imagery. The map in panel (j) contains public sector information licensed under the Open Government Licence v3.0. https://www.nationalarchives.gov.uk/doc/open-government-licence/version/3/. The detailed information for these regions can be found in Supplementary Tab. 1.

Full size image

Branching networks with recursive structures are not unique to tidal areas; they represent a universal pattern observed across various natural systems, from dendritic rivers that carve through landscapes to the intricate nerves and blood vessels that sustain human life8,9. These networks are strongly reminiscent of fractals10, patterns that are widely recognized as nature’s optimized strategy for filling space and distributing resources11. Researchers have long applied power-law scaling relationships to describe the order and universality within branching networks12,13,14,15. A prime example in fluvial studies is Hack’s law16, which relates the length of the main (longest) channel of a given river network (\({L}_{m}\)) to its drainage area (\({A}_{d}\)) using a power function: \({L}_{m}\propto {A}_{d}^{h}\). The exponent \(h\) is typically 0.6 across diverse basin sizes and geographic contexts17,18,19,20. Power-law functions are appropriate for describing relationships that remain the same from local to global levels, which is also called scale invariance21. The prevailing conceptual models of branching networks, including those based on a random topology22,23 and principles of optimization24,25, are considered credible because they replicate this scale invariance.

However, tidal networks do not exhibit the scale invariance that forms the basis for understanding branching networks1. These channel networks defy conventional simplifications using power-law functions, particularly in the relationship between \({L}_{m}\) and \({A}_{d}\) (i.e., Hack’s law)1 and in the probability distribution of the unchanneled length (\({L}_{u}\))26,27—a metric representing the flow distance across a marsh before a channel is reached. Researchers have generally linked this observed divergence to the complex interactions among ocean dynamics (e.g., tides, waves and sea level rise)28,29,30, sediment transport31,32,33, and biotic processes (e.g., vegetation colonization)34,35,36,37,38. As a result, extensive studies have investigated these factors2,6,39, identifying various mechanisms that can affect tidal network morphology, such as the role of vegetation in disturbing flow and stabilizing the sediment bed40, and the nonlinear responses of hydrodynamics and coastal ecosystems to sea level rise41. However, beyond these case-specific mechanisms, an alternative possibility remains unexplored: that the anomalous behavior of tidal networks might represent a more holistic, fundamental principle not captured by existing theoretical models.

In this study, we explored this hypothesis by investigating the ideal branching network that should result from such a general principle. First, by analyzing the common features of tidal networks collected worldwide (Fig. 1), we derived that the ideal branching network is a space-filling fractal, although its specific shape remains undetermined. We then formulated a simple theoretical model, generalized from the morphodynamics of tidal networks, suggesting that channels occupying larger drainage areas are more likely to extend. The model produced a regular space-filling fractal pattern, supporting our hypothesis regarding the ideal branching network geometry. A comparison between the modeled and field networks revealed an inherent “laziness”, where tidal networks in settings more conducive to channel formation exhibit more pronounced deviations from the ideal space-filling branching fractal. This result explains the commonly observed anomalous scaling behavior in tidal networks. Furthermore, the identification of “laziness” opens up the possibility for controlling the morphogenesis of tidal networks by strategically manipulating environmental factors. Our model also potentially provides an alternative theoretical explanation for the emergence of branching networks in nature.

Results

Conceptualization of the ideal branching network

We begin our analysis by investigating the scaling properties of real tidal networks, leveraging an extensive dataset including more than 300 networks from 21 sites worldwide (Fig. 1). We geometrically delineated every watershed by assigning each point within the drainage platform to its nearest channel (Fig. 2a). The distance from each point to its corresponding channel approximately represents the unchanneled length (\({L}_{u}\))34. To assess the scaling of \({L}_{m}\) to \({A}_{d}\) (i.e., Hack’s law), we decomposed each network level by level, resulting in a collection of subnetworks (with branches) and terminal channels (without branches) spanning multiple scales (Fig. 2a). By fitting the exponent \(h\) of the scaling law \({L}_{m}\propto {A}_{d}^{h}\), we observed significant variations across sites and scales (Fig. 2b). Moreover, the exceedance probability distributions of unchanneled length \(P\left({L}_{u}\ge {l}_{u}\right)\) of all collected tidal networks (Fig. 2c-d) exhibited a cluster of curves that defy conventional empirical models based on a power-law or exponential fits26,27,42,43. Here \({L}_{u}\) denotes the set of unchanneled lengths extracted from a tidal network (Fig. 2a), and \({l}_{u}\) represents a specific value of the unchanneled length.

Fig. 2: Scaling relationships of tidal networks that defy conventional simplification through simple empirical fits.
figure 2

a Data processing workflow in this study, from image data to the distribution of the unchanneled length (\({L}_{u}\)), followed by level-by-level decomposition of the drainage basin, which results in an array of subnetworks and terminal channels. The randomly colored patches with white borderlines represent the extracted subnetworks. b Discontinuity in the scaling relationships between the main channel length (\({L}_{m}\)) and drainage area (\({A}_{d}\)). c, d Divergence of the exceedance probability distributions of the geometric unchanneled length, \(P\left({L}_{u}\ge {l}_{u}\right)\), represented in log-log and semi-log coordinates. e Three sets of length–area relationships resulting from isolating the data for unbranched terminal channels, including \(\Sigma L\) versus \({A}_{d}\) (in blue) and \({L}_{m}\) versus \({A}_{d}\) (in red) for all subnetworks and \({L}_{m}\) versus \({A}_{d}\) (in green) for all terminal channels (\(\Sigma L={L}_{m}\) in this case). f Violin plots summarizing the variability in the three sets of length–area scaling exponents for all collected tidal networks. g Exceedance probability distributions of the nondimensional geometric unchanneled length \({L}_{u}^{*}\) for all collected tidal networks. The blue-to-white color code denotes the spatial density of these curves. The dashed line denotes the linear function \(P\left({L}_{u}^{*}\ge {l}_{u}^{*}\right)=1-{l}_{u}^{*}\), which is derived based on a simple scenario of a straight channel bisecting a rectangular basin.

Full size image

Our objective is to identify the common factors among these various patterns and deduce the characteristics of an ideal branching network. To achieve this objective, we introduced the total channel length (\(\Sigma L\)) of a network as a comparative metric. Notably, \(\Sigma L={L}_{m}\) in the case of an unbranched terminal channel (where only one channel exists) and \(\Sigma L \, > \, {L}_{m}\) in a subnetwork with branches. We then isolated the data for unbranched terminal channels, resulting in three sets of length–area scaling relationships as shown in Fig. 2e (refer to Supplementary Fig. 1 for the scaling exponents for each site). For all collected data (Fig. 2f), subnetworks with branches show that \(\Sigma L\) scales linearly with \({A}_{d}\), denoted as \(\Sigma L\propto {A}_{d}\) (in blue); while \({L}_{m}\) scales with \({A}_{d}\), with the exponent fluctuating around 0.6 (in red). The unbranched terminal channels indicate intermediate scaling exponent values (in green). This transition suggests that the widely-observed deviation from Hack’s law (e.g., Fig. 2b) reflect a compromise between \({L}_{m}\propto {A}_{d}^{h}\) and \(\Sigma L\propto {A}_{d}\), constrained by the inability of real networks to branch infinitely. Note that in each network of Fig. 2b, subnetworks and terminal channels are separated, creating a scale break similar to that in Fig. 2e.

We then derived the exceedance probability of unchanneled length \(P\left({L}_{u}\ge {l}_{u}\right)\) for the ideal branching network. As \(\sum L\propto {A}_{d}\) indicates that each channel segment occupies an equal proportion of the drainage area, a straightforward scenario that fulfills this criterion is a straight channel bisecting a rectangular basin (Fig. 2g). Under this scenario, we can directly derive that \(P\left({L}_{u}\ge {l}_{u}\right)=1-2{D}_{d}{l}_{u}\), where \({D}_{d}=\Sigma L/{A}_{d}\) is the drainage density. By introducing the dimensionless parameter \({L}^{*}_{u}=2{L}_{u}{D}_{d}\), we obtain:

$$P\left({L}_{u}^{*}\ge {l}_{u}^{*}\right)=1-{l}_{u}^{*}$$
(1)

To verify this ideal form, we nondimensionalized the geometric unchanneled lengths of all collected tidal networks and integrated their exceedance probability distributions within a linear coordinate system (Fig. 2g). After these operations, the curves exhibit denser clustering than their dimensional counterparts (Fig. 2c, d). In addition, these curves closely align to the right of the black dashed line representing Eq. (1). This result suggests that the ideal branching network may represent a limiting condition that real-world networks strive toward but cannot fully achieve.

Intriguingly, \(\varSigma L\propto {A}_{d}\) also presents a fascinating scenario where one-dimensional lines exhibit the same dimensionality as a two-dimensional plane. This phenomenon readily evokes associations with space-filling fractals10. Such fractals can fill the space they occupy completely, as their fractal dimension equals the topological dimension of the space, such as the Peano curve44. Given this relevance, we hypothesized that the ideal branching network might be a standard space-filling fractal network, even though its exact geometry remains undetermined.

Deviation from an ideal branching network

The above analysis motivated us to characterize the network morphology based on the deviation from the hypothetical space-filling model (i.e., the ideal branching network), instead of conventional power-law descriptions. To measure this deviation, we defined a metric named “space-filling deviation index”:

$${D}_{i}={\int_{0}^{\infty } }P\left({L}_{u}^{*}\ge {l}_{u}^{*}\right){{{\rm{d}}}}{l}_{u}^{*}-0.5$$
(2)

The graphical representation of this index is the area enclosed between a probability curve in Fig. 2g and the black dashed line that denotes the hypothetical baseline case (i.e., Eq. 1). A value of 0 indicate no deviation as we can see by substituting Eq. (1) into Eq. (2). A larger \({D}_{i}\) value indicates greater deviation.

We then examined the morphological variability in tidal networks through a diagram of the drainage density \({D}_{d}\) versus the index \({D}_{i}\) (Fig. 3a). Here, \({D}_{d}\) is a measure of the progress in network growth, while \({D}_{i}\) denotes the degree to which branching structures occupy the drainage platform. Note that by definition, a larger \({D}_{i}\) value indicates reduced effectiveness in filling the drainage platform space.

Fig. 3: Deviation of tidal networks from the theoretical space-filling model.
figure 3

a Diagram of the drainage density \({D}_{d}\) versus the space-filling deviation index \({D}_{i}\) for different tidal networks worldwide. The size of the circle denotes the normalized number of branches \({N}_{b}^{*}\), and the color of the circle denotes the sinuosity \({\sigma }_{s}\). b, c Variability in \({D}_{i}\) in relation to \({N}_{b}^{*}\) and \({\sigma }_{s}\). The solid line in (c) denote the linear regression results visualizing the trend. d–f Examples of networks corresponding to \({D}_{d}\)–\({D}_{i}\) data pairs transitioning from the bottom-left to top-right corners in the diagram in (a). The three networks are indicated by the black squares in (a). The color scales in (d–f) represent the spatial patterns of the geometric unchanneled length \({L}_{u}\).

Full size image

The data points exhibit a diverse distribution, with \({D}_{d}\) spanning approximately from 0.01 to 0.1 \({{{{\rm{m}}}}}^{-1}\) and \({D}_{i}\) ranging between 0.05 and 0.35 (Fig. 3a). While the pattern seems random with no clear trend, it unveils several counterintuitive findings. Some networks develop a dense distribution of channels (large \({D}_{d}\)) yet deviating markedly from the ideal space-filling condition (large \({D}_{i}\)). Conversely, some achieve higher spatial occupation (small \({D}_{i}\)) despite much sparser channel dissection (small \({D}_{d}\)). This contrast underscores the significant variability among these networks in their efficiency at utilizing the available space.

To explore the dynamics behind this variability, we examined three networks that span the entire distribution (Fig. 3d-f and the black boxes in Fig. 3a). In these three examples, networks with larger \({D}_{i}\) values exhibit small branches and fewer sinuous extensions. Inspired by this observation, we introduced two normalized parameters, namely, \({N}_{b}^{*}\) and \({\sigma }_{s}\), to quantify the number of branches and the sinuosity of the network. These two parameters are represented by the sizes and colors of the bubbles, respectively, in Fig. 3a. We found that \({D}_{i}\) did not exhibit a clear correlation with \({\sigma }_{s}\) (Fig. 3b, Spearman’s rank correlation coefficient \(\rho \approx 0\), with a p value > 0.05). In contrast, \({D}_{i}\) significantly increased with \({N}_{b}^{*}\) (Fig. 3c, \(\rho \approx 0.7\), with a p value < 0.05). This result suggests system-level “laziness”—where increased ease of new channel formation (i.e., more channel branches) paradoxically leads to decreased effectiveness in occupying the available space (i.e., larger \({D}_{i}\) value).

Constructing the ideal branching network

Based on the positive correlation between channel length and drainage area (Fig. 2b), we proposed a simple rule: channel edges, including channel heads and side banks, fed by larger drainage areas have a higher probability of longitudinal extension or lateral tributary development. We developed a theoretical model to implement this rule. The model operates on a two-dimensional lattice, where each cell can occur in one of three states: platform, channel, or edge between the platform and channel (Fig. 4a). In each iteration step, the model computes the erosion potential (\(E\)) of each edge cell, which increases monotonically with drainage area. Subsequently, the model randomly transforms an edge cell to a channel cell with probability \({p}_{i}={E}_{i}/\Sigma {E}_{i}\), where the subscript \(i\) denotes the edge cell index. This selection procedure represents competition among all edge cells within a catchment.

Fig. 4: Branching networks based on the competition principle.
figure 4

a Sketch of the model lattice. Each cell can exhibit one of three states: platform, channel, or edge between the platform and channel. b Variability in the erosion potential \(E\) calculated with different activation drainage areas \({A}_{{act}}\) (Eq. 3). Here, \(E\) is normalized by the erosion potential calculated as \({A}_{d}=200\) m2, which allows clear comparison between these curves. c, d Modeled network growth within a \(1200{{{\rm{m}}}}\times 1200{{{\rm{m}}}}\) lattice using small and large \({A}_{{act}}\) values. e Network growth modeled by restricting channel extension toward the edge cells that occupy the largest drainage area (equivalent to\(\,{A}_{{act}}\to \infty\)).

Full size image

The erosion potential is calculated as

$$E\propto {A}_{d}\exp \left(-\frac{{A}_{{act}}}{{A}_{d}}\right)$$
(3)

where \({A}_{{act}}\) is the activation drainage area above which erosion is significant. \({A}_{{act}}\) is conceptually similar to the energy barrier that must be overcome for a chemical reaction to proceed at an appreciable rate, as expressed in the Arrhenius equation45. Within the present context, an edge cell with drainage area \({A}_{d} \, < \, {A}_{{act}}\) exhibits a substantially reduced potential for transformation into a channel. Correspondingly, a high \({A}_{{act}}\) value skews the competition by amplifying the advantage of edge cells occupying larger drainage areas, as shown in the graphical representation of Eq. (3) (Fig. 4b). Equation (3) follows a well-established formula for assessing sediment erosion by flow46. While the drainage area is directly employed as erosion power in Eq. (3), the underlying rationale is that a larger drainage area supplies more water, thereby enhancing flow erosion and promoting channel formation47.

We assessed the model using different basin shapes and initial conditions (Supplementary Fig. 2). Here, we describe scenarios run with a classic square lattice initialized by a channel cell at the bottom-left corner (e.g., Fig. 4a). As shown in Fig. 4c, d, branching networks spontaneously developed in all the simulations (see also Supplementary Movie 1). A larger \({A}_{{act}}\) value yields increased proximity to the hypothesized ideal space-filling state characterized by a reduced deviation index \({D}_{i}\). Pursuing this observation, we explored the limit state of the model by adopting \({A}_{{act}}\to \infty\), which restricts channel extension to the edge cells that occupy the largest drainage area. Strikingly, this limiting scenario produced a regular network (Fig. 4e and Supplementary Movie 2), which is identical to the space-filling fractal tree mathematically constructed based on a recursive scheme (Supplementary Fig. 3)48. This consistency directly elucidates the exact shape of the ideal branching network we hypothesized above, and shows how it can emerge from a localized competition principle. We denoted the ideal branching network as \({T}_{{square}}\).

Laziness in the morphogenesis of branching networks

As shown in Fig. 4c, d, the modeled networks exhibit less efficient spatial occupation (i.e., large \({D}_{i}\)), when channel extension activation requires less effort (i.e., smaller drainage area \({A}_{{act}}\)). This variation mirrors our observations regarding laziness in real-world tidal networks (Fig. 3c). The model explains this laziness as competitive dynamics: a reduction in \({A}_{{act}}\), while facilitating channel formation, causes an increase in the frequency of inefficient extension events at occupying the available space. Notably, inefficient channel formations manifest as small branches in both the simulated networks (Fig. 4c) and real-world networks (Fig. 3f).

However, the above explanation relies heavily on the parameter \({A}_{{act}}\), which cannot be directly measured in natural systems. Thus, the following uncertainty remains: does the laziness observed in natural systems originate from the same competitive dynamics assumed in the model? To investigate this issue, we introduced the metric \({A}_{{tip}}\), defined as the drainage area feeding each end tip of the branching network (illustrated in Fig. 5a). The rationale behind this parameter is that the evolution of tidal networks with well-developed branching structures is generally slow49. This suggests that the drainage area occupied by each channel tip should be close to the threshold required for its extension. Since the drainage density and basin size of real-world networks vary, we normalized \({A}_{{tip}}\) to create a dimensionless metric \({A}_{{tip}}^{*}={A}_{{tip}}{D}_{d}^{2}\). A higher \({A}_{{tip}}^{*}\) value indicates that, relative to the basin size, a larger drainage area is needed to generate sufficient erosive force for channel formation. Therefore, \({A}_{{tip}}^{*}\) parallels \({A}_{{act}}\) by representing the difficulty of channel formation. However, \({A}_{{tip}}^{*}\) can be directly measured from imagery data (Fig. 5a, b), allowing a more empirical approach to examine the proposed competitive dynamics.

Fig. 5: Illustration of the laziness of branching morphogenesis utilizing the drainage areas of channel tips.
figure 5

a An example from Saeftinghe, Netherlands, depicting the spatial distribution of drainage areas at all channel tips across a drainage platform (i.e., \({A}_{{tip}}\)). The colored blocks denote the catchment of each tidal network, with the gray lines representing channel lines. The white patches with thin gray borders indicate the drainage areas of the channel tips. b Frequency distribution of \({A}_{{tip}}^{*}\) (\(={A}_{{tip}}{D}_{d}^{2}\)) calculated from the Saeftinghe drainage platform, with the median value of \({A}_{{tip}}^{*}\) also noted. c Variation in the space-filling deviation index \({D}_{i}\) with the median of \({A}_{{tip}}^{*}\). The gray points denote the field data collected in this study. The small, gradient-colored dots indicate the results of more than 6000 model simulations with varying \({A}_{{act}}\) values substituted into Eq. (3). The pink dashed lines delineate the characteristic thresholds derived from \({T}_{{square}}\). d Detailed view of the dashed box in the diagram of (c). The light gray squares correspond to the model results (consistent with the gradient-colored dots in (c)). The circles of varying sizes and colors denote the field data, detailed with mean tidal range and vegetation coverage information. The green line indicates the linear regression of the field data from vegetated flats, with a transparent band showing 95% confidence intervals.

Full size image

We used the median value of \({A}_{{tip}}^{*}\) to measure the overall difficulty of channel formation within a given drainage platform (Fig. 5a, b). We then compared the laziness between the natural and modeled networks through a diagram of \({D}_{i}\) versus the median of \({A}_{{tip}}^{*}\) (Fig. 5c). We analyzed a comprehensive suite of 6000 simulations with parameter \({A}_{{act}}\) ranging from \(0\) to \(\exp \left(9\right)\) (\(\approx 8103\)) \({{{{\rm{m}}}}}^{2}\). The modeled networks exhibited a linear decreasing trend in \({D}_{i}\) with \({A}_{{tip}}^{*}\), as shown by the colored scatter points in Fig. 5c. The field data acquired from sites worldwide (the gray dots in Fig. 5c), while indicating higher variance, largely overlapped with the modeled scenarios. Furthermore, the analytical solutions derived from \({T}_{{square}}\) yielded \({A}_{{tip}}^{*}=4/9\) and \({D}_{i}\approx 0.037\). These values closely conform with the extremities of both the modeled and natural networks, as indicated by the two intersecting pink dashed lines in Fig. 5c.

We further considered the dashed box in Fig. 5c to describe the pattern of the real-world networks (Fig. 5d). Vegetated flats (green circles) notably agreed with the model predictions, in contrast to the more scattered distribution of the bare-flat data (brown circles). Additionally, the bare-flat data generally appeared to be located in the upper-right portion relative to the vegetated flat data. This suggests that calculating \({A}_{{tip}}^{*}\) solely based on geometric distance could overestimate the actual drainage area feeding the channel tip in bare flats. These discrepancies can be attributed to the higher overall bed elevation and the mature, stable developmental phase of vegetated flats, which leads to a stronger concentration of flow within the channels1,34,50.

Therefore, we performed a linear regression exclusively for the vegetated flats, yielding a linearly decreasing trend, specifically, \({D}_{i}=-0.8617\times {A}_{{tip}}^{*}+0.4445\) (\({R}^{2}=70\%\)). This result supports that the “laziness” observed in natural systems and the model is rooted in the same competitive dynamics. However, it is also worth noting that \({A}_{{tip}}\) is not the only factor influencing the spatial occupation of real-world tidal networks. Other factors, such as channel meandering and bed sediment heterogeneity, could reduce the space-filling efficiency, resulting in a larger deviation index \({D}_{i}\) than anticipated. Therefore, the data points derived from the idealized theoretical model tend to fall below the fitted line for data from real-world networks as shown in Fig. 5c.

Discussion

The competition principle proposed in this study incorporates a deterministic increase in erosion potential with drainage area, and a random selection procedure. This configuration resonates with the two ingredients of nature’s design—necessity and chance51. By eliminating randomness in morphological evolution, the model physically produces a regular space-filling tree (\({T}_{{square}}\)), which is identical to the synthetic fractal generated mathematically48. This consistency indicates a clear delineation of the respective roles of necessity and chance in branching morphogenesis: necessity establishes the underlying fractality and order, while chance introduces variability and diversity, leading to a myriad of possible outcomes that share consistent statistical properties.

\({T}_{{square}}\) serves as a baseline for analyzing tidal network morphology by enabling the identification of characteristic thresholds through analytical solutions (e.g., \({D}_{i}\approx 0.037\)). Furthermore, the model revealed a systematic deviation from \({T}_{{square}}\) with decreasing difficulty in channel formation. These findings lay the groundwork for a universal framework to assess the effects of various dynamics on tidal network morphology. For instance, recent studies have shown that the presence of vegetation leads to more efficient drainage patterns34,36; notably, our model supports this observation, as vegetation amplifies the erosion resistance, thereby increasing the difficulty of channel formation. On the other hand, we also observed that a single factor, whether the tidal range or vegetation coverage, does not exclusively dictate the space-filling efficiency of the network (Fig. 5d). This complexity indicates the need for a more integrated approach in future studies to address the multifaceted dynamics that shape network morphology.

Our model also suggests a potential paradigm shift in capturing the morphogenesis of branching networks. Established theories conventionally view branching morphogenesis as a goal-oriented process, actively evolving toward states that enhance efficiency, minimize energy dissipation, or optimize other thermodynamic properties24,52,53. Our model proposes an alternative explanation by presenting a passive, shortsighted evolutionary process. The model does not pursue any optimized state; instead, it evolves following a straightforward, localized competition principle. This principle does not force the channels to branch, but branching occurs spontaneously as channels expand to occupy space. This behavior highlights the simplicity and elegance with which nature organizes itself without any deliberate guiding hand.

This interpretation complements rather than contradicts the prevailing understanding of branching networks. It emphasizes the dual nature in the growth of branching networks, where the observed “laziness” acts as a counterpoint to the well-recognized optimization properties. Indeed, several studies have argued that natural systems are not perfectly optimized in terms of defined state variables54, such as entropy53, energy dissipation55, and the index \({D}_{i}\) constructed in this study. Here, tidal networks offer an ideal empirical subject for exploring this imperfection, owing to their pronounced morphological diversity (e.g., Fig. 1). This diversity enables us to identify another key facet of branching morphogenesis in contrast to optimization, namely, the laziness phenomenon. This laziness reconciles the divergence and deviation from power-law scaling relationships (e.g., Hack’s law and the probability distribution of the unchanneled length illustrated in Fig. 2), bridging the gap in optimization theories.

The results of this study contribute to the conservation and restoration efforts of coastal wetlands—a pressing issue that has received widespread attention56. Our morphometric analysis introduces innovative metrics with theoretical benchmarks for quantifying the drainage conditions in coastal wetlands (e.g., \({D}_{i}\) and \({A}_{{tip}}^{*}\)). These metrics could provide researchers and wetland managers with tools to evaluate the functional capacity of wetlands to sustain biodiversity and ecosystem services, such as flood mitigation and water filtration. In addition, the inherent laziness of tidal networks suggests effective ways to adjust the drainage patterns of constructed wetlands in a nature-driven manner57. For example, network efficiency could be enhanced through targeted interventions, such as planting specific types of vegetation38 and controlling tidal forces58. Rather than serving as a definitive guideline, these findings represent a fundamental step towards more refined design and management of coastal wetlands. Further research is needed to explore the effects of tidal network drainage patterns on key ecological processes, including species distribution, habitat connectivity, and nutrient cycling. This will be critical for applying the “laziness” concept across diverse coastal wetland ecosystems and integrating it into practical management.

The impact of this study could extend beyond tidal channels, as it introduces a fundamental approach for controlling the space-filling procedure of branching networks. In our theoretical model, the concept of “drainage area” that drives the competition could be generalized to represent other critical resources for network growth across various contexts, such as sunlight for plant growth59, chemical gradients or mechanical tension for lung branching60, and high voltage potential for electrical breakdowns61. Exploring this competitive mechanism in different disciplines could open up various potential applications, including the optimization of agricultural irrigation systems62, biomimetic designs in engineering63, astrophysical studies of planetary surface morphology and erosion processes64, and enhancement in search algorithms and data mining techniques within complex databases65.

Methods

Extracting tidal networks from satellite imagery

We employed ArcMap software to manually trace the axial lines of tidal networks worldwide from Google Earth satellite imagery. The digitized data resulting from this effort were integrated with preexisting datasets to constitute the dataset of this study1,26,42. The visibility of tidal networks in satellite imagery is often obscured by vegetation, clouds, and floodwaters. Therefore, in our selection criteria, we prioritized regions with distinctly visible tidal networks, exemplified by regions such as Mokpo in South Korea and the Drowned Land of Saeftinghe in the Netherlands. Despite these constraints, the collected morphological data herein are the most extensive to date. In addition, the selected regions encompass diverse tidal dynamics and vegetation cover conditions (Table S1), facilitating the derivation of general patterns beyond case-specific mechanisms.

We used MATLAB to process the extracted tidal networks. The segmentation of the tidal network, as shown in Fig. 2a, follows the classic Hack method16: the longest channel reaching the network outlet is designated the main channel of order 1, and all tributaries joining an order-\(n\) channel are assigned order \(n+1\). This level-by-level hierarchical ordering process is shown in Supplementary Fig. 4. For instance, by excluding the main (i.e., the longest) channel of order 1, one can identify all subnetworks of order 2. The subsequent steps involved calculating the main channel length (\({L}_{m}\)), total channel length (\(\varSigma L\)), and associated drainage area (\({A}_{d}\)) of these decomposed subnetworks, thus allowing power-law fitting of the length–area scaling relationships.

We delineated the drainage area of each subnetwork by geometrically assigning each point on the drainage platform to its nearest channel, which differs from the common downslope direction approach employed for terrestrial river networks. This difference occurs because the drainage platforms of tidal networks are generally shallow and flat. In such cases, water typically flows into the nearest channel to minimize the frictional energy loss26. This delineation technique has been validated in previous morphometric studies34,42,43. In practice, we discretized the drainage platform into small square cells and assigned each cell to the channel line closest to the center of the cell at a straight-line distance (Supplementary Fig. 5). This process also aimed to determine the geometric unchanneled length (\({l}_{u}\)). We adopted different sizes of discrete cells based on the size of the drainage platform to ensure a magnitude of approximately \({10}^{6}\) cells, which usually satisfies the accuracy requirements.

To compute the smallest distances between the centers of the cells and the channel lines, the calculation of all distances and subsequent selection of the minimum value comprises a time-consuming process. Here, we employed the K-nearest neighbors algorithm embedded in MATLAB to efficiently determine these distances66.

Metrics derived from the network morphology

Here, we explain the calculation of the metrics that represent the branching morphology of tidal networks. We directly calculated the drainage density (\({D}_{d}=\Sigma L/{A}_{d}\)) and the space-filling deviation index (\({D}_{i}={\int }_{0}^{\infty }P({L}_{U}^{*}\ge {l}_{u}^{*}){{{\rm{d}}}}{l}_{u}^{*}-0.5\)) based on their definitions. The power-law scaling exponents of the channel lengths (\(\Sigma L\) and \({L}_{m}\)) versus the drainage area (\({A}_{d}\)) were derived through regression analysis utilizing the nonlinear least squares method, supported by the Curve Fitting Toolbox in MATLAB.

The sinuosity of a curvilinear channel can generally be quantified by the \(L/{L}_{E}\) ratio67, where \(L\) is the curvilinear channel length and \({L}_{E}\) is the straight-line distance between the start and end points. Here, we extended this definition and proposed the following equation for the average sinuosity of the channel network:

$${\sigma }_{s}=\frac{\Sigma L}{\Sigma {L}_{E}}$$
(4)

where \(\Sigma L\) is the total channel length and \(\Sigma {L}_{E}\) is the sum of the straight-line distances of all the channels ordered by Hack’s method (Supplementary Fig. 6).

The number of branches (\({N}_{b}\)) was determined by counting the number of the end tips of the channel network. The drainage area adjoining each end tip (\({A}_{{tip}}\)) was computed using the same technique as the watershed delineation described above (Fig. 5a). These two metrics were normalized to facilitate comparison between channel networks of different sizes and drainage densities. To achieve this normalization, we introduced a scaling factor \(\alpha\) to adjust the drainage density to 1. The key is that if the total length of the channel network is scaled by a factor of \(\alpha\), the corresponding drainage area is scaled by a factor of \({\alpha }^{2}\). Therefore, the factor can be determined as follows:

$${D}_{d}=\frac{\alpha \Sigma L}{{\alpha }^{2}{A}_{d}}=1$$
(5)

Solving this equation yields \(\alpha=\Sigma L/{A}_{d}\), which is exactly the drainage density \({D}_{d}\). Consequently, we defined the normalized tip area as \({A}_{{tip}}^{*}={D}_{d}^{2}{A}_{{tip}}\) and the normalized number of branches as \({N}_{b}^{*}={N}_{b}/({D}_{d}^{2}{A}_{d})\), which denotes the number of branches per unit area.

Algorithms of the theoretical model

In this section, we provide a detailed description of the theoretical model established in this research. Model coding and simulation were conducted using MATLAB. The competition principle, two-dimensional lattice, cell attributes (i.e., channel, platform and edge), and iteration process of the model have been described in the main text and are not repeated here.

The functional form of the erosion potential (Eq. 3) originated from the well-established sediment erosion model, which takes the following form46:

$${E}^{*}=\frac{{\varphi }_{s}}{{T}_{B}^{*}}{D}^{*1.5}\theta \exp \left[-\frac{8}{{\left(\theta /{\theta }_{c}\right)}^{0.7}}\right]$$
(6)

where \({E}^{*}\) is the sediment erosion rate, \({\varphi }_{s}\) is the ratio of the dry bulk density to the grain density of the sediment, \({T}_{B}^{*}\) is the nondimensional bursting period, \({D}^{*}\) is the nondimensional particle diameter, \(\theta\) is the Shields parameter, and \({\theta }_{c}\) is the critical Shields parameter. Equation (3) simplifies the complex empirical parameters in Eq. (6), while retaining the basic form of the product of the linear and exponential terms.

Here, we describe the algorithm for computing the drainage area that contributes to each edge cell in each iteration step in the model simulation process. Tidal networks develop in intertidal zones, where tidal movements are typically restricted to shallow water depths. This results in a flow field dominated by friction1, which allows the approximation of the flow paths across the drainage platform by straight lines (Supplementary Fig. 7a)34. We assumed that each straight flow line starts from a platform cell and directly extends to its nearest channel cell. Thus, these flow lines cross the edge cells and enter the channel cells from different angles. These angles can be categorized into eight main directions, namely, north, south, west, east, northeast, northwest, southeast, and southwest. Each direction corresponds to one of the eight neighboring cells surrounding a channel cell. We assessed the number of platform cells that contributed to each edge cell by counting the number of flow lines passing through each edge cell (Supplementary Fig. 7b). Subsequently, we sequentially calculated the drainage areas of the edge cells, the erosion potentials of the edge cells, and the probabilities that the edge cells are transformed into channel cells. The computation of the smallest distances between the flat and channel cells also entailed the K-nearest neighbors algorithm embedded in MATLAB66.

Scaling properties of \({{{{\boldsymbol{T}}}}}_{{{{\boldsymbol{square}}}}}\)

In this section, we explain the derivation of various scaling properties of \({T}_{{square}}\), which serves as the baseline network in our model and is mathematically equivalent to a regular space-filling fractal tree48. The iterative algorithm of \({T}_{{square}}\) can be understood as a recursive self-replicating procedure. For simplicity, we defined \({{{{\mathcal{T}}}}}_{{square}}\) as an infinite sequence of trees, i.e., \(\left\{{T}_{1},\,{T}_{2},\,{T}_{3},\cdots \right\}\). \({T}_{1}\) denotes a channel segment linking the bottom-left corner and the center of a \(1\times 1\) square (Supplementary Fig. 3a). To generate \({T}_{2}\), we first created four symmetric copies of \({T}_{1}\) by mirroring it horizontally, vertically and diagonally in sequence. Then, we extended the channel segment to the bottom-left corner of the newly formed \(2\times 2\) square. With the use of the same procedure, we generated \({T}_{3}\) from \({T}_{2}\), etc. Thus, the main channel length \({L}_{m}\), total channel length \(\Sigma L\) and drainage area \({A}_{d}\) of \({T}_{i}\) were derived as follows:

$${L}_{m}=\frac{\sqrt{2}}{2}\left({2}^{i}-1\right)$$
(7)
$$\Sigma L=\frac{\sqrt{2}}{6}\left({4}^{i}-1\right)$$
(8)
$${A}_{d}={4}^{i-1}$$
(9)

Adopting \(i\to \infty\), we can derive \({L}_{m}\propto {2}^{i}\), \(\Sigma L\propto {4}^{i}\) and \({A}_{d}\propto {4}^{i}\). This derivation demonstrates the scaling relationships of \({L}_{m}\propto {A}_{d}^{0.5}\) and \(\Sigma L\propto {A}_{d}\) for \({T}_{{square}}\). In addition, when \(i\to \infty\), the drainage density is \({D}_{d}=\Sigma L/{A}_{d}=2\sqrt{2}/3\), and the drainage area \({A}_{{tip}}\) adjacent to each channel tip is \(1/2\) as shown in Supplementary Fig. 3b. Thus, the normalized metric is as follows: \({A}_{{tip}}^{*}={A}_{{tip}}{D}_{d}^{2}=4/9\).

To determine the unchanneled length \({l}_{u}\), we decomposed the drainage basin into \(1\times 1\) square units. As shown in Supplementary Fig. 3a, these units can be categorized into two types: Unit A, in which the channel line links the center and corner; and Unit B, in which the channel line runs diagonally. The exceedance probability distributions of \({l}_{u}\) for these two units can be derived as follows:

$${P}_{A}\left({L}_{U}\ge {l}_{u}\right)=\left\{\begin{array}{c}{\left(\frac{\sqrt{2}}{2}-{l}_{u}\right)}^{2}+\dfrac{1}{2}-\dfrac{\pi }{2}{l}_{u}^{2},{l}_{u}\le \dfrac{1}{2} \hfill \\ {\left(\frac{\sqrt{2}}{2}-{l}_{u}\right)}^{2}+\dfrac{1}{2}-\dfrac{\pi }{2}{l}_{u}^{2}+{2l}_{u}^{2}\arccos \frac{1}{2{l}_{u}}-\sqrt{{l}_{u}^{2}-\frac{1}{4}},\frac{1}{2} < {l}_{u} \, \le \, \frac{\sqrt{2}\,}{2}\end{array}\right.$$
(10)
$${P}_{B}\left({L}_{U}\ge {l}_{u}\right)={\left(1-\sqrt{2}{l}_{u}\right)}^{2}$$
(11)

Regarding \({T}_{i}\), the values of Units A and B are \({4}^{i}/6+1/3\) and \({4}^{i-1}/3-1/3\), respectively. For \(i\to \infty\), these quantities converge to stable proportions, at \(2/3\) for Unit A and \(1/3\) for Unit B. Therefore, the exceedance probability of \({l}_{u}\) over the entire drainage basin is \(2{P}_{A}/3+{P}_{B}/3\) (Supplementary Fig. 3d). The associated space-filling deviation index \({T}_{{square}}\) is \({D}_{i}\approx 0.037\) based on Eq. (2).