Inverse streamflow routing

The process whereby the spatially distributed runoff (generated through saturation/infiltration excesses, subsurface flow, etc.) travels over the hillslope and river net- work and becomes streamflow is generally referred to as "routing". In short, routing is a runoff-to-streamflow process, and the streamflow in rivers is the response to runoff inte- grated in both time and space. Here we develop a method- ology to invert the routing process, i.e., to derive the spa- tially distributed runoff from streamflow (e.g., measured at gauge stations) by inverting an arbitrary linear routing model using fixed interval smoothing. We refer to this streamflow- to-runoff process as "inverse routing". Inversion experiments are performed using both synthetically generated and real streamflow measurements over the Ohio River basin. Re- sults show that inverse routing can effectively reproduce the spatial field of runoff and its temporal dynamics from suf- ficiently dense gauge measurements, and the inversion per- formance can also be strongly affected by low gauge density and poor data quality. The runoff field is the only component in the terrestrial water budget that cannot be directly measured, and all pre- vious studies used streamflow measurements in its place. Consequently, such studies are limited to scales where the spatial and temporal difference between the two can be ig- nored. Inverse routing provides a more sophisticated tool than traditional methods to bridge this gap and infer fine- scale (in both time and space) details of runoff from aggre- gated measurements. Improved handling of this final gap in terrestrial water budget analysis may potentially help us to use space-borne altimetry-based surface water measurements for cross-validating, cross-correcting, and assimilation with other space-borne water cycle observations.


Introduction
Runoff is a very important component in the terrestrial water budget (precipitation, evapotranspiration, runoff, and soil/snow water storage) in terms of both its magnitude and temporal variability (Hagemann and Dumenil, 1998;Pan et al., 2012).And runoff is also the only component in the terrestrial water budget that cannot be measured directly at the time and location it occurs.When precipitation is measured by rain gauges, radars, or satellite sensors, the measured value is validated at the same time and location it rains or snows, and so is evapotranspiration by towers/satellites and soil moisture by probes/microwave sensors.But so far there seems to be no way of measuring the spatial field of runoff as it occurs.Therefore, previous studies used the streamflow measurements in place of runoff (Sahoo et al., 2011;Sheffield et al., 2009).However, as streamflow can be measured at river gauges and will be measured at large scales by space-borne altimetry sensors through Manning's equation (Alsdorf and Lettenmaier, 2003), for example, the planned Surface Water and Ocean Topography (SWOT) mission (Alsdorf et al., 2011), it is not fully equivalent to runoff.So all such studies are limited to the situations/scales where the temporal and spatial difference between the two can be either ignored or somehow accounted for.For example, the river travel time may be ignored at long-term timescales.Hydrologically, streamflow differs from runoff by one process called "routing".
The process by which the spatially distributed runoff generated through various mechanisms (e.g., saturation excess, infiltration excess, and subsurface flow) travels over the hillslope and river network and becomes streamflow is referred to as routing.During the routing process, the streamflow at a particular location in the river channel is a collective result of runoff from different locations and times.In other words, Published by Copernicus Publications on behalf of the European Geosciences Union.
the streamflow is the response to the runoff field integrated in both time and space.Routing essentially provides a runoffto-streamflow conversion, and it is a quite well studied process in hydrology.Routing models have been developed to parameterize this runoff-to-streamflow process and predict the streamflow at desired gauging locations given rainfall or runoff inputs (Brutsaert, 1994;Lohmann et al., 1998).
Now the question for our study is how to bridge the gap between streamflow and runoff in both time and space, i.e. to derive the spatial fields of runoff from streamflow measurements at gauging points, such that our water budget analysis or any related studies are no longer limited by the gap between the two.Obviously, this requires us to invert the routing process and invent a way to realize the streamflow-torunoff conversion.We refer to such a streamflow-to-runoff process as "inverse routing".Note that solving inverse problems is nothing new in hydrology.For example, inverse problems are frequently studied in groundwater hydrology for parameter estimation purposes (McLaughlin and Townley, 1996).While the general methods for solving inverse problems are no different from any other optimal estimation problems like data assimilation (McLaughlin, 2002;Reichle, 2008), different problems may require very different methodological considerations.For example, parameterestimation-related inverse problems usually solve for static (time-invariant) unknowns.Thus complicated and computationally intensive methods may be used to invert subtly behaved nonlinear models with non-Gaussian errors.The inverse routing problem needs to solve for dynamic fields of runoff repeatedly in time, and thus it requires a higher computational efficiency.Also, the streamflow values are always correlated in time as a result of the time integration nature of the routing process, and that implies the unknown runoff fields across multiple time steps need to be solved together, which dramatically increases the size of the estimation problem (number of simultaneous unknowns).
For these above reasons, we look for a linear routing model to invert such that the most efficient methods for linear systems like Kalman filters/smoothers (Anderson and Moore, 1979;Kalman, 1960) can be applied.Also, inverse routing is not the only way to estimate runoff -hydrological models infer it from rainfall and other inputs, and the basin-specific discharge has always been calculated from gauge measurements and water balance relationship between gauge basins and their contributing sub-basins (though the fine-scale variability below sub-basins is ignored).So the potential contribution from inverse routing should also be evaluated in terms of its "added value" to the traditional methods.
In Sect.2, we will first introduce the routing model to use and show how to invert it using a special type of data assimilation technique called fixed interval smoothing.Then in Sect. 3 an inverse routing experiment will be performed using synthetically generated runoff/streamflow data where the inversion errors and related performance issues can be investigated against the synthetic truth.Later in Sect. 3 an inverse routing experiment will be performed using real river gauge measurements from the United States Geological Survey (USGS) to evaluate the inversion performance in real world applications.And finally the various limitations of inverse routing will be studied in Sect. 4.

A linear routing model
Here we choose the University of Washington (UW) routing model (Lohmann et al., 1996(Lohmann et al., , 1998)), which is a relatively simple linear routing model developed for coupling with land surface models (LSMs), and it has been calibrated, implemented and validated in many large-scale streamflow studies (Mitchell et al., 2004;Nijssen et al., 1997).The inputs to the UW routing model are runoff fields defined on a rectangular computing grid -the format used by most LSMs.The UW model routes the runoff water in two stages: first the runoff water drains from within the grid pixel (over the hillslope) to a conceptual "outlet" of the pixel following a known unit hydrograph function (UHF) u(t), and the pixel outflow o(t) is the convolution between the UHF u(t) and pixel runoff r(t): (1) Then the water travels in channels between pixels following the 1-D diffusive wave equation: Here q = q(x, t) is the streamflow generated by the pixel outflow o(t) at distance x downstream from the pixel.C and D are referred to as channel wave velocity and diffusivity parameters following our upstream literature (Lohmann et al., 1996(Lohmann et al., , 1998)).Note that a more precise name for C is in fact the wave celerity (Beven, 1979), which better distinguishes it from the fluid velocity measured in river channels.The 1-D diffusive wave equation is a standard advection-diffusion equation for transport, and it is always linear as long as neither C nor D is a function of q.In other words, the wave velocity C and diffusivity D can change in both time (t) (e.g., from summer to winter) and space (x) (e.g., from flat areas to mountains) as long as the values can be prescribed and are independent of q.Retention effects (e.g., lakes and reservoirs) are not considered.The equation is solved analytically using the convolution between the impulse response function (IRF) and pixel outflow: where i(x, t) is called the IRF.Note that mathematically UHF is identical to IRF in their functional roles, and the two convolutions can be combined because the convolution operations here are associative.Define a combined IRF h(x, t) as the convolution between u(t) and i(x, t): The combined IRF h(x, t) is the "overall" hydrograph function in response to a unit runoff input from one pixel.And the two-stage routing is solved at once using h(x, t): At a given gauge location g, we calculate the streamflow value Q(g, t) by integrating (summing up) the contributions from all upstream pixels (i.e., the entire sub-basin that drains to g, noted as basin (g)): Equation ( 7) fully defines the integration of runoff field in space and time into the streamflow at a gauge location.As the routing model always runs in discretized time steps, the integration Eq. ( 7) is implemented as summations: The UW routing model is fully contained in Eq. ( 8), and, in short, the streamflow at a gauge point is nothing but the sum of runoff from all contributing pixels in all possible lag times weighted by the overall IRF.This routing model is linear and simple, though the number of runoff inputs for streamflow calculation is very large, making the inverse problem challenging.

Inversion through fixed interval smoothing
In dynamic system analysis and related estimation theories, the prediction model is mostly written in a "state space" form: the observations are written as a function of input states in a vector/matrix form with some model error term ε like y = Hx + ε.Now we rewrite Eq. ( 8) in this form.Say we have m gauge locations and n runoff computing pixels in the study area, and we define the streamflow observation vector y t and runoff state vector x t as the collection of all m gauge measurements Q 1 , Q 2 , . . ., Q m and runoff states in all n pixels r 1 , r 2 , . . ., r n at time t: As a result of the time integration in Eq. ( 8), the calculation of y t requires not only x t but also runoff fields at previous time steps x t−1 , x t−2 , . . ., x t−k .Physically, all the lag times within the longest travel time to the gauges should be included (i.e., until the last bit of runoff from the farthest pixel in the basin passes the most downstream gauge).Say the maximum travel time is k+1 time steps, and the observation equation is H 0 , H 1 , . .., H k are the measurement operator matrices for different lag times and each has the size m×n.The entries in the operator reflect how much of the runoff from one specific pixel contributes to one specific gauge at a specific lag time.
All of them are calculated from the combined IRF h(x, t) according to the downstream travel distance and lag time.Again, because of the time integration, direct inversion of Eq. ( 10) is impossible and incomplete because streamflow at future time steps also contains information about the runoff at current time step, and the time series of streamflow are highly correlated.This means the inverse estimation must be done for multiple time steps at once, and observation and state vectors need to be "augmented" to include multiple time steps.Writing Eq. (10) for all s + 1 time steps in the time interval [t − s, t], we have And define the time-augmented streamflow/runoff/error vectors as , and where y t has the size m(s + 1) and x t has the size n(s + 1).
Then we can write the augmented observation equation as In the above, the augmented observation operator matrix H is H is mostly empty except the upper diagonal belt is filled with H 0 , H 1 , . . ., H k , and it has the size m(s + 1) × n(s + 1).The augmented observation equation has one extra term compared to the un-augmented one: L has the same size as H : m(s +1)×n(s +1).It is empty for the first m(s + 1 − k) rows and n(s + 1 − k) columns, and the remaining mk × nk block at the lower right corner is filled with H 1 , H 2 , . . ., H k in the lower diagonal.The extra term L x t−k provides the routing model the initial conditions -the water stored in river channels contributed by runoff in the k time steps prior to time t −s (i.e., interval The time interval of augmentation needs to be larger than the maximum travel time (i.e., s +1 > k +1).This is because the streamflow measured at the outlet of the basin k + 1 time steps later still contains information about the runoff generated at the most upstream pixel of the basin at the current time step.To enable the inversion to use all possible streamflow information from all gauges to update the most upstream pixel, there must be s + 1 > k + 1.Now we can apply the Kalman-type methods for the inversion.Since we always have a fewer number of gauges than runoff pixels (i.e., m < n), the inverse problem is underconstrained.Therefore, an initial guess of runoff fields, noted as x t , is usually needed to represent the prior information we have on this variable.This initial guess could just be a uniform field of long-term mean runoff (or simply zeros), which is referred to as the "null" initial guess, or LSM estimates forced with some baseline rainfall (better choice) if possible.Note that the "null" initial guess is equivalent to having no initial guess at all.Given the initial guess x t and streamflow measurements y t , the Kalman filter equation gives an updated estimate x t as And the Kalman gain K t is calculated as K t has the size n(s + 1) × m(s + 1), and P t is the error covariance matrix of the initial guess of runoff and has the size n(s + 1) × n(s + 1).P t can simply be a diagonal matrix of the long-term mean runoff error variance, or with diagonal entries proportional to the squared runoff if the initial guess is not uniform.R t is the error covariance matrix of the gauge measurements and has the size m(s + 1) × m(s + 1).R t can be looked up from river gauge documentations or empirically estimated from instrumentation type, flow rate, channel morphology, etc.However, here we want to force the updated runoff field to reproduce the streamflow observed at all the gauges exactly, i.e., to make the updated runoff field x t to satisfy Eq. ( 13) with no errors: This exact match is achieved by treating the gauge measurements as error-free (i.e., R t ≡ 0), and Eq. ( 16) will push all the streamflow errors in the initial guess y t − H x t − L x t−k back to the runoff guess and effectively force Eq. ( 18) to be exactly satisfied.Caution is needed here that treating the gauge measurements as perfect does not mean we believe there are no errors in them.There are definitely errors, and such treatment is merely a way to force the exact reproduction of streamflow, and all of the gauge errors will be carried into the inverted runoff.Also, such a setting maximizes the correction Eq. ( 16) can impose onto the initial runoff guess (Pan and Wood, 2010).This is a same measure as taken by the constrained data assimilation procedures (Pan and Wood, 2006;Pan et al., 2012).Now the Kalman gain becomes Note that the above update procedures are no different than a regular data assimilation when R t is not particularly chosen.And in fact, many studies have been devoted to the assimilation of streamflow or water altimetry measurements (Andreadis et al., 2007;Biancamaria et al., 2011;Durand et al., 2008).We would like to call the runoff estimation with the particular setting of R t ≡ 0 as "inverse routing", in order to differentiate it from the general practice of streamflow assimilation.Also, since the inversion involves multiple time steps, the procedure is no longer a filtering operation but a smoothing operation or, more precisely, a s + 1 step fixed interval smoothing.
During the Kalman gain calculation in Eq. ( 19), the matrix H P t H T to invert has the size m(s + 1) × m(s + 1).As matrix inversion has its computational complexity grow cubically against matrix size, the interval size s+1 and number of gauges m to use cannot be too large.Per discussion on the time augmentation, s + 1 has to be larger than k + 1.Note that even with s + 1 > k + 1, updating the runoff in the last k + 1 time steps in the [t − s, t] smoothing interval (i.e., not receive all possible streamflow information.To eliminate such an "edge effect" of fixed interval smoothing, the smoothing will be done interval by interval sequentially, and consecutive intervals will overlap by k + 1 steps.In other words, only the first s − k steps in the s + 1 smoothing interval are usable, and the last k + 1 steps will be re-updated in the next smoothing interval.

Study area and general setups
We choose the Ohio River basin in United States for the inversion experiments.Figure 1 shows the definition of the Ohio River basin, which also includes the Tennessee River in the south and covers an area of 490 600 km 2 .Rivers in the area are very well monitored by the United States Geological Survey (USGS), and we select 75 USGS river gauge stations out of all available ones (i.e., m = 75).Gauge stations that are too close to each other were selectively removed to reduce redundancy.The computing grid for the UW routing model is set up at 0.125 • and the flow network on this grid is derived from 30 arcsec digital elevation model (DEM) data, as shown in Fig. 2. The computing grid consists of 3681 0.125 • pixels (i.e., n = 3681).All routing model parameters are identical to those used in the National Land Data Assimilation (NLDAS) project (Lohmann et al., 2004, Mitchell et al., 2004) over the same area, where the same routing model has been calibrated and validated against USGS-observed streamflow.The channel wave velocity C = 1.4 m s −1 and diffusivity D = 0 m 2 s −1 are constant all over the basin.The time step of the routing model is 1 day, and because the 0.125 • pixel is small enough for any runoff to flow out of the pixel within 1 day, the outflow UHF u(t) = 1 when t = 1 and u(t) = 0 when t > 1.The resulted runoff water travel time to the basin outlet can be found in Fig. 2, and the maximum travel time is 16 days (i.e., k + 1 = 16).The study period is the entire year (365 days) of 2009.Two types of inversion experiments will be performed where the streamflow (y t in Eq. 16) to be inverted is generated differently: 1. Inversion experiment with synthetically generated streamflow, or in short synthetic experiment.
In this experiment, we assume the "true" values of runoff and streamflow are known, and the synthetically "true" streamflow values will be inverted (i.e., be assimilated into the initial guess of runoff) ( x t in Eq. 16).Then errors in the inverted runoff ( x t in Eq. 16), and initial guess as well, will be calculated against the synthetically "true" runoff.To do this, the synthetically "true" runoff fields will first be created by running the variable infiltration capacity (VIC) LSM (Liang et al., 1994(Liang et al., , 1996) ) forced with the 0.125 • NLDAS meteorological data set (Cosgrove et al., 2003).NLDAS rainfall combines hourly WSR-88D radar analyses and daily gauge reports (∼13 000 per day) and is considered the best available surface forcing over United States.Given the NLDASderived "true" runoff, the synthetically "true" streamflow is created using the UW routing model.Then with an initial guess of runoff, the inversion is performed, and the errors are calculated against the NLDAS-derived synthetic truth.The benefit of a synthetic experiment is that the performance of the inversion method can be well evaluated using the synthetic truth.Also, as the model-derived streamflow is assimilated, the complications caused by errors/biases in the routing model are avoided.
2. Inversion experiment with real streamflow measurements, or in short real experiment.
This experiment differs from the synthetic one only in that the real USGS streamflow measurements will be inverted.All routing models have errors, and errors arise from simplifying assumptions of the model, imperfect model parameters, model inputs, and so on.Figure 3 shows the comparison between the NLDAS-derived synthetic "true" streamflow and actual measurements from USGS at four gauge stations with drainage area ranging from 160 511 km 2 (Fig. 3a) to 2928 km 2 (Fig. 3d).As the routing model reasonably reproduces the gauge measurements, both timing and magnitude errors can be seen, especially for larger sub-basins.Such errors will corrupt the performance of inverse routing, and the performance assessment from the synthetic experiment will be discounted when real gauge measurements are used.The purpose of the real experiment is to find out how much the performance degradation is and how that affects the usefulness of inverse routing.Besides routing model errors, LSM and its parameters/input data also have errors, and that makes the NLDAS-derived runoff not really a "truth".In the real experiment, errors will still be calculated against NLDASderived runoff, even though the latter is no longer a "truth".
In both the synthetic and real experiments, the smoothing interval for the inversion is 70 days (i.e., s + 1 = 70).The routing model spins up during the first 16 days of 2009 to provide initial conditions for later routing, and the inversion starts on day 17 of 2009.

Experiment with synthetically generated streamflow
Here we further consider two different ways of creating the initial guess of runoff ( x t in Eq. 16).The first is the "null" initial guess mentioned in the methodology section.This is to assume we have no alternative source of runoff information at all, and the initial guess is simply a uniform field of constant value for all the time steps.Here we set it to a longterm mean value of 1.474 mm day −1 .Figure 4 presents the inversion results for day 75 of 2009 using the null initial guess of runoff.The inverted runoff field shown in Fig. 4c shows a very similar spatial pattern to the NLDAS-derived synthetic truth in Fig. 4b.The two rainfall/runoff belts in the synthetic truth, one of which goes through the northeastern and northwestern tips of the basin and the other in the southern half (mostly over Tennessee River), have been well recovered in the inverted runoff.Note that the initial guess has no spatial variability at all (Fig. 4a).The empty area between the two belts is also fairly well cleared in the inverted runoff.The inversion increment (Fig. 4d), i.e., the difference between the inverted and the initial guess or as in Eq. ( 16), shows where the runoff water has been added to or removed from the initial guess.The spatial pattern in the inverted runoff is mildly patchy, and the shape of patches follows the boundaries of sub-basins that drain to the input gauge locations.Besides the sub-basin shaped patchiness, banding along the equal travel time contours (Fig. 2) also exists within sub-basins because it is impossible for the inversion to distinguish the water arriving at the same time but from different points on the band of equal travel time.The close similarity between the inverted and synthetic truth indicates that the inversion procedure developed here is potentially able to recover the runoff patterns without any prior information, provided that sufficient streamflow data are available within the area.Inversion results for 4 more randomly selected days (day 37, 107, 177, and 317 of 2009) are shown in Fig. 5.In this figure, the inversion recovers a very reasonable spatial pattern for all days, with some days (e.g., day 317) performing better than others (e.g., day 177).In other words, the inversion also recovers a reasonable timing of the runoff.Figure 5 (bottom row) also shows the 7-day mean antecedent precipitation for those 4 days, and most of the spatial patterns in the antecedent precipitation are reasonably reflected in the inverted runoff fields.The inversion problem from streamflow to runoff is extremely under-constrained here (m = 75 versus n = 3681), and such results suggest the inversion method has a very strong capability.The ability to work under null initial guess (i.e., no initial guess at all) has a critical meaning for our study.This is because in this case the streamflow measurements are the only input to the estimation problem and that means the runoff fields derived from the observed streamflow are also a purely "observationally based" quantity with no influence from an LSM.Another way to create the initial guess is to use the LSM to calculate runoff from some baseline rainfall inputs that are considered always available for all locations and all times.This is supposedly a better initial guess than the null guess.It should provide a better assessment of the inverse routing than the null initial guess experiment since it is not realistic that we know nothing about the rainfall all the time.Here we force VIC LSM with the satellite rainfall product TRMM Multi-satellite Precipitation Analysis (TMPA) version 3B42RT (Huffman et al., 2007) to obtain an initial guess of runoff.TMPA is available globally between 60 • S and 60 • N every 3 h at 0.25 • resolution.Though much less accurate than the ground-based NLDAS, it relies only on satellites and thus is available almost everywhere.The 0.25 • data are interpolated to 0.125 • in order to force VIC simulations at 0.125 • (Pan et al., 2010).Figure 6 shows the inversion results using the TMPA-derived initial guess of runoff for the same day 75 of 2009 as in Fig. 4. The inverted runoff here is similar to the null guess case in Fig. 4 but with a slightly better definition of the rainfall/runoff plumes.The transition from wet to dry areas is also smoother (i.e., less gauge basinshaped patchiness).This suggests a good initial guess that can reasonably represent the spatial and temporal dynamics of rainfall will help improve the quality of the inverted runoff.
Figure 7 shows the time series of basin mean bias and root mean squared errors from the synthetic experiment with TMPA-derived runoff as the initial guess.In Fig. 7a, the inverted runoff (red line) shows a consistently and significantly lower basin mean bias than the TMPA-derived initial guess (blue line).The time average of absolute bias is 0.145 mm day −1 for the inverted runoff and 0.553 mm day −1 for the initial guess, and the relative bias reduction is about 74 %.However, smaller bias in the basin mean does not necessarily imply small errors in the pixel-to-pixel comparisons since positive and negative errors on the same map can average out.Figure 7b shows the time series of basin mean root mean squared errors (RMSEs), and the inverted runoff still has a consistently lower RMSE than the initial guess but to a lesser degree.The time average of RMSE is 1.283 mm day −1 for the inverted runoff and 1.962 mm day −1 for the initial guess (about 35 % RMSE reduction).
Figure 8 shows the time series of streamflow calculated from the synthetic truth runoff (NLDAS-derived), initial guess runoff (TMPA-derived), and inverted runoff for the same four USGS gauge stations as in Fig. 3.The difference between the synthetic truth (thick green line) and initial guess (blue line) of streamflow (i.e., y t − H x t − L x t−k in Eq. 16) is referred to as "innovation" in data assimilation literature, and it basically drives the update of the initial guess.For all the stations shown here, the innovation (difference between thick green and blue lines) is considerably large compared to the magnitude of streamflow itself, suggesting that the inversion delivers a significant amount of information to the inverted runoff.Note that the streamflow time series reconstructed from the inverted runoff (red line) lies exactly on top of the synthetic truth (thick green line) nearly all the time.This verifies the fact that our setting of R t ≡ 0 will force the inverted runoff to reproduce the input streamflow exactly (Eq.18).Perfect reproduction of input streamflow is an important purpose of our inversion design, and it distinguishes the "inversion" from the streamflow data assimilation in a general sense.There are a few occasions, for example, during the few days before day 350 of 2009 at USGS gauge station 03216600 (Fig. 8a), where the inverted runoff leads to a slightly higher streamflow than the synthetic truth.This is because the inversion update may, under some rare occasions, produce negative runoff values, and such negative values are reset to zero.

Experiment with real streamflow measurements
The real experiment using USGS gauge measurements is carried out only with the TMPA-derived runoff as the initial guess.Some of 75 USGS gauge stations used here have missing data, but a vast majority (67 stations) of them are fairly complete for the year 2009 (available more than 95 % of the time).Figure 9 shows the inverted runoff fields versus the NLDAS-derived synthetic truth (no longer treated as a "truth" here though) for the same 4 days as in Fig. 5.Here the inverted runoff fields have larger discrepancies against the NLDAS-derived values than those in the synthetic experiment (Fig. 5).The inverted runoff can capture some large spatial features in the NLDAS-derived map, but a lot of detailed spatial features differ from the NLDAS-driven synthetic truth.Figure 10 plots the same time series of basin mean bias and RMSE as Fig. 7.In Fig. 10a, the inverted runoff (red line) still shows a consistently lower basin mean bias than the initial guess (blue line), though not as significant as in Fig. 7.The time average of the absolute bias is reduced from the 0.553 mm day −1 for the initial guess to 0.392 mm day −1 for the inverted runoff.The relative bias reduction is 29 % (compared to 74 % in the synthetic experiment in Fig. 7).The basin mean RMSE (Fig. 10b) for the inverted runoff (red line) has lower peaks than the initial guess (red line) but is often higher than the initial guess elsewhere.The time average of RMSE is 1.970 mm day −1 , and this is even slightly higher than the TMPA-derived initial guess (1.962 mm day −1 ).This suggests it is more difficult to make significant improvement to the initial guess using real gauge measurements, especially when the initial guess is already very reasonable.Large biases can be easily corrected, but small spatial details are much more difficult to recover.Many factors contribute to this degraded inversion performance in the real experiment.Generally speaking, it is because of the routing model errors (Fig. 3).For example, model assumptions like runoff water flows in 0.125 • digitized stream channels can be a reason, and model parameters (constant wave velocity/diffusivity everywhere) are far from perfect as well.Another very important reason is that many water regulation structures (dams, reservoirs, etc.) operate in this area, and the USGS-measured streamflow is not the natural flow.Moreover, while the inversion treats the streamflow inputs as perfect observations, the USGS measurements are never perfect, and errors exist in both the instruments and the stage-flow regressions.Such errors are entirely carried into the inverted runoff because of this treatment.Figure 3 also shows that the NLDAS-derived streamflow compares better to the USGS measurements at gauges of smaller drainage basins (Fig. 3c, d) than large ones (Fig. 3a, b).A possible reason is that smaller basins are less affected by flow regulations.Unfortunately, dam/reservoir operations are mostly nonlinear (i.e., to cut flood peak or retain water for dry season release), and the best way to reduce their impact is to perform streamflow naturalization (Wurbs, 2006) separately before the inversion or to avoid using gauges of heavily regulated large basins.Finally, NLDAS and VIC LSM have errors too, and the synthetic truth being compared is not an exact truth.

Limitations of inverse streamflow routing
Some limitations of the inverse routing method have been discussed.For example, a well-performing forward routing model is an absolute prerequisite for the inverse routing to work.Also, no inversion or data assimilation techniques can create new information from nothing.That means the inverted runoff fields will not contain any more information (derived or direct) than what is already in the inputs, i.e., the streamflow observations (all those within the same contributing basin and maximum travel time) and the physical knowledge (routing equations/parameters) we have on how the water may flow.The inversion serves only to "infer" fine-scale details from aggregated measurements and not to "create" fine-scale information.That said, the inversion will degrade and fail as all these inputs degrade, for example, gauges too sparse, observations too sporadic, large errors in gauge measurements, errors in routing parameters, etc.The inversion experiments in the previous section are conducted over a heavily gauged basin with good-quality data, and those findings may not apply to many other places.Now we study how the performance of inverse routing will degrade by repeating the synthetic experiment using TMPAderived runoff as initial guess (Figs. 6 and 7) and adding errors and missing data to various inputs.The same two performance metrics, the time-averaged absolute bias and RMSE in the inverted runoff, are evaluated, and the degradation is tested with respect to four factors: (1) number of gauges available (from full 75 down to 4); (2) amount of missing streamflow observations, as measured by the missing percentage (0-100 %); (3) gauge observation errors, as mea-sured by the log standard deviation of lognormal errors (0-30 %); and (4) wave velocity errors in the routing model, as measured by the standard deviation of normal errors (0-2 m s −1 ).All the errors added are uncorrelated in time or space.Figure 11 shows how the bias and RMSE grow as more errors and missing data are added.In the case of fewer gauges and missing observations (Fig. 11a, b, e, f), both the bias and RMSE grow slowly at first, then faster, and finally approach the level of initial guess (dashed lines in Fig. 11) as the availability is reduced to zero.Note that performance worse than the initial guess means that the inversion is completely useless.Errors (lognormal) in gauge observations have a large impact on the runoff RMSE (Fig. 11g) -20 % lognormal errors can consume all the potential RMSE improvement while the bias is less affected (Fig. 11c).The wave velocity errors in the routing model also matter but to a lesser extent (Fig. 11d, h), and this parameter can usually be well estimated through careful calibrations.No joint sensitivity test is carried out although multiple factors will exacerbate the performance degradation.The above findings suggest that some reasonable gauge density (e.g., more than 10 gauges per 10 6 km 2 ) and data quality (e.g., lower than 25 % errors) are required for the inverse routing to provide improvements to TMPA-derived runoff fields.Such requirements may be lowered at places where TMPA performs less well.For example, TMPA is much less well calibrated in Africa than in United States because the rain gauge data there are very sparse.

Conclusions
We propose the concept of inverse routing as the process to estimate the spatial fields of runoff from point measurements of streamflow and develop the methodology to achieve it by inverting a linear routing model using fixed interval smoothing.In theory, the inversion method introduced here applies to any linear routing models.
The synthetic experiment shows that the inversion method is able to reproduce reasonably the spatial and temporal dynamics of the synthetically true runoff fields from point measurements of streamflow without any meaningful initial guess provided that sufficient gauge data exist.Besides the routing model and its parameters, the only input required by the inversion is streamflow.So inverse routing is always possible as long as sufficiently dense streamflow data are available.If a reasonable initial guess of runoff exists (e.g., from LSM), such an initial guess can help improve the quality of the inverted runoff.
The real experiment illustrates how the inversion performance will degrade when real river gauge measurements are used.The difference between the real and synthetic streamflow data is basically caused by the routing model errors.Such errors could be due to imperfect model design or parameters, but a large part is due to the human regulation of flow -an effect unaccounted for in the routing model.In short, the inverse routing can work well only if the (forward) routing model works well, and that requires efforts in routing model calibration, streamflow naturalization, etc.
The limitations of inverse routing arise primarily from two aspects -(1) the performance of the forward routing model and its failure to account for nonlinear flow behavior like retention, and (2) gauge data availability and quality.Low gauge density and poor data quality can strongly affect the usefulness of inverse routing.
Inverse routing may be performed to infer runoff at any fine scales where the routing model is parameterized, though the quality of fine-scale information so derived is determined by the density and quality of input gauge data.Historically, runoff has not been an observationally based variable, and streamflow measurements are used in its place and such studies are limited to the occasions where the mismatch between the two in time and space can be ignored.Now inverse routing provides a more sophisticated tool than traditional methods to resolve the scale mismatch and infer fine-scale details (in both time and space) of runoff from aggregated measurements.Improved handling of this final gap in terrestrial water budget analysis may potentially help us to use spaceborne altimetry-based surface water measurements for crossvalidating and cross-correcting other space-borne water cycle observations.For example, given the strong tie between the rainfall and runoff, as shown in Fig. 5, runoff fields inverted from the future SWOT mission can be used to identify and correct missing or overestimated precipitation estimates from the Global Precipitation Measurement (GPM) mission (Tapiador et al., 2012) through the water balance relationship and with the help of a hydrological model to resolve the rainfall-to-runoff process.It also makes it more convenient to assimilate the surface water measurements into other water cycle observations without worries about scale mismatch.The inverse routing is also a good tool to disaggregate streamflow information in time and space and provide more continuous and better river information for water resources management.Without satellite altimetry measurements, runoff fields derived from ground river gauges can help us study the long-term terrestrial water budget at a much higher spatial and temporal resolution.

Fig. 1 .
Fig. 1.The Ohio River basin (shaded area) and 75 USGS river gauges in use.

Fig. 2 .
Fig. 2. Flow paths over the 0.125 • grid and runoff travel time to the basin outlet.

Fig. 3 .
Fig. 3.Streamflow predictions from the routing model using the synthetically "true" NLDAS rainfall (green line) versus USGS measurements (black dots) over four gauge stations.The USGS gauge station number and drainage area are noted in the panel title.

Fig. 4 .
Fig. 4. Runoff estimates for day 75 of 2009 from the synthetic experiment using the null initial guess of runoff (constant field of 1.474 mm day −1 ): (a) initial guess, (b) synthetic truth, (c) inverted fields, and (d) the difference between the inverted and initial guess (i.e., inversion increment during the smoothing update).

Fig. 7 .Fig. 8 .
Fig. 7. Time series of basin mean bias (a) and root mean squared errors (b) from the synthetic experiment using TMPA-derived runoff as initial guess.Blue lines are for the initial guess of runoff (TMPA-derived) and red lines for the inverted runoff.All error measures are calculated against NLDAS-derived synthetic truth runoff.

Fig. 10 .Fig. 11 .
Fig. 10.Time series of basin mean bias (a) and root mean squared errors (b) from the real experiment using TMPA-derived runoff field as initial guess and real USGS gauge measurements.All lines are plotted the same way as in Fig. 7.