Sensitivity equations for hyperbolic conservation law-based flow models

Advances in Water Resources 30 (2007) 1943–1961 www.elsevier.com/locate/advwatres

Sensitivity equations for hyperbolic conservation law-based flow models Vincent Guinot *, Matthieu Leme´nager, Bernard Cappelaere Hydrosciences Montpellier UMR 5569 (CNRS, IRD, UM1, UM2), Maison des Sciences de l’Eau, Universite´ Montpellier 2, 34095 Montpellier Cedex 5, France Received 19 September 2006; received in revised form 13 March 2007; accepted 15 March 2007 Available online 2 April 2007

Abstract The present paper focuses on the governing equations for the sensitivity of the variables to the parameters in flow models that can be described by one-dimensional scalar, hyperbolic conservation laws. The sensitivity is shown to obey a hyperbolic, scalar conservation law. The sensitivity is a conserved scalar except in the case of discontinuous flow solutions, where an extra, point source term must be added to the equations in order to enforce conservation. The propagation speed of the sensitivity waves being identical to that of the conserved variable in the original conservation law, the system of conservation laws formed by the original hyperbolic equation and the equation satisfied by the sensitivity is linearly degenerate. A consequence on the solution of the Riemann problem is that rarefaction waves for the variable of the original equation result in vacuum regions for the sensitivity. The numerical solution of the hyperbolic conservation law for the sensitivity by finite volume methods requires the implementation of a specific shock detection procedure. A set of necessary conditions is defined for the discretisation of the source term in the sensitivity equation. An application to the one-dimensional kinematic wave equation shows that the proposed numerical technique allows analytical solutions to be reproduced correctly. The computational examples show that first-order numerical schemes do not yield satisfactory numerical solutions in the neighbourhood of moving shocks and that higher-order schemes, such as the MUSCL scheme, should be used for sharp transients.  2007 Elsevier Ltd. All rights reserved. Keywords: Hyperbolic conservation laws; Model sensitivity; Riemann problem; Godunov-type schemes; Shock detection; Kinematic wave

1. Introduction Model calibration methodologies and sampling network design (that is, the definition of the optimal location of one or several sensors so as to increase the predictive capacity of a model) are crucial issues in hydrological modelling [4,5,19,31,32]. In the theory of optimal experiment design [2,24], the sampling network should be designed in such a way that the parameters of the model are determined with optimal accuracy. Conversely, calibration procedures should be conducted in such a way that only the parameters to which the model output is sensitive are modified while the parameters to which output sensitivity is small

*

Corresponding author. Tel.: +33 0 4 67 14 90 56; fax: +33 0 4 67 14 47

74. E-mail address: [email protected] (V. Guinot). 0309-1708/$ - see front matter  2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2007.03.004

should not be subject to calibration, otherwise leading to failure and possibly decreasing or destroying the predictive capability of the model [16]. Consequently, the determination of the sensitivity of model output to model parameters may be seen as a key step in the optimization of model performance. In what follows, a model is understood as a set of partial differential equations (PDEs) with associated parameters, initial and boundary conditions describing one or several physical processes, without any application to a specific region of space and time. The present paper focuses on the derivation of sensitivity equations for hydrological and hydrodynamic models that can be described using one-dimensional scalar, hyperbolic conservation laws with source terms and non-uniform parameters. Sensitivity equations are classically used in optimisation problems such as optimal structure design [20], optimal flow control [17,18] and design [11], data assimilation [11] or experiment

1944

V. Guinot et al. / Advances in Water Resources 30 (2007) 1943–1961

design [1]. The sensitivity of a model output to the parameters can be carried out in two possible ways: either (i) the equations that describe the evolution of the state variables of the model are differentiated with respect to the model parameters, or (ii) the differentiation is performed on the numerical and informatic implementation of the model (i.e. the computational code). The former approach is known as the continuous approach, while the latter is known as the discrete approach [13,20]. Although well-suited to practical and real-world situations, the discrete approach provides no theoretical insight into the general behaviour of a given model and must be repeated for any new situation. Therefore, it does not allow general guidelines to be derived for the calibration of a given model. Moreover, the discrete approach uses the numerical solution provided by the model, that is, an approximate solution of the equations embedded in the model. Consequently, sensitivities computed using the discrete approach suffer from the inevitable inaccuracies generated by the numerical scheme used to solve the PDEs in the original model. Conversely, the continuous approach, based on first- or second-order perturbation approaches, may become inaccurate when strongly nonlinear PDEs are to be solved (see e.g. [9,10] for recent examples of results given by the two approaches via an adjoint formulation). The present paper focuses on the application of the continuous sensitivity approach to flow models described by one-dimensional scalar, hyperbolic conservation laws. The purpose is to investigate the analytical properties of the sensitivity of the flow variables to the model parameters, in particular when the flow becomes discontinuous. The possible discontinuous character of the flow is an essential feature of non-linear hyperbolic PDEs. It triggers a number of issues in terms of existence and uniqueness of the solutions. As shown in the present paper, the system of PDEs satisfied by the flow variable and the sensitivity is linearly degenerate, which leads to a very unusual behaviour for the sensitivity, such as the discontinuous character of the sensitivity profile a the boundaries of simple waves and the presence of Dirac source terms across shock waves. Similar problems have been pointed out in [17], albeit focusing on the sensitivity to model initial conditions rather than model parameters. In contrast with the solutions presented in [17], the sensitivity equations presented hereafter do not require that the analytical solution of the flow problem be known. Section 2 states the problem and presents the perturbation method used to derive the sensitivity equations. The need for an extra source term in the sensitivity equations is acknowledged and the expression for this source term is derived. Section 3 focuses on the properties of the solution of the Riemann problem. The results of Sections 2 and 3 are used in Section 4 to propose a discretization of the flow and sensitivity equations in the framework of Godunov-type schemes. Section 5 demonstrates the validity of the approach using various computational examples. Section 6 is devoted to concluding remarks.

2. Governing equations 2.1. Definitions The present paper focuses on one-dimensional, scalar hyperbolic conservation laws, that is, partial differential equations (PDEs) that can be written in conservation form as oU oF þ ¼ S; ot ox

ð1Þ

where U is the conserved variable, F is the flux, S is the source term, and t and x are the time and space coordinates, respectively. F and S are assumed to be functions of U and a set of parameters /k (k = 1, . . . , m) that may be variable in space F ¼ F ðU ; /1 ; . . . ; /m Þ; S ¼ SðU ; /1 ; . . . ; /m Þ:

ð2Þ

Eq. (1) can be rewritten in nonconservation form by expanding the derivative oF/ox and using the chain rule oU oU þk ¼ S0; ot ox

ð3Þ

where k and S 0 are defined as follows: k¼

oF ; oU

S0 ¼ S 

m X oF o/m : o/m ox k¼1

ð4Þ

The PDE (1) is said to be hyperbolic if k is a real number. The quantity k, called the wave celerity or wave propagation speed, is the speed at which small variations in the quantity U propagate in the solution domain. Eq. (1) can also be written in the so-called characteristic form dU ¼ S0 dt

for

dx ¼ k: dt

ð5Þ

A well-known example of such a conservation law is the kinematic wave equation for free-surface flow [23,25] oh o þ ðKI 1=2 h5=3 Þ ¼ 0; ot ox

ð6Þ

where h is the water depth, I is the local ground slope and K is Strickler’s friction coefficient. The kinematic wave equation is a simplified version [23] of the 2 · 2 hyperbolic system of conservation laws known as the Saint–Venant equations [8], that describe the behaviour of long waves in shallow water [27]. Eq. (6) may be written in the form (1) and (2) by defining U, F, S and / as U ¼ h; 1=2

F ¼ /1 /2 h5=3 ; /1 ¼ K; /2 ¼ I; S ¼ 0:

ð7Þ

V. Guinot et al. / Advances in Water Resources 30 (2007) 1943–1961

The celerity k is obtained by substituting the definitions (7) into the definition (4) 5 5 k ¼ KI 1=2 h2=3 ¼ u; 3 3

ð8Þ

where u = F/U is the flow velocity.

The sensitivity equations being developed for one parameter at a time, the subscript k for the parameters /k is omitted without loss of generality in what follows. The relationships (2) are simplified into S ¼ SðU ; /Þ:

ð9Þ

The sensitivity equations are obtained by assuming that the parameter / is known certainly over the entire simulation domain, except over a region X, where it may be subject to adjustment for the purpose e.g. of calibration. No particular assumption is made about the nature of X in that it may be restricted to a subdomain of the solution domain as well as extend over the entire solution domain, depending on the problem to be solved. The adjustment of the parameter may be seen as a perturbation of the original value by a perturbation field / 0 (x) equal to zero everywhere except in the region X. Since most models are calibrated using piecewise constant values for the parameters in practice, the perturbation is assumed to be uniform, equal to /0 over X in what follows. The parameter /(x) then becomes (/ + / 0 )(x), with  0 for x 62 X; /0 ðxÞ ¼ /0 vðxÞ; vðxÞ ¼ ð10Þ 1 for x 2 X; where v(x) is the so-called characteristic function of X. The perturbation / 0 yields a perturbation U 0 in the solution U of the equation. The sensitivity s of U to / is defined as sðx; tÞ 

oU U 0 ðx; tÞ ¼ lim : o/0 /0 !0 /0

ð11Þ

The governing PDE for the sensitivity s is obtained by stating that the solution U + U 0 with parameter / + / 0 should also verify Eq. (1). Therefore, the following equality holds: o o ðU þ U 0 Þ þ F ðU þ U 0 ; / þ /0 Þ ot ox ¼ SðU þ U 0 ; / þ /0 Þ:

ð12Þ

As shown in Appendix A, a first-order Taylor series expansion is sufficient to derive the sensitivity equation under the assumption of an infinitesimal perturbation /0 (see Appendix A.1 for the details of the derivation) oF 0 oF 0 U þ /; oU o/ oS 0 oS 0 U þ /: SðU þ U 0 ; / þ /0 Þ  SðU ; /Þ þ oU o/

Substituting Eqs. (13) into Eq. (12), using the definition (4) for k and subtracting Eq. (1) yields   oU 0 o oS 0 oS 0 o 0 oF 0 U þ /  / : ð14Þ þ ðkU Þ ¼ ox oU o/ ox o/ ot Dividing by /0 and using the definition (11) yields the following equation:

2.2. Development of the equations for sensitivity

F ¼ F ðU ; /Þ;

1945

F ðU þ U 0 ; / þ /0 Þ  F ðU ; /Þ þ

ð13Þ

os o þ ðksÞ ¼ Q; ot ox

ð15Þ

where the source term Q is given by the right-hand side of Eq. (14) divided by /0. Substituting the definition (10) for / 0 (x) yields   oS oS o oF sþ vðxÞ  vðxÞ : ð16Þ QðxÞ ¼ oU o/ ox o/ 2.3. Weak solutions Weak solutions of Eq. (1) are solutions of the weak form of Eq. (1)  Z t2 Z x2  oU oF þ  S dx dt ¼ 0; ð17Þ ot ox t1 x1 where [t1, t2] and [x1, x2] are the integration bounds over space and time. Eq. (17) can also be rewritten as Z

x2

½U ðt2 Þ  U ðt1 Þ dx þ

x1

¼

Z

t2

Z

t2

½F ðx2 Þ  F ðx1 Þ dt t1

x2

S dx dt: t1

Z

ð18Þ

x1

Eq. (1) is a particular case of Eq. (18) when the size of the integration intervals [t1, t2] and [x1, x2] tends to zero. In contrast with the solutions of Eq. (1), the solutions of Eq. (18) do not need to be class C1 over space and time. Therefore, the class of weak solutions of Eq. (1) includes both the set of solutions of Eq. (1) and the additional set of discontinuous and non-differentiable solutions of Eq. (18). Discontinuous solutions may arise from three types of configurations: (i) the initial profile U(x, 0) is discontinuous, (ii) the boundary condition U(0, t) is discontinuous with respect to time, or (iii) the flux function is non-linear in U, which means that k is a function of U. In the latter case, situations where ok/ox is negative usually result in the derivative oU/ox reaching infinite values after a time given by [22]   1 ok tinf ¼  min : ð19Þ x ox Note that the three configurations above are not mutually exclusive. There exist two types of discontinuities, namely contact discontinuities and shock waves. A contact discontinuity verifies the following conditions: U L 6¼ U R ; kL ¼ kR ;

ð20Þ

1946

V. Guinot et al. / Advances in Water Resources 30 (2007) 1943–1961

where the subscripts L and R indicate that the values of the variable are taken on the left- and right-hand side of the discontinuity, respectively. A shock verifies the following criteria U L 6¼ U R ;

ð21Þ

kL > kR :

Note that this condition is valid regardless of the sign of kL and kR. In both cases, the discontinuity moves at a speed cs given by the so-called jump relationship, or Rankine– Hugoniot condition (see e.g. [7,15,21]) cs ¼

FL  FR : UL  UR

ð22Þ

The jump relationship is obtained from a mass balance over a control volume of infinitesimal width that contains the discontinuity (Fig. 1), as shown in Appendix A.2. Such a relationship, however, is not applicable to Eq. (15), as pointed out in [17,3]. Note that Eq. (22) may also become invalid when the source term S is non-regular, see [29]. If Eq. (22) was valid for the sensitivity, setting U = s and F = ks would yield the following equality: cs ¼

kL s L  kR s R : sL  sR

ð23Þ

Ns X os o þ ðksÞ ¼ Q þ dk q k ; ot ox k¼1

v ¼ Const;

ð24Þ

where both Eqs. (22) and (23) yield the same expression cs ¼ v:

ð25Þ

s

s1

dk ðxÞ ¼ dðx  xs;k Þ;

x2

ð27Þ

where xs,k is the x-coordinate of the kth discontinuity. The expression for the amplitude qk of the source term is determined as follows. Consider a single shock moving at the speed cs with left and right values s1 and s2, respectively. A control volume [x1, x2] is defined, in which the discontinuity is contained (Fig. 1). Integrating Eq. (26) over the control volume yields the following equality: Z x2 Z t2 ½sðx; t2 Þ  sðx; t1 Þ dx þ ½ðksÞðx2 ; tÞ  ðksÞðx1 ; tÞ dt ¼

Z t1

t2

Z

t1 x2

Qðx; tÞ dx dt þ x1

Z

t2

qðtÞ dt:

When both the size of the control volume and the integration time tend to zero, Eq. (28) tends to the following expression: ðt2  t1 Þq ¼ ðs1  s2 Þðt2  t1 Þcs þ ðt2  t1 Þðk2 s2  k1 s1 Þ  ðt2  t1 Þ½ðxs  x1 ÞQ1 þ ðx2  xs ÞQ2 ;

x

dx/dt = λ t2

x

Fig. 1. Control volume across a shock in the physical space (top) and the phase space (bottom).

ð30Þ

If the flow is continuous or if the wave is a contact discontinuity, k1 = k2 = cs, consequently q = 0. Therefore, Eqs. (26) and (30) are applicable to both continuous and discontinuous solutions. Generalising Eq. (30) to multiple discontinuities yields the final form of Eq. (15) for the sensitivity Ns X os o þ ðksÞ ¼ Q þ dk qk ; ot ox k¼1

t1 x2

ð29Þ

where Q1 and Q2 denote the value of the source term on the left- and right-hand side of the discontinuity, respectively, and xs is the abscissa of the discontinuity. Dividing by (t2–t1) and considering that both x1 and x2 tend to xs leads to the following expression for q ¼ ðcs  k1 Þs1 þ ðk2  cs Þs2 :

t

x1

ð28Þ

t1

q ¼ ðs1  s2 Þcs þ k2 s2  k1 s1

s2

x1 xs

ð26Þ

where Ns is the number of discontinuities in the solution domain, qk is the amplitude of the correction term for the kth discontinuity, and dk is defined as Dirac’s distribution for the kth discontinuity:

x1

Since U and s are independent variables, the values sL and sR are independent from UL and UR in the general case and Eqs. (22) and (23) cannot be satisfied simultaneously. The complete expression of the jump relationship for the sensitivity, given in Appendix A.3, is much more complex than Eq. (23). A trivial exception is that of the linear advection equation with a uniform advection velocity v F ¼ vU ;

Bearing in mind that cs is known from Eq. (22), it is proposed that the governing equation (15) should be modified so as to include an extra source term that takes effect only across discontinuities. The proposed formulation is the following

ð31Þ

qk ¼ ðcs;k  k1;k Þs1;k þ ðk2;k  cs;k Þs2;k : The numerical solution of Eq. (31) involves a specific procedure for shock detection, as discussed in Section 4.

V. Guinot et al. / Advances in Water Resources 30 (2007) 1943–1961

1947

3. The Riemann problem

3.2. Contact discontinuities and shock waves

3.1. Problem and wave patterns

When the conservation law is concave or convex, a contact discontinuity appears in the solution U if the wave celerities kL and kR associated with the left and right states UL and UR are identical

The Riemann problem is an essential tool in the finite volume solution of hyperbolic conservation laws. The purpose of the present section is to investigate the properties of its solution. The Riemann problem is an initial-value problem defined as follows oU oF þ ¼ 0; ot ox os o þ ðksÞ ¼ q; ot ox  U L for x < x0 ; U ðx; 0Þ ¼ U R for x > x0 ;  sL for x < x0 ; sðx; 0Þ ¼ sR for x > x0 ;

kL ¼ kR :

ð33Þ

A shock wave appears if the wave celerity of the left state is larger than the celerity of the right state ð34Þ

kL > kR :

ð32Þ

/ðxÞ ¼ const; where sL and sR are the (constant) left and right states of the Riemann problem in s respectively, UL and UR are the (constant) left and right states for the Riemann problem in s and U, respectively, and x0 is the location of the initial discontinuity. The parameter / in the flux function F(/, U) is assumed to be uniform over the entire solution domain, that is, v = 0. The source term S is assumed to be zero because only the hyperbolic part of the conservation law is considered. Therefore, the sensitivity source term Q is also zero. The source term q is the source term arising from the possible shocks in the solution U. For the sake of clarity, the summation operator has been omitted in the notation. The behaviour of s cannot be analysed independently from that of U because the flux ks is a function of U via the first definition (4). Note that the system formed by the PDEs in U and s is not a hyperbolic system of conservation laws strictly speaking because its two eigenvalues are identical (it is recalled that a m · m hyperbolic system of conservation laws has m real, distinct eigenvalues). The degenerate character of the system of conservation laws PDEs has consequences on the behaviour of the solution of the Riemann problem in s when the solution U contains rarefaction waves (see Section 3.3). The analysis of the behaviour of s is based on the behaviour of U. Four types of wave may be present in the solution U: contact discontinuities, shock waves, rarefaction waves and mixed (or compound) waves. These various patterns are analysed in the following subsections. The behaviour of contact discontinuities and shock waves is discussed under the assumption of a monotonous (i.e. convex or concave) conservation law in Section 3.2. The behaviour of rarefaction waves is analysed under the assumption of a monotonous conservation law in Section 3.3. The issue of non-convex conservation laws that may lead to compound waves is discussed in Section 3.4.

In both cases (see Fig. 2) the solution U is given by  U L for x < x0 þ ct; U ðx; tÞ ¼ ð35Þ U R for x > x0 þ ct; where c is the speed at which the discontinuity moves. In the case of a contact discontinuity, c is equal to the wave celerities kL and kR. In the case of a shock wave, c is equal to the shock speed cs defined as in Eq. (22). The parameter / being uniform over the solution domain, the wave celerity is also piecewise constant  kL for x < x0 þ ct; kðx; tÞ ¼ ð36Þ kR for x > x0 þ ct: And therefore the equation in s reduces to os os þk ¼0 ot ox

ð37Þ

that is ds ¼0 dt

for

dx ¼ k: dt

ð38Þ

Consequently, the solution s is also piecewise constant (Fig. 2)  sL for x < x0 þ ct; sðx; tÞ ¼ ð39Þ sR for x > x0 þ ct:

U,s

sL

U,s

sL sR

UL

sR UL

UR x0

UR

x

t

x0

x

x0

x

t

dx/dt = λ L

dx/dt = λ L

dx/dt = λ R

dx/dt = λ R x0

x

Fig. 2. Solution of the Riemann problem. Behaviour of the solution in the physical space (top) and the phase space (bottom) in the case of a contact discontinuity (left) and a shock wave (right).

1948

V. Guinot et al. / Advances in Water Resources 30 (2007) 1943–1961

3.3. Rarefaction waves A rarefaction wave appears if the wave celerity of the left state is smaller than the celerity of the right state (Fig. 3) ð40Þ

kL < kR :

In this case the solution U is shown to be self-similar and to depend only on the ratio (x  x0)/t [21] 8 for x < x0 þ kL t; > < UL e ð41Þ U ðx; tÞ ¼ U ðx; tÞ for x0 þ kL t < x < x0 þ kR t; > : UR for x > x0 þ kR t; e ðx; tÞ in the rarefaction wave verifies the where the profile U following condition x  x0 e ðx; tÞÞ: ¼ kð U ð42Þ t Eq. (42) can be solved uniquely for U when the flux function is convex or concave. Its implicit character most often requires that an iterative procedure be used when the flux function is nonlinear in U. In practice the shape of the proe ðx; tÞ is more easily determined by expressing x as a file U function of U and t rather than expressing U as a function of x and t as done usually (see e.g. [30] for an illustration of the method) e ; tÞ ¼ x0 þ kð U e Þt: xð U

ð43Þ

The wave celerity k is therefore given by 8 > < kL for x < x0 þ kL t; kðx; tÞ ¼ xxt 0 for x0 þ kL t < x < x0 þ kR t; > : kR for x > x0 þ kR t:

sL

U,s UR sR

UL x0 dx/dt = λ L

ð44Þ

x dx/dt = (x – x0)/t

t

dx/dt = λ R

x0

x

Fig. 3. Solution of the Riemann problem. Behaviour of the solution in the physical space (top) and the phase space (bottom) in the case of a rarefaction wave.

The solution s therefore becomes 8 for x < x0 þ kL t; > < sL sðx; tÞ ¼ ~sðx; tÞ for x0 þ kL t < x < x0 þ kR t; > : sR for x > x0 þ kR t;

ð45Þ

where the expression of ~sðx; tÞ remains to be determined. This is done by expanding the equation in s as follows: o~s o~s ok þ k ¼  ~s: ð46Þ ot ox ox Eq. (46) can be rewritten as d~s ok dx ¼  ~s for ¼ k: ð47Þ dt ox dt From the expression (44) for k in the rarefaction wave it follows that ok 1 dx ¼ for kL 6 6 kR : ð48Þ ox t dt Substituting Eq. (48) into Eq. (47) yields ~s d~s dx dx ¼ ¼ k; kL 6 6 kR : for ð49Þ dt dt dt t The solution of this ordinary differential equation is dð~stÞ ¼ 0

for

dx ¼ k; dt

kL 6

dx 6 kR : dt

ð50Þ

Therefore, the product ~sðx; tÞt is constant along a characteristic curve in the rarefaction wave. Since the rarefaction wave originates from t = 0, the product ~sðx; tÞt is equal to zero and remains so at further times. Consequently, ~sðx; tÞ is uniformly equal to zero in the rarefaction wave. Eq. (45) becomes (see Fig. 3) 8 > < sL for x < x0 þ kL t; sðx; tÞ ¼

> :

0 sR

for x0 þ kL t < x < x0 þ kR t; for x > x0 þ kR t:

ð51Þ

This unusual, discontinuous profile results directly from the degenerate character of the system of conservation laws in U and s. A similar behaviour has been observed for another degenerate, 2 · 2 system of conservation laws that describe traffic flow [15]. Note that this result can also be obtained via a balance in s over a variable control volume located between two abscissae x1 and x2 defined as  x 1 ¼ x 0 þ k1 t ð52Þ ; kL 6 k1 6 k2 6 kR : x 2 ¼ x 0 þ k2 t In other words, the control volume is contained in the rarefaction wave (Fig. 4). A balance over this control volume yields Z x1 d sðx; tÞ dx ¼ f1  f2 ; ð53Þ dt x2 where f1 and f2 are the fluxes across the left-and right-hand boundaries of the control volume, respectively. The sensitivity being advected at the speed k, f1 and f2 are given by the product of the local values of s and the local values of k relative to the boundaries of the control volume

V. Guinot et al. / Advances in Water Resources 30 (2007) 1943–1961

compound wave travelling to the right, the solution is given by 8 for x < x0 þ kL t; > < UL e ð56Þ U ðx; tÞ ¼ U ðx; tÞ for x0 þ kL t < x < x0 þ cs t; > : UR for x > x0 þ cs t;

s sL

λ1

λ2 sR x1

x2

e ðx; tÞ is given by Eq. (42) or Eq. (43). The propawhere U gation speed cs of the shock, that verifies the Rankine– Hugoniot condition, is also equal to the wave celerity immediately behind the shock (see [30] for more details)

x

t dx/dt = λ 2 dx/dt = λ L

cs ¼ kðU s Þ ¼ dx/dt = λ 2

F ðU s Þ  F ðU R Þ ; Us  UR

ð57Þ

where Us is the value of U immediately behind the shock. Reproducing the reasoning of Sections 3.2 and 3.3 across the shock and in the rarefaction wave, respectively, leads to the following profile for s

dx/dt = λ R x Fig. 4. Solution of the Riemann problem. Balance over a variable control volume extending over the rarefaction wave. Definition sketch in the physical space (top) and the phase space (bottom).

  dx1 f1 ¼ k1  s1 ; dt   dx2 f2 ¼ k2  s2 : dt

1949

ð54Þ

From the definition (52), f1 and f2 are both equal to zero. Consequently, Eq. (53) reduces to Z x1 sðx; tÞ dx ¼ const: ð55Þ x2

At t = 0, x1 = x2 = x0, therefore the integral in Eq. (55) is equal to zero. Since this equality holds for any combination of x1 and x2, bringing x1 and x2 arbitrarily close to each other by bringing k1 arbitrarily close to k2 yields the immediate consequence that the point value of s is equal to 0 everywhere within the rarefaction wave. Note that this particular behaviour of the sensitivity in the rarefaction wave stems from the assumption (that is inherent to the definition of the Riemann problem) of a zero sensitivity source term Q in the governing equation (32). If S depends on U and/or /, or if F depends on / and / is variable in space, or if v is nonzero, Q is nonzero (see Eq. (16)) and the rarefaction wave is not a vacuum region. 3.4. Compound waves Compound waves may arise when the flux function F is non-convex, that is, the wave celerity k is not a monotonic function of U. There exists at least one value Uex of U for which the celerity k reaches an extremum (this extremum may be local or global). If UL and UR are not located on the same side with respect to Uex, a compound wave appears in the solution. This wave is made of a rarefaction wave connected to a shock wave (Fig. 5). In the case of a

8 > < sL sðx; tÞ ¼ 0 > : sR

for x < x0 þ kL t; for x0 þ kL t < x < x0 þ cs t;

ð58Þ

for x > x0 þ cs t:

Conversely, a compound wave heading to the left yields the following solution for U 8 for x < x0 þ cs t; > < UL e U ðx; tÞ ¼ U ðx; tÞ for x0 þ cs t < x < x0 þ kR t; ð59Þ > : UR for x > x0 þ kR t; where cs is given by cs ¼ kðU s Þ ¼

F ðU s Þ  F ðU L Þ : Us  UL

ð60Þ

U,s UR

sL

sR

UL x0 dx/dt = λ L

x dx/dt = (x – x0)/t

t

dx/dt = λ R

x0

x

Fig. 5. Solution of the Riemann problem. Behaviour of the solution in the physical space (top) and the phase space (bottom) in the case of a compound wave.

1950

V. Guinot et al. / Advances in Water Resources 30 (2007) 1943–1961

The solution for s is given by 8 > < sL for x < x0 þ cs t; sðx; tÞ ¼ 0 for x0 þ cs t < x < x0 þ kR t; > : s for x > x0 þ kR t:

ð61Þ

4. Numerical solution 4.1. Finite volume discretisation A finite volume discretisation is proposed for Eqs. (1) and (31). The solution domain is discretised into computational cells (Fig. 6) over which the average value of the sought variable is to be computed every time step using the following formulae: Dt nþ1=2 nþ1=2 nþ1=2 U inþ1 ¼ U ni þ ðF  F iþ1=2 Þ þ S i Dt; Dxi i1=2 h i Dt nþ1=2 nþ1=2 nþ1=2 nþ1=2 sinþ1 ¼ sni þ ðksÞi1=2  ðksÞiþ1=2 þ qi1=2;i þ qiþ1=2;i Dxi nþ1=2

þ Qi

Dt;

ð62Þ where U ni and sni are the average values of U and s, respecnþ1=2 nþ1=2 tively, over the cell i at the time level n, S i and Qi are the average values of the source terms S and Q, respectively, over the cell i between the time levels n and nþ1=2 nþ1=2 n þ 1; F iþ1=2 , ðksÞiþ1=2 are the average values of the fluxes F and ks respectively at the interface i + 1/2 between the nþ1=2 cells i and i þ 1; qiþ1=2;i is the source term arising in the cell

i due to the possible presence of a shock appearing at the interface i + 1/2, Dxi is the width of the cell i, and Dt is the computational time step. Note that the finite volume formalism used in Eq. (62) implies that the parameter / should be estimated at the interfaces between the computational cells so as to allow the fluxes to be computed. nþ1=2 The flux F iþ1=2 between the cells i and i + 1 is computed by solving the following Riemann problem: oU oF þ ¼ 0; ot ox ( U iþ1=2;L n U ðx; t Þ ¼ U iþ1=2;R

for x < xiþ1=2 ;

ð63Þ

for x > xiþ1=2 ;

F ¼ F ðU ; /iþ1=2 Þ; where /i+1/2 is the value of / for the interface i + 1/2, and Ui+1/2,L and Ui+1/2,R are the left and right states of the Riemann problem, respectively. There are many ways of defining the left and right states of the Riemann problem. In Godunov’s scheme [12], the left and right states are defined as the average cell values U ni and U niþ1 , respectively, while in the MUSCL [28] or higher-order methods [6,14], Ui+1/2,L and Ui+1/2,R are defined as linear combinations of the average value of U over the three to five closest cells to the interface. Solving the Riemann problem (63) yields the value nþ1=2 U iþ1=2 at the interface i + 1/2 over the computational time step. The fluxes in U and s are estimated as nþ1=2

nþ1=2

F iþ1=2 ¼ F ðU iþ1=2 Þ; 8 nþ1=2 k siþ1=2;R > > > iþ1=2 < nþ1=2 ðksÞiþ1=2 ¼ 0 > > > : nþ1=2 kiþ1=2 siþ1=2;L nþ1=2

nþ1=2

if kiþ1=2 < 0; nþ1=2

if kiþ1=2 ¼ 0;

ð64Þ

nþ1=2

if kiþ1=2 > 0;

nþ1=2

kiþ1=2 ¼ kðU iþ1=2 Þ: Note that when the first-order Godunov method is used, the left and right states Ui1/2,L and Ui1/2,R used in Eq. (63) and the states si1/2,L and si1/2,R used in Eq. (64) are defined as (the discretisation is given here only for U but is strictly identical for s) U i1=2;L ¼ U ni1 ; U i1=2;R ¼ U ni :

ð65Þ

When the second-order MUSCL scheme [28] is used, the variable U is reconstructed using a linear profile in the cell i U ðx; tn Þ ¼ U ni þ ðx  xi ÞðU x Þni ;

ð66Þ

where xi is the abscissa of the centre of the cell i. The slope n ðU x Þi of the variable U over the cell i is initially defined as Fig. 6. Finite volume discretisation of the governing equations for U (a) and s (b). The bold lines in (b) indicate the propagation paths of possible shocks arising from the interfaces i  1/2 and i + 1/2.

n

ðU x Þi ¼ 2

U niþ1  U ni1 : Dxi1 þ 2Dxi þ Dxiþ1

ð67Þ

V. Guinot et al. / Advances in Water Resources 30 (2007) 1943–1961

It is then limited as follows 8 h i > n > > > max ðU x Þi ; li1=2 ; liþ1=2 > > > < n ðU x Þi  0 > > > h i > > n > > : min ðU x Þi ; li1=2 ; liþ1=2

( if

li1=2 < 0; liþ1=2 < 0;

if li1=2 liþ1=2 6 0; ( li1=2 > 0; if liþ1=2 > 0; ð68Þ

where the slope li1/2 is the maximum permissible slope for the variable U between the cells i  1 and i. It is defined so as to prevent a local extremum from occurring in the reconstructed profile between the cells i and i  1 li1=2 ¼ 2

U ni  U ni1 : Dxi1 þ Dxi

ð69Þ

This leads to the following definition of the left and right states of the Riemann problem nþ1=2

U i1=2;L ¼ U ni1 þ U i1=2;R ¼ U ni þ

1  Cri1=2;L

2 nþ1=2 Cri1=2;R  1 2

ðU x Þni1 ;

ð70Þ

n ðU x Þi ; nþ1=2

nþ1=2

where the Courant numbers Criþ1=2;L and Criþ1=2;R are defined as (see [26] for a detailed derivation) ! nþ1=2 kiþ1=2 Dt nþ1=2 ;0 ; Criþ1=2;L ¼ max Dxi1 ð71Þ ! nþ1=2 kiþ1=2 Dt nþ1=2 Criþ1=2;R ¼ min ;0 : Dxi 4.2. Discretisation of the ‘regular’ source terms S and Q The source term S is discretised using an explicit estimate over the cell i S ni ¼ SðU ni Þ:

ð72Þ

The source term Q is discretised as follows. Local explicit estimates are used for the average cell values of oS/oUs and v oS/o/ over the cell i  n oS oS ðU n ; / Þ; ¼ oU i oU i i ð73Þ  n oS oS n v ¼ ðU i ; /i Þvni : o/ i o/ The average value of the derivative o(v oF/o/)/ox over the cell i is estimated using an upwind approach. The term is evaluated as the sum of two contributions from the two interfaces i  1/2 and i + 1/2 of the cell i #  n "  n  n o oF o oF o oF v ¼ v v þ ; ð74Þ ox o/ i ox o/ i;i1=2 ox o/ i;iþ1=2 where the subscripts (i, i  1/2) and (i, i + 1/2) represent the contributions from the interfaces i  1/2 and i + 1/2 to the

1951

evaluation of the derivative within the cell i, respectively. In the upwind approach the contribution from the interface i  1/2 is zero if the wave celerity at the interface is negative (in other words if the characteristic passing at the interface is leaving the cell i). Otherwise it is assessed using an upwind estimate of the space derivative. The resulting estimate is 8 nþ1=2  n <0 if ki1=2 < 0; oF n n v ¼ ð75Þ v  oF v nþ1=2 o/ i;i1=2 : 2 ðoF o/ Þi ðo/ Þi1 if k P 0: Dxi1 þDxi

i1=2

Conversely, the contribution from the interface i + 1/2 is equal to zero if the wave celerity at the interface is positive, thus leaving the cell. Otherwise, it is again assessed using an upwind estimate of the derivative. The resulting estimate is 8 oF n oF n  n < 2 ðo/vÞiþ1 ðo/vÞi if knþ1=2 < 0; oF iþ1=2 Dxi þDxiþ1 v ð76Þ ¼ nþ1=2 o/ i;iþ1=2 : 0 if kiþ1=2 P 0: A necessary condition for the discretisation of the source term Q is that the sensitivity should remain equal to zero in any cell where v is zero, regardless of the value of s, v or Q in the neighbouring cells. The proposed discretisation fulfills this condition provided that oF/o/ is discretised as follows. Consider the following situation: vi1 ¼ 1; vi ¼ 0;

ð77Þ

kni1=2 > 0; kniþ1=2 > 0: Then the following conditions should hold  sni ¼ 0 8sni1 : snþ1 ¼0 i

ð78Þ

Substituting Eqs. (77) and (78) into Eq. (62) and simplifying yields  n 1 nþ1=2 n 2 oF 0¼ ki1=2 si1  : ð79Þ Dxi Dxi1 þ Dxi o/ i1 The proper estimate of oF/o/ satisfying Eq. (78) is therefore  n oF Dxi1 þ Dxi nþ1=2 n ¼ ki1=2 si1 : ð80Þ o/ i1 2Dxi For a regular cell width, Eq. (80) can be simplified into  n  nþ1=2  n oF oF oU nþ1=2 n ¼ ki1=2 si1 ¼ : ð81Þ o/ i1 oU i1=2 o/ i1 4.3. Discretisation of the point source term q nþ1=2

The point source term qi1=2;i is non-zero only if a shock is detected at the interface i + 1/2 and if this shock enters the cell i (see Section 4.4 for the details of the shock detecnþ1=2 tion procedure). In this case, qi1=2;i is given by an explicit

1952

V. Guinot et al. / Advances in Water Resources 30 (2007) 1943–1961

formulation directly derived from the application of Eq. (30) to the Riemann problem with left and right states si1/2,L and si1/2,R, respectively 8 if cs 6 0; > <0 n nþ1=2 ð82Þ qi1=2;i ¼ ðcs  ki1=2;L Þsi1=2;L > : þðkn  c Þs if c > 0 s i1=2;R s i1=2;R

h

and the shock speed cs is computed as

t

F ðU i1=2;L ; /i1=2 Þ  F ðU i1=2;R ; /i1=2 Þ cs ¼ U i1=2;L  U i1=2;R

where the left and right states and the shock speed are obtained by applying Eq. (83) to the interface i + 1/2 F ðU iþ1=2;L ; /iþ1=2 Þ  F ðU iþ1=2;R ; /iþ1=2 Þ : U iþ1=2;L  U iþ1=2;R

ð85Þ

4.4. Shock detection The purpose of the present subsection is to define a set of criteria that allows shocks possibly arising from the Riemann problem (63) to be identified. This issue is trivial when / is uniform over the solution domain. When / is variable in space, however, the problem becomes more complex, as illustrated by the following example. Consider the kinematic wave equation (6) with a uniform slope I and a piecewise constant Strickler coefficient  K 1 for x < xiþ1=2 ; KðxÞ ¼ ð86Þ K 2 for x > xiþ1=2 : Assume that steady state is reached. Mass conservation imposes that the following condition should be verified by the analytical solution Kh5=3 ¼

q I

1=2

¼ Const:;

ð87Þ

where q is the unit discharge over the solution domain. Therefore, the analytical solution h yields the following equality: ( hðx; tÞ ¼

h1 ¼ q3=5 ðI 1=2 K 1 Þ3=5

for x < xiþ1=2 ;

3=5

for x > xiþ1=2 :

h2 ¼ q

3=5

ðI

1=2

K 2Þ

K = K2 hi+1/2,L = h2

hi+1/2,L = h1

x dx/dt = λ

ð83Þ

A similar reasoning leads to the following expression for nþ1=2 qiþ1=2;i 8 n > < ðcs  kiþ1=2;L Þsiþ1=2;L nþ1=2 n qiþ1=2;i ¼ þðkiþ1=2;R  cs Þsiþ1=2;R if cs < 0 ð84Þ > : 0 if cs P 0;

cs ¼

K = K1

ð88Þ

Assume that K1 < K2. Then h1 > h2 and it stems directly from Eq. (8) that k1 < k2, as illustrated by Fig. 7. Therefore, there is no shock at the interface xi+1/2 because the wave celerity is larger on the right-hand side of the interface than on the left-hand side. Assume now that Eq. (6) is to be solved numerically, using Eq. (88) as an initial state. In this

xi+1/2

x

Fig. 7. Example of a steady state configuration for the kinematic wave equation with a piecewise constant Strickler coefficient increasing downstream. Schematisation in the physical space (top) and in the phase space (bottom).

case the following Riemann problem is defined at the interface i + 1/2 oh o þ ðK iþ1=2 I 1=2 h5=3 Þ ¼ 0; ot ox  h1 for x < xiþ1=2 ; n hðx; t Þ ¼ h2 for x > xiþ1=2 ;

ð89Þ

where Ki+1/2 is a consistent estimate of K at the interface (e.g. the arithmetic mean of K1 and K2). The numerical estimates of k on the left- and right-hand sides of the interface are given by 5 2=3 kiþ1=2;L ¼ K iþ1=2 I 1=2 h1 ; 3 5 2=3 kiþ1=2;R ¼ K iþ1=2 I 1=2 h2 3

ð90Þ

and since h1 > h2, ki+1/2,L is identified to be larger than ki+1/2,R, which results in a shock being detected while there is none in reality. This example shows that a shock identification procedure based only on the left and right states of the Riemann problems at the cell interfaces can be misleading when the parameters of the flux function are variable in space. The proposed shock detection procedure is based on the following consideration: if a shock crosses the interface i + 1/2 from right to left during the time step between the time levels n and n + 1 (Fig. 8a), the wave celerity is smaller on the right-hand side of the shock than on the left-hand side of the shock. Therefore, the wave celerity at the interface is smaller after the shock has crossed the interface and the solution of the Riemann problem at the interface at the time level n + 1 gives a smaller value of k than the solution of the Riemann problem at the time level n. Consequently, the following condition is verified   nþ3=2 nþ1=2 k U iþ1=2 ; /iþ1=2 < k U iþ1=2 ; /iþ1=2 < 0; ð91Þ

V. Guinot et al. / Advances in Water Resources 30 (2007) 1943–1961

t

λ (x, t n+1)

t xi+1/2

1953

λ (x, t n+1) λ (x, t ) n

x

λ (x, t n)

t

xi+1/2

x

t

xi+1/2

x

λ (x, t n+1)

xi+1/2

t

λ (x, t n) λ (x, t n+1)

λ (x, t n) xi+1/2

x

x

xi+1/2

x

t

t

xi+1/2

x

xi+1/2

x

Fig. 8. Configurations where a shock should be detected by the procedure: negative celerities on both sides of the shock (a), positive celerities on both sides (b), positive celerity on the left-hand side, negative celerity on the right-hand side with a mobile (c) or stationary (d) shock. nþ1=2

where U iþ1=2 is the solution of the Riemann problem at the interface i + 1/2 between the time levels n and n + 1. Conversely, if a shock crosses the interface from left to right (Fig. 8b), the wave celerity increases from the time level n to the time level n + 1, that is nþ1=2

nþ3=2

0 < kðU iþ1=2 ; /iþ1=2 Þ < kðU iþ1=2 ; /iþ1=2 Þ

ð92Þ

A third case is encountered when the wave celerities on the left- and right-hand sides of the interface have opposite signs kðU iþ1=2;L ; /iþ1=2 Þ > 0; kðU iþ1=2;R ; /iþ1=2 Þ < 0:

ð93Þ

The shock may be mobile (Fig. 8c) or immobile, located on the interface (cs = 0, Fig. 8d). If any of the conditions (91), (92) or (93) is fulfilled, a shock is detected at the interface i + 1/2 between the time level n and the time level n + 1. The shock celerity is computed using Eq. (85). In the partic-

ular case where cs is found to be zero, the shock must be ‘assigned’ arbitrarily to either the cell i or i + 1. If the shock is assigned to the cell i, the value of k used for the computation of the flux ks at the interface i + 1/2 must be computed using Ui+1/2,R. Conversely, if the shock is assigned to the cell i + 1, the value of k must be computed using Ui+1/2,L. 5. Computational examples 5.1. Test 1. Propagation of a sensitivity wave over a domain with constant parameters The objective of this purely numerical test case is to assess the performance of the shock detection algorithm proposed in Section 4.4 under transient conditions. The sensitivity approach is applied to the kinematic wave equation presented in Section 2.1. The expressions for the flux, the source term and the wave celerity are recalled

1954

V. Guinot et al. / Advances in Water Resources 30 (2007) 1943–1961

U ¼ h;

s 1=2 5=3

F ¼ KI h ; 5 k ¼ KI 1=2 h2=3 ; 3 S ¼ 0:

sus

ð94Þ

The test case consists in simulating the flow and sensitivity pattern resulting from sudden variations in the water depth at the upstream boundary of the domain. The Strickler coefficient and the bottom slope are assumed uniform all over the domain. The initial water depth is constant, equal to h0 and the sensitivity is initially zero. At t = 0 s, the water depth rises instantaneously to a fixed value hus at the upstream boundary of the domain. This corresponds to an injection at a constant unit discharge. The injection stops instantaneously at the time ts. The sensitivity s at the upstream end of the domain is kept to a constant value sus all throughout the simulation. The parameters of the test case are given in Table 1. The analytical profile for the water depth is the following. The instantaneous rise to a depth hus yields a shock wave (thick line in Fig. 9) propagating at the constant speed cs given by Eq. (22). The instantaneous decrease to a zero depth at t = tus yields a rarefaction wave, the head of which propagates at a speed kus = k(hus). The rarefaction wave eventually catches up the shock wave at a time tc = kusts/(kus  cs). Before this time, the two waves are connected by a region of constant state with h=hus. To summarise, the profile for h is the following 8 ~ > < hðx; tÞ for x 6 ðt  ts Þkus ; hðx; tÞ ¼ hus ð95Þ for ðt  ts Þkus 6 x < cs t; > : h0 for x > cs t; where kus = k (hus) and the profile ~ hðx; tÞ in the rarefaction wave is given by hðx; tÞÞ: xð~ hðx; tÞÞ ¼ ðt  ts Þkð~ The analytical 8 > <0 sðx; tÞ ¼ sus > : 0

ð96Þ

sensitivity profile is piecewise constant for x 6 ðt  ts Þkus ; forðt  ts Þkus 6 x < cs t;

ð97Þ

for x > cs t:

Table 1 Parameters for Test 1 Symbol

Meaning

Value

h0 hus I K s0 sus Tmax ts Dt Dx

Initial water depth Upstream water depth Bottom slope Strickler’s friction coefficient Initial sensitivity Upstream boundary condition for the sensitivity Length of the simulation Injection end time Computational time step Cell width

0m 0.01 m 102 40 m1/3 s1 0 m2/3 s 1 m2/3 s 300 s 200 s 1s 1m

x h hus

x t

ts

h = hus h=0 x

Fig. 9. Test 1. Behaviour of the analytical solutions for h and s in the physical space (top, middle) and in the phase space (bottom).

Fig. 10 shows the water depth and the sensitivity profiles computed using Eqs. (62) in combination with the first-order Godunov method (65) and the MUSCL scheme (70) at t = 300 s. The water depth profile computed by the Godunov scheme exhibits the classical front smearing associated with highly dissipative schemes. Such smearing is reduced to a large extent when the MUSCL scheme is used. The effects of numerical diffusion on the computed sensitivity profile are even more dramatic. Owing to the large numerical diffusion induced by the small Courant number in the rarefaction wave, the tail of the sensitivity wave is smeared excessively in the solution computed by the Godunov scheme. The sensitivity profile computed by the MUSCL scheme is much closer to the analytical solution. Note that the width of the shock identified by the detection algorithm at t = 300 s is 16 cells for the Godunov scheme, against seven cells for the MUSCL scheme. This is illustrated by Fig. 11, where the ‘shock region’ is represented by black dots in the (discretised) phase space. 5.2. Test 2. Propagation of a sensitivity wave over a domain with piecewise constant parameters The objective of this test is to assess the performance of the shock detection algorithm over a domain with piecewise constant parameters. The sensitivity approach is applied to the kinematic wave equation presented in Section 2.1. The test case (see Table 2) consists in simulating the behaviour of a water and sensitivity front propagating

V. Guinot et al. / Advances in Water Resources 30 (2007) 1943–1961

Analytical

h (m) 0.01

Analytical

s (m2/3 s)

Numerical

Numerical

1.0

x (m)

0.00

1955

x (m) 0.0

0 h (m)

50

0

Analytical

s (m2/3 s)

Numerical

0.01

x (m)

0.00 0

50 Analytical Numerical

1.0

x (m) 0.0 0

50

50

Fig. 10. Test 1. Comparison of the analytical and numerical profiles for the water depth (left) and the sensitivity (right) at t = 300 s. Results given by the original Godunov scheme (top) and the MUSCL scheme (bottom).

300

300

t (s)

0 0

50

x (m)

t (s)

0 0

50

x (m)

Fig. 11. Test 1. Regions where shocks are identified by the proposed detection algorithm. The shaded area indicates the location of the shocks in the discretised phase space for the Godunov scheme (left) and the MUSCL scheme (right).

into an initially dry domain where the Strickler coefficient is equal to a constant value K1 everywhere, except in a subdomain delimited by abscissae x1 and x2 where it is equal to a smaller value K2. The initial water depth and sensitivity are both equal to zero. At t = 0 s, a constant water depth hus and a constant sensitivity sus are prescribed at the upstream boundary of the computational domain. This triggers the propagation of a shock wave into the domain. When the wave reaches the subdomain where the Strickler coefficient is smaller (thus indicating high friction), continuity imposes that the water level should increase because the flux F given by the second equation (94) is uniform behind the shock. The water level h2 in the high friction subdomain verifies 5=3

5=3 ¼ K 2 I 1=2 h2 K 1 I 1=2 hus

and therefore  3=5 h2 K1 ¼ : hus K2

ð98Þ

ð99Þ

The wave celerity k2 in the high friction subdomain is related to the celerity k1 in the rest of the domain via the third equation (94)  2=3  3=5 k2 K 2 h 2 K2 ¼ ¼ : ð100Þ k1 K 1 h us K1 The sensitivity s2 is related to the sensitivity s1 = sus in the rest of the domain by the following relationship:  3=5 s2 K1 ¼ : ð101Þ s1 K2 Fig. 12 shows the analytical solution (solid line) for the water depth and the sensitivity at t = 400 s. The numerical solutions computed using the Godunov and the MUSCL schemes are indicated by the dots. Fig. 13 shows the regions of the discretised phase space where a shock has been detected by the detection procedure. At t = 400 s, the width of the shock is 20 cells for the Godunov scheme, against 7 cells for the MUSCL scheme. The reduced wave celerity in

1956

V. Guinot et al. / Advances in Water Resources 30 (2007) 1943–1961

Table 2 Parameters for Test 2 Symbol

Meaning

Value

h0 hus I K1 K2

Initial water depth Upstream water depth Bottom slope Strickler’s friction coefficient in the domain Strickler’s friction coefficient in the high friction subdomain Initial sensitivity Upstream boundary condition for the sensitivity Length of the simulation Abscissa of the upstream boundary of the high friction subdomain Abscissa of the downstream boundary of the high friction subdomain Computational time step Cell width

0m 0.01 m 102 40 m1/3 s1 20 m1/3 s1

s0 sus Tmax x1 x2 Dt Dx

0 m2/3 s 1 m2/3 s

5.3. Test 3. Rainfall-runoff model with infiltration

400 s 20 m

The following model is considered (see Table 3)

40 m

U ¼ h;

1s 1m

F ¼ KI 1=2 h5=3 ; 5 k ¼ KI 1=2 h2=3 ; 3 S ¼ Rðx; tÞ  k inf h;

the high friction subdomain leads to a reduced Courant number, which triggers an increased numerical diffusion and an increased smearing of the front. The Godunov scheme does not allow the initial shock width to be recovered after the wave has left the subdomain, while the MUSCL scheme allows the initial shock width to be retrieved. This contrast in the behaviour of the Godunov and the MUSCL schemes triggers again a large contrast in the behaviour of the sensitivity profiles computed by the two schemes. While the MUSCL scheme yields a solution in close agreement with the analytical solution, the Godunov scheme fails to compute a correct value for the sensitivity in the neighbourhood of the shock. This is because, due to a

Analytical

h (m)

0.02

large smearing of the sensitivity profile, the magnitude of the source term q is overestimated when the shock wave crosses the region [x1, x2]. As a consequence, too much sensitivity is lost. Apart from this problem with front resolution when the Godunov scheme is used, the behaviour of the analytical solution is accounted for correctly by the numerical solution. Moreover, the shock detection procedure allows the location of the shock to be identified correctly.

ð102Þ

where R(x, t) is time- and space-dependent rainfall intensity and kinf is an infiltration constant. This simple kinematic wave-based model is classical in rainfall-runoff and catchment modelling [23,25]. The solution domain is the interval [0, L]. The behaviour of the sensitivity s of h with respect to the Strickler coefficient K is investigated. The initial and boundary conditions of the simulation are the following: hðx; 0Þ ¼ 0

8x;

sðx; 0Þ ¼ 0 8x; hð0; tÞ ¼ 0 8t; sð0; tÞ ¼ 0

ð103Þ

8t:

s (m2/3 s)

2.0

Analytical

Numerical

Numerical

0.01

0.0

0.00 0

50

Analytical

h (m)

0.02

x (m )

0

50 s (m2/3 s)

2.0

x (m )

Analytical

Numerical

Numerical

1.0

0.01

0.0

0.00

0

50

x (m )

0

50

x (m )

Fig. 12. Test 2. Comparison of the analytical and numerical profiles for the water depth (left) and the sensitivity (right) at t = 300 s. Results given by the original Godunov scheme (top) and the MUSCL scheme (bottom).

V. Guinot et al. / Advances in Water Resources 30 (2007) 1943–1961

400

t (s)

400

0 0

x (m)

50

1957

t (s)

0 0

50

x (m)

Fig. 13. Test 2. Regions where shocks are identified by the proposed detection algorithm. The shaded areas represent the location of the shocks in the discretized phase space for the Godunov scheme (left) and the MUSCL scheme (right). The dashed lines indicate the boundaries of the subdomain where the Strickler coefficient is reduced.

Such conditions describe the runoff generated over a hillsope from initially dry conditions with no inflow from the upstream part of the slope. The following simple parameter configuration is assumed KðxÞ ¼ K 0

8x;

k inf ðxÞ ¼ 0 8x; IðxÞ ¼ I 0 8x; Rðx; tÞ ¼ R0

ð104Þ

8ðx; tÞ:

where K0, I0 and R0 are constants. Although very simple (and seemingly oversimplified), the assumption of uniform (or piecewise constant) parameters is not unusual in the field of catchment modelling, because large areas are usually involved, with very little opportunities and/or financial means for parameter estimation. In that sense the distribution (104) is legitimised by engineering practice. Moreover, Eqs. (103) and (104) allow a steady-state solution to be determined for h and s. The numerical values adopted for the parameters of the test are given in Table 3. The following calibration problem is simulated: while I0 is assumed to be known precisely, the value of K0 has been estimated very roughly over a subdomain [x1, x2] of the solution domain [0, L]. At the initial stage of the calibration process it is assumed to be equal to the average value Table 3 Parameters for Test 3 Symbol

Meaning

Value

h0 hus I K1 L s0 sus

Initial water depth Upstream water depth Bottom slope Strickler’s friction coefficient in the domain Length of the computational domain Initial sensitivity Upstream boundary condition for the sensitivity Length of the simulation Abscissa of the upstream boundary of the calibration subdomain Abscissa of the downstream boundary of the calibration subdomain Computational time step Cell width

0m 0m 102 40 m1/3 s1 100 m 0 m2/3 s 0 m2/3 s

Lmax x1

x2 Dt Dx

400 s 10 m

K0 but it could be subject to change in order to improve the calibration. In order to refine the estimate of K0 over [x1, x2] it is necessary to measure the depth h at a number of locations in the model under a precipitation event, the (uniform and constant) intensity R0 of which is known. Owing to financial and equipment constraints, the number of locations where the measurements are to be taken is to remain limited in practice. The optimal measurement location should therefore be that where the sensitivity is the highest in the model. Note that the parameter K being modified over the subdomain [x1, x2] yields the following distribution for the support function v(x)  vðxÞ ¼

1

if x1 6 x 6 x2 ;

0

otherwise:

ð105Þ

Solving Eqs. (1) and (31) with the parameter distribution given by Eq. (104), the initial and boundary conditions (103) and the parameter perturbation (105) under the assumption of steady state reached after an infinite time yields the following ordinary differential equations in h and s dh ¼ R0 ; dx   d o F ðksÞ ¼  vðxÞ : dx ox K

k

ð106Þ

Note that from the definition (102), the derivative of the source term in ks is   o F R0 1=2 5=3 1=2 5=3 vðxÞ ¼ I 0 h1 dðx  x1 Þ  I 0 h2 dðx  x2 Þ þ vðxÞ; ox K K ð107Þ where h1 and h2 are the water depths at x = x1 and x = x2, respectively and d(x) is Dirac’s distribution. The analytical solution of the system (106) is 

3=5 R0 x ; KI 1=2 8 > <0 R0 k1 s1  kK ðx  x1 Þ sðxÞ ¼ kðxÞ > : 0 hðxÞ ¼

80 m 1s 1m

for x < x1 ; for x1 < x < x2 ; for x > x2 ;

ð108Þ

1958

V. Guinot et al. / Advances in Water Resources 30 (2007) 1943–1961 3.E-02

s (m2/3 s)

Analytical

x (m )

0.E+00

Numerical

0

50

100

h (m)

x (m ) 0.E+00 0

50

100

-5.E-04

Analytical Numerical

Fig. 14. Test 3. Steady state water depth and sensitivity profiles for the parameters in Table 3.

where s1 is the value of s at x ¼ xþ 1 (the + superscript indicates that the abscissa is approached from the right), and k1 is the value of k at x = x1. Eqs. (102), (107) and (108) lead to the following expressions for k(x) and s1 5 1=2 3=5 2=5 kðxÞ ¼ ðK 0 I 0 Þ ðR0 xÞ ; 3 1=2 5=3 I h s1 ¼  0 1 : k1

ð109Þ

Note that it follows immediately from the expression of s in Eq. (108) that s is a negative, monotonically decreasing function of x over the interval [x1, x2], in other words djsj > 0 for x1 < x < x2 : dx

ð110Þ

Consequently, the optimal point for a measurement of h in view of the calibration of K over the subdomain [x1, x2] is x = x2 because this is the point where the absolute value of the sensitivity s of h with respect to K is maximum. The water depth and the sensitivity profiles computed by the first-order Godunov scheme are compared with the analytical solutions (108) and (109) in Fig. 14. A very good agreement is observed with the analytical solution. 6. Discussion and conclusions The main results of the present study may be summarised as follows: (i) the governing PDEs for the sensitivity s of one-dimensional scalar, hyperbolic conservation laws are hyperbolic PDEs, (ii) the system formed by the original conservation law and the PDE for the sensitivity is not strictly hyperbolic but linearly degenerate, (iii) a specific source term in the form of a Dirac function must be added to the sensitivity PDE when the solution of the original conservation law becomes discontinuous, (iv) an analytical solution is available for the Riemann problem in regions where the parameters are known with certainty (i.e. where the support function v is equal to zero), (v) owing to the linearly degenerate character of the system, a rarefaction wave for the solution U of the original conservation law in regions where v = 0 coincides with a vacuum region for the sensitivity s (that is, s = 0), (vi) efficient shock-capturing schemes can be proposed to solve the governing equations numerically, provided that necessary conditions for source term discretisation are satisfied and (vii) a specific shock detection pro-

cedure is available, and (viii) steady-state sensitivity solutions are correctly computed by the first-order Godunov scheme but the use of higher-order schemes such as the MUSCL scheme is strongly advised for the solution of transient problems involving shock propagation. The vacuum region in the sensitivity profile is one of the most striking features of the solution of the Riemann problem when the source term Q is zero. This feature is directly related to the linearly degenerate character of the system formed by the original conservation law and the governing PDE for the sensitivity. The equality between the wave propagation speed of the solution U of the original conservation law and the sensitivity s indeed prevents the sensitivity from crossing the boundary of the rarefaction wave, thus preserving the initial zero integral of the sensitivity within the rarefaction wave. It is not clear yet whether this specific property is preserved in the extension of the sensitivity approach to hyperbolic systems of conservation laws. The properties of the solutions of the sensitivity Riemann problem, the extension of the proposed discretization techniques and the necessary generalisation of the shock detection algorithm to flow models described by hyperbolic systems of conservation laws (such as the Saint–Venant or the two-dimensional shallow water equations) are the subject of ongoing research. Acknowledgements This work was supported by Research Contract ANR05-PGCU-004, ‘‘RIVES’’ funded by the French Agence Nationale de la Recherche/National Research Agency. The contribution and support of the following partners is gratefully acknowledged: Cemagref, CETE Me´diterrane´e, CETMEF, LMFA UMR 5509 (CNRS, ECL, INSA Lyon, UCB), INSAVALOR, LCPC, LRPC Bordeaux, LRPC Clermont-Ferrand, Sogreah, URGC (INSA Lyon), UTC. The constructive comments and discussion with Dr J.A. Cunge are gratefully acknowledged. Appendix A. Derivation of the sensitivity equations A.1. Continuous case The present subsection is devoted to the detailed derivation of Eq. (15). Eq. (12) is recalled

V. Guinot et al. / Advances in Water Resources 30 (2007) 1943–1961

o o ðU þ U 0 Þ þ F ðU þ U 0 ; / þ /0 Þ ¼ SðU þ U 0 ; / þ /0 Þ ot ox with the function / 0 (x) given by Eq. (10)  0 for x 62 X; 0 / ðxÞ ¼ /0 vðxÞ; vðxÞ ¼ 1 for x 2 X:

Dividing by /0 yields the following equation   os o oS oS o oF þ ðksÞ ¼ sþ v v þ /0 aðv; s; wÞ; ot ox oU o/ ox o/ aðv; s; wÞ ¼ RS  RF ;

0

0

A second-order Taylor series expansion in U and / leads to U0 

1959

oU o2 U /0 þ 2 /20 ; o/ o/

oF 0 oF 0 U þ / oU o/ o2 F U 02 o2 F /02 o2 F U 0 /0 ; þ 2 þ þ 2 oU o/ oU 2 o/ 2 oS 0 oS 0 U þ / SðU þ U 0 ; / þ /0 Þ  SðU ; /Þ þ oU o/ 2 o2 S U 0 o2 S /02 o2 S U 0 /0 : þ þ þ oU o/ oU 2 2 o/2 2 ðA:1Þ

F ðU þ U 0 ; / þ /0 Þ  F ðU ; /Þ þ

Substituting Eq. (10) and the definition (11) of the sensitivity into Eq. (A.1) and retaining only the terms up to the second order leads to   oF /20 0 0 s/0 þ w F ðU þ U ; / þ / Þ  F ðU ; /Þ þ oU 2 oF 1 o2 F 2 /0 v þ ðs/0 Þ o/ 2 oU 2 2 o2 F ð/ vÞ o2 F vs/20 ; þ þ 2 0 oU o/ 2 o/   oS /20 0 0 s/0 þ w SðU þ U ; / þ / Þ  SðU ; /Þ þ oU 2

ðA:5Þ where a(v, s, w) is a second-degree polynomial in v and s. It is easy to check that the first equation (A.5) becomes equivalent to Eq. (16) when the amplitude /0 of the perturbation tends to zero. This is because /0 a(v, s, w) tends to zero when the amplitude /0 of the perturbation tends to zero. Therefore all the terms in /0 a(v, s, w), that stem from the second-order terms in the Taylor series expansions for F and S, vanish in the limit of vanishing /0 and only the first-order terms of the Taylor series expansion remain in the final equation. A.2. Jump relationships for the flow variable The present subsection is devoted to the details of the derivation of the jump relationship. Consider a shock with left and right states UL and UR, respectively. Let x1 and x2 be the abscissae of the shock at times t1 and t2, respectively (Fig. 1). A balance over the control volume [x1, x2] between t1 and t2 yields the following equation: Z x2 Z x2 U ðx; t2 Þ dx  U ðx; t1 Þ dx x1

þ

¼

x1

Z

t2

F ðx1 ; tÞ dt 

t1

ðA:2Þ

oS 1 o2 S 2 /0 v þ ðs/0 Þ o/ 2 oU 2 2 o2 S ð/ vÞ o2 S vs/20 ; þ þ 2 0 oU o/ 2 o/ þ

where w is the second-order derivative of U with respect to /0, or the first-order derivative of s with respect to /0. Substituting Eqs. (A.2) into Eq. (12) and subtracting Eq. (1) leads to   os o oF oF /0 þ s/ þ /v ot ox oU 0 o/ 0 oS oS s/ þ / v þ ðRS  RF Þ/20 ; ¼ ðA:3Þ oU 0 o/ 0 where RF and RS contain the remaining terms in the Taylor series expansions for F and S, respectively   o oF w s2 o2 F o2 F v2 o2 F þ vs ; þ RF ¼ þ ox oU 2 2 oU 2 o/2 2 oU o/ ðA:4Þ oS w s2 o2 S o2 S v2 o2 S RS ¼ þ vs: þ þ oU 2 2 oU 2 o/2 2 oU o/

Z t1

t2

F ðx2 ; tÞ dt þ

Z

t2

Z

x2

Sðx; tÞ dx dt: t1

x1

ðA:6Þ The profile is approximated using a Taylor series expansion on both sides of the shock 8 < U L þ ðx  x1 ÞU ðxÞ for x < x1 ; L þ HOT1 ðx  x1 Þ U ðx; t1 Þ ¼ : ðxÞ U R þ ðx  x2 ÞU R þ HOT2 ðx  x2 Þ for x > x1 ; 8 ðxÞ ðtÞ > U L þ ðx  x1 ÞU L þ ðt2  t1 ÞU L > > > < þHOT ðx  x ; t  t Þ 3 1 2 1 U ðx; t2 Þ ¼ ðxÞ ðtÞ > > U R þ ðx  x2 ÞU R þ ðt2  t1 ÞU R > > : þHOT4 ðx  x2 ; t2  t1 Þ

ðA:7aÞ

for x < x2 ; for x > x2 ðA:7bÞ

ðxÞ UL

ðxÞ UR

where and are the derivatives of U with respect to ðtÞ ðtÞ x at x = x1 and x = x2, respectively, and U L and U R are the derivatives of U with respect to time. The integrals of U over the control volume at the times t1 and t2 are therefore given by Z x2 U ðx; t1 Þ dx x1 h i x2  x1 ðxÞ ¼ ðx2  x1 Þ U R þ U R þ HOT5 ðx2  x1 Þ ; 2 ðA:8aÞ

1960

Z

V. Guinot et al. / Advances in Water Resources 30 (2007) 1943–1961

x2

U ðx; t2 Þ dx x1 h x2  x1 ðxÞ ðtÞ U þ ðt2  t1 ÞU L ¼ ðx2  x1 Þ U L þ 2 i L þHOT6 ðx  x1 ; t2  t1 Þ :

ðA:8bÞ

UL and UR, respectively. This also causes a variation c0s in the speed cs of the shock. Since the shock moves at a different speed, its abscissa xs is subject to a change x0s . Therefore, the perturbed values of U and F on the left- and right-hand sides of the shock are given by ðxÞ

U L;p ¼ U L þ U 0L þ x0s U L ;

Consequently, Z x2 Z x2 U ðx; t2 Þ dx  U ðx; t1 Þ dx x1 x1 h x2  x1 ðxÞ ðxÞ ¼ ðx2  x1 Þ U L  U R þ ðU L  U R Þ 2 i

ðxÞ

U R;p ¼ U R þ U 0R þ x0s U R ; ðxÞ

F R;p ¼ F R þ F 0R þ x0s F R ;

ðtÞ

þðt2  t1 ÞU L þ HOT7 ðx  x1 ; t2  t1 Þ :

ðA:9Þ

When t2  t1 tends to zero, x1 tends to x2 and the following equivalence is obtained by neglecting the higher-order terms Z x2 Z x2 U ðx; t2 Þ dx  U ðx; t1 Þ dx  ðx2  x1 ÞðU L  U R Þ: x1

x2 x1 !0

x1

ðA:10Þ A similar reasoning leads to the following equivalence for the integrals of F and S Z t2 Z t2 F ðx1 ; tÞ dt  F ðx2 ; tÞ dt  ðt2  t1 ÞðF L  F R Þ; Z

t1 t2

Z

Sðx; tÞ dx dt t1

t2 t1 !0

t1 x2

x1

 ðt2  t1 Þðx2  x1 Þ

t2 t1 !0 x2 x1 !0

SL þ SR : 2 ðA:11Þ

Substituting Eqs. (A.10) and (A.11) into Eq. (A.6) and dividing by (t2  t1) leads to ðU L  U R Þ

x2  x1 SL þ SR : ¼ ðF L  F R Þ þ ðx2  x1 Þ t2  t1 2 ðA:12Þ

By definition, (x2  x1)/(t2  t1) is the speed cs of the discontinuity. Noticing again that x2 tends to x1 when t2 tends to t1, the term (x2  x1)(SL + SR)/2 is negligible compared to the other terms in Eq. (A.11) and the final jump relationship is obtained ðU L  U R Þcs ¼ F L  F R :

ðA:13Þ

A.3. Jump relationships for the sensitivity As mentioned in [3], the jump relationship (A.13) is not valid for the sensitivity. Note that if this was the case, the following relationship would hold ðsL  sR Þcs ¼ kL sL  kR sR :

ðA:15Þ

ðxÞ

F L;p ¼ F L þ F 0L þ x0s F L ;

ðA:14Þ

However, Eq. (A.14) is not valid, as shown hereafter. This is shown easily by considering a small perturbation / 0 in / as given by Eq. (10) in Section 2.2. At a given point x, the perturbation / 0 yields variations U 0L and U 0R in the value of

where the subscript p denotes the perturbed value and the superscript (x) denotes the derivative with respect to space of the variable of concern. The jump relationship (A.13) therefore becomes ðA:16Þ ðU L;p  U R;p Þðcs þ c0s Þ ¼ F L;p  F R;p : Substituting Eq. (A.15) into Eq. (A.16) gives the following relationship: h i ðxÞ ðxÞ U L þ U 0L þ x0s U L  ðU R;p U R þ U 0R þ x0s U R Þ ðcs þ c0s Þ ðxÞ

ðxÞ

¼ F L þ F 0L þ x0s F L  ðF R þ F 0R þ x0s F R Þ:

ðA:17Þ

Subtracting Eq. (A.13), retaining only the first-order terms, noticing that U 0 ¼ /0 s; oF v/ ; F 0 ¼ /0 ks þ ðA:18Þ o/ 0 oxs 0 / xs ¼ o/ 0 and simplifying by /0, the following relationship is obtained   oxs ðxÞ oxs ðxÞ dcs U L  sR  U R cs þ ðU L  U R Þ sL þ o/ o/ d/ oF L oxs ðxÞ oF R oxs ðxÞ v þ F  v  F : ¼ kL sL  kR sR þ o/ L o/ L o/ R o/ R ðA:19Þ This rather complex expression can be simplified in the case of the solution of the Riemann problem, where the derivatives of U and F are zero on both sides of the shock, and where the support function v is assumed to be zero everywhere. In this case, (A.19) can be simplified into dcs ¼ kL s L  kR s R : ðA:20Þ ðsL  sR Þcs þ ðU L  U R Þ d/ Even in this simple case, the jump relationship is different from the intuitive (and incorrect) formulation (A.14). A.4. Particular case: the Burgers equation The purpose of the present subsection is to relate the general jump relationship (A.19) to the result given in [3] for the inviscid Burgers equation. In [3], the influence of a perturbation in the initial condition for the Burgers equation is examined, therefore the flux F, the parameter / and the characteristic function v of X are given by

V. Guinot et al. / Advances in Water Resources 30 (2007) 1943–1961

U2 ; 2 / ¼ U ðx; 0Þ;  1 for t ¼ 0; v¼ 0 for t > 0:

F ðU Þ ¼

ðA:21Þ

In this case, the celerity k and the shock speed cs are given by k ¼ U; ðA:22Þ UL þ UR cs ¼ : 2 Substituting Eqs. (A.21) and (A.22) into Eq. (A.19) yields the following equation:   oxs ðxÞ oxs ðxÞ U L þ U R dcs sL þ U  sR  U þ ðU L  U R Þ o/ L o/ R 2 d/ oxs oxs ðxÞ ðxÞ U LU L  U RU R : ¼ U L sL  U R sR þ ðA:23Þ o/ o/ Rearranging leads to ðU L  U R Þ

dcs U L  U R UL  UR ¼ ðsL þ sR Þ þ d/ 2 2 oxs ðxÞ ðxÞ  ðU þ U R Þ: o/ L

ðA:24Þ

Simplifying by UL  UR leads to ðxÞ

ðxÞ

dcs sL þ sR U L þ U R oxs ¼ þ d/ 2 2 o/

ðA:25Þ

which is the relationship (3b) derived in [3]. References [1] Anderson ML, Bangerth W, Carey GF. Analysis of parameter sensitivity and experimental design for a class of nonlinear partial differential equations. Int J Numer Meth Fluids 2005;48:583–605. [2] Atkinson AC, Donev AN. Optimum experimental design. Oxford: Clarendon Press; 1992. [3] Bardos C, Pironneau O. A formalism for the differentiation of conservation laws. C R Acad Sci 2002;I335:839–45. [4] Beven K. Prophecy, reality and uncertainty in distributed hydrological modelling. Adv Water Resour 1993;16(1):41–51. [5] Beven K. Rainfall-runoff modelling: the primer. Chichester: Wiley; 2001. [6] Colella P, Woodward PR. The piecewise parabolic method (PPM) for gas-dynamical simulations. J Comput Phys 1984;54:174–201. [7] Courant R, Friedrichs KO. Supersonic flow and shock waves. New York: Interscience Publishers; 1948. [8] Cunge JA, Holly Jr FM, Verwey A. Practical aspects of computational river hydraulics. Pitman Publishing; 1980.

1961

[9] Fraidenraich A, Jacovkis PM, De Lima FRA. Sensitivity computations using first and second orders perturbative methods for the advection-diffusion-reaction model of pollutant transport. J Brazilian Soc Mech Sci Eng 2003;25(1):23–9. doi:10.1590/S167858782003000100004. [10] Fraidenraich A, Jacovkis PM, Lima FRA. Sensitivity analysis of shallow water problems via perturbative methods. J Brazilian Soc Mech Sci Eng 2006;28(3):286–92. [11] Gejadze IY, Copeland GJM. Adjoint sensitivity analysis for fluid flow with free surface. Int J Numer Meth Fluids 2005;47:1027–34. [12] Godunov. A difference method for the calculation of discontinuous solutions of hydrodynamics. Matematicheski Sbornik 1959;47:271–300. [13] Griewank A. Evaluating derivatives: principles and techniques of algorithmic differentiation. Philadelphia (PA): SIAM; 2000. [14] Guinot V. The discontinuous profile method for simulating twophase flow in pipes using the single component approximation. Int J Numer Meth Fluids 2001;37:341–59. [15] Guinot V. Godunov-type schemes. an introduction for engineers. Elsevier; 2003. [16] Guinot V, Gourbesville P. Calibration of physically-based models: back to basics? J Hydroinformatics 2003;5:233–44. [17] Gunzburger MD. Sensitivities, adjoints and flow optimization. Int J Numer Meth Fluids 1999;31:53–78. [18] Gunzburger MD. Perspectives in flow control and optimization. Philadelphia (PA): SIAM; 2002. [19] Gupta HV, Beven KJ, Wagener T. Model calibration and Uncertainty estimation. In: Anderson MG, Bates PD, editors. Encyclopedia of hydrological sciences. New York: Wiley; 2005 [Chapter 131]. [20] Kleiber M, Antu´nez H, Hien TD, Kowalczyk P. Parameter sensitivity in nonlinear mechanics. New York: Wiley; 1997. [21] Lax PD. Systems of conservation laws II. Comm Pure Appl Math 1957;10:537–66. [22] LeVeque RJ. Numerical methods for conservation laws. Basel: Birkhau¨ser; 1992. [23] Lighthill MJ, Whitham GB. On kinematic waves. I. Flood movement in long rivers. Proc R Soc A 1955;229:281–316. [24] Pa´zman A. Foundations of optimum experimental design. Dordrecht: Reidel Publishing; 1986. [25] Singh VH. Kinematic-wave modelling in water resources. New York: Wiley; 1997. [26] Soares-Fraza˜o S, Guinot V. An eigenvector-based linear reconstruction scheme for the shallow-water equations on two-dimensional unstructured meshes. Int J Numer Meth Fluids 2007;53:23–55. [27] Stoker JJR. Water waves. Interscience; 1957. [28] Van Leer B. Toward the ultimate conservative difference scheme. IV. A new approach to numerical convection. J Comput Phys 1977;23: 276–99. [29] Vasseur A. Well-posedness of scalar conservation laws with singular sources. Meth Appl Anal 2002;9:291–312. [30] Whitham GB. Linear and nonlinear waves. New York: Wiley; 1974. [31] Wu J, Zheng C, Chien CC, Zheng L. A comparative study of Monte Carlo simple genetic algorithm and noisy genetic algorithm for costeffective sampling network design under uncertainty. Adv Water Resour 2006;29:899–911. [32] Zhang D, Beven K, Mermoud A. A comparison of non-linear least square and GLUE for model calibration and uncertainty estimation for pesticide transport in soils. Adv Water Resour 2006;29:1924–33.