Consistent Initial Conditions for the Saint-Venant Equations in River Network Modeling

1 Initial conditions for flows and depths (cross-sectional areas) throughout a river network are 2 required for any time-marching (unsteady) solution of the one-dimensional (1D) hydrodynamic 3 Saint-Venant equations. For a river network modeled with several Strahler orders of tributaries, 4 comprehensive and consistent synoptic data are typically lacking and synthetic starting con5 ditions are needed. Because of underlying nonlinearity, poorly-defined or inconsistent initial 6 conditions can lead to convergence problems and long spin-up times in an unsteady solver. Two 7 new approaches are defined and demonstrated herein for computing flows and cross-sectional ar8 eas (or depths). These methods can produce an initial condition data set that is consistent with 9 modeled landscape runoff and river geometry boundary conditions at the initial time. These new 10 methods are: (1) the Pseudo-Time-Marching Method (PTM) that iterates toward a steady-state 11 initial condition using an unsteady Saint-Venant solver, and (2) the Steady-Solution Method 12


Motivation
Setting initial conditions for unsteady simulations of the Saint-Venant equations (SVEs) across large river networks can be challenging.Every element of the river network must be given initial values of flow and depth, and these values should be consistent with the inflow boundary conditions (e.g., from a land surface model) at the starting time to prevent instabilities.This issue has not been previously addressed in the literature, arguably because (i) adequate initial conditions are fairly trivial for small SVE reaches, (ii) hydrological models with large river networks often do not use the full SVE (e.g., Beighley et al., 2009) or use it over a smaller set of reaches (e.g., Paiva et al., 2012), and (iii) timemarching models only consider results after spin-up time is complete (i.e., after the effects of the initial conditions have been washed out of the system solution), which implies the initial conditions are entirely irrelevant in analyzing the model results.However, for SVE solutions of river networks, the initial conditions can dramatically affect both spin-up time and convergence.Indeed, in our experience, näive initial conditions can cause the spin-up time to be longer than the time-marching period of interest.In extreme cases, this can result in divergence and complete failure of an SVE solver.Consistency between the initial conditions and the boundary conditions appears to be necessary for short spin-up times, and thus is the focus of this study.Note that the need for consistency in an SVE model initialization is due to the coupling of nonlinearity and the water surface slope in the momentum equation, so common reduced-physics models (not Published by Copernicus Publications on behalf of the European Geosciences Union.discussed herein) may not be as sensitive to consistent initial conditions.
Saint-Venant equation modeling arguably dates from Preissmann's seminal work (Preissmann, 1961;Preissmann and Cunge, 1961), followed by decades of advances in techniques and applications (Cunge, 1974;Ponce et al., 1978;Cunge et al., 1980;Abbott et al., 1986;Zhao et al., 1996;Sanders, 2001;Pramanik et al., 2010).These models focused on hydraulics of short river reaches or main stem rivers that are easy to initialize for flow and depth.It is only recently that the solvers for large river networks have become practical (Hodges, 2013;Liu and Hodges, 2014), and it is with large networks that initial conditions are problematic.Indeed, initial conditions and associated spin-up problems have been recently acknowledged and investigated for hydrological models (e.g., Ajami et al., 2014) but without consideration of a separate river network model.Work by Seck et al. (2015) and Rahman and Lu (2015) show that hydrological model spin-up computational times could be significant and were dominated by the selected initial hydrological conditions.
Our experience with Saint-Venant river network modeling is that simple approaches to initial conditions often cause localized numerical instabilities, slow convergence of the timemarching numerical solution, and long model spin-up times.Herein, we investigate the initial condition problem for a Saint-Venant river network model for a given set of inflows from a hydrological model.

Synthetic vs. observed initial conditions
Model spin-up time is completed when the effects of initial conditions cannot be observed in the model results.Thus, by definition, the initial conditions cannot affect the unsteady solution beyond spin-up time.It follows that initial conditions are irrelevant to the quality of the time-marching simulation and are only important in how they affect the spin-up duration.We can imagine a "perfect" set of initial conditions with zero spin-up time, which would require initial flows and depths consistent with (i) the actual unsteady behavior prior to the model start time and (ii) the model boundary conditions; the latter includes both the bathymetric model for the river channels and the coupled hydrological model providing runoff and baseflows.Such perfect initial conditions are practically unattainable due to the sparsity of synoptic flow/depth data as well as unavoidable uncertainty and errors in both bathymetric and hydrological models.Another way to think of this is that perfect initial conditions also require perfect boundary conditions (perfect bathymetry and hydrology), or else some spin-up time is required to wash out inconsistencies.In general, the spin-up time will be affected by how far the initial conditions are from the theoretical perfect conditions.
The key point is that the exact observed river initial conditions (if such were available throughout a network) will not eliminate or necessarily reduce spin-up time if the observed data are inconsistent with the model boundary conditions.Similarly, interpolations of sparse synoptic data will not be a priori consistent with the boundary conditions and thus cannot eliminate spin-up time.Inconsistency is a critical concept: the mismatches between the initial conditions and the boundary conditions can lead to unrealistic destabilizing impulses in time marching the SVE solution.Such impulses can require extensive spin-up time to damp their effects.An extreme example is a high runoff rate into an almost dry stream that can cause a Gibbs phenomenon at a wave front and negative values for the computed cross-sectional area (Lax, 2006;Kvočka et al., 2015;Yang et al., 2012).Although several studies show such numerical discontinuities can be resolved (Kazolea and Delis, 2013;Caleffi et al., 2003;Liang et al., 2006), the high computational cost of damping or resolving is a burden (Kvočka et al., 2015) that seems unnecessary during spin-up since it cannot affect the time-marching results.
We argue that the primary goal of initial conditions is providing consistency with the boundary conditions to allow smooth, convergent spin-up of an unsteady solver.This task can be accomplished with synthetic initial conditions that are independent of observations.The only practical discriminators between using observed and synthetic initial conditions are (i) the effort required to prepare the initial condition data and (ii) the length of spin-up time.
As demonstrated herein, consistent synthetic initial conditions can be readily generated for even large complex networks -a task that is daunting for interpolation/extrapolation of sparse synoptic observations.Indeed, developing consistent synthetic initial conditions only requires the channel geometry and hydrological model that are used for the SVE time marching, supplemented by a steady-state SVE solver (Sect.2.3).In contrast, interpolation/extrapolation of synoptic data requires analysis of the data locations and model geometry, which is likely to require customization for each river network.
Unfortunately, we cannot definitively prove our second discriminator; the wide range of possible hydrological models, channel geometry models, and synoptic observation systems makes it impossible to prove that initial conditions based upon observations will always require longer spin-up times.However, no one has proposed an approach for interpolation/extrapolation of observed synoptic data that provides simultaneous consistency with a hydrological model, a river geometry model, and the SVE.It follows that we are on firm ground in stating that any existing approach using observed data as SVE initial conditions will result in inconsistencies.Finally, as our work demonstrates that inconsistent initial conditions cause long spin-up times in SVE models of sufficiently large river networks, it can be argued that observed data should be deprecated as initial conditions until a consistent approach is demonstrated and the spin-up can be shown to be shorter than provided by synthetic initial conditions.

Initial condition approaches
Approaches for specifying initial conditions for the SVE can be grouped into three main categories: (i) a "synoptic start" applying an interpolated/extrapolated set of sparse observational data, (ii) a "cold start" with initial flow rates and flow depths prescribed either as zero (e.g., Chau and Lee, 1991) or from some analytical values, e.g., mean annual flows and depths, and (iii) a "steady-state" start, which we describe herein.The metric for evaluating initial conditions is not how well they reflect available real-world observations but how effective they are in efficiently providing a consistent set of initial conditions.
Based on our discussion above, the first approach (synoptic start) is unlikely to be efficient for SVE initial conditions in a large river network due to inconsistencies between observations and model boundary conditions as well as inconsistencies caused by interpolating/extrapolating sparse observations throughout a network.There are no proven approaches to analyzing consistency and melding observations to hydrological model runoff, so the river network model spin-up will be subject to random inconsistencies and instabilities that can delay or prevent convergence.
The second approach, a cold start, provides innumerable possible ways to create initial conditions.For example, mean annual flows and depths (e.g., from the NHDPlus data in the US) can provide a smooth and consistent set of flows and elevations throughout a network.Although such cold start initial conditions can be internally consistent, they may be far from the flows/depths implied by the initial hydrological forcing.For example, a river network model that is started with mean annual values would be substantially in error if the initial hydrological inflows were from the monsoon season.As a result of inconsistencies between the selected cold start values and the hydrological inflows, a cold start can require extensive spin-up time to dilute or wash out the error.Indeed, the spin-up time dominated the computational time for the large SVE networks that we previously modeled in Liu and Hodges (2014) when we used a cold start with mean annual values.It might be possible to design a cold start approach that is consistent with the hydrological inflows; however, we suspect that any such approach is likely to be merely a variant of the steady-state approach discussed herein.
Herein, we investigate the third approach, steady-state initial conditions, as a preferred method for initializing a large river network model.With this idea, a set of consistent initial conditions is one that satisfies both the t = 0 hydrological forcing and the steady-state Saint-Venant equations at t = 0.That is, we know that we cannot match the unsteady SVE at t = 0 as we cannot perfectly know the time-varying nonlinear effects before the model starts (unless, of course, we run a model for the prior period, which would be simply time shifting the initial condition problem).The steady-state approach has the advantage of providing flows and depths that are consistent across the entire network with all the boundary conditions (inflows and channel geometry) as well as consistent with the nonlinear governing equations.These consistencies eliminate destabilizing impulses otherwise caused by mismatches between the flow/depth in a river reach and the runoff, so subsequent time marching of the unsteady solution is smooth.Furthermore, the steady-state solution is the closest available proxy to the unknown unsteady solution at t = 0, so this approach should minimize the spin-up time required to reach an unsteady time march that is independent of the initial conditions.

Overview
Herein, we present an efficient approach to establishing a set of steady-state conditions that provides a consistent and smooth starting point for time marching an unsteady Saint-Venant simulation.A full model initialization problem has two parts: (i) determining a set of flows and water surface elevations that are consistent steady solutions of the SVE for starting an unsteady solver and (ii) determining the spinup time needed to ensure errors in the initial conditions are washed out of the unsteady solution.The second problem is highly dependent on the network characteristics and the particular flow and boundary conditions during spin-up, so for brevity, this work deals quantitatively with solving the first problem and then illustrates the effects on the second problem.

Saint-Venant equations
The Saint-Venant equations for temporal (t) evolution of flow and water surface elevation along one spatial dimension (x) following a river channel are generally derived using the hydrostatic and Boussinesq approximations applied to the incompressible Navier-Stokes and continuity equations (de Saint-Venant, 1871).Cross-section averaging to obtain the 1-D equations is considered reasonable where crosssectional gradients are smaller than along-channel gradients.However, the equations are widely used even where such assumptions are questionable (e.g., near a bridge with multiple immersed piers), with the effects of significant crosssection gradients or non-hydrostatic behavior being represented as empirical energy losses.A number of conservative and non-conservative equation forms have been used, with different advantages and disadvantages (Hodges and Liu, 2014).Herein, we follow Liu and Hodges (2014) in using cross-sectional area (A) and flow rate (Q) as principle solution variables of the numerical system and the local water depth (h) and friction slope (S f ) as a secondary variables (i.e., variables that depend on A and Q through auxiliary re-C.-W.Yu et al.: Initial conditions for Saint-Venant equations lationships).The equation set can be written as where boundary conditions are the local channel bottom slope (S 0 ) and the local lateral net inflow (q l ), the latter representing both inflows from the landscape and outflows to groundwater.Auxiliary equations for h = h(A) are derived from river cross-section data.The Chezy-Manning equation can be used to provide the friction slope as where n is the standard Manning's n roughness coefficient and F is a convenient equivalent friction geometry (Liu and Hodges, 2014), which subsumes the conventional hydraulic radius (R h ) using a definition of where P = P (A) is the wetted perimeter and R h = AP −1 .Note that Eq. ( 4) fixes a typographical error in Eq. ( 10) of Liu and Hodges (2014) and Eq.(3.55) of Hodges and Liu (2014).
Required boundary conditions for the unsteady Saint-Venant solution are q l (t) for each stream segment, Q bc (t) at the furthest upstream node (headwater) in river branches with a Strahler order of 1, and h with an h(A) relationship at the downstream boundary (assumed subcritical).The timemarching unsteady solution requires initial conditions for (Q, A), which can also be given as (Q, h) with A = A(h).
Implementation details of the unsteady solver used herein can be found in Liu and Hodges (2014) and Liu (2014).

Pseudo time-marching approach
The most obvious approach for finding steady-state initial conditions is to time march an unsteady solver until a steady state is achieved.That is, we apply the unsteady solver with time-invariant boundary conditions of q l (t) = q l (0) and Q bc (t) = Q bc (0) for t 0 ≤ t < 0 where t 0 is our pseudo time start and t = 0 is the time for which we want a set of initial conditions.We call this the pseudo time-marching method (PTM).The initial condition for PTM is a set of Q(t 0 ) and A(t 0 ) for each stream segment (e.g., some cold start method as described above).At first glance, the logic here might seem circular: we are trying to solve for initial condition set {Q(0), A(0)} of the unsteady model and PTM requires specifying {Q(t 0 ), A(t 0 )}.This begs the question as to why PTM should be used rather than simply applying a cold start of the unsteady solver with The answer is that the key difference between the PTM using Q(t 0 ) and A(t 0 ) and a cold start of the unsteady solver with the same values is that the former has time-invariant boundary conditions while the latter's are time varying.Thus, an unsteady solver with time-varying boundary conditions is trying to take an inconsistent starting condition and converge it to a moving target.In contrast, the PTM takes the inconsistent starting conditions and attempts to converge them to a time-invariant target, which is more likely to be successful.However, the PTM does not a priori ensure consistency between the Q(t 0 ) and A(t 0 ) starting conditions and the t = 0 boundary conditions.It follows that PTM performance can be subject to the same type of problems as a cold start depending on the choice of Q(t 0 ) and A(t 0 ) and the skill of the modeler in their selection.In Sect.4.5, we show that PTM typically has problems for large river systems with complex geometry because the complexity of selecting a reasonable set of {Q(t 0 ), A(t 0 )} to ensure convergence.
The PTM is outlined as Algorithm 1.A user-selected parameter ( ) is used as a threshold tolerance value for declaring convergence to the steady state.A typical choice of the tolerance is the square root of the computer hardware tolerance.For example, on a 64-bit Intel architecture, the hardware tolerance for a double precision floating point floating number is 2.2204 × 10 −16 , which means a good choice of is 1.4901 × 10 −8 .As a practical matter, of 10 −6 or even 10 −4 is likely to be sufficient for initial conditions; that is, as further spin-up time is still required to dilute initial condition errors, the convergence needs only to be sufficient for consistency across the network.The method can use a timestep size that is either constant or varying, with an automatic reduction in step size when convergence is not achieved in a given time step (Liu and Hodges, 2014).To avoid infinite runtimes for non-convergent behavior (e.g., due to instabilities developed with inconsistent starting conditions), the solution is terminated (failure to converge) in Algorithm 1 after the user-selected N max iterations.The starting conditions for {Q(t 0 ), A(t 0 )} are discussed in Appendix A.

Steady-solution method
The PTM approach (above) results in a steady solution of the unsteady Saint-Venant equations that satisfies both momentum and continuity for time-invariant q l (0) and Q bc (0) boundary conditions in the unsteady solver.However, we can achieve a similar effect more directly by writing a steadystate version of the Saint-Venant equations as A key point, implied by Eq. ( 5), is that the spatial gradient of steady-state Q over a stream segment is entirely due to the lateral inflow (q l ) without any influence of A. It follows that for steady q l and Q bc boundary conditions, the flow in the ith river segment (Q i ) that has Strahler order S i must be the Algorithm 1 Pseudo time-marching method Solve SVE at time point t i using unsteady method 8: Compute error: if e < then 10: return Success 11: end if 12: t i+1 ← t i + t i 13: i ← i + 1 14: end for 15: return Failure 16: end procedure sum of all the Q j for all the j connected reaches of Strahler order S j < S i .That is, the steady flow at any point is simply the sum of all the connected upstream t = 0 boundary conditions.The corresponding A (and hence depth h) can then be computed with a numerical partial differential equation (PDE) solution of Eq. ( 6) for known Q values.Note that for large river networks, the natural downstream boundary condition is subcritical, which requires specification of h and the corresponding A as the starting point.We call this a "steadysolution method" (SSM).To look at this from another viewpoint, if A is uniform, then Eq. ( 6) devolves the fundamental equation of gradually varying flow, dE/dx = S 0 − S f , where E is the specific energy.Thus, the SSM corresponds to using the steady-state flow based on all boundary conditions and solving for surface elevations with a gradually varying flow solution for non-uniform cross sections.
To efficiently compute the conservative initial Q throughout the river network, it is useful to apply graph theory as discussed in Hodges and Liu (2014).A river network can be classified as a "direct acyclic graph" (DAG) as a river may split upstream or downstream at a junction, but the flow cannot loop back to a starting point.The connectivity of a DAG can be efficiently computed by applying existing graph methods, such as depth-first search (DFS) or breadth-first search (BFS), which provide simple and efficient approaches to computing Q(0) for each stream segment over an entire network.Note that these methods were designed and named by computer scientists, so "depth" in DFS and "breadth" in BFS do not refer to hydraulics or river geometry but instead are jargon referring to the graph network characteristics.For simplicity in the present work, we confine ourselves to the subset of DAG systems that are simply connected trees, i.e., where there is never more than a single downstream reach from any junction (as shown in Fig. 1) so that there are no uncertainties in flow directions or magnitudes.Extending the method to geometry with multiple downstream reaches (e.g., braiding, canals, deltas) requires additional rules for downstream splitting of flows that are beyond the scope of the present work.
A simple DFS traversal (Cormen et al., 2001) for Q is shown in Algorithm 2. From each headwater node (Q j ), the inflow boundary condition is propagated downstream by adding the value to the downstream node and including any lateral q l from the upstream reach (stored in Q k ).For river networks, the DFS traversal is highly efficient and requires negligible computational time for river networks of 10 5 computational nodes (e.g., Liu and Hodges, 2014).Based on our experience, the DFS computational costs should be essentially trivial for even continental-scale systems of 10 7 nodes.
Algorithm 2 DFS traversal for Q 1: procedure QTRAVERSAL 2: for all i do {initialization} 3: Q i ← 0 4: end for 5: for each headwater node j with BC Q j (t) do 6: k ← downstream node of node j 8: while k is not empty do 9: k ← downstream node of node k 11: end while 12: end for 13: return 14: end procedure After the steady Q i for each stream segment is computed, Eq. ( 6) can be solved for the corresponding A i .We discretize this equation with the Preissmann scheme, similar to the approach used for the unsteady Saint-Venant equation in Liu and Hodges (2014).The value and derivative for any term are approximated as where subscripts indicate a node in the discrete system.Using j + 1/2 to represent geometric data that are logically between nodes (i.e., roughness n and S 0 ), Eq. ( 6) becomes These nonlinear equations are similar to the unsteady discrete equations, except that Q for each computational node is known from the DFS traversal.Newton's method is used to solve this system for A without linearization, similar to the approach in Liu and Hodges (2014).The SSM requires a starting guess for A to solve the steady-state problem.Herein, we use a bisection method with the Chezy-Manning equation for normal depth conditions (discussed in Appendix A).The overall algorithm for SSM is illustrated in Algorithm 3.

Overview
The performance of PTM and SSM are examined with a series of test cases ranging from simple uniform cross sections over short river reaches to 15 000 km of a real river network.To demonstrate the robustness and performance of the SSM, we conduct tests from three perspectives: (i) effects of different cross-section geometries; (ii) scalability with an increasing number of computational nodes; and (iii) realworld river networks.Two different computers are used: the cross-section and scalability tests are run on a computer with 2.00 GHz Intel Xeon D-1540 CPUs and 64 GB of RAM, while the large network tests are run on a computer with 2.52 GHz Intel i7-870 CPUs and 8 GB of RAM.In both cases, Ubuntu Linux is the operating system and GNU C++ compiler is used.

Effects of cross-section geometry
Test cases for cross-section geometry effects were conducted for synthetic geometry of simple river reaches without tributaries.Cases included rectangular, parabolic, trapezoidal, and non-uniform cross sections, with a range of channel lengths, widths, and computational nodes, as provided in Table 1.

Scalability
To demonstrate the scalability as the number of computational nodes increases, we use the geometry and flow conditions of Case 4 in Table 1 and generate synthetic test cases with increasing numbers of nodes from a few hundred to over a million in the set: { 560, 2800, 5600, 11 200, 22 400, 44 800, 89 600, 179 200, 358 400, 716 800, 1 433 600 }.

Large river networks
To examine the robustness of PTM and SSM for more realistic conditions over both small and large scales, we use a section of Waller Creek (Texas, USA) as well as the entire watershed of the San Antonio and Guadalupe river basins (Texas, USA).The former is a small urban watershed for which dense cross-section survey data are available, whereas the latter is a large river basin that has been previously modeled with the RAPID Muskingum routing model (David et al., 2011) and the SPRNT Saint-Venant model (Liu and Hodges, 2014).
The Waller Creek study includes two stream reaches and the catchment area illustrated in Fig. 2. The total stream length is 11.6 km, which drains an area of 14.3 km 2 .The layout of Waller Creek is shown in Fig. 2a, and parts of the bathymetry surveyed data from City of Austin are shown in Fig. 2b for clarity.Two different model geometries were considered, which are designated as WCA and WCB.For WCA, the stream is discretized by 373 computational nodes based on separation of the surveyed cross sections.WCA neglects the minor tributary of Waller Creek and includes the full complexity of the surveyed cross sections shown in Fig. 2b.In contrast, WCB includes both tributaries but uses wider computational node separation with only 30 of the 373 surveyed cross sections.
To test the initial condition approach for a large river network, we use the San Antonio and Guadalupe river basins (Fig. 3), which have a combined total stream length of 12 728 km (excluding some minor first-order segments).The model herein uses 63 777 computational nodes, 59 594 segments, and 2643 junctions but is otherwise similar to the model setup with 1.3 × 10 5 nodes used in Liu and Hodges (2014)   synthetic inflow data set for the headwater inflows.The synthetic flow at each headwater stream was computed based on a downstream peak flow rate distributed uniformly across all the headwater reaches.We used the peak flow rate recorded on the main stem of Guadalupe River at Victoria (Texas) on 19 January 2010 by USGS gauge 08176500.As this gauge does not include the San Antonio River flows, we divided the peak flow rate (453 m 3 s −1 ) by the total number of headwater streams in the Guadalupe River (815) to get a single inflow value that was applied to each headwater reach (0.55 m 3 s −1 ).The same flow rate was used for the 725 headwater reaches of the San Antonio River network.This approach ensures that there is flow in every branch in the river network.
As is often the case in large river networks, comprehensive cross-section geometry data were not available for the San Antonio and Guadalupe rivers.Indeed, Hodges (2013) noted data availability and our ability to effectively use syn-thetic geometry as one of seven fundamental challenges to continental river dynamics modeling.
Because the geometry affects both PTM and SSM solutions, we tested four different estimation approaches for synthesizing geometry (Cases A, B, C, and D).Case A uses synthetic trapezoidal cross sections using the approach applied in Liu and Hodges (2014) based on Western et al. (1997).In this method, trapezoidal widths (W ) are computed from mean annual flows (Q m ) from the NHDPlus data set as W = αQ 0.5 m with α = 1.5.For the side slope of the trapezoidal cross section, an identical sidewall slope (45 • ) is used throughout the river network.Case B channels were similar to Case A but included some minor changes to Manning's n, inflow boundary conditions, and channel bottom slopes in reaches where instabilities occurred, which was necessary to provide convergence for the PTM (see Sect. 4.5).Case C channels were based on work of Santibanez (2015), who used USGS streamflow measurements in the San Anto- nio and Guadalupe river network along with the at-a-station hydraulic geometry approach (Rhodes, 1977) to find the best trapezoidal cross-section approximation for the drainage area.Using this approach, the bed width (b 0 ) is an exponential function of cumulative drainage area (A D ) as where b 0 is meters, A D is km 2 , and the coefficients are γ = 12.59 and λ = 0.382.The Santibanez (2015) approach provides reasonable values for trapezoidal channel sidewall slopes over most of the basin but fails in many of the firstorder streams with small drainage areas (< 25 km 2 ) where the computed sidewall slopes are near zero.For simplicity in the present test cases, a uniform value of 45 • is used for the sidewall slopes throughout the river network.Case D uses channel bathymetry data generated from Zheng (2016), which use a height above nearest drainage (HAND) analysis (Nobre et al., 2011) applied to the National Elevation Dataset (NED) to provide an automated approach for estimating trapezoid-based composite cross sections.

Comparison metrics
The overall algorithm efficiency is evaluated by the number of Newton iterations required for convergence to steady state.The number of Newton iterations reflects the difficulty in converging the nonlinear solution and is proportional to the simulation runtime.As this metric is independent of computer architecture, it provides a universal measure of algorithm performance.For SSM, we use the number of iterations to converge the area (A) solution of Eq. ( 6), which is the dominant computational cost (i.e., the non-iterative graphtraversal solution for Q is negligible in comparison).For PTM, we use the cumulative sum of Newton iterations for the (Q, A) solution over all pseudo time steps.Where converged solutions of PTM and SSM both exist, comparisons (not shown) indicate the resulting (Q, A) steady-state results are identical within the convergence tolerance ( = 10 −6 ).

Effects of cross-section geometry
Table 2 provides a comparison of Newton iterations for the test cases of Table 1 for single reaches with different channel cross sections.The SSM converges quickly across all cases, whereas the performance of the PTM is always substantially slower than that of the SSM.The performance of the PTM appears somewhat erratic, which is likely because the overall number of pseudo time steps depends on how far the starting guess is from the converged answer and the size of the time step used in the PTM pseudo time march.By comparing the geometric data from Table 1 with the results in Table 2, it can be seen that the largest discrepancies between PTM and SSM performance (Cases 6, 7) are with non-uniform cross sections.In both of these, the SSM performs O(10 3 ) times better, compared to O(10) to O(10 2 ) improvements for simple geometry cases.This result is consistent with the idea that the performance of PTM depends on how close the starting guess for {Q, A} is to the steadystate solution.With non-uniform cross-section geometry, the starting guess is generally quite far from the steady-state condition, as it is difficult to a priori estimate gradients of the water surface that match the nonlinear acceleration associated with cross-section variability.In contrast, the benchmark tests with simple cross-section geometry (Cases 1-5) show more modest speed-up by SSM, which is consistent with the steady-state solution for PTM with simple geometry being closer to the starting guess.For short reaches with simple geometry and only a few computational nodes (Cases 2, 3), the speed-up by SSM is essentially irrelevant.

Scalability
Computing initial conditions using models with varying numbers of computational nodes for Case 4 in Table 1 provides the speed-up results shown in Fig. 4.These tests use simple trapezoidal cross sections and, consistent with the results above, the speed-up advantage of the SSM is relatively modest with less than 10 3 nodes.However, beyond this point, the effective speed-up with SSM is quite dramatic.It appears that the SSM method becomes more effective than PTM both with increasing complexity of the cross-sectional geometry and the increasing number of computational nodes.

Waller Creek test cases
The results of initial condition convergence for two Waller Creek simulations are shown in Table 3.The SSM method dramatically reduces the total number of iterations to convergence, which is also reflected in reducing the computer runtime by 99 and 92 % for WCA and WCB, respectively.Although the absolute runtime for this small system is trivial for either PTM or SSM, the disparity provides insight into the performance that is confirmed with the more complicated river network (discussed below).

San Antonio and Guadalupe river basins
The results of the full river network computations are provided in Table 4, which show the SSM was successful and used a relatively small number of Newton iterations despite the complexity of the system.Although results in Sect.4.2 and Sect.4.4 indicate that PTM method is acceptable for small systems, results for large-scale river network simulations are less promising.In the San Antonio and Guadalupe river network configurations, three out of four PTM solutions failed; that is, the method diverged from any selected starting condition and finally caused convergence failure.Configurations A, C, and D all showed evidence of numerical instabilities leading to divergence.For example, configuration D in PTM method failed at convergence after 497 iterations; the maximum L2 convergence norm of the matrix reached up to 1.51e 235 , which is unquestionably divergent.Our inability to converge PTM with any A, C, or D cases led to development of the ad hoc B setup that provided the only PTM solution on the full river network.To obtain the B configuration, we identified reaches where instabilities developed and made minor ad hoc adjustments for Manning's n, inflow boundary conditions, and channel bottom slopes until the model showed reasonable convergence behavior (see below).Note that the modeler's time to tune the system for the PTM method to   4.
successfully converge is not included in the comparisons of Table 4.The convergence behavior of the PTM for configuration B is shown in Fig. 5.It can be seen that for several hundred time-marching steps the solution was oscillating rather dramatically but eventually settled down to a slow, smooth behavior.We believe this is evidence of the PTM trying to overcome inconsistencies between the {Q(t 0 ), A(t 0 )} starting conditions and the boundary conditions in the network.Note that PTM for B was not converged to the same = 10 −6 tolerance used for SSM.Instead, the solution was manually terminated after more than 9 h, when the convergence norm reached 1.6×10 −4 and was sufficiently smooth so that it was clear that the method would eventually converge.

Effects on spin-up
As alluded to in the introduction, obtaining an effective model initial condition is only one step in the initialization of an unsteady model.A second step is understanding at what time the model results are independent of any errors or inconsistencies in the initial conditions -i.e., the spin-up time.Some model spin-up time is generally unavoidable as we never have exactly the correct spatially distributed initial conditions that are exactly consistent with spatially distributed boundary conditions.In effect, eliminating spin-up time requires a set of initial conditions that are not only consistent with the boundary conditions at t = 0 but also consistent with the boundary conditions for t m < t < 0, where t m represents the system "memory" (or the time interval to wash out a transient impulse).
As an illustration of the scale of the spin-up problem compared to the initial condition problem, we have run the SPRNT unsteady SVE model (Liu and Hodges, 2014) for the San Antonio and Guadalupe river network using over 30 000 data points of unsteady lateral inflows for 14 days in January 2010.These boundary condition data were generated from the North American Land Data Assimilation System (NLDAS).The initial conditions were generated using SSM, as described above.The initial conditions were then perturbed by ±20 % in every first-order reach, which provides two slightly different initial condition data sets to compare to the baseline.In Fig. 6, the time-marching results for the perturbed and baseline initial condition cases reach the same state throughout the network (0.001 % threshold value) at 152 and 154 h of simulation time, respectively.Thus, approximately 160 h represents a conservative estimate of the expected time for errors in first-order streams to be diluted in the higher-order (larger) river branches.Note that it only takes 3.8 s of CPU time to compute initial conditions using SSM and an additional 5 min of CPU time to compute the time marching during the spin-up interval with the SPRNT unsteady model.This is 2 orders of magnitude faster than the 9 h or more of CPU time required just to compute initial conditions using PTM for the same system.

Model performance
In general, the PTM performed poorly except on very simple systems.As the river network complexity increases, the PTM changes from being somewhat slower than SSM to being non-convergent.Indeed, the PTM has only one advantage over the SSM in providing initial conditions to an unsteady SVE solver: specifically, no new code is needed as PTM uses the same unsteady SVE code.However, using PTM for large Hydrol.Earth Syst.Sci., 21, 4959-4972, 2017 www.hydrol-earth-syst-sci.net/21/4959/2017/ Figure 6.Spin-up for the San Antonio and Guadalupe river network with the SPRNT unsteady SVE model initialized using the SSM approach.The positive and negative 20 % perturbations are for the Q initial conditions in first-order reaches.
systems requires a frustrating trial and error approach to tuning the system to obtain convergence.In contrast, the SSM provides a rapid solution to the initial condition problem because Q is computed from simple graph traversal (once through the network), and the subsequent computation of A is "local" (in the sense there is no coupling between distant computational nodes).Note that in contrast to PTM, the SSM does not require the modeler to select a set of starting conditions; i.e., the SSM starting condition (Q) is defined solely by the t = 0 boundary conditions.Thus, different modelers will produce exactly the same Q(0) and A(0) with the same number of iterations when using the SSM on identical geometry and boundary conditions.This same cannot be said for PTM as modelers must select their {Q(t 0 ), A(t 0 )} starting condition and may need to resort to custom model tuning to obtain convergence.Herein, we only tested two methods for initial conditions, both based on finding the steady-state {Q, A} that are consistent with the boundary conditions.However, we can also argue that the cold start and synoptic start (see Sect. 1.3) would likely perform as bad or worse than PTM.That the cold start would perform poorly follows from the fact that it has the exact same problem as the PTM (converging over time from inconsistent starting data) but increases the difficulty by trying to converge to the unsteady boundary conditions.A cold start effectively turns the initial condition problem into a spin-up problem.For a cold start model performing similarly to the PTM for the San Antonio and Guadalupe river network, we can expect spin-up to require more than 10 4 time steps of the unsteady solver.
Although it is possible that a synoptic start could perform better than PTM or a cold start, it seems likely that any approach to interpolating/extrapolating sparse observational data across a larger river network will necessarily result in inconsistencies between the initial {Q, A} and the boundary conditions.If such inconsistencies result in model instabilities (a difficult thing to predict), the overall model spin-up time could be extensive.The key problem for the synoptic start is that it requires judgment as to how to best interpolate/extrapolate observational data for initial conditions, which is contrasted to the SSM approach of simply using the actual Q(0) boundary conditions and the steady solver without any further choices by the modeler.
Note that the poor performance of the PTM cannot be attributed to inefficiencies in the unsteady solver.As discussed in Sect.5.1, the unsteady solver computed 150 h of unsteady simulation from the SSM about 5 or 1800× faster than real time.This implies that the unsteady solver computed 114 min of real-world unsteady flows in the 3.8 s required for the SSM to compute the steady-state initial condition.Thus, the SPRNT unsteady solver is quite efficient -except when given an inconsistent set of initial and boundary conditions as in the PTM solution.

Limitations
For simplicity, the SSM algorithms presented herein are for common river networks with only a single reach downstream of each junction.Such tree networks are a subset of the more general DAG river networks.Extension to more complex DAG geometry requires the definition of splitting rules to uniquely define partitioning of Q where the flow splits into multiple downstream branches.Search algorithms for more complicated DAG forms are available in the graph theory lit-erature (e.g., Cormen et al., 2001) and could be readily extended to the SSM.
Fundamental to the SSM approach is the presumption that there are no places of net flow accumulation or loss throughout the network.All reservoirs and hydraulic structures are treated as pass through so that all the upstream flow is propagated through the downstream reaches.However, if there are known locations of accumulation or loss, the river network could be divided into subnetworks with separate SSM solutions.For example, upstream of a reservoir an SSM solution can be used to obtain the inflow to the reservoir, which can be used with operating information from the reservoir to provide the correct outflow Q for use in the SSM solution for downstream reaches.Hydraulic structures that affect h (or dh/dx) as a function of Q can be readily included in the SSM.To do so merely requires changing the momentum equation in a reach to model the physics of the structure (as is similarly done in unsteady SVE models).

Conclusions
It is demonstrated that inconsistencies between initial conditions and boundary conditions for a large river network solver of the Saint-Venant equations can lead to long spin-up times or solution divergence.We note that synthetic initial conditions are preferred over observed synoptic initial conditions due to the ability of the former to provide smooth and consistent spin-up.Two methods to compute synthetic initial conditions for flow (Q) and cross-sectional area (A) for an unsteady Saint-Venant river network model have been presented.Both approaches use the steady-state solution for t = 0, which provides initial conditions that are smooth and globally consistent with the boundary conditions for the model start time.The PTM is likely similar to undocumented approaches previously used; i.e., application of an unsteady model with constant boundary conditions to achieve a steady solution consistent with initial land surface inflows.For large river networks, the PTM method is slow and inconsistent, arguably depending on the quality of the first guesses for Q and A and the size of the time step required for a stable pseudo time march.A new SSM is developed to address these issues.The SSM computes the initial condition Q in each reach from the inflow boundary conditions of the entire network at t = 0 by applying a mass-conservative graph traversal technique.The initial condition A in each reach is found from the solution of the steady-state 1-D momentum equation with known Q.The first-guess problem for A is solved as a normal-flow problem with the Chezy-Manning equation and the Q from graph traversal.Although SSM requires writing an additional numerical solver rather than relying on an existing unsteady solver (as used in PTM) our numerical experiments show that SSM is more robust and consistently faster than PTM.The code for both initial condition solvers is publicly available at GitHub (Liu, 2014).
Data availability.Cross-section geometry test cases can be found in the open-source repository mentioned in Liu (2014).Test cases and files for scalability, Waller Creek, San Antonio and Guadalupe River network are uploaded to a public repository under Texas ScholarWorks (http://hdl.handle.net/2152/61757).

Figure 1 .
Figure 1.Propagation of flow rate Q at a junction.

Figure 2
Figure 2. (a) Waller Creek and catchment in Austin (Texas, USA).(b) Surveyed cross sections of main channel for Waller Creek (Texas).Only 149 of 327 cross sections are shown for clarity.Elevations are relative to mean sea level (data courtesy of the City of Austin).

Figure 3 .
Figure 3. San Antonio and Guadalupe river network from an NHDPlus V2 flowline.

Figure 4 .
Figure 4. Speed-up multiplier of SSM compared to PTM for Case 4 as a function of the number of computational nodes.

Figure 5 .
Figure 5. Convergence of the L2 norm between consecutive pseudo time-marching solutions for the PTM with configuration B of the San Antonio and Guadalupe river network.Note the above figure uses the number of time-marching steps as compared to the larger number of Newton iterations provided in Table4.
. Although the unsteady SPRNT model is typically run by coupling with a land surface model for headwater and lateral inflows, for the present steady-state tests we used a . Earth Syst.Sci., 21, 4959-4972, 2017 www.hydrol-earth-syst-sci.net/21/4959/2017/

Table 1 .
Cross -section geometry test cases.W B and S sw represent bottom width and sidewall slope, respectively; f represents the focal length of parabolic shape.

Table 2 .
Newton iterations required to achieve convergence for benchmark geometry test cases.The converged results are identical for both methods.

Table 3 .
Total Newton's iterations required to achieve convergence of the Waller Creek test case.