Flood hydrograph estimation using an adjoint shallow-water model

+

MODEL

Available online at www.sciencedirect.com

ScienceDirect Journal of Hydro-environment Research xx (2014) 1e12 www.elsevier.com/locate/jher

Research paper

Flood hydrograph estimation using an adjoint shallow-water model Keisuke Yoshida a,*, Tadaharu Ishikawa b b

a Graduate School of Environmental and Life Science, Okayama University, Tsushima-naka 3-1-1, Kita-ku, Okayama, 700-8530, Japan Department of Environmental Science and Technology, Tokyo Institute of Technology, 4259-G5-3, Nagatsuta-cho, Midori-ku, Yokohama, Japan

Abstract This paper describes the estimation of the flood discharge hydrograph for a compound open channel floodplain using data assimilation. The hydrograph was produced iteratively using an adjoint shallow-water model with time-series data of observed water levels. Data assimilation was applied to a flood event that occurred in mid-September 1998 along a 20 km stretch of the lower Tone River in Japan. Both the estimated hydrograph and the flow field simulated using the hydrograph were verified by comparison with the results of aerial photographic analysis. Results show that the hydrograph based on the stageedischarge was overestimated, but it was reasonably corrected by the data assimilation technique used for this study. In addition, the resulting relation during the flood event exhibited a hysteresis loop characteristic, which is typical of compound open channel flows during flooding. © 2014 International Association for Hydro-environment Engineering and Research, Asia Pacific Division. Published by Elsevier B.V. All rights reserved.

Keywords: Data assimilation; Discharge estimation; Adjoint model; Shallow-water equation; Compound open channel

1. Introduction A discharge hydrograph includes crucially important data that are necessary for river engineering planning and management tasks. A hydrograph is usually obtained using the Manning formula or some sort of empirical formulas based on a single-valued stageedischarge (HeQ) relation or a rating curve. Generally, the relation is developed for a specific river location. It is uniquely determined for the entire flow range including extreme flows such as flooding. For this relation, measurements of river discharges are routinely conducted for several stages (water levels): for normal flows, direct measurements are conducted using a velocity meter such as a Price meter, but for high water flows they are typically estimated from the traveling time of a floating rod dropped from bridges between two streamwise locations. It is true that this workable engineering practice for river discharge determination enables monitoring of the current state of flow condition, and enables

* Corresponding author. Tel.: þ81 86 251 8149. E-mail address: [email protected] (K. Yoshida).

the prediction of subsequent hydraulic phenomena. However, it is necessary that the HeQ relation be valid only for quasisteady flows in a stable cross-section of the river channel. Previous studies have demonstrated that the relation might change for unsteady flows because of flood wave propagation (Chow, 1959; Fread, 1975), or in compound open channels when the flood plain is inundated (Fukuoka et al., 2005). Consequently, the discharge estimated using such a simple relation might include some errors for high-water events. For instance, the practical method based on the Manning formula is reportedly difficult in the accurate estimation of river discharge because of inappropriate modelling of energy slope in the compound channel during flooding (Bousmar and Zech, 1999). Replaced by the conventional method using the HeQ relation of limited accuracy, a new instrument using a nonintrusive sensor has been developed by several researchers for river discharge measurements. Gordon (1989) presented a new moving-boat method for the rapid measurement of a river discharge using an acoustic Doppler current profiler (ADCP). Although the entire system was not ready for easy deployment in practical measurements of river flows at early times of the history of ADCPs, improvements including weight reduction of

http://dx.doi.org/10.1016/j.jher.2014.12.003 1570-6443/© 2014 International Association for Hydro-environment Engineering and Research, Asia Pacific Division. Published by Elsevier B.V. All rights reserved.

Please cite this article in press as: Yoshida, K., Ishikawa, T., Flood hydrograph estimation using an adjoint shallow-water model, Journal of Hydro-environment Research (2014), http://dx.doi.org/10.1016/j.jher.2014.12.003

+ 2

MODEL

K. Yoshida, T. Ishikawa / Journal of Hydro-environment Research xx (2014) 1e12

instruments, simple operation, and upgraded performance have recently enabled researchers to conduct continuous measurements more accurately during flooding (Nihei and Kimizu, 2008). However, measurements using updated ADCPs still present several challenges for practical use: (1) signal decay because of high turbidity during flooding; (2) difficulty in measurements for floods with piles of drift wood and a great deal of litter in water; and (3) wash-out of equipment for measurements during severe high water events, and so on. Fujita et al. (1998), based on the idea of PIV for laboratory flumes, developed large-scale particle image velocimetry (LSPIV) for the instantaneous measurement of surface velocity fields in rivers using a protocol for image pattern recognition (Adrian, 1991). This image velocimetry technique has been improved for practical use by several researchers including Bradley et al. (2002), Fujita and Hino (2003), and Muste et al. (2004). They used the traceless LSPIV to obtain valuable data clarifying scientific hypotheses related to river processes including flood wave propagation. Nevertheless, this technique still exhibits difficulties in continuous measurements of river discharge for some cases: (1) during normal water conditions because the free surface has a mirror-like appearance; (2) when lighting is limited, for example, at night; and (3) when strong winds blow over a targeted river reach. Recently, flooding in small and medium-sized rivers has occurred frequently in Japan because of local intensive rainfall, which is reportedly induced by abnormal weather. As a result, tremendous damage from such flooding has increasingly occurred in many regions throughout Japan. It is acknowledged that updated observation methods offer valuable data for river engineering tasks. However, direct measurements of river flows using ADCPs and LSPIV during high water events threaten the safety of engaged personnel and engender the loss of expensive equipment and related missing data in severe natural conditions. Therefore, diversity in observation methods must be accepted for the development of a stable low-cost observation system for practical work. Flexible choices must be made among several observation methods considering strong and weak points of the methods, including total cost, robustness, user-friendliness, and accuracy might contribute to river engineering work not only in Japan but also in locations world-wide. This paper therefore presents a new methodology of using data assimilation for improvement of the conventional HeQ relation. An inverse estimation of the flood discharge hydrograph was developed for compound open channels in actual rivers using an adjoint shallow-water (backward) model. The method used here is based on a variational approach, which is a type of data assimilation technique. Previous studies (e.g., Atanov et al., 1998; Ding and Wang, 2006; Sanders and Katopodes, 1999) have used identical twin experiments to demonstrate the applicability of such an approach to hydraulic problems. The identical twin experiment is a numerical procedure by which synthetic data can be generated by the forward model to which data assimilation is applied. This numerical experiment assures us of feasibility and creditability. These studies emphasized one-dimensional solutions, which cannot be applied easily to horizontal hydraulic phenomena in natural

rivers where the cross-section varies irregularly in the longitudinal direction or where the floodplains exhibit a range of bed roughness values. Additionally, it is important to use real data for data assimilation and to check the feasibility and applicability in river engineering, although such a numerical experiment is effective for designing an observation system in rivers. This paper is organized as follows. The governing equations for both two-dimensional (2-D) shallow water-model and corresponding adjoint model are described using numerical procedures to produce a river discharge hydrograph. This study applied data assimilation methods to a flood event that occurred in midSeptember 1998 along a 20 km stretch of the lower Tone River in Japan. The assimilated data comprised a time series of observed water levels. In Japan, observation data are usually difficult to obtain during flooding, except for data of water levels at established hydraulic stations. Measurements of river discharge using ADCPs and LSPIV have been conducted for experimental purposes. Practical work conducted using updated instruments is not prevalent in Japan today because of observation costs and a scarcity of sufficiently skilled engineers. Alternatively, for a flood event, aerial photographs were taken over the targeted river reach during around half an hour after the peak period of the flood wave. Therefore, the data assimilation results were verified by the surface flow field obtained from analyses of a series of aerial photographs. Consequently, discussion is conducted mainly for river discharge data and simulated data of velocity profiles of the main flow across the river during that period. Additionally, the potential of this work is presented in light of upgrading the HeQ relation considering 2-D shallow-water hydrodynamics. 2. Governing equations We used the shallow-water equations as a forward model in a boundary-fitted coordinate system as follows (Nagata et al., 2000).


 v h v M v N fh ¼ þ þ ¼0 vt J vx J vh J

ð1Þ




 v M v M2 v MN fM ¼ þ þ vt J vx Jh vh Jh

 m vxx M vxx N n vxy M vxy N þ þ   J vx h vh h J vx h vh h ! x2x þx2y vH xx hx þxy hy vH þ þgh vh J vx J 1 0 2 x2x vth xx hx vth xy vjh þ þ C B C J vh J vx tx B J vx C¼0 þ B B rJ @ x h vjh 2x x vfh x h þx h vfh C A x y y x y x y þ þ þ x vh J vh J vx J ð2Þ

Please cite this article in press as: Yoshida, K., Ishikawa, T., Flood hydrograph estimation using an adjoint shallow-water model, Journal of Hydro-environment Research (2014), http://dx.doi.org/10.1016/j.jher.2014.12.003

+

MODEL

K. Yoshida, T. Ishikawa / Journal of Hydro-environment Research xx (2014) 1e12




 v N v MN v N2 fN ¼ þ þ vt J vx Jh vh Jh

 m vhx M vhx N n vhy M vhy N þ þ   J vx h vh h J vx h vh h
 xx hx þxy hy vH h2x þh2y vH þgh þ vx J J vh 0 1 2 xx hx vth hx vth xy hy vjh B J vx þ J vh þ J vx C C th B C¼0 þ B B rJ @ h2 vjh x h þx h vfh 2h h vfh C A x y y x y x y þ þ þ vx J vh J J vh

subcritical flows can be achieved by minimizing the following objective function, which is defined as the integration of the square errors between simulated and observed data (Atanov et al., 1998). 1 F¼ 2

ZT ZLx ZLh X Nk 0

0

ð3Þ fm ¼ m  xx M  xh N ¼ 0; fn ¼ n  yx M  yh N ¼ 0
 v v m 2 þ k¼0 ft ¼ t  2ε xx þ hx vx vh h 3
 v v n 2 fj ¼ j  2ε xy þ hy þ k¼0 vx vh h 3

  v v m v v n  ε xx þ hx ¼0 f f ¼ f  ε xy þ h y vx vh h vx vh h

3

ð4Þ ð5Þ ð6Þ

 2 Ki di Hi  Hi dxdhdt

ð9Þ

i¼1

0

In that equation, T represents the computation period, Lx and Lh respectively denote the longitudinal and cross-sectional spatial lengths in the computational domain, i is an index that distinguishes the observation station, Nk stands for the number of stations for which the data are assimilated, Hi is the observed water level at the ith station, di is Dirac's delta, andKi is a weighting factor of the data measured at the station. In this study, Ki ¼ 1.0, independent of the station. The objective function varies with the discharge hydrograph. Therefore, minimization of the function results in optimal estimation of the discharge hydrograph. 3.2. Variational approach

ð7Þ

Therein, functions fh, fM, and fN are differential operators for the mass balance and momentum conservations. The other f symbols are operators used to define the corresponding dependent variables m, n, t, j, and f, where t is time, x and y are Cartesian coordinates, x and h are boundary-fitted coordinates, h is the local water depth, and H is the water level. In addition, m (¼uh) and n (¼vh) are x- and y-components of the discharge flux vector, u and v are x- and y-components of depth-averaged velocity vector, M (¼xxm þ xyn) and N (¼hxm þ hyn) are contravariant components of the discharge flux vector, xx, hx, xy, and hy are metrics (e.g., xx ¼ vx/vx), tx and th are contravariant components of the bed shear stress vector, g is the acceleration of gravity, n* is Manning's roughness coefficient, t, j and f are depth-averaged Reynolds stress tensors, J is the Jacobian defined as J≡ðxx yh  xh yx Þ1 , ε is the horizontal eddy viscosity evaluated by Richardson's power law, and k is the depth-averaged turbulent kinetic energy evaluated using a quasi-empirical formula (Nezu and Nakagawa, 1993). The physical meaning of the contravariant components of discharge flux vector is the discharge per cell spacing, which is defined perpendicular to the cell boundary. In this study, the Reynolds stress terms are calculated using a linear zero-equation model. The bed roughness is evaluated using Manning's formula as presented below. pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi gn2 M m2 þ n2 gn2 N m2 þ n2 tx ¼ * 7=3 ; th ¼ * 7=3 ð8Þ h J h J 3. Adjoint analysis 3.1. Objective function The framework of the inverse problem of estimating the river discharge or upstream boundary condition in the

The objective function is reconstructed into a functional (Sasaki, 1970) to minimize the objective function and to ensure that the water level satisfies the shallow-water model (1)e(7). L¼Fþ

ZT ZLx ZLh X 0

0

ð10Þ

lq fq dxdhdt

q

0

Therein, L stands for an augmented objective function, q is a dependent flow variable, and lq represents Lagrangian multipliers with respect to the corresponding operators. According to the variational principle, the first variation in L must be zero to minimize the function with the constraints, which are expressed as the forward model. After integrating Equation (10) by parts along the longitudinal direction x, the following equation is derived by taking the variation of the functional with respect to all flow variables as shown below. dL ¼

ZT ZLx ZLh X 0

0

0

dqfq* dxdhdt

ZT ZLh 

q


  1 2M N lh þ lM þ lN  u¼ J Jh Jh x¼0

udM0 dhdt 0

ð11Þ

0

ð12Þ

In those equations, dq represents variations of the flow variables, fq* represents differential operators with respect to the Lagrangian multipliers, and M0 is the value of M at the upstream end (x ¼ 0). In the process of the integration by parts, the terms multiplied by dM0 remain while the others vanish because variable M0 is related to the river discharge at the upstream end, as described later. According to the extremum conditions, all terms multiplied by the variations of the flow variables in Equation (11) can be

Please cite this article in press as: Yoshida, K., Ishikawa, T., Flood hydrograph estimation using an adjoint shallow-water model, Journal of Hydro-environment Research (2014), http://dx.doi.org/10.1016/j.jher.2014.12.003

+ 4

MODEL

K. Yoshida, T. Ishikawa / Journal of Hydro-environment Research xx (2014) 1e12

set to zero to establish the optimal conditions, except for the term related to variableM0. Consequently, the adjoint equations (Yoshida, 2010) and the sensitivity of the objective function with respect to variable M0 are derived as follows.

piecewise profile. We then calculated a discrete value of the discharge at regular time intervals of 1 h. Therefore, a vector consisting of an array of sensitivity coefficients is defined as follows.

adjoint equations : fq* ¼ 0

g¼

ð13Þ

dL dM0j

ð15Þ

dL ¼ u ð14Þ dM0 The complete form of the adjoint equations is shown in Appendix I.

Therein, j is a time index and dL=dM0j is the sensitivity coefficient at discrete time j. River discharge Q at the upstream end is computed using the following equation.

4. Numerical procedure

Q¼

4.1. Discretization scheme

The sensitivity of the discharge at the discrete time is readily calculated from that of the discharge flux at the upstream end according to Equation (16). The nonlinearity of the problem requires an iterative procedure to minimize the objective function and to estimate the optimal discharge. Fig. 1 portrays a flowchart for finding the optimal solution. We used a quasi-Newton method to minimize the function. Specifically, we used a typical algorithm using a BroydeneFletchereGoldfarbeShanno (BFGS) formula for solving the present unconstrained nonlinear optimization problem. The Newton direction was specified with an inverse Hessian matrix and the steepest descent direction, which can be computed efficiently using Equation (15). The initial estimate of the hydrograph was the observed one. In addition, a convergence condition was imposed to make the relative error between the previous and the updated discharge in the iteration less than a certain tolerance b as follows.   X Qjlþ1  Qjl  b ð17Þ Qjl j

sensitivity :

ZLh

M0 dh J

ð16Þ

0

Equations (1)e(7) were discretized using the finite difference method. Time development of the flow variables was computed using the explicit Euler scheme. The spatial derivatives were discretized using the second-order central scheme except for the advection terms, which were discretized by the first-order upwind scheme. The water edge varies over time during a flood. Therefore, a threshold water depth of 0.35 m was set to determine drying and wetting. The computation nodes in the dry bed were assigned a velocity of zero. In addition, the slip-velocity condition was applied at the water's edge. The explicit Euler scheme might be a low-order one to solve unsteady flow systems. However, this scheme and several explicit schemes were used in earlier technical reports related to the numerical investigation of past flood flows that occurred in the lower reach of the Tone River. In addition, data assimilation has been applied only rarely in the field of river engineering. Therefore, we used the Euler explicit scheme because we specifically examine a fundamental element for the application of data assimilation in river management. Improvement of the numerical schemes in both data assimilation and flood flow simulations will be examined in our future studies. In this study, the backward model was derived directly from the linearization of the discretized forward model, which indicates that the discretized forms of the present adjoint equations were derived from transposition of a tangential form of the discretized shallow-water equations. Even though the discretized form of the continuous adjoint equation is not always identical to the discrete adjoint equation because of the numerical schemes that were used, previous studies (e.g., Lardner et al., 1993) have shown that the exact sensitivity can be computed using the discrete adjoint equations. Moreover, the discretized boundary condition of the adjoint equations can be derived using the same method as that of the adjoint equations. The derivation process for a simplified shallowwater model is presented in Appendix II. 4.2. Estimation of river discharge hydrograph To facilitate the task of estimating the hydrograph numerically, we assumed that the discharge hydrograph had a linear

In that expression, l is the iteration number; j stands for the time index defined in Equation (15). The river discharge per unit width, M0, at the first iteration was estimated assuming a uniform flow at the upstream end using both the observed data of river discharge Q and the riverbed bathymetry. 5. Application 5.1. Site description Fig. 2 shows a plan view of the land use in the target domain, the 20-km-long lower reach of the Tone River. The kilo post (KP) numbers in this figure show the longitudinal distance from the river mouth. Fig. 2 is based both on aerial photographs taken during normal water conditions and on a field survey of bed materials and vegetation. The numbers in parentheses in this figure represent values of the distributed Manning's roughness coefficients used in the numerical simulations. These values were found using tables and figures of earlier papers and technical reports on the bed roughness in the lower part of the Tone River (Fukuoka et al., 2005, 2006). The sand particles were not distributed densely along the main

Please cite this article in press as: Yoshida, K., Ishikawa, T., Flood hydrograph estimation using an adjoint shallow-water model, Journal of Hydro-environment Research (2014), http://dx.doi.org/10.1016/j.jher.2014.12.003

+

MODEL

K. Yoshida, T. Ishikawa / Journal of Hydro-environment Research xx (2014) 1e12

5

Fig. 3. Plan view and bathymetry contour map of the lower Tone River.

Fig. 1. Flowchart for finding the optimal discharge hydrograph using the adjoint method described as a result of this study.

channel in this targeted domain. Therefore, we specified the uniform roughness coefficient for the numerical simulation. The effect of irregular cross-sections within the domain on the flow resistance was considered implicitly using the momentum transfer attributable to the convection terms and the Reynolds stress terms in the Equations. (2) and (3). Fig. 3 shows the plan view and bathymetry contour map of the target domain, the 20-km-long lower reaches of the Tone River. In this figure, Yedogawa Peil (Y.P.) is the datum level for level data in the Edogawa River in Japan. The channel bathymetry was obtained by direct survey every 500 m along the river. The bed profile was estimated over the complete target domain using linear interpolation based on these survey data. Three hydraulic stations (Fukawa, Suga, and Shinkawa) operate in this area. They record water levels only every hour with accuracy of less than an order of magnitude of 1 cm. The river discharge hydrograph is estimated using the stageedischarge relation at Fukawa Station. The data margin for river discharge using the HeQ relation is reportedly estimated as ± 20% in Japan, compared with data of ADCPs. The water level data measured at these stations were used for the data assimilation. 5.2. Observation data Fig. 4 shows the hourly recorded water levels at three hydraulic stations during the flood in September, 1998. Heavy

rainfall during Typhoon No. 5, which passed over the south Kanto region in mid-September, caused a sharp rise in the water levels. Table 1 shows the river discharge measured by floating rods at several sections in the lower Tone River at 8:00 a.m. on September 17. The discharge measured at Fukawa Station differed considerably from data obtained at other stations. It might have been overestimated for some reason. Japanese river discharge is usually estimated using a crosssectional integral of the flow velocity, which is calculated using the travel time of floating rods between two longitudinal sections. The rods are dropped manually from bridges. Therefore, their behavior depends on the local flow around the bridge piers. Generally, upwelling flows downstream of the bridge piers generate parallel helical flows (Kinoshita, 1984) between the piers, producing a series of high-speed and lowspeed flow bands that result in convergence in the highspeed region and divergence in the low-speed region. Measurements based on floating rods tend to approach the values for the high-speed region, leading to overestimation of the discharge. The narrowness of the river near Fukawa Station might also contribute to the lower accuracy of observations made there. The availability of observation data was limited for several reasons during the flood. Sequential measurements of accurate hydraulic data in many areas along the river channel were difficult, except for some measurements of water levels. Consequently, few observational data were available to verify the inverse estimation of the river discharge. Fortunately, a series of aerial photographs with an overlap of 80% was taken in the area on a sunny day during the flood, as portrayed in Fig. 5(a). Originally, aerial photographs were taken under the instructions of a local governmental office for river management to gain a quick overview of disaster sites in rivers. Fig. 5(a) presents an example of aerial photographs taken of the downstream sections in this river. Flow tracers such as foam and suspended loads on the free surface appear as clear shading patterns in the photographs. For this study, flow field data were extracted from these photographs based on a hybrid PIV algorithm for the digitalized aerial photographs (Minoura et al., 2008). Their method fundamentally determines a surface velocity field based on the

Fig. 2. Plan view of land use and corresponding Manning's roughness coefficient along a 20-km stretch of the lower Tone River. Please cite this article in press as: Yoshida, K., Ishikawa, T., Flood hydrograph estimation using an adjoint shallow-water model, Journal of Hydro-environment Research (2014), http://dx.doi.org/10.1016/j.jher.2014.12.003

+ 6

MODEL

K. Yoshida, T. Ishikawa / Journal of Hydro-environment Research xx (2014) 1e12

Fig. 4. Hourly observation data of water level at hydraulic stations.

PIV method by tracking movement of the tracers between two times when two photographs were taken in sequence. Furthermore, recognizable errors in the results of surface velocity vectors are manually revised appropriately with support of stereo image analysis based on Cameron's effect (Cameron, 1952). Consequently, the surface velocity field was determined as presented in Fig. 5(b). Additionally, the river discharge was estimated from the cross-sectional profile of depth-averaged velocity in conjunction with a cross-sectional bed profile and water level data at the time the photographs were taken. For the estimation process, we assumed that the depth-averaged velocity was related to the surface velocity by a velocity index, which was set to 0.85 for this study. These datasets were used to verify the estimated discharge. It has been commonly acknowledged to date that accurate observation data of flow fields has been rarely obtained during flooding with large discharge, considering the difficulty of river discharge observation in extreme floods, as described both in this section and in the Introduction. Therefore, river discharge data estimated based on the flow field data using accurate aerial photographic analysis are valuable as validation data during high water levels, although they might be limited in terms of the time they are obtained. For this study, comparison of data assimilation results with estimated data was conducted. The comparison is meaningful for modern river engineering. 5.3. Computational conditions Numerical simulations of the forward and backward models were conducted as follows. The time increment Dt was 0.5 s. The computational mesh consisted of Table 1 River discharge measured by floating rods on the lower Tone River during the period when photographs were taken. Station

River discharge (m3/s)

Toride (85 KP) Fukawa (76.5 KP) Ryuudai (61 KP) Kanaetsu (54 KP) Kanzaki (49 KP) Sawara (40 KP)

6641 7678 6438 6334 6730 6314

(Nx,Nh) ¼ 504  24 cells, with 504 cross-sections in the longitudinal direction and 24 nodes in each cross-section. The time increment was adjusted to allow for numerical stability so that the Courant number was less than 0.1. The mesh size was 20e50 m. The computation period was set at 48 h, from 5:00 p.m. on September 16 to 5:00 p.m. on September 18. The average CPU time for each computation was about 48 h for each case. The water level measured in the downstream section at Shinkawa Station (Fig. 2) and the river discharge estimated in the upstream section at 77 KP near Fukawa Station were used as boundary conditions of the numerical simulations of the forward model. The initial conditions of the forward model were prepared assuming a non-uniform flow, whereas those of the backward model were assigned a transversality condition. The adjoint equations were calculated backward in time. Therefore, all Lagrange multipliers were zero at the final time of 5:00 p.m. on September 18. The assimilated data that were given were prepared by interpolating the water levels measured hourly at the observation stations. Table 2 presents the computational conditions of the observation stations where the water level data were assimilated for the backward model throughout the computation period. Data of only one kind were used for Cases 2 and 3, whereas data of two kinds were assimilated simultaneously for Case 1. However, data assimilation (backward simulation) was not conducted for Case 4, for which the observed discharge was simply used for the simulation of the forward model with no modifications. 5.4. Results and discussion Fig. 6 presents the iteration histories of the objective function values, where F0 denotes the initial value of the function in the iteration. After 10 iterations in each case, the resulting discharge was regarded as the optimal solution because the function was minimized. According to previous studies (e.g., Ding et al., 2004), both the iteration number and convergence behavior depend on the numerical scheme adopted for optimization. The value of F0 in Case 1 is almost 3% of the initial objective function. That rate in Case 1 is the greatest of all cases. The results of data assimilation invariably contain various errors: measurement error, model error, numerical error, and so on (Khatibi et al., 1997). Additionally, real data of the initial flow field cannot be provided for the forward model such as a spatial variation of flow velocities and water levels. Therefore, the converged value of the function was greater in Case 1 because the model outputs were not simultaneously consistent with the data at two stations, and data assimilation gives the solution with the least squared error. Even in Cases 2 and 3 with just one kind of data assimilated, the model output did not match the observed data completely because the initial condition in the forward model was not identical to the actual condition. Fig. 7 shows the estimated and observed discharges, with the results of aerial photographic analysis added for comparison. The observed discharge denotes the river discharge, as

Please cite this article in press as: Yoshida, K., Ishikawa, T., Flood hydrograph estimation using an adjoint shallow-water model, Journal of Hydro-environment Research (2014), http://dx.doi.org/10.1016/j.jher.2014.12.003

+

MODEL

K. Yoshida, T. Ishikawa / Journal of Hydro-environment Research xx (2014) 1e12

7

Fig. 5. Example of aerial photographic analysis in the downstream part of the river showing (a) aerial photographs taken at 8:06 a.m. on September 17 and (b) the corresponding surface velocity field.

estimated empirically using a single-valued stageedischarge relation with observed data of water levels measured at Fukawa Station. This relation was based on measurement data of river discharges at the upstream end and water levels at Fukawa Station. Cases 1, 2, and 3 refer to the inverse estimation of the river discharge using observation data of water levels at the station(s) described in Table 2, whereas Case 4 presents the forward (normal) simulation of the shallow water model. Because the aerial photographs were taken at different times while an airplane flew over the river reach, the data shown for discharge were reasonably derived from simulated results, considering the time adjustment. As Fig. 7 shows, the observed discharge was 20%e30% greater than the estimates. The result presented in Case 1 was consistent with the results of the aerial photographic analysis, indicating that more data improved the estimation if the assimilated data were accurate. Moreover, Fig. 7(b) shows that the estimated results based on data assimilation show good agreement with aerial photographic analysis data, except for Case 2, although data obtained using the aerial photographic analysis were limited. For Case 2, where the data of water levels at Suga Station were assimilated, the channel has a mildly curved section near Suga Station. Additionally, detailed data of the bed level and bed roughness within the golf course around 66 KP were not fully incorporated into the numerical simulation. Therefore, simulated data of water levels using the two-dimensional (2-D) shallow water equation do not fully reproduce the measured data of water levels, which indicates both the importance of data quality for data assimilation and limitations of the 2-D shallow-water model to the flood flow in this actual river reach. Fig. 8 shows the relation between the water level and the river discharge at Fukawa Station during the computation

period. It is natural that the observed discharge is correlated uniquely with the water level data because the discharge is calculated from the single-valued stageedischarge relation. The estimated discharges agreed closely with the results of aerial photographic analysis at the time the photographs were taken. They also show the hysteresis loop characteristic. This type of stageedischarge relation is quite characteristic of flows in compound open channels during a flood. Therefore, the method used here appears to yield reliable results. For practical works of river discharge measurements, the rating curve using a single-valued stageedischarge relation is generally prevalent because of its easy implementation and low cost at present. However, using the data assimilation technique facilitates accurate estimation of river discharge hydrographs. Moreover, the results are useful as fundamental data for river management tasks because they include 2-D shallow-water hydrodynamics including the hysteresis loop characteristics observed during flooding. Fig. 9 presents examples of numerical results for the crosssectional profiles of the main flow simulated with the estimated discharge every 3 km longitudinal distance, compared with the results of aerial photographic analysis. In Case 4, a numerical simulation was conducted using the observed

Table 2 Computational conditions for data assimilation showing station(s) where the water level data were assimilated. Case

Station

1 2 3 4

Fukawa, Suga Suga Fukawa None

Fig. 6. Iteration histories of objective function F where F0 is the initial value of F.

Please cite this article in press as: Yoshida, K., Ishikawa, T., Flood hydrograph estimation using an adjoint shallow-water model, Journal of Hydro-environment Research (2014), http://dx.doi.org/10.1016/j.jher.2014.12.003

+ 8

MODEL

K. Yoshida, T. Ishikawa / Journal of Hydro-environment Research xx (2014) 1e12

Fig. 7. Simulated and observed data of river discharge Q (a) for hydrographs at the upstream end, and (b) at three hydraulic stations, to which the results of aerial photographic analysis have been added.

discharge based on the stageedischarge relation as the upstream boundary condition. Data from aerial photograph analysis represent instantaneous surface velocities, whereas the numerical results show the Reynolds-averaged and depthaveraged velocities. Unsteady horizontal eddies might be observed on the river surface along the main channel edges onsite, according to their appearance in aerial photographs. However, coarse mesh was used in numerical simulations for this study because data assimilation is costly in terms of both computer memory and simulation time. Therefore, these eddies were not reproduced for the present numerical simulations. Alternatively, the eddy viscosity based on Richardson's power law was included for the Reynolds stress terms in the present numerical simulations, considering the effects of the momentum transfer because of the horizontal eddies. These figures show that the magnitude of the depth-averaged velocity in Case 4 was consistent with the surface velocity estimated from aerial photograph analysis, although the surface velocity was affected strongly by organized large-scale motions induced in the shear layer between the main channel and the

Fig. 8. Relation between water level H and river discharge Q at the upstream end.

floodplain, as shown in Fig. 5. This result implies that the observed discharge might be overestimated because the depthaveraged velocity predicted with the observed discharge are expected to be less than the surface velocity estimated from aerial photographic analysis. 6. Conclusions This study used a data assimilation technique that is applicable to actual flood events to examine the improvement of prevalent stageedischarge relations for river discharge hydrograph estimation in practical works. Specifically the study demonstrated inverse estimation of the discharge hydrograph in the lower part of the Tone River during a flood that occurred in mid-September 1998. The adjoint shallowwater model developed by Yoshida (2010) was used for inverse estimation of the discharge hydrograph. The assimilated data consisted of a time series of the water levels, which were measured easily at three established hydraulic stations along the 20 km compound open channel during flooding. The estimated discharge was verified by comparison with data from the aerial photographic analysis. Results showed that the observed discharge was overestimated but reasonably corrected using the method described here, although insufficient verification data are available to estimate the discharge hydrograph exactly. In addition, the estimated hydrograph of the water levels and discharges showed a hysteresis loop, which is typical of the hydrodynamics in compound open channels during flood events. For simplicity, because large bed deformation did not occur during the flooding, this study did not consider sand deposits or erosion. Bed deformation can cause time variation in the cross-sectional flow area and can therefore strongly affect the water level. Additionally, few validation data for comparison were available for this study, although the methodology used for this study is universally applicable. Despite these limitations, the present work can provide a low-cost, reasonable tool to estimate river discharge hydrographs considering 2-D shallow-water hydrodynamics in actual rivers, and also to contribute to reasonable solutions for river management tasks.

Please cite this article in press as: Yoshida, K., Ishikawa, T., Flood hydrograph estimation using an adjoint shallow-water model, Journal of Hydro-environment Research (2014), http://dx.doi.org/10.1016/j.jher.2014.12.003

+

MODEL

K. Yoshida, T. Ishikawa / Journal of Hydro-environment Research xx (2014) 1e12

9

Fig. 9. Cross-sectional profile of main flow velocity U from numerical simulation and aerial photographic analysis (KP numbers denote the longitudinal distance in kilometers from the river mouth).

Future work will include the collection of sufficient verification data during a flood as well as the design of numerical simulations that include bed deformation and related phenomena. Acknowledgments This work was supported in part by JSPS KAKENHI Grant Number 26820203. The first author thanks Dr. Minoura for his assistance to address the data of the surface flow field from aerial photographic analysis. Valuable comments by anonymous reviewers are also gratefully acknowledged for supporting the improvement of this work.

Notation The following symbols were used herein. fq operators for the shallow-water model fq* operators for the adjoint shallow-water model F objective function F0 initial value of objective function g acceleration of gravity g vector of sensitivity coefficients defined in Equation (15) h local water depth H water level H i observed water level at the ith station i index identifying an observation station j time index in Equation (15) J Jacobian defined as J≡ðxx yh  xh yx Þ1 k depth-averaged turbulent kinematic energy Ki weighting factor at the ith station

L Lx Lh l M

functional longitudinal spatial length in the x direction longitudinal spatial length in the h direction iteration number in quasi-Newton method contravariant component of discharge flux vector in the x direction M0 value of M at the upstream end m x- component of discharge flux vector (¼uh) N contravariant component of discharge flux vector in the h direction Nk number of observation stations at which data were assimilated Nx,Nh number of numerical grid cells in x and h directions: n y-component of discharge flux vector (¼vh) n* Manning's roughness coefficient Q river discharge t time Dt time increment T computation period u x-component of depth-averaged velocity vector U main flow velocity v y-component of depth-averaged velocity vector x,y Cartesian coordinates ε horizontal eddy viscosity di Dirac's delta function at the ith station dq variation of the flow variable q lq Lagrangian multiplier with respect to fq x,h boundary-fitted coordinates xx ; hx ; xy ; hy metrics (e.g.,xx ¼ vx/vx) q arbitrary flow variable t; j; f depth-averaged Reynolds stress tensors tx contravariant component of bed shear stress vector in the x direction

Please cite this article in press as: Yoshida, K., Ishikawa, T., Flood hydrograph estimation using an adjoint shallow-water model, Journal of Hydro-environment Research (2014), http://dx.doi.org/10.1016/j.jher.2014.12.003

+ 10

MODEL

K. Yoshida, T. Ishikawa / Journal of Hydro-environment Research xx (2014) 1e12

th contravariant components of bed shear stress vector in the h direction u coefficient defined in Equation (12)



  v lM 1 vlh 2M vlM N vlM N vlN  þ   vt J J vx hJ vx hJ vh hJ vx  nma na3 o mb1 nb3 1  lM þ þ  lN hJ hJ hJ hJ p ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi gn2 m2 þ n2 lM  xx lm  yx ln þ * 7=3 h J ¼0

f*M ¼ 

Appendix I. Complete form of adjoint shallow-water equations fq* Based on the deviation process described in a body text of this paper, the adjoint shallow-water equations are given in detail as presented below. Equations I.1 to I.8 represent relations among eight Lagrange multipliers (lh,lM,lN,lm,ln,lt,lj, and lf). Equation (I.1) only has the source term signifying discrepancies between observed and simulated water levels. These multipliers are computed backward in time. Using values of these multipliers obtained at the end of these computations, we calculate the sensitivity of the objective function, as shown in Equation (14). Based on any algorithm for nonlinear optimization problems including the BFGS formula for this study, the discharge values of M at the upstream end are iteratively updated such that the objective function can be minimized further. Consequently, the best estimate of river discharge hydrograph is acquired based on Equation (16) after the iteration process is terminated satisfying a prescribed tolerance.

2  v lh M vlM MN vlM MN vlN N 2 vlN f*h ¼  þ 2 þ 2 þ 2 þ 2 vt J h J vh h J vx h J vh h J vx

  m a1 M a2 N n a3 M a4 N þ þ þ þ lM J h2 h2 J h2 h2

  m b1 M b2 N n b3 M b4 N þ 2 þ þ 2 þ lN J h2 h J h2 h
 vghG1 lM vghG2 lM vghU1 lN vghU2 lN þ þ þ  vx vh vx vh  

vH vH vH vH þ glN U1 þ U2 þ glM G1 þ G2 vx vh vx vh

 vG3 lM vU3 lN vG4 lM vU4 lN vU7 lN þ þ þ þ txx þ vx vx vh vh vx
 vG7 lM vG8 lM vU6 lN vG5 lM vU5 lN þ þ þ þ txy þ vx vh vh vx vx pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  2 2 2 vU4 lN vG6 lM 7 gn* m þ n þ ðMlM þ NlN Þ þ tyy  3 vh vh h10=3 J
  14 ggn2* 2 2εm vxx lt vhx lt 2  þ þ n þ l Þ  m ðl t j 9 h10=3 h2 vx vh

  vh l vh l 2εn vxy lj εm vxy lf y j y f þ þ  2  2 h h vx vh vx vh


 εn vxx lf vhx lf þ  2 þ Ki di Hi  H i h vx vh ¼0

ðI:2Þ f*N



  v lN 1 vlh M vlM M vlN 2N vlN   þ ¼  vt J J vh hJ vh hJ vx hJ vh  nma o na4 mb2 nb4 2 þ þ  lM  lN hJ hJ hJ hJ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 gn m þ n lN  xh lm  yh ln þ * 7=3 h J ¼0 ðI:3Þ  

lM a1 M a2 N lN b1 M b2 N þ þ  h h h h J J
 2 gn* m 2ε vxx lt vhx lt pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðMlM þ NlN Þ þ þ þ h vx vh h7=3 J m2 þ n2
 2 ε vxy lf vhy lf 4 ggn* þ þ mðlt þ lj Þ þ h vx 3 h7=3 vh

f*m ¼ lm 

¼0 ðI:4Þ f*n

 

lM a3 M a4 N lN b3 M b4 N þ þ ¼ ln   h h h h J J
 2 gn* n 2ε vxy lj vhy lj pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðMlM þ NlN Þ þ þ þ h vx vh h7=3 J m2 þ n2
 2 ε vxx lf vhx lf 4 ggn* þ þ nðlt þ lj Þ þ h vx 3 h7=3 vh ¼0 ðI:5Þ

f*t ¼ lt þ h



 vG3 lM vG4 lM vU3 lN vU4 lN þ þ þh ¼0 vx vh vx vh ðI:6Þ

f*j ¼ lj þ h

ðI:7Þ f*f ¼ lf þ h

ðI:1Þ



 vG5 lM vG6 lM vU5 lN vU6 lN þ þ þh ¼0 vx vh vx vh



 vG7 lM vG8 lM vU7 lN vU8 lN þ þ þh ¼0 vx vh vx vh ðI:8Þ

Please cite this article in press as: Yoshida, K., Ishikawa, T., Flood hydrograph estimation using an adjoint shallow-water model, Journal of Hydro-environment Research (2014), http://dx.doi.org/10.1016/j.jher.2014.12.003

+

MODEL

K. Yoshida, T. Ishikawa / Journal of Hydro-environment Research xx (2014) 1e12

vxy vxy vxx vx vh ; a2 ¼ x ; a3 ¼ ; a4 ¼ ; b1 ¼ x ; b2 vx vh vx vh vx vh vh vh y y ; b4 ¼ ¼ x ; b3 ¼ vh vx vh ðI:9Þ

a1 ¼

G1 ¼

x2x þ x2y xx hx þ xy hy x2 xh ; G2 ¼ ; G3 ¼ x ; G4 ¼ x x J J J J

Therein, sis the streamwise direction, q is the flow rate per unit width, Ls is the reach length, q is the discharge per unit width at an upstream end, and H is the water level at a downstream end. The discretized forward model is obtainable using a finite differential approximation with a uniform timeespace resolution in a staggered numerical grid as shown below.

ðI:10Þ xy hy xx hy þ xy hx 2xx xy G5 ¼ ; G6 ¼ ; G7 ¼ ; G8 ¼ J J J J x2y

U1 ¼

U5 ¼

In principle, adjoint equations for the shallow-water model can be derived on both discretized and continuous forms. The complete form of adjoint equations are shown on a continuous form in Appendix I. For this study, the adjoint equations are derived directly from a discretized form of shallow-water equations, as demonstrated by Das and Lardner (1992). For illustration, we consider a simplified one-dimensional flow system in a straight open channel having a rectangular cross section for the subcritical flow condition herein to derive a discretized form of adjoint equations and relevant boundary conditions for a shallow-water flow. In such a flow, the governing equations and boundary conditions in the Cartesian coordinate system can be described using a simplified shallow water model as follows.  vh vq Continuity equation ¼ 0  t  T; 0  s  Ls vt vs ðII:1Þ Momentum equation

 vq vH ¼ gh 0  t  T; 0  s vt vs  Ls ðII:2Þ



Boundary conditions



q ¼ q at s ¼ 0; H ¼ H at s ¼ Ls ðII:3Þ

 hwþ1  hwp p Dt

qwpþ1  qwp 0  w  Nw ; 0  p  Lp Ds

Discretized momentum equation

Discretized boundary conditions

ðII:4Þ

qwþ1  qwp p

Dt w hwp þ hwp1 Hpw  Hp1 0  w  Nw ; 0  p  Lp ¼ g 2 Ds



xy hy xx hy þ xy hx 2hx hy ; U6 ¼ ; U7 ¼ ; U8 ¼ J J J J h2y

Appendix II. Discretized form of adjoint equations and boundary conditions for a simplified shallow-water flow model





h2x þ h2y xx hx þ xy hy xh h2 ; U2 ¼ ; U3 ¼ x x ; U4 ¼ x J J J J ðI:12Þ

ðI:13Þ



Discretized continuity equation ¼

ðI:11Þ

11



ðII:5Þ

qwp ¼ q at p ¼ 0; Hpw

¼ H at p ¼ Lp ðII:6Þ In those equations, p is the integer suffix signifying grid points along the s direction, w is the integer suffix indicating the discrete time levels, Nw denotes the terminal number which corresponds to the simulation period (w ¼ Nw at t ¼ T ), and Lp represents the terminal number for the stream reach ( p ¼ Lp at s ¼ Ls). In addition,Ds is the grid interval. The variable q is discretized as qwp at the grid point ( p,w), whereas the variables h and H are discretized respectively as hwp and Hpw at the grid point ( p þ 1/2,w). According to work reported by Lardner et al. (1993), the discretized form of adjoint equations and boundary conditions can be derived from transposing a tangential form of these equations and conditions. For the present simplified model, the results are given as follows.
w  w ðlh Þwp  ðlh Þwþ1 Hp  Hp1 hwp þ hwp1 wþ1 p ¼ g þ lq p Dt 2Ds 2Ds
w  Hpþ1  Hpw hwpþ1 þ hwp wþ1  g lq pþ1 2Ds 2Ds ðII:7Þ

w wþ1 lq p  lq p Dt

ðlh Þp

wþ1

¼

 ðlh Þp1 Ds

wþ1

w Dt wþ1 w lq 0 ¼ ðlh Þ0 ; ðlh ÞLp Ds
w  HLp  HLwp 1 hwLp þ hwLp 1 wþ1 þ lq Lp ¼ gDt 2Ds 2Ds

ðII:8Þ

ðII:9Þ

Therein, lq is the adjoint variable (or the Lagrange multiplier) for the variable q. The algebraic equations above are derived using the same numerical schemes as described in the

Please cite this article in press as: Yoshida, K., Ishikawa, T., Flood hydrograph estimation using an adjoint shallow-water model, Journal of Hydro-environment Research (2014), http://dx.doi.org/10.1016/j.jher.2014.12.003

+ 12

MODEL

K. Yoshida, T. Ishikawa / Journal of Hydro-environment Research xx (2014) 1e12

main text. Equation (II.9) presents the discretized boundary conditions for the adjoint variables lq and lh. As shown in Equations (II.7) and (II.8), the adjoint equations are computed backward in time from the future to the past. The initial conditions of the adjoint variables are usually assigned zero at the terminal time level (w ¼ Nw). The derivation process shown here was used and extended to obtain the discretized form of adjoint equations and relevant boundary conditions used for this study. References Adrian, R.J., 1991. Particle-imaging techniques for experimental fluid mechanics. Annu. Rev. Fluid Mech. 23, 261e304. Atanov, G.A., Evseeva, E.G., Work, P.A., 1998. Variational problem of waterelevel stabilization in open channels. J. Hydraul. Eng. 124 (1), 50e54. Bousmar, D., Zech, Y., 1999. Momentum transfer for practical flow computation in compound channels. J. Hydraul. Eng. 125 (7), 696e706. Bradley, A.A., Kruger, A., Meselhe, E., Muste, M., 2002. Flow measurement in streams using video imagery. Water Resour. Res. 38 (12), 1315. Cameron, H.L., 1952. The measurement of water current velocity by parallax methods. Photogramm. Eng. 18 (1), 99e104. Chow, V.T., 1959. Open-channel Hydraulics. McGraw-Hill, New York, pp. 604e605. Das, S.K., Lardner, R.W., 1992. Variational parameter estimation for a twodimensional numerical tidal model. Int. J. Numer. Methods Fluid. 15, 313e327. Ding, Y., Wang, S.S.Y., 2006. Optimal control of openechannel flow using adjoint sensitivity analysis. J. Hydraul. Eng. 132 (11), 1215e1228. Ding, Y., Jia, Y., Wang, S.S.Y., 2004. Identification of Manning's roughness coefficients in shallow water flows. J. Hydraul. Eng. 130 (6), 501e510. Fread, D.L., 1975. Computation of stageedischarge relationships affected by unsteady flows. JAWRA J. Am. Water Resour. Assoc. 11, 213e227. Fujita, I., Hino, T., 2003. Unseeded and seeded PIV measurements of river flows videotaped from a helicopter. J. Vis. 6 (3), 245e252.

Fujita, I., Muste, M., Kruger, A., 1998. Large-scale particle image velocimetry for flow analysis in hydraulic applications. J. Hydraul. Res. 36 (3), 397e414. Fukuoka, S., Nagai, S., Sato, H., 2005. Evaluation of hydrograph of flood discharge and flood water storage in the main river and tributary e in the case of the Tone River and the Watarase River. Annu. J. Hydraul. Eng. JSCE 49, 625e630 (in Japanese). Fukuoka, S., Watanabe, A., Tabata, K., Kazama, S., Gocho, H., 2006. Evaluation of hydrograph of flood discharge and roughness coefficient at the branched section in the Tone River and the Edo River. Annu. J. Hydraul. Eng. JSCE 50, 1165e1170 (in Japanese). Gordon, R. Lee, 1989. Acoustic measurement of river discharge. J. Hydraul. Eng. 115 (7), 925e936. Khatibi, R.H., Williams, J.J.R., Wormleaton, P.R., 1997. Identification problem of open-channel friction parameters. J. Hydraul. Eng. 123 (12), 1078e1088. Kinoshita, R., 1984. Present status and future prospects of river flow analysis by the aerial photograph. J. JSCE 345 (II-1), 1e19 (in Japanese). Lardner, R.W., Al-Rabeh, A.H., Gunay, N., 1993. Optimal estimation of parameters for a two-dimensional hydrodynamical model of the Arabian gulf. J. Geophys. Res. 98 (C10), 18229e18242. Minoura, Y., Ishikawa, T., Yoshida, K., 2008. Flood flow analysis of aerial photographs by image correlation method with support of stereo visualization. In: Proceeding of 16th APD-IAHR, Nanjing, pp. 999e1004. Muste, M., Xiong, Z., Sch€one, J., Li, Z., 2004. Flow diagnostic in hydraulic modeling using image velocimetry. J. Hydraul. Eng. 130 (3), 175e185. Nagata, N., Hosoda, T., Muramoto, Y., 2000. Numerical analysis of river channel processes with bank erosion. J. Hydraul. Eng. 126 (4), 243e252. Nezu, I., Nakagawa, H., 1993. Turbulence in Open-channel Flows. IAHR Monogragh, Balkema, The Netherlands, pp. 53e56. Nihei, Y., Kimizu, A., 2008. A new monitoring system for river discharge with H-ADCP measurements and river-flow simulation. Water Resour. Res. 44, W00D20. Sanders, B.F., Katopodes, N.D., 1999. Control of canal flow by adjoint sensitivity method. J. Irrig. Drain. Eng. 125 (5), 287e297. Sasaki, Y., 1970. Some basic formalisms numerical variational analysis. Mon. Weather Rev. AMS 98, 875e883. Yoshida, K., 2010. Inverse estimation of timeeseries variation of bed roughness coefficients using adjoint shallow-water model. J. Appl. Mech. JSCE 14, I_533eI_540 (in Japanese).

Please cite this article in press as: Yoshida, K., Ishikawa, T., Flood hydrograph estimation using an adjoint shallow-water model, Journal of Hydro-environment Research (2014), http://dx.doi.org/10.1016/j.jher.2014.12.003