A minimalist model for co-evolving supply and drainage networks

Numerous complex systems, both natural and artificial, are characterized by the presence of intertwined supply and/or drainage networks. Here we present a minimalist model of such co-evolving networks in a spatially continuous domain, where the obtained networks can be interpreted as a part of either the counter-flowing drainage or co-flowing supply and drainage mechanisms. The model consists of three coupled, nonlinear partial differential equations that describe spatial density patterns of input and output materials by modifying a mediating scalar field, on which supply and drainage networks are carved. In the 2-dimensional case, the scalar field can be viewed as the elevation of a hypothetical landscape, of which supply 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 supply and drainage networks occurs above a critical 'erosion' strength. The steady-state solutions are 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 congested extreme regimes, within which the resulting networks are characterized based on the absolute as well as the relative values of two non-dimensional indices.


Introduction
Many natural and man-made systems consist of materials being conveyed in and/or out of the domain through preferred routes, which result in the evolution of supply and/or drainage networks. In some biological systems, motile cells regulate their movement based on the affinity or aversion to specific environmental factors (temperature, chemical/biological signal) [1,2,3,4]. Two co-existing materials, moving up and down a signal gradient, * I am corresponding author Email addresses: skanand@princeton.edu (Shashank Kumar Anand), hooshyar@princeton.edu (Milad Hooshyar), jan.nordbotten@math.uib.no (Jan Martin Nordbotten), aporpora@princeton.edu (Amilcare Porporato ) drive the formation of the competing networks. In other systems, the material is supplied throughout a domain and gets collected once it has been utilized (and often also transformed), resulting in the formation of co-flowing supply and drainage networks. Examples include the cardiovascular network of blood and nutrients in animals, the supply-chain network of a commodity from the manufacturer to the customer and the related disposal, the aqueduct, and waste-flow network in urban water systems [5,6,7,8,9,10,11]. In all these systems, 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, the cost and modes of transportation for supply and drainage material, etc.
A great deal of research has explored the quantitative laws that explain the structure of networks in different disciplines, but this has been typically done considering either the supply or the drainage network separately [12,13,14,15,16]. In many cases, the general framework for studying such systems has been a static cost optimization problem typical of optimal transport theory [17,18,19,20]. As a result, the topology of the underlying supply or drainage network depends on the definition of the cost, including minimum energy dissipation, geometrical constraints, etc. [21,22,23,24]. Recently, there have been efforts to provide the interpretation of this static principle as the result of a dynamic evolution based on partial differential equation (PDE) [25,26,27,28].
Less efforts have been devoted to analyze the co-evolution of supply and/or drainage transport systems within a continuous domain, which is complicated by the presence of common and individual factors that affect transportation for both materials, including shape and 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 model that captures the essential interactions between two materials being conveyed in a continuous domain, where the system can be interpreted either as a counter-flowing drainage system or a co-flowing supply and drainage system. The model can be generalized to incorporate multi-species interplay; however, we keep the discussion up to two-species interactions.
The conceptual framework developed here stems from observing the complex ridges and valleys patterns in topographic landscapes, and the related work in the fields of image processing, geomorphology, and hydrology to formalize the duality between the interlocking network of ridges and valleys [29,30,31]. For mathematical formulization, we draw inspiration from landscape evolution models (LEMs) which have been successful in describing the 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 essential parameterization. We introduce a scalar field in a continuous domain that mediates two competing 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 boundary conditions on the obtained steady-state network patterns. The two channelization indices for both materials are obtained by the non-dimensionalization of the coupled PDEs, which allow us to define the role of various common and individual factors on the extent and patterns of the formed networks.
The paper is structured as follows. In Section 2, we first present the conceptual framework for both viewpoints of the model. We construct the 3-field mathematical model and define non-dimensional indices to describe the relative importance of various factors that alters the characteristics of the coupled networks. We also show that for unity value of exponents, the proposed model can be re-written as a 2-field model. The steady-state closed-form solutions for non-channelized flows in 2 and 3 dimensions are derived in Section 3. In Section 4, the numerical simulation results for the 2-dimensional and 3-dimensional cases are presented and the spatial patterns are analyzed for different levels of complexity and branching. Conclusions and future research directions are discussed in Section 5. In Appendix A, we discuss the 2-field equivalent formulation for the proposed PDE model and the complexity in the boundary conditions that emerges from this model reduction.

Conceptual model inspired by the ridge and valley duality
We begin by considering the geometry of a topographic field (Figure 1a), visualized as a scalar field (h) in 3-dimensional Euclidean space, where the vertical direction points in the direction of gravity. Avoiding maxima (summits) and minima (pits), the curves of particular significance here are the ridges and valleys, which provide a skeleton for the structure of the drainage network [35,36]. With the assumption of negligible inertial effects, a fluid present in the domain flows under gravity over the scalar field along the direction of the steepest descent, 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 the stream network and flows out of the system at the boundary.
Inverting artificially the initial topography, as shown in Figure 1c, the duality between ridges and valleys is apparent, as ridges become valleys and valleys become ridges. The interlocked network of ridge-lines and valley-lines extracted from the original topography is shown in Figure 1e. Based on this duality, and similarly to the density of a − for the drained material, one can imagine another flow with density a + in this inverted topography, which is produced within the domain and gets drained by the ridge network. The field of material density (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 the ridge/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 be formulated as a co-flowing supply-drainage problem (Problem II), where a + represents the density 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 the boundary through the ridge network, a + enters from the boundary where the ridge forms peak, and flows along the ridge network following the topographic steepest descent. Figure  1f 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 dominated by either material. One can envision the supplied material (density, a + ) enter the area at the boundary concentrated at the ridges of the scalar field (h), flowing and getting distributed over the hillslopes as it gets exhausted. In turn, the consumption of the supplied material produces the drained material (density, a − ), which moves under the scalar field potential, and gets discharged out of the domain preferentially via the valleys.

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 two materials (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 For simplicity, we assume that the velocity fields, v + and v − , follow the positive and negative gradient of h, respectively, with unit speed as 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 and source terms due to the feedback from both materials as where D is the diffusion coefficient, K > 0, m ± > 0 and n ± > 0 are model parameters and r + /r − indicate the production rates for the respective material. The coupled nonlinear Equations (1) and (3) form a closed system for the interaction of counter-flowing drainage mechanisms by modifying the scalar field (h) with appropriate initial and boundary conditions for h, a + and a − . In this work, we consider 2-dimensional and 3-dimensional domains in the shape of a rectangle 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 remaining side edges/faces (Ω s ) follow zero Neumann boundary conditions in h, which provide closed boundary conditions in a ± . The proposed arrangement induces a directionality to the movement of the two materials in the domain with top (Ω t ) and bottom (Ω b ) edges/faces functioning as the exit boundaries for a + and a − , respectively. Assuming that the densities of the two materials are negligible at their upstream domain boundaries, the boundary conditions become simple and time-independent as a + (Ω b ) = a − (Ω t ) = 0. Under such boundary conditions and the assumption of spatially uniform production rates (r + and r − ), the governing equations compute the counter flow of the materials across the domain (including at the 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 drained out 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 input material supplied at the boundaries, that moves with velocity v + and gets consumed at the unitary rate, as where both velocity fields, v + and v − , follow the negative gradient of h with unit speed as Equations (3) and (4) form a general minimalist model for the interaction of two underlying mechanisms of supply and drainage in a spatially continuous domain by modifying the scalar field (h), as apparent from Equation (3), where r + now represents the consumption rate of the supplied material (a + ). For parsimony, we here assume that the supply is consumed uniformly in space at constant rate, which is immediately disposed giving rise to a uniform and constant source of material that gets drained. More complicated patterns of supply 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 consider the viewpoint of supply and drainage mechanisms (Problem II) from this point for the interpretation of the solutions.
The sink and source terms mathematically formalize the conceptual framework shown in 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 an hypothetical landscape over which input and output materials move following Equation (5). As indicated by Equation (3), the accumulation of drainage material decreases the elevation that results in the formation of valleys (sink term). Conversely, the aggregation of the supply material 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 is concentrated in valleys.
In the 3-dimensional case, the scalar field can be interpreted as the strength of a chemical signal 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 opposite to its gradient. Vascularization of the supply and drainage networks takes place in the domain with high material density of the supply material in high-valued scalar field region and high material density of the drainage material in low-valued scalar field region due to the feedback of sink and source terms in Equation (3). The mathematical structure of the proposed model bears some resemblance with more complex models of vasculogenesis and chemotaxis [3,37]. Specifically, the core component of the model resembles minimalist versions of the well-known KellerSegel model for chemotaxis under the negligible diffusion of 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 ) and bottom (Ω 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 supply material is conveyed out of Ω b , the boundary conditions of the material densities remain simple 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 the domain (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 − ). c.
a. 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 as − ∇ · a * ∇h |∇h| = −1.
Equation (3) then can be re-written using Equations (6) and (7) as where 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 of the new spatial field (a * ) for the reduced model require the knowledge of a + and a − in advance to obtain the same solution given by the 3-field model with the time-independent boundary conditions. The reader is referred to Appendix A, where this problem of the boundary condition for the 2-field model is discussed in detail for simulation results in the 2-dimensional case.

Non-dimensionalization
For a typical value H of the scalar field and a typical length scale of the domain L, the following dimensionless quantities are established: Using these quantities, Equations (3) and (4) can be written in dimensionless form, where This shows that the overall behavior of the system can be described by the two 'channelization indices', C I + and C I − . For a constant value of exponents in Equation (9), an increase in 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 the feedback of the sink term, keeping all other factors the same. This mechanism results in a correlation of the density of the two materials to the value of the scalar field at steady-state which can be visualized by looking at the level set L c (h) of the scalar field h for a constant value c. High density of input/output material accruing on the different level sets of the scalar field is shown in the steady-state solutions of the 2-dimensional and 3-dimensional cases (Section 4).

Closed-form solution
At steady-state, the closed-form solution can be obtained for the case where diffusion in Equation (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 top edge, which is at a fixed higher elevation (H) compared to the bottom edge, is separated by the distance L from the bottom edge. This situation is analogous to the flow of two materials 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 compared to the fixed zero chemical signals strength at the bottom face, which drives the flow of the materials.
Smooth profiles using Equation (13) and the corresponding slope variations following Equation (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 the surface attains a higher profile compared to the case for C I + = 0.

Numerical solutions
Numerical experiments are started for the 2-dimensional and 3-dimensional cases with a linear initial condition containing a small amount of random spatial noise. We limit our discussion to the case with unity exponents (m ± = n ± = 1), and utilize the efficient algorithm presented in [41] to update the scalar field h over the entire domain until the steady-state is reached. The fundamental concept behind this algorithm is inspired from the notion of a flow network of the material over the entire scalar field, which is traversed in a way to make the matrix system upper/lower triangular for the efficient implicit computation. The accuracy of the numerical algorithm has been carefully tested for the case of drainage-network evolution model for the natural landscape against analytical solutions in non-channelized/vascularized conditions, as well as against analytical results of the onset of linear stability analysis [11,41] and with exact mean field solutions obtained in the condition of fully channelized/vascularized regime [41]. We refer to these references for further details.

2-dimensional case 4.1.1. Code Verification
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 ± to verify the implemented code. The first channelization in the domain occurs at C I ± = 3.5. The closed-form solution is applicable for the cases when the field h is smooth enough (no channelization). The mean surface profile starts deviating from the closed-form solution for C I ± ≥ 3.5 (Figure 4(a)) due to channel formation. This is apparent from the accumulation plot of 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 onset of first channelization) and 12.5, the interface becomes a curve due to the emergence of channels in supply and drainage networks (Figure 4 (b,c,d)).

Effect of individual factors
In this numerical experiment, we focus on the impact of individual factors on the couplednetwork formations. Varying the relative rate of consumption (r + ) of the supplied material versus the generation rate (r − ) of the drained material affects the feedback on the scalar field, which in turn affects the structure of the coupled networks. Having all other parts parameters the same, the variation in these rates can be expressed as changing the values of non-dimensional indices C I + and C I − . We study the extent and spatial patterns of a − and a + for 55 cases with C I ± ∈ [50, 500] and C I + ≤ C I − . Figure 5 presents simulation results, where the supply network is represented in red (highlighting high density region of a + , i.e., a + > a − ) and the drainage network is shown in blue color (accentuating high density region of a − , i.e., a − > a + ) with the white curve representing the interface a + = a − . These spatial networks that evolve for various values 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 where the number of channels of the supply and drainage network is high, with mostly straight channels and very little branching. Panels (c,e,g) present the plots where comparatively less number of channels are observed with more branching . Panels (d,f,h) show the plots where maximum branching is observed with curved channels compared to previous other cases. c. d. 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 the 3-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 in the formation of shallow ridge and valley patterns. Panels (b) and (d) show the surface plot for the branched region of C I + = 100, C I − = 400 and C I + = 50, C I − = 250 where high value of C I − compared to C I + results in branched channels of the networks with wide valleys and thin ridges at the steady-state.
The interplay between model parameters and boundary conditions becomes apparent when considering panels f and g. These cases have the same total C I + + C I − for the 2field 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.
This shows that the non-trivial boundary conditions allowed by the three-field model (as opposed two-field model) influences the solution throughout the domain, both quantitatively and qualitatively. A full discussion of this is included in Appendix A.
a. 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 the surface 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, which increases 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 together to form branched patterns. For C I − = 500, changing C I + = 50 to C I + = 500 increases the length 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 redcolored region in Figure 8(A). Conversely, for disproportionate values of C I ± , the interface length (L i ) is relatively smaller as shown in the blue-colored region in Figure 8(A). We define a quantity N c which refers to the maximum number of either main supply or drainage channels of length greater than the half of the width of the domain (50 in this case) originating from the boundaries of the domain. N c is plotted for 55 cases of various values of C I ± as Figure 8(B), which looks similar to the plot of L i as expected. High values of N c occur for large and similar values of C I ± again shown as the red-colored region. More branching results in less number of main channels for C I + << C I − , as indicated by the blue-colored in Figure 8(B). The scatter plot of L i versus N c with best-fit line having correlation coefficient r = 0.988 reconfirms the close relationship between number of main channels and the interface length (Figure 8(C)).
The simulation results shown in Figure 5 can be mapped to different regions in the colorplot of the contour length L c . Panels (a,b) in Figure 5 belong to the red-colored region in Figure 8(A), where a high density of nearly unswerving main channels with a few offshoots is observed. For this reason, we refer to it as the congested region. Panels (d,f,h) in Figure 5 exhibit the plots for the blue (branched) area in Figure 8(A), with heavily branched channels. Panels (c,e,g) in Figure 5 display the plots from the yellow/green (transient) region, which lies within these two extremes.

3-dimensional case
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 chemical signal's strength). There are fixed boundary conditions for two faces (h(x, 0, z) = H = 10 and h(x, 80, z) = 0) and zero Neumann boundary conditions at the remaining faces. a − is zero at h(x, 0, z) = H = 10 and a + is zero at h(x, 80, z) = 0, with closed boundary conditions in 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 plots for the field h are drawn on the two side faces (closed boundary conditions in a ± ) along with contour 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 witnessed in the 2-dimensional results. As shown in Figure 9(a), a large number of red-colored contour curves (high density region of h) in cross-section near the face h(x, 0, z) remain in crosssection near the face h(x, 80, z), which is dominated by the blue-colored contour curves (low density region of h). This pattern for C I ± = 1000 resembles the equal strength of shallow ridges and valleys in the 2-dimensional case. We display the largest drainage conduit from the steady-state solution in Figure 9(c,d), where green-colored haze in panel c indicates the points in the domain from which the flow is collected in the given conduit. c. Similarly, the contour patterns for C I + = 200 and C I − = 1000 on the faces resemble the thin ridges and wide valleys obtained in the 2-dimensional branched case when the values of C I + and C I − are disproportionate (Figure 9(b)). This is apparent as the tiny red-colored contour curves (high density region of h) in the cross-section near the face h(x, 0, z) vanish in the cross-section near the face h(x, 80, z) dominated by the blue-colored contour curves (low density region of h). This parallels the thin ridges that start from the fixed elevation end (h = H) and disappear near the fixed elevation end (h = 0) surrounded by deep and wide valleys, leading to the branched supply and drainage networks.

Conclusions
The minimalist model developed in this work leads to the formation of spatial patterns of combined supply and drainage networks in a continuous domain, whereby the corrugations of a 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 intensity of the branching of these networks and the instability in the profile of h. The crucial role of the boundary condition for these coupled PDEs is particularly evident when reducing the presented 3-field model to a 2-field model for unit values of the exponents in source and sink terms, as the achieved simplification in the number of equations entails a complication in the boundary conditions, which is necessary to solve the same co-existing supply and drainage networks of the 3-field model.
While we limit our discussion here to unit exponents of the source and sink terms in Equation (9), the solutions for non-unitary values of the exponents have qualitatively similar features and reflect an analogous spectrum of branched versus congested regime after the first channelization for a different range of C I ± . However, the specific results do depend on these nonlinearities and the source and sink terms: future work will be devoted to adjusting 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 complexity of the implicit solver by making the matrix system upper/lower triangular. This has been a vital improvement in 2-dimensional cases, where the space complexity of the algorithm is not an issue [41]. However, the simulations in the 3-dimensional domain require a large amount of memory space compared to the 2-dimensional cases due to the increased input size of nodes 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 numerical solver so that the coupled patterns for a 3-dimensional domain can be analyzed in more depth.
equations 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 channelization indices that are defined for the 3-field model. We discuss here the dependency of complex boundary conditions of a * on the solution of spatial fields a + and a − by presenting steady-state solutions for a 2-dimensional square domain 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 −3 , K = 10 −5 and (r + + r − ) = 5, two cases were simulated as r + = 1, r − = 4 (C I + = 100, C I − = 400), and r ± = 2.5 (C I ± = 250). Figures A.10(a,b) show the plots of steady-state supply and drainage material densities (a + and a − ) for the two cases. The difference in the obtained supply and drainage networks can be interpreted as the role of different values of C I ± in the 3-field model. However, the two cases correspond to the same value of K * = 5 × 10 −5 for the 2-field model, which indicates the crucial role of time-dependent boundary condition of a * on the obtained supply and drainage networks.
For 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 condition for a − is fixed in time throughout the simulations with a − = 0 at the top edge of the square domain. For the 2-field model, the boundary condition for a * in the two cases is different and is defined by specifying individual values of r + and r − initially, as shown in Figure  A.10(c,d). The blue curves for both cases, representing a * = 0, vary in time, indicating the contribution of time-dependent boundary condition of a * on the simulation results. This dependency is further shown in Figure A.11, where the value of a * at top and bottom edges of the square domain at different time-steps are displayed for both cases. The different values of a * in time at domain boundaries indicate that, the steady-state solutions with the same value of parameters (K * = 5 × 10 −5 ) are different because of the distinct time-varying boundary conditions for a * . Therefore, the model can be simulated using the two fields of supply and drainage density with simple boundary conditions for the densities of the materials. This way, the results can be interpreted in simple terms as the interplay of 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, the obtained co-existing networks can only be constructed from each of the two fields from which the sum a * originates.