Coding, Web, Hydrology and more.

System Response in Hydrological Model


0X00 What is System Response?

From my own perspective, System Response is system’s inner structures and components’ feedback to a change in single input or combination of multiple inputs. If you are familiar with model sensitivity analysis, it is well-known that this analysis is based on analyzing model outputs’ differences on given inputs’ variation, and this differences, in a sense, sound very similar to System Response, the only difference between them is that one is about model and the other one is about system. Since that we often use models to simulate the real world systems like ecosystem, hydrological system, financial system and etc, if these models would be regarded as approximation to real world systems, and we can presume that  errors in approximation could be ignored, then aforementioned model outputs’ differences can be called as System Response or, to be more precisely, simulation of System Response. 

0x01 What is System Response In Hydrological Model?

Given a system called y=f(x) , in this system, y is system’s output , x is system’s input and f represents system’s response to input. If we are talking about a hydrological system, this formula can be described more specifically as Q=f(P,E,θ), in which Q is streamflow, P is precipitation, E is evapotranspiration and θ is system’s parameters.

We can use Δy=f(x+Δx)-f(x) to denote the system response, if we use Taylor’s Formula to expand this, we can have

Δy=f(x+Δx)-f(x)=\frac{\partial f(x)}{\partial x} \Delta x+o(\Delta x)

If Δx, the input’s variation, is small enough, we can ignore o(\Delta x) , which is higher order infinitesimal of Δx, here we assume that it can be ignored.

Moreover, if we consider y is a continuous system’s output, then aforementioned formula can be described as

Δy(t)=f(x+Δx,t)-f(x,t)=\frac{\partial f(x,t)}{\partial x} \Delta x

Here we have the basis of system response curve, and Δy(t) is the System Response , note that x, Δx and \frac{\partial f(x,t)}{\partial x} can be vector, especially in hydrological system, latter I will explain why.

Now another question arises, how can we determined the specific output of f(x,t) when we are talking about complex system like hydrological system? As a matter of fact, we always use hydrological model to approximate that function’s output, and it is pretty easy for us to obtain hydrological model output, then to consider it as f(x,t) will solve the problem.

0x02 What System Response Can Do In Hydrological Model

You may found that system response is similar to unit hydrograph, but they’re different things in practice, unit hydrograph represents the watershed’s average runoff characteristics, it should be used along with net rain process (precipitation after deducting losses such as evapotranspiration, plant retention, recharge soil moisture and etc), but system response, on the other hand, can be utilized with many hydrological factors, including net rain, runoff and even evapotranspiration, because all these factors can be regarded as system’s inputs and can be used to generate system response. Then here we come to the ultimate question, what system response can do in hydrological model?

As mentioned before, hydrological system’s formula can be describe as Q=f(P,E,θ), any variation of input will create a system response, in other words, if we change precipitation, evapotranspiration and any kind of input, the difference in system’s output can be utilized to do marvelous things. They can be real-time flood forecast error correction and model parameters calibration.

Before we go any further, please allow me to clarify one point, as I mentioned before, in formula

Δy(t)=f(x+Δx,t)-f(x,t)=\frac{\partial f(x,t)}{\partial x} \Delta x,

Δx can be vector, especially in hydrological system. Because in hydrological system, any fluctuation of input at a certain time will have impact on a series of model’s output, unit hydrograph is a great representation of this phenomenon. In other words, model’s output at a certain time is affected by a series of inputs within a time range. Therefore, if we take this into account, the overall system response Δy(t) can be decomposed into a set of system responses to inputs within a time range, here we can call this System Response Curve (SRC). In this case, aforementioned formula should be converted to:

\Delta y(t)=f(x+\Delta x, t)-f(x, t)=\left\{\begin{array}{l}\frac{\partial f_{t=1}}{\partial x_{1}} \Delta x_{1}+\frac{\partial f_{t=1}}{\partial x_{2}} \Delta x_{2}+\cdots+\frac{\partial f_{t=1}}{\partial x_{N}} \Delta x_{N} \\\frac{\partial f_{t=2}}{\partial x_{1}} \Delta x_{1}+\frac{\partial f_{t=2}}{\partial x_{2}} \Delta x_{2}+\cdots+\frac{\partial f_{t=2}}{\partial x_{N}} \Delta x_{N} \\\vdots \\\frac{\partial f_{t=L}}{\partial x_{1}} \Delta x_{1}+\frac{\partial f_{t=L}}{\partial x_{2}} \Delta x_{2}+\cdots+\frac{\partial f_{t=L}}{\partial x_{N}} \Delta x_{N}\end{array}\right.,

in this formula, \operatorname{SRC}(n)=\left[\begin{array}{l}\frac{\partial f_{t=1}}{\partial x_{n}} \\\frac{\partial f_{t=2}}{\partial x_{n}} \\\vdots \\\frac{\partial f_{t=L}}{\partial x_{n}}\end{array}\right], \mathrm{n}=1,2, \ldots, \mathrm{N}.

And if we describe it in matrix, it will be like:

\Delta \mathrm{Y}=\left[\begin{array}{ccc}\frac{\partial f_{1}}{\partial x_{1}} & \cdots & \frac{\partial f_{1}}{\partial x_{N}} \\\vdots & \cdots & \vdots \\\frac{\partial f_{L}}{\partial x_{1}} & \cdots & \frac{\partial f_{L}}{\partial x_{N}}\end{array}\right]\left[\begin{array}{c}\Delta x_{1} \\\vdots \\\Delta x_{N}\end{array}\right]=U \Delta X

in this formula, U is called as System Response Matrix (SRM), which is crucial for us to perform parameters calibration and forecast errors correction on hydrological model.

Given one simple hydrological model


in this formula, I is variable model’s input like precipitation, θ is the only model parameter.

Because model’s output will usually be different to real system output, we call the difference as simulation error, and in hydrological model, that would be error of simulated streamflow to observed one:

\Delta Q(t)=Q_{o b s}(t)-Q_{c a l}(t),

in this formula, \Delta Q(t) is simulation error, Q_{o b s}(t) is observed streamflow and Q_{c a l}(t) is simulated streamflow.

Now the question is how we can correct simulation error. The answer is simple, now we know that the SRC is determined by SRM and input’s variation, and SRM is static for a hydrological model at a specific time point, so if we manage to find an input’s variation to generate System Reponse which is similar to simulation error, we can add this variation to model’s inputs, and finally generate model outputs which are as near to observed features as possible. Based on that, the equation we need to solve is:

\Delta Q(t)=\Delta Y=U \Delta X

Based on least squares estimation theory, the optimal estimation of \Delta X is

\overline{\Delta \mathrm{X}}=\left(\mathrm{U}^{\mathrm{T}} \mathrm{U}\right)^{-1} U^{T} \Delta Y,

note that SRM U should be invertible, otherwise this equation cannot be solved.

In practice, sometimes we may found U is invertible, there are some resorts which help to solve this problem, the most simple one is to transform previous formula into

\begin{aligned}\overline{\Delta \mathrm{X}}=&\left(\mathrm{M}^{\mathrm{T}} \mathrm{M}\right)^{-1} M^{T} \Delta Y \\& \mathrm{M}=\mathrm{U}+\alpha \mathrm{I}\end{aligned},

in this formula, I is Identity Matrix, and \alpha is called as ridge coefficient, this method introduce Ridge Estimation into matrix equation solving, for more details, check out reference at the end of post.

Now everything is clear, if we are able to calculate SRM and obtain simulation error, we can find a way to calculate model’s input variation, and this is going to be used for generation of corrected model output.

\mathrm{Y}_{\mathrm{new}}=f(X+\overline{\Delta X})

0x03 System Response on Parameters Calibration

When referring to hydrological model’s inputs, you may think about precipitation, evapotranspiration. But it shouldn’t be neglected that model’s parameters are also treat as a typical kind of model’s input, though they are static compared to other inputs. As a kind of model inputs, model’s parameters can be used to generate system response by aforementioned method. Now it’s time to solve an important question, how can we use system response to calibrate parameters?

Base on aforementioned theory, now it is easy for us to derive the formula to calibrate parameters of hydrological model

\begin{array}{l}\overline{\Delta \theta}=\left(\mathrm{U}^{\mathrm{T}} \mathrm{U}\right)^{-1} U^{T} \Delta Q \\\mathrm{U}=\left[\begin{array}{ccc}\frac{\partial f_{1}}{\partial \theta_{1}} & \cdots & \frac{\partial f_{1}}{\partial \theta_{N}} \\\vdots & \ddots & \vdots \\\frac{\partial f_{L}}{\partial \theta_{1}} & \cdots & \frac{\partial f_{L}}{\partial \theta_{N}}\end{array}\right] \\\Delta Q=\left[\begin{array}{c}Q_{o b s_{1}}-Q_{c a l_{1}} \\\vdots \\Q_{o b s_{L}}-Q_{c a l_{L}}\end{array}\right]\end{array}

0x04 System Response on Precipitation Correction

Precipitation is the most important factor which contribute to the generation of runoff. But according to relevant research, the measurement of precipitation always has errors, and this kind of errors will be propagated through hydrological model, along with model’s structure, parameters errors, and finally contribute to the output errors of hydrological model. Moreover, due to lack of precipitation gauges, in practical usage, we often take precipitation gauge’s observation records, which is point observation, as surface observation, this kind of approximation will introduce errors as well. Therefore, it is plausible to use System Response on precipitation correction.

Just like what we do on parameters calibration, the correspondent formula for precipitation correction will be like

\begin{array}{l}\overline{\Delta \mathrm{P}}=\left(\mathrm{U}^{\mathrm{T}} \mathrm{U}\right)^{-1} U^{T} \Delta Q \\\mathrm{U}=\left[\begin{array}{ccc}\frac{\partial f_{1}}{\partial P_{1}} & \cdots & \frac{\partial f_{1}}{\partial P_{N}} \\\vdots & \ddots & \vdots \\\frac{\partial f_{L}}{\partial P_{1}} & \cdots & \frac{\partial f_{L}}{\partial P_{N}}\end{array}\right] \\\Delta Q=\left[\begin{array}{c}Q_{o b s_{1}}-Q_{c a l_{1}} \\\vdots \\Q_{o b s_{L}}-Q_{c a l_{L}}\end{array}\right]\end{array}

0x05 Reference

  • [1]司伟,包为民,瞿思敏,石朋.基于面平均雨量误差修正的实时洪水预报修正方法[J].湖泊科学,2018,30(02):533-541.
  • [2]司伟,包为民,瞿思敏.洪水预报产流误差的动态系统响应曲线修正方法[J].水科学进展,2013,24(04):497-503.
  • [3]张小琴,刘可新,包为民,李佳佳,赖善证.产流误差比例系数的系统响应修正方法[J].水科学进展,2014,25(06):789-796.
  • [4]刘可新,张小琴,包为民,赵丽平,束慧连,李佳佳.产流误差平稳矩阵的系统响应修正方法[J].水利学报,2015,46(08):960-966.
  • [5]Si W, Gupta H V, Bao W, et al. Improved Dynamic System Response Curve Method for Real‐Time Flood Forecast Updating[J]. Water Resources Research, 2019, 55(9): 7493-7519.
  • [6]Si W, Bao W, Gupta H V. Updating real‐time flood forecasts via the dynamic system response curve method[J]. Water resources research, 2015, 51(7): 5128-5144.
  • [7]Sun Y, Bao W, Jiang P, et al. Development of Multivariable Dynamic System Response Curve Method for Real‐Time Flood Forecasting Correction[J]. Water Resources Research, 2018, 54(7): 4730-4749.

About the author


[stay(d) for d in ('determined','diligent','devoted')]

By EDLinus
Coding, Web, Hydrology and more.