A minimalist model for co-evolving supply and drainage networks
Shashank Kumar Anand, Milad Hooshyar, Jan Martin Nordbotten, Amilcare Porporato
AA minimalist model for co-evolving supply and drainage networks
Shashank Kumar Anand a , Milad Hooshyar b,c , Jan Martin Nordbotten d , AmilcarePorporato a,b, ∗ a Department of Civil and Environmental Engineering, Princeton University, USA b Princeton Environmental Institute, Princeton University, Princeton, USA c Princeton Institute for International and Regional Studies, Princeton University, USA d Department of Mathematics, University of Bergen, Bergen, Norway
Abstract
Numerous complex systems, both natural and artificial, are characterized by the presence ofintertwined supply and/or drainage networks. Here we present a minimalist model of suchco-evolving networks in a spatially continuous domain, where the obtained networks can beinterpreted as a part of either the counter-flowing drainage or co-flowing supply and drainagemechanisms. The model consists of three coupled, nonlinear partial differential equationsthat describe spatial density patterns of input and output materials by modifying a medi-ating scalar field, on which supply and drainage networks are carved. In the 2-dimensionalcase, the scalar field can be viewed as the elevation of a hypothetical landscape, of whichsupply and drainage networks are ridges and valleys, respectively. In the 3-dimensional case,the scalar field serves as the chemical signal strength, in which vascularization of the supplyand drainage networks occurs above a critical ‘erosion’ strength. The steady-state solutionsare presented as a function of non-dimensional channelization indices for both materials.The spatial patterns of the emerging networks are classified within the branched and con-gested extreme regimes, within which the resulting networks are characterized based on theabsolute as well as the relative values of two non-dimensional indices.
Keywords: optimal transport, coupled networks, pattern formation, nonlinearity
1. Introduction
Many natural and man-made systems consist of materials being conveyed in and/orout of the domain through preferred routes, which result in the evolution of supply and/ordrainage networks. In some biological systems, motile cells regulate their movement based onthe affinity or aversion to specific environmental factors (temperature, chemical/biologicalsignal) [1, 2, 3, 4]. Two co-existing materials, moving up and down a signal gradient, ∗ I am corresponding author
Email addresses: [email protected] (Shashank Kumar Anand), [email protected] (Milad Hooshyar), [email protected] (Jan Martin Nordbotten), [email protected] (Amilcare Porporato ) a r X i v : . [ n li n . AO ] A ug rive the formation of the competing networks. In other systems, the material is suppliedthroughout a domain and gets collected once it has been utilized (and often also trans-formed), resulting in the formation of co-flowing supply and drainage networks. Examplesinclude the cardiovascular network of blood and nutrients in animals, the supply-chain net-work of a commodity from the manufacturer to the customer and the related disposal, theaqueduct, and waste-flow network in urban water systems [5, 6, 7, 8, 9, 10, 11]. In all thesesystems, the co-existing networks must evolve or be designed in a way that is coordinated,depending on different constraints, such as the configuration of the distribution region, thecost and modes of transportation for supply and drainage material, etc.A great deal of research has explored the quantitative laws that explain the structureof networks in different disciplines, but this has been typically done considering either thesupply or the drainage network separately [12, 13, 14, 15, 16]. In many cases, the generalframework for studying such systems has been a static cost optimization problem typicalof optimal transport theory [17, 18, 19, 20]. As a result, the topology of the underlyingsupply or drainage network depends on the definition of the cost, including minimum energydissipation, geometrical constraints, etc. [21, 22, 23, 24]. Recently, there have been effortsto provide the interpretation of this static principle as the result of a dynamic evolutionbased on partial differential equation (PDE) [25, 26, 27, 28].Less efforts have been devoted to analyze the co-evolution of supply and/or drainagetransport systems within a continuous domain, which is complicated by the presence ofcommon and individual factors that affect transportation for both materials, including shapeand size of the region, extra scalar/vector field, production/consumption rate, velocity, etc.As a step in this direction, this study aims at formulating and analyzing a minimalist modelthat captures the essential interactions between two materials being conveyed in a continuousdomain, where the system can be interpreted either as a counter-flowing drainage systemor a co-flowing supply and drainage system. The model can be generalized to incorporatemulti-species interplay; however, we keep the discussion up to two-species interactions.The conceptual framework developed here stems from observing the complex ridges andvalleys patterns in topographic landscapes, and the related work in the fields of image pro-cessing, geomorphology, and hydrology to formalize the duality between the interlockingnetwork of ridges and valleys [29, 30, 31]. For mathematical formulization, we draw inspi-ration from landscape evolution models (LEMs) which have been successful in describingthe formation of river and stream networks [32, 33, 11, 34]. Generalizing these models,we develop a simple system consisting of three nonlinear coupled PDEs with the essentialparameterization. We introduce a scalar field in a continuous domain that mediates twocompeting mechanisms of two counter-flowing drainage or co-flowing supply and drainage.We show the influence of rules of production and/or consumption as well as the boundaryconditions on the obtained steady-state network patterns. The two channelization indicesfor both materials are obtained by the non-dimensionalization of the coupled PDEs, whichallow us to define the role of various common and individual factors on the extent andpatterns of the formed networks.The paper is structured as follows. In Section 2, we first present the conceptual frame-work for both viewpoints of the model. We construct the 3-field mathematical model and2efine non-dimensional indices to describe the relative importance of various factors thatalters the characteristics of the coupled networks. We also show that for unity value ofexponents, the proposed model can be re-written as a 2-field model. The steady-stateclosed-form solutions for non-channelized flows in 2 and 3 dimensions are derived in Section3. In Section 4, the numerical simulation results for the 2-dimensional and 3-dimensionalcases are presented and the spatial patterns are analyzed for different levels of complexityand branching. Conclusions and future research directions are discussed in Section 5. InAppendix A, we discuss the 2-field equivalent formulation for the proposed PDE model andthe complexity in the boundary conditions that emerges from this model reduction.
2. Mathematical model
We begin by considering the geometry of a topographic field (Figure 1a), visualized as ascalar field ( h ) in 3-dimensional Euclidean space, where the vertical direction points in thedirection of gravity. Avoiding maxima (summits) and minima (pits), the curves of particularsignificance here are the ridges and valleys, which provide a skeleton for the structure of thedrainage network [35, 36]. With the assumption of negligible inertial effects, a fluid presentin the domain flows under gravity over the scalar field along the direction of the steepestdescent, resulting in the distribution of material density, say a − , as shown in Figure 1b,highlighting the drainage network for the topography. This density, a − , is drained by thestream network and flows out of the system at the boundary.Inverting artificially the initial topography, as shown in Figure 1c, the duality betweenridges and valleys is apparent, as ridges become valleys and valleys become ridges. Theinterlocked network of ridge-lines and valley-lines extracted from the original topography isshown in Figure 1e. Based on this duality, and similarly to the density of a − for the drainedmaterial, one can imagine another flow with density a + in this inverted topography, whichis produced within the domain and gets drained by the ridge network. The field of materialdensity ( a + ), marking the drainage network for the flipped topography, is shown in Figure 1d,where the main courses of flow follow the ridge-lines of the original topography. Therefore,the flow of a + /a − moving up/down the slope of the topographic field to be drained by theridge/valley network becomes the counter-flow problem (we will refer to this as Problem I).Reversing the flow direction of the density a + in above scenario, the problem can beformulated as a co-flowing supply-drainage problem (Problem II), where a + represents thedensity of the supplied material that flows down the slope similar to the drained material( a − ). Instead of having a distributed source through the domain and exiting from theboundary through the ridge network, a + enters from the boundary where the ridge formspeak, and flows along the ridge network following the topographic steepest descent. Figure1f displays the accumulation of a + along ridges (red-colored region) and accumulation of a − along the valleys (blue-colored region), with the white curve subdividing regions dominatedby either material. One can envision the supplied material (density, a + ) enter the area at theboundary concentrated at the ridges of the scalar field ( h ), flowing and getting distributedover the hillslopes as it gets exhausted. In turn, the consumption of the supplied material3roduces the drained material (density, a − ), which moves under the scalar field potential,and gets discharged out of the domain preferentially via the valleys. x [𝑚] y [𝑚] 0 500 − ℎ 𝐻 y [𝑚] x [𝑚] y [𝑚] ℎ 𝐻 l o g 𝑎 − l o g 𝑎 + a. b.c. d.e. f. x [𝑚]0500 y [𝑚]0500 y [𝑚]0500 y [𝑚] x [𝑚]0 x [𝑚]0 x [𝑚] Figure 1: Conceptualization of supply and drainage networks using the dual ridge and valley networks in atopographic landscape. (a): 3-dimensional surface ( h ) for the selected topography. (b): Drainage networkfor a − following the negative gradient of h . (c): The inverted 3D surface for the original topography, where H is the maximum elevation in the domain. (d): Drainage network for a + following the negative gradientof the inverted h . (e) Interlocked planar ridge and valley networks with prominent ridge-lines (red) andvalley-lines (blue). (f) The white curve represents the interface a + = a − and separates the red-coloredregion representing high accumulation of a + ( a + > a − ) from the blue-colored region showing aggregationof a − ( a − > a + ). The selected topography is from the Calhoun Critical Zone landscape in South Carolina(obtained from the the OpenTopography facility). .2. Governing equations The 2-dimensional illustration presented above can be formalized and extended to an n -dimensional space ( R n ), considering a scalar field h : R n → R , defined inside a domain Ωalong with two scalar fields a + and a − playing the role of the densities of materials.For the counter-flow drainage problem (Problem I), the continuity equation for the twomaterials ( a + and a − ), that are produced at a unitary rate and flow with opposite velocity v + and v − , respectively, can be written under the assumption of quasi steady-state as ∇ · ( a ± v ± ) = 1 , (1)For simplicity, we assume that the velocity fields, v + and v − , follow the positive andnegative gradient of h , respectively, with unit speed as v ± = ± ∇ h |∇ h | . (2)The scalar field h is assumed to co-evolve with the fields of both materials. Specifically,the temporal evolution of h consists of a diffusion term and nonlinear, nonlocal sink andsource terms due to the feedback from both materials as ∂h∂t = D ∇ h + K ( r + a + ) m + |∇ h | n + − K ( r − a − ) m − |∇ h | n − , (3)where D is the diffusion coefficient, K > m ± > n ± > r + / r − indicate the production rates for the respective material. The coupled nonlinearEquations (1) and (3) form a closed system for the interaction of counter-flowing drainagemechanisms by modifying the scalar field ( h ) with appropriate initial and boundary condi-tions for h , a + and a − .In this work, we consider 2-dimensional and 3-dimensional domains in the shape of arectangle or parallelepiped, respectively, with the top edge/face (Ω t ) at a fixed higher value( h = H ) compared to the bottom edge/face (Ω b ) at a fixed lower value ( h = 0). The re-maining side edges/faces (Ω s ) follow zero Neumann boundary conditions in h , which provideclosed boundary conditions in a ± . The proposed arrangement induces a directionality to themovement of the two materials in the domain with top (Ω t ) and bottom (Ω b ) edges/facesfunctioning as the exit boundaries for a + and a − , respectively. Assuming that the densitiesof the two materials are negligible at their upstream domain boundaries, the boundary con-ditions become simple and time-independent as a + (Ω b ) = a − (Ω t ) = 0. Under such boundaryconditions and the assumption of spatially uniform production rates ( r + and r − ), the gov-erning equations compute the counter flow of the materials across the domain (including atthe boundaries where the densities are not specified i.e., a + (Ω t ) and a − (Ω b )).From the viewpoint of the co-flowing supply and drainage mechanism (Problem II), a + represents the density of the input material in the domain, which is utilized and drainedout of the domain as the output material with density a − . The continuity equation for a − ,therefore, remains the same, with the modification in the continuity equation for the input5aterial supplied at the boundaries, that moves with velocity v + and gets consumed at theunitary rate, as ∇ · ( a ± v ± ) = ∓ , (4)where both velocity fields, v + and v − , follow the negative gradient of h with unit speed as v ± = − ∇ h |∇ h | . (5)Equations (3) and (4) form a general minimalist model for the interaction of two un-derlying mechanisms of supply and drainage in a spatially continuous domain by modifyingthe scalar field ( h ), as apparent from Equation (3), where r + now represents the consump-tion rate of the supplied material ( a + ). For parsimony, we here assume that the supply isconsumed uniformly in space at constant rate, which is immediately disposed giving rise toa uniform and constant source of material that gets drained. More complicated patterns ofsupply and drainage are certainly of interest and will be investigated in the future work.The model can be analyzed from either of two discussed formulations; however, we con-sider the viewpoint of supply and drainage mechanisms (Problem II) from this point for theinterpretation of the solutions.The sink and source terms mathematically formalize the conceptual framework shownin Figure 1, where the movement of materials carves out the preferential paths. Thus,for the 2-dimensional case the scalar field ( h ) may be viewed as an elevation field of anhypothetical landscape over which input and output materials move following Equation (5).As indicated by Equation (3), the accumulation of drainage material decreases the elevationthat results in the formation of valleys (sink term). Conversely, the aggregation of the supplymaterial increases the surface elevation that leads to the formation of ridges (source term).Consequently, the input material is accumulated on ridges, while the output material isconcentrated in valleys.In the 3-dimensional case, the scalar field can be interpreted as the strength of a chemicalsignal that drives the movement of the materials (chemotaxis). As Equation (5) indicates,the concentration of the chemical signal ( h ) stimulates the migration of the materials op-posite to its gradient. Vascularization of the supply and drainage networks takes place inthe domain with high material density of the supply material in high-valued scalar fieldregion and high material density of the drainage material in low-valued scalar field regiondue to the feedback of sink and source terms in Equation (3). The mathematical structureof the proposed model bears some resemblance with more complex models of vasculogenesisand chemotaxis [3, 37]. Specifically, the core component of the model resembles minimalistversions of the well-known KellerSegel model for chemotaxis under the negligible diffusionof biological cells [38, 39, 40, 4].The boundary conditions play an important role in the proposed model. For Problem II,the boundary conditions for h are the same as discussed for Problem I, with top (Ω t ) andbottom (Ω b ) edges/faces functioning now as the entry and exit boundaries for the domain,respectively. Under the assumption that no drainage material exits from Ω t and no supplymaterial is conveyed out of Ω b , the boundary conditions of the material densities remain6imple and time-independent as a + (Ω b ) = a − (Ω t ) = 0. Under such boundary conditions,the governing equations determine the flow of supplied and drained materials across thedomain (including at the boundaries where the densities are not specified i.e., a + (Ω t ) and a − (Ω b )) under the assumption of spatially uniform consumption ( r + ) and production rates( r − ). Ω 𝑠 𝑎 ± = Ω 𝑡 𝑎 − = 0Ω 𝑏 𝑎 + = 0 Ω 𝑠 𝑎 ± = l o g 𝑎 − Ω 𝑏 d.c.a. ℎ𝐻 l o g 𝑎 + Ω 𝑡 b. Figure 2: Schematic representation of the boundary conditions used in the model. (a): Surface profile ofthe scalar field ( h ) in a portion of the rectangular domain near the first channel instability, where H isthe maximum elevation value in the domain (see Section 4.1 for details). Three (shallow) ridges and two(shallow) valleys can be observed in the plotted profile. (c): Boundary conditions of a + and a − counter-flowdrainage problem (Problem I) and co-flow supply and drainage problem (Problem II). Two (red-colored)channels of a + and three (blue-colored) channels of a − corresponding to the ridges and valleys in panel (a)are observed. The white curve, representing the interface a + = a − , separates the regions dominated by theeither material. Four contour lines of the scalar field ( h ) are plotted along with black-colored streamlineswhich indicate the flow direction of the materials. (b,d): Obtained signals of a + and a − at Ω t and Ω b ,respectively, with peaks indicating channel formation at the corresponding domain boundaries. It is interesting to observe that the model can be reduced to a 2-field system for the case m ± = n ± = 1. Multiplying the continuity equations for a + and a − (Equation (4)) with r + and r − , and subtracting, one can write the single equation for a new spatial field a ∗ = r + a + − r − a − r + + r − , (6)as − ∇ · (cid:18) a ∗ ∇ h |∇ h | (cid:19) = − . (7)Equation (3) then can be re-written using Equations (6) and (7) as ∂h∂t = D ∇ h + K ∗ a ∗ |∇ h | , (8)7here K ∗ = K ( r + + r − ). Equations (7) and (8) form a 2-field equivalent formulation ( a ∗ , h )to the proposed 3-field model ( a + , a − , h ) for unit values of the exponents in Equation (3).The achieved simplification is, in practice, only apparent as the boundary conditions ofthe new spatial field ( a ∗ ) for the reduced model require the knowledge of a + and a − inadvance to obtain the same solution given by the 3-field model with the time-independentboundary conditions. The reader is referred to Appendix A, where this problem of theboundary condition for the 2-field model is discussed in detail for simulation results in the2-dimensional case. For a typical value H of the scalar field and a typical length scale of the domain L , thefollowing dimensionless quantities are established: ˆ h = hH , ˆ a + = a + L , ˆ a − = a − L , ˆ t = L D , ˆ x = xL and ˆ y = yL . Using these quantities, Equations (3) and (4) can be written in dimensionlessform, ∂ ˆ h∂ ˆ t = ˆ ∇ ˆ h + C I + ˆ a m + + | ˆ ∇ ˆ h | n + − C I − ˆ a m − − | ˆ ∇ ˆ h | n − , (9) − ˆ ∇ · (cid:32) ˆ a ± ˆ ∇ ˆ h | ˆ ∇ ˆ h | (cid:33) = ∓ , (10)where C I + = Kr m + + L m + − n + DH − n + , C I − = Kr m − − L m − − n − DH − n − . (11)This shows that the overall behavior of the system can be described by the two ‘channel-ization indices’, C I + and C I − . For a constant value of exponents in Equation (9), an increasein the value of C I + by high consumption rate, r + , enhances the feedback of the source term.On the other hand, a rise in the value of C I − by high production rate, r − , strengthens thefeedback of the sink term, keeping all other factors the same. This mechanism results in acorrelation of the density of the two materials to the value of the scalar field at steady-statewhich can be visualized by looking at the level set L c ( h ) of the scalar field h for a constantvalue c . High density of input/output material accruing on the different level sets of thescalar field is shown in the steady-state solutions of the 2-dimensional and 3-dimensionalcases (Section 4).
3. Closed-form solution
At steady-state, the closed-form solution can be obtained for the case where diffusion inEquation (9) inhibits the instability formation in the scalar field. In the 2-dimensional case,it can be visualized as the smooth elevation field in a semi-infinite domain where the topedge, which is at a fixed higher elevation ( H ) compared to the bottom edge, is separatedby the distance L from the bottom edge. This situation is analogous to the flow of twomaterials before vascularization across two infinite parallel plates placed at a finite distance L in 3-dimension, with a fixed high chemical signal’s strength ( H ) at the top face compared8o the fixed zero chemical signals strength at the bottom face, which drives the flow of thematerials.Assuming that the scalar field (ˆ h ) decreases monotonically in the 1D transect, Equation(10) can be solved with the boundary conditions ˆ a + (ˆ y = 1) = ˆ a − (ˆ y = 0) = 0 to obtainˆ a + = (1 − ˆ y ) and ˆ a − = ˆ y . For the case of m ± = n ± = 1, substituting the expressions for ˆ a + and ˆ a − in Equation (9) at steady-state can be written asˆ h (cid:48)(cid:48) + C I − ˆ yh (cid:48) − C I + (1 − ˆ y )ˆ h (cid:48) = 0 . (12)Solving Equation (12) givesˆ h = erf (cid:16) C I− √ C I− + C I + ) (cid:17) − erf (cid:16) C I− ˆ y −C I + (1 − ˆ y ) √ C I + + C I− ) (cid:17) erf (cid:16) C I− √ C I− + C I + ) (cid:17) + erf (cid:16) C I + √ C I− + C I + ) (cid:17) (13) | ˆ h (cid:48) | = e − ( CI− ˆ y −CI +(1 − ˆ y ))22( CI− + CI +) erf (cid:16) C I− √ C I− + C I + ) (cid:17) + erf (cid:16) C I + √ C I− + C I + ) (cid:17) × (cid:114) C I − + C I + ) π . (14) |ℎ ′ | ℎ ො𝑦 → a. b. 𝒞 𝔩 − = 𝒞 𝔩 + = 25 𝒞 𝔩 + = 0 𝒞 𝔩 − = 0 ො𝑦 → Figure 3: (a,b): Steady-state solutions given by Equation (13) and (14) for three cases of C I + = 0, C I − = 0and C I ± = 25. Smooth profiles using Equation (13) and the corresponding slope variations followingEquation (14) for C I + = 0, C I − = 0 and C I + = C I − = 25 are displayed in Figure 3 (a,b).As expected, for C I − = 0, the contribution from the nonlinear sink term goes away and thesurface attains a higher profile compared to the case for C I + = 0.
4. Numerical solutions
Numerical experiments are started for the 2-dimensional and 3-dimensional cases witha linear initial condition containing a small amount of random spatial noise. We limit9ur discussion to the case with unity exponents ( m ± = n ± = 1), and utilize the efficientalgorithm presented in [41] to update the scalar field h over the entire domain until thesteady-state is reached. The fundamental concept behind this algorithm is inspired fromthe notion of a flow network of the material over the entire scalar field, which is traversedin a way to make the matrix system upper/lower triangular for the efficient implicit com-putation. The accuracy of the numerical algorithm has been carefully tested for the case ofdrainage-network evolution model for the natural landscape against analytical solutions innon-channelized/vascularized conditions, as well as against analytical results of the onset oflinear stability analysis [11, 41] and with exact mean field solutions obtained in the conditionof fully channelized/vascularized regime [41]. We refer to these references for further details. We first simulate a rectangular domain with high aspect ratio (length = 500, width =100) and compare the mean elevation profile along the length for varying values of C I ± toverify the implemented code. The first channelization in the domain occurs at C I ± = 3 . h is smooth enough (nochannelization). The mean surface profile starts deviating from the closed-form solution for C I ± ≥ . a + and a − , where the white curve represents the interface for a − = a + . For C I ± = 1,the interface is a straight line, while for C I ± = 3 . .
5, the interface becomes a curve due to the emergence of channels in supply and drainagenetworks (Figure 4 (b,c,d)). d. 𝒞 𝔩 ± = 𝑥 →𝑥 → Supply network
Drainage network a. ො𝑦 → 𝒞 𝔩 ± = 1𝒞 𝔩 ± = 3.5 𝒞 𝔩 ± = 12.5 ℎ , ഥ ℎ 𝐻 → b. 𝒞 𝔩 ± = 1 c. 𝒞 𝔩 ± = 𝑦 ↓𝑦 ↓ 𝑦 ↓ Figure 4: (a): Solid lines represent the computed mean surface profile (¯ h ) along the length for a rectangulardomain (width = 100, length = 500) for C I ± = 1, 3 . . a + (red-colored) and a − (blue-colored) for C I ± = 1, 3 . .
5. The white curve represents the interface a + = a − , separating the two regions dominated by the either material. .1.2. Effect of individual factors In this numerical experiment, we focus on the impact of individual factors on the coupled-network formations. Varying the relative rate of consumption ( r + ) of the supplied materialversus the generation rate ( r − ) of the drained material affects the feedback on the scalarfield, which in turn affects the structure of the coupled networks. Having all other partsparameters the same, the variation in these rates can be expressed as changing the valuesof non-dimensional indices C I + and C I − . b. 𝒞 𝔩 − = 500, 𝒞 𝔩 + = 350 a. 𝒞 𝔩 − = 500, 𝒞 𝔩 + = 500 c. 𝒞 𝔩 − = 500, 𝒞 𝔩 + = 200 d. 𝒞 𝔩 − = 500, 𝒞 𝔩 + = 50 e. 𝒞 𝔩 − = 400, 𝒞 𝔩 + = 300 f. 𝒞 𝔩 − = 400, 𝒞 𝔩 + = 100 g. 𝒞 𝔩 − = 250, 𝒞 𝔩 + = 250 h. 𝒞 𝔩 − = 250, 𝒞 𝔩 + = 50 Supply network Drainage network Figure 5: Simulation results for various values of C I + and C I − in a rectangular domain (width = 100,length = 500). The accumulation of the input material a + is represented in red (highlighting a + > a − ),while the accumulation of the output material a − is shown in blue (highlighting a + < a − ). The white curverepresents the interface a + = a − and separates region dominated by the input material from the outputmaterial. The grey-colored boxes on the panels e, f, g and h indicate the regions of the domain for which3-dimensional surface profiles of the elevation field ( h ) are shown in Figure 6. We study the extent and spatial patterns of a − and a + for 55 cases with C I ± ∈ [50 , C I + ≤ C I − . Figure 5 presents simulation results, where the supply network is representedin red (highlighting high density region of a + , i.e., a + > a − ) and the drainage network isshown in blue color (accentuating high density region of a − , i.e., a − > a + ) with the whitecurve representing the interface a + = a − . These spatial networks that evolve for various11alues of C I ± are quite distinctive, indicating the role of absolute as well as relative values of C I + and C I − on the overall pattern formation. Panels (a,b) of Figure 5 display the plots wherethe number of channels of the supply and drainage network is high, with mostly straightchannels and very little branching. Panels (c,e,g) present the plots where comparatively lessnumber of channels are observed with more branching. Panels (d,f,h) show the plots wheremaximum branching is observed with curved channels compared to previous other cases. a. b. c. d. Figure 6: 3-dimensional surface profiles of h for the regions highlighted in Figure (5). (a): C I − = 400, C I + = 300. (b): C I − = 400, C I + = 100. (c): C I − = 250, C I + = 250. (d): C I − = 250, C I + = 50. Comparablevalues of C I + and C I − result in the formation of shallow ridges and valleys with less branching, while thedisproportionate values of C I + and C I − result in wide and more branched valleys with sharp ridges. The relative strength of C I + and C I − affecting the shape of the surface ( h ) and hence,the spatial patterns of both networks is apparent from Figure 6. The figure displays the3-dimensional surface plots of h from the selected regions in Figure 5. Panels (a) and (c)display the surface plots where the comparable opposing strength of a + and a − results inthe formation of shallow ridge and valley patterns. Panels (b) and (d) show the surface plotfor the branched region of C I + = 100, C I − = 400 and C I + = 50, C I − = 250 where high valueof C I − compared to C I + results in branched channels of the networks with wide valleys andthin ridges at the steady-state.The interplay between model parameters and boundary conditions becomes apparentwhen considering panels f and g. These cases have the same total C I + + C I − for the 2-field model, and thus, both solutions satisfy the same differential equations (7) and (8).Nevertheless, as is apparent in Figure 5, resulting branched structures are vastly different.12his shows that the non-trivial boundary conditions allowed by the three-field model (asopposed two-field model) influences the solution throughout the domain, both quantitativelyand qualitatively. A full discussion of this is included in Appendix A. a. y↓
0 x → b.
0 x → y↓ Figure 7: Plot of the interface a + = a − for ( C I + = 500, C I − = 500) and ( C I + = 50, C I − = 500). The variety of patterns in above cases can be explained by the structure of Equation(9) for the 3-field model. Increasing values of C I ± hike the tendency of a + / a − to mold thesurface as ridge and valley respectively. For very high and comparable values of C I − and C I + , the primary channels get stuck, can not coalesce to form branched channels. Therefore,the number of main channels (channels originating from the boundaries) increases, whichincreases the length of the interface a + = a − . Reducing the value of C I + with respect to C I − results in more branching as channels of a − dominate the space and coalesce togetherto form branched patterns. For C I − = 500, changing C I + = 50 to C I + = 500 increases thelength of interface ( a + = a − ) as shown in Figure 7.We plot the length of the interface a + = a − (denoted by L i ) for various values of C I ± .High values of L i occur for large and comparable values of C I ± which is shown as the red-colored region in Figure 8(A). Conversely, for disproportionate values of C I ± , the interfacelength ( L i ) is relatively smaller as shown in the blue-colored region in Figure 8(A). We definea quantity N c which refers to the maximum number of either main supply or drainage chan-nels of length greater than the half of the width of the domain (50 in this case) originatingfrom the boundaries of the domain. N c is plotted for 55 cases of various values of C I ± asFigure 8(B), which looks similar to the plot of L i as expected. High values of N c occurfor large and similar values of C I ± again shown as the red-colored region. More branchingresults in less number of main channels for C I + << C I − , as indicated by the blue-coloredin Figure 8(B). The scatter plot of L i versus N c with best-fit line having correlation coeffi-cient r = 0 .
988 reconfirms the close relationship between number of main channels and theinterface length (Figure 8(C)).The simulation results shown in Figure 5 can be mapped to different regions in the color-plot of the contour length L c . Panels (a,b) in Figure 5 belong to the red-colored region inFigure 8(A), where a high density of nearly unswerving main channels with a few offshootsis observed. For this reason, we refer to it as the congested region. Panels (d,f,h) in Figure 5exhibit the plots for the blue (branched) area in Figure 8(A), with heavily branched channels.13anels (c,e,g) in Figure 5 display the plots from the yellow/green (transient) region, whichlies within these two extremes. C. 𝑁 𝑐 𝐿 𝑖
30 40 50 60 70 𝒞 𝔩 − 𝒞 𝔩 + 𝒞 𝔩 − 𝑁 𝑐 B. 𝒞 𝔩 + 𝐿 𝑖 d c ba ef gh A. Figure 8: (A): Color-plot of the interface length, L i , for various values of C I ± ∈ (50 , , , ... N c , for variousvalues of C I + and C I − ∈ (50 , , , ... L c vs N c with the best-fit line (Correlationcoefficient r = 0 . We apply the proposed model to a 3-dimensional domain for a parallelepiped ( x =50 , y = 80 , z = 60), where h now refers to a density field (can be viewed as a chemicalsignal’s strength). There are fixed boundary conditions for two faces ( h ( x, , z ) = H = 10and h ( x, , z ) = 0) and zero Neumann boundary conditions at the remaining faces. a − iszero at h ( x, , z ) = H = 10 and a + is zero at h ( x, , z ) = 0, with closed boundary conditionsin a ± for the remaining four faces. We explore two cases keeping C I − = 1000, while changing C I − from 1000 to 200. The simulation results are shown in Figure 9, where the contour plotsfor the field h are drawn on the two side faces (closed boundary conditions in a ± ) along withcontour line plots for the two cross-sections near the faces along y -axis.The steady-state solutions for the 3-dimensional cases agree with the patterns witnessedin the 2-dimensional results. As shown in Figure 9(a), a large number of red-colored contour14urves (high density region of h ) in cross-section near the face h ( x, , z ) remain in cross-section near the face h ( x, , z ), which is dominated by the blue-colored contour curves (lowdensity region of h ). This pattern for C I ± = 1000 resembles the equal strength of shallowridges and valleys in the 2-dimensional case. We display the largest drainage conduit fromthe steady-state solution in Figure 9(c,d), where green-colored haze in panel c indicates thepoints in the domain from which the flow is collected in the given conduit. ← x ←y ℎ𝐻 a.b. y → x → l o g 𝑎 − d. c. Figure 9: Simulation results for the 3-dimensional domain ( x = 50 , y = 80 , z = 60). Contour plot presentingthe scalar h at the side boundaries and at the two cross-sections along the y -axis near the respective ends(top and bottom face) for (a): C I + = 1000 and C I − = 1000. (b): C I + = 200 and C I − = 1000. (c,d): Thelargest drainage conduit extracted in the domain for the case of C I ± = 1000 with light green colored haze(in panel c) indicating the points in the domain from which the flow of drainage material is received in theconduit. Similarly, the contour patterns for C I + = 200 and C I − = 1000 on the faces resemble thethin ridges and wide valleys obtained in the 2-dimensional branched case when the valuesof C I + and C I − are disproportionate (Figure 9(b)). This is apparent as the tiny red-coloredcontour curves (high density region of h ) in the cross-section near the face h ( x, , z ) vanishin the cross-section near the face h ( x, , z ) dominated by the blue-colored contour curves(low density region of h ). This parallels the thin ridges that start from the fixed elevationend ( h = H ) and disappear near the fixed elevation end ( h = 0) surrounded by deep andwide valleys, leading to the branched supply and drainage networks.15 . Conclusions The minimalist model developed in this work leads to the formation of spatial patterns ofcombined supply and drainage networks in a continuous domain, whereby the corrugations ofa mediating scalar field, h , cleave these competing networks in the same continuous domain.A channelization index ( C I ± ) corresponding to each material governs the relative intensityof the branching of these networks and the instability in the profile of h . The crucial role ofthe boundary condition for these coupled PDEs is particularly evident when reducing thepresented 3-field model to a 2-field model for unit values of the exponents in source and sinkterms, as the achieved simplification in the number of equations entails a complication in theboundary conditions, which is necessary to solve the same co-existing supply and drainagenetworks of the 3-field model.While we limit our discussion here to unit exponents of the source and sink terms inEquation (9), the solutions for non-unitary values of the exponents have qualitatively similarfeatures and reflect an analogous spectrum of branched versus congested regime after thefirst channelization for a different range of C I ± . However, the specific results do dependon these nonlinearities and the source and sink terms: future work will be devoted toadjusting them to cater to specific applications, as has been done in various other models,such as the minimalist versions of the well-known KellerSegel model for chemotaxis [4, 38],mechanochemical models of angiogenesis and vasculogenesis [2, 37].From the numerical point of view, the employed algorithm decreases the time complexityof the implicit solver by making the matrix system upper/lower triangular. This has been avital improvement in 2-dimensional cases, where the space complexity of the algorithm is notan issue [41]. However, the simulations in the 3-dimensional domain require a large amountof memory space compared to the 2-dimensional cases due to the increased input size ofnodes and the corresponding auxiliary space utilized by the algorithm during the execution.This increases the overall computational cost of the simulations in the 3-dimensional cases.A part of the future work, therefore, is to reduce the space complexity of the numericalsolver so that the coupled patterns for a 3-dimensional domain can be analyzed in moredepth. Acknowledgment
The authors acknowledge support from the US National Science Foundation (NSF) grantsEAR-1331846 and EAR-1338694, and BP through the Carbon Mitigation Initiative (CMI)at Princeton University. A.P. and M.H. also acknowledge the support from the PrincetonInstitute for International and Regional Studies (PIIRS) and the Princeton EnvironmentalInstitute (PEI). J. M. N. was supported in part through Norwegian Research Council grantnumber 250223.The authors are pleased to acknowledge that the simulations presented in this article wereperformed on computational resources managed and supported by Princeton Research Com-puting, a consortium of groups including the Princeton Institute for Computational Scienceand Engineering (PICSciE) and the Office of Information Technology’s High PerformanceComputing Center and Visualization Laboratory at Princeton University.16he code used for the simulations is available at https://github.com/ShashankAnand1996/Supply-Drainage . References [1] E. M. Hedgecock, R. L. Russell, Normal and mutant thermotaxis in the nematode caenorhabditiselegans, Proceedings of the National Academy of Sciences 72 (10) (1975) 4061–4065. doi:10.1073/pnas.72.10.4061 .[2] D. Manoussaki, S. Lubkin, R. Vemon, J. Murray, A mechanical model for the formation of vascularnetworks in vitro, Acta biotheoretica 44 (3-4) (1996) 271–282. doi:10.1007/BF00046533 .[3] M. A. Chaplain, Mathematical modelling of angiogenesis, Journal of neuro-oncology 50 (1-2) (2000)37–51. doi:10.1023/A:1006446020377 .[4] T. Hillen, K. J. Painter, A users guide to pde models for chemotaxis, Journal of mathematical biology58 (1-2) (2009) 183. doi:10.1007/s00285-008-0201-3 .[5] A. Rinaldo, I. Rodriguez-Iturbe, R. Rigon, E. Ijjasz-Vasquez, R. Bras, Self-organized fractal rivernetworks, Physical review letters 70 (6) (1993) 822. doi:10.1103/PhysRevLett.70.822 .[6] T. Sun, P. Meakin, T. Jøssang, Minimum energy dissipation river networks with fractal boundaries,Physical Review E 51 (6) (1995) 5353. doi:10.1103/PhysRevE.51.5353 .[7] G. C. Dandy, A. R. Simpson, L. J. Murphy, An improved genetic algorithm for pipe network optimiza-tion, Water resources research 32 (2) (1996) 449–458. doi:10.1029/95WR02917 .[8] S. Abdinnour-Helm, Network design in supply chain management, International Journal of Agile Man-agement Systems doi:10.1108/14654659910280929 .[9] S. Suweis, M. Konar, C. Dalin, N. Hanasaki, A. Rinaldo, I. Rodriguez-Iturbe, Structure and controlsof the global virtual water trade network, Geophysical Research Letters 38 (10). doi:2011GL046837 .[10] D. Mahlke, A. Martin, S. Moritz, A simulated annealing algorithm for transient optimization ingas networks, Mathematical Methods of Operations Research 66 (1) (2007) 99–115. doi:10.1007/s00186-006-0142-9 .[11] S. Bonetti, M. Hooshyar, C. Camporeale, A. Porporato, Channelization cascade in landscape evolution,Proceedings of the National Academy of Sciences doi:10.1073/pnas.1911817117 .[12] J. R. Banavar, A. Maritan, A. Rinaldo, Size and form in efficient transportation networks, Nature399 (6732) (1999) 130–132. doi:10.1038/20144 .[13] M. Hooshyar, A. Singh, D. Wang, Hydrologic controls on junction angle of river networks, WaterResources Research 53 (5) (2017) 4073–4083. doi:10.1002/2016WR020267 .[14] M. Minoux, Networks synthesis and optimum network design problems: Models, solution methods andapplications, Networks 19 (3) (1989) 313–360. doi:10.1002/net.3230190305 .[15] J. R. Banavar, F. Colaiori, A. Flammini, A. Maritan, A. Rinaldo, Topology of the fittest transportationnetwork, Physical Review Letters 84 (20) (2000) 4745. doi:PhysRevLett.84.4745 .[16] H. Ronellenfitsch, E. Katifori, Global optimization, local adaptation, and the role of growth in distribu-tion networks, Physical review letters 117 (13) (2016) 138301. doi:10.1103/PhysRevLett.117.138301 .[17] S. N. Dorogovtsev, J. F. Mendes, Evolution of networks, Advances in physics 51 (4) (2002) 1079–1187. doi:10.1080/00018730110112519 .[18] S. Bohn, M. O. Magnasco, Structure, scaling, and phase transition in the optimal transport network,Physical review letters 98 (8) (2007) 088702. doi:10.1103/PhysRevLett.98.088702 .[19] B. Danila, Y. Yu, J. A. Marsh, K. E. Bassler, Optimal transport on complex networks, Physical ReviewE 74 (4) (2006) 046106. doi:10.1103/PhysRevE.74.046106 .[20] M. Durand, Architecture of optimal transport networks, Physical Review E 73 (1) (2006) 016116. doi:10.1103/PhysRevE.73.016116 .[21] K. Chen, Simple learning algorithm for the traveling salesman problem, Physical Review E 55 (6)(1997) 7809. doi:10.1103/PhysRevE.55.7809 .[22] I. Rodr´ıguez-Iturbe, A. Rinaldo, Fractal river basins: chance and self-organization, Cambridge Univer-sity Press, 2001.
23] C. D. Murray, The physiological principle of minimum work: I. the vascular system and the cost ofblood volume, Proceedings of the National Academy of Sciences of the United States of America 12 (3)(1926) 207. doi:10.1073/pnas.12.3.207 .[24] S. T. Rachev, The monge–kantorovich mass transference problem and its stochastic applications, Theoryof Probability & Its Applications 29 (4) (1985) 647–676. doi:10.1137/1129093 .[25] L. C. Evans, Partial differential equations and monge-kantorovich mass transfer, Current developmentsin mathematics 1997 (1) (1997) 65–126.[26] E. Facca, F. Cardin, M. Putti, Towards a stationary monge–kantorovich dynamics: The physarumpolycephalum experience, SIAM Journal on Applied Mathematics 78 (2) (2018) 651–676. doi:10.1137/16M1098383 .[27] F. Cardin, J. R. Banavar, A. Maritan, Optimal transport from a point-like source, Continuum Mechan-ics and Thermodynamics (2019) 1–11 doi:10.1007/s00161-019-00844-5 .[28] M. Hooshyar, S. Anand, A. Porporato, Variational analysis of landscape elevation and drainage net-works, Proceedings of the Royal Society A 476 (2239) (2020) 20190775. doi:10.1098/rspa.2019.0775 .[29] J. J. Koenderink, A. J. van Doorn, Local features of smooth shapes: Ridges and courses, in: GeometricMethods in Computer Vision II, Vol. 2031, International Society for Optics and Photonics, 1993, pp.2–13. doi:10.1117/12.146617 .[30] C. Werner, Several duality theorems for interlocking ridge and channel networks, Water ResourcesResearch 27 (12) (1991) 3237–3247. doi:10.1029/91WR02322 .[31] W. R. Dawes, D. Short, The significance of topology for modeling the surface hydrology of fluviallandscapes, Water Resources Research 30 (4) (1994) 1045–1055. doi:10.1029/93WR02479 .[32] A. Fowler, Mathematical geoscience, Vol. 36, Springer Science & Business Media, 2011. doi:10.1007/978-0-85729-721-1 .[33] E. Istanbulluoglu, R. L. Bras, Vegetation-modulated landscape evolution: Effects of vegetation onlandscape processes, drainage density, and topography, Journal of Geophysical Research: Earth Surface110 (F2). doi:10.1029/2004jf000249 .[34] J. T. Perron, P. W. Richardson, K. L. Ferrier, M. Lapˆotre, The root of branching river networks, Nature492 (7427) (2012) 100. doi:10.1038/nature11672 .[35] J. J. Koenderink, A. J. Van Doorn, The structure of relief, in: Advances in imaging and electronphysics, Vol. 103, Elsevier, 1998, pp. 65–150. doi:10.1016/S1076-5670(08)70015-6 .[36] S. Bonetti, A. Bragg, A. Porporato, On the theory of drainage area for regular and non-regular points,Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 474 (2211) (2018)20170693. doi:10.1098/rspa.2017.0693 .[37] D. Manoussaki, A mechanochemical model of angiogenesis and vasculogenesis, ESAIM: MathematicalModelling and Numerical Analysis 37 (4) (2003) 581–599. doi:10.1051/m2an:2003046 .[38] E. F. Keller, L. A. Segel, Initiation of slime mold aggregation viewed as an instability, Journal oftheoretical biology 26 (3) (1970) 399–415. doi:10.1016/0022-5193(70)90092-5 .[39] T. Hillen, K. Painter, C. Schmeiser, Global existence for chemotaxis with finite sampling radius, Discrete& Continuous Dynamical Systems-B 7 (1) (2007) 125. doi:10.3934/dcdsb.2007.7.125 .[40] M. Alber, N. Chen, T. Glimm, P. M. Lushnikov, Multiscale dynamics of biological cells with chemotacticinteractions: from a discrete stochastic model to a continuous description, Physical Review E 73 (5)(2006) 051901. doi:10.1103/PhysRevE.73.051901 .[41] S. K. Anand, M. Hooshyar, A. Porporato, Linear layout of multiple flow-direction networks forlandscape-evolution simulations, Environmental Modelling & Software (2020) In press. Preprint avail-able at https://arxiv.org/abs/1909.03176 . Appendix A. 2-field model
In Section 2.2, we show that the original 3-field model can be reduced to a 2-field modelconsisting of Equations (7) and (8) for the new spatial field, a ∗ , and the scalar field, h ,under the assumption of unit exponents of the source and sink term in Equation (3). These18quations form a closed system where the dynamics of h depends on the parameter K ∗ ,which is determined by the summation of r + and r − only, instead of the two channelizationindices that are defined for the 3-field model. b.a. ↓ d.c. 𝑥 → 𝑥 →𝑦 ↓ Figure A.10: (a,b): Steady-state solutions for C I + = 100, C I − = 400, and C I + = 250, C I − = 250, respectively.The accumulation of a + is in red ( a + > a − ), the accumulation of a − is in blue ( a + < a − ) and the whitecurve is interface a + = a − . (c,d): a ∗ = 0 at t = 0 (blue), t = intermediate (green) and t = steady state(orange) for cases (a) and (b), respectively. We discuss here the dependency of complex boundary conditions of a ∗ on the solutionof spatial fields a + and a − by presenting steady-state solutions for a 2-dimensional squaredomain with top edge ( y = 0) at fixed high elevation ( H = 10) and bottom edge ( y = L =100) at fixed zero elevation, with zero Neumann boundary conditions on the side edges.With the same values of D = 10 − , K = 10 − and ( r + + r − ) = 5, two cases were simulatedas r + = 1, r − = 4 ( C I + = 100, C I − = 400), and r ± = 2 . C I ± = 250). Figures A.10(a,b)show the plots of steady-state supply and drainage material densities ( a + and a − ) for thetwo cases. The difference in the obtained supply and drainage networks can be interpretedas the role of different values of C I ± in the 3-field model. However, the two cases correspondto the same value of K ∗ = 5 × − for the 2-field model, which indicates the crucial role oftime-dependent boundary condition of a ∗ on the obtained supply and drainage networks.19or the 3-field model, time-independent boundary condition for a + is well defined, with a + = 0 at the bottom edge ( h = 0) for the both cases. Similarly, the boundary conditionfor a − is fixed in time throughout the simulations with a − = 0 at the top edge of the squaredomain. For the 2-field model, the boundary condition for a ∗ in the two cases is differentand is defined by specifying individual values of r + and r − initially, as shown in FigureA.10(c,d). The blue curves for both cases, representing a ∗ = 0, vary in time, indicating thecontribution of time-dependent boundary condition of a ∗ on the simulation results. d.c. b.a. 𝑙 𝑜 𝑔 𝑎 ∗ 𝑙 𝑜 𝑔 𝑎 ∗ 𝑙 𝑜 𝑔 𝑎 ∗ 𝑙 𝑜 𝑔 𝑎 ∗ Figure A.11: Blue ( t = 0), green ( t = intermediate) and orange ( t = steady state) curves represent thevalue of a ∗ at the domain boundaries at different time-steps. Panels (a,b): a ∗ at the top edge for C I + = 100, C I − = 400, and C I + = 250, C I − = 250, respectively. Panels (c,d): a ∗ at the bottom edge for C I + = 100, C I − = 400, and C I + = 250, C I − = 250, respectively. This dependency is further shown in Figure A.11, where the value of a ∗ at top andbottom edges of the square domain at different time-steps are displayed for both cases. Thedifferent values of a ∗ in time at domain boundaries indicate that, the steady-state solutionswith the same value of parameters ( K ∗ = 5 × − ) are different because of the distincttime-varying boundary conditions for a ∗ . Therefore, the model can be simulated using the20wo fields of supply and drainage density with simple boundary conditions for the densitiesof the materials. This way, the results can be interpreted in simple terms as the interplayof two indices of the supply and drainage density fields. If the 2-field model is employed,the time-dependent boundary conditions for a ∗ are extremely complex, and in practice, theobtained co-existing networks can only be constructed from each of the two fields from whichthe sum a ∗∗