Home

User's Guide - Applicaties Helpdesk Water

1. 154 5 4 5 Use of the Hessian sssssssennnee 154 54 0 Concl siOns onec due e p e eden 155 5 5 THE USE OF CURRENT VELOCITY MEASUREMENTS 156 5 6 FRICTION FORMULATIONS 0cccccceseeeeeeeeeeseeeeeeeeeeeees 157 5 6 1 INTOGUCTION 23 4 21 a aono eed 157 5 6 2 Theoretical background ee 158 5 6 3 Implementation esee 160 5 7 FINITE DIFFERENCE CALCULATIONS eenene 162 5 8 VELOCITY AND DISCHARGE BOUNDARIES 164 6 165 Z APPENDIGES vin oh ono e D E dada 167 7 1 167 1 2 REFERENCES dts etel deti rdeqn 168 7 3 LIST OF SYMBOLS cccccccccccceeeeceeseeeeeeseaeeseessaeeeeeees 170 7 4 ADVECTION IN THE CASE OF BOUNDERIES AND DRYFALL 174 7 4 1 Closed 668 174 7 4 2 Open Boundaries 192 7 4 3 Spherical and curvilinear 5 193 Version 1 0 October 25 1996 Introduction Introduction It is generally recognized that mathematical models based on shallow water equations are useful for solving hydraulic problems in civil engineering Such mode
2. seseeee 129 5 1 4 Validation of the calculation 55 131 5 1 5 nione a daa 133 5 1 6 Penalty for deviation from prior estimates 133 5 31 6 14 dntroduction inira 133 5 1 6 2 Theoretical 135 5 1 6 3 Practical US irie ee edet rete tege 139 5 1 7 Concl siOorn 5 esha e ed i dois 140 5 2 THE USE OF WEIGHTING FUNCTIONS FOR PARAMETRIZATIQON indier n HER ERE RRRRRXEYERYAR ENTER 141 Version 1 0 October 25 1996 V Technical documentation of WAQAD 922 15 Introduction ee eee 141 5 2 2 Theoretical background in the case of triangles 142 5 2 3 Computation of the gradient in the case of triangles 143 5 2 4 Implementation in the case of triangles 143 5 2 5 Conclusion deg 144 5 3 COUPLING HARMONIC CONSTITUENTS eee 145 5 3 1 Introd ction 2 145 5 3 2 Theoretical background e 145 5 3 3 Implementation iani eee eee 149 5 4 ANALYSIS OF THE GRADIENT IN TIME AND SPACE 150 5 41 Introd ctlon 5 eterno ei ed 150 5 4 2 Theoretical Background e 151 5 4 3 e ddaa adaa aa Tid 153 5 4 4 Practical use and interpretation of results
3. AU as 3S Vs AUr 4544 FEC Vinn AUR nt FSC Vmn m n mni 4 m n which equals m n amp AVing Aut S Vina c In this expression certain derivatives of discrete model variables are found namely Av and Au So this cross advective term uv has influences on the adjoint momentum equations which are solved in the first and second stages First it contributes to the adjoint v momentum equation in the first stage corresponding to Avi Therefore we have to examine the v momentum equations at the locations m 2 n m 1 n m 1 n and m 2 n U m 2 n 4U m 1 n A teks re 1 wy 1 es Via va Ees v os 3 98 v y AU ntn v y eo n ane REDIT A e Version 1 0 October 25 1996 The 2D adjoint model Second uv contributes to the adjoint u momentum equation in the second stage corresponding to Therefore we have to examine the v momentum equations at the locations m n m 1 n m 1 n 1 and m n 1 Aun vv dan ES Vin n FSC Vann 3 99 vy Vena 1 n 1 aol ME vy 184 Vinni 2 53 Technical documentation of WAQAD Summary In summary the adjoint advective terms are ignoring the boundary treatment and the degenerations caused by impeded points e Inthe adjoint umomentum equation associated with mn n Um 1 n 7 Um 1 n 7 Um 4n 1 1 T 3 3 100 Vu Jmn
4. afl Contribution of the continuity equation in the first stage in The 2D adjoint model m n 1 1 vna i 2 mn Kenn Kl net Ene An Contribution of the WAQUA boundary condition vc Joni By adding the contributions and rearranging the variables the adjoint boundary condition is described by 1 LE vv ant ve dma Caran 5 1 pum Kan Jmn t Ve mmnet Am n Bnet Diagonal open boundary Some models have an open boundary which is diagonal In this case the boundary condition is a combination of two boundary conditions For example in the case of the Continental Shelf Model the boundary condition is a combination of the left hand and upper boundary conditions In this case the boundary condition becomes Po 1 vcn E kc acc RM Um M vul vu 2 Uma veta A amp Kien Eas 1 g t frvwl m n 1 L vy ae Vv inci ve mn Ean 2 1 Ve mnt t nsa Kri n 1 Ehn Remarks e Only if uncertain parameters of the WAQUA boundary conditions are estimated is it necessary to calculate the corresponding adjoint variables at the boundary e The adjoint boundary conditions are solved after solving the adjoint model equations This is allowed since the Version 1 0 October 25 1996 67 Technical documentation of WAQAD boundary condition is not defined by state variables at interior points of the model and since it i
5. 24 05 y fruwl mn 1 36 g Una Um vna C In10 CZ H E M 2 4 i 2 frvwl n 1 36 g Vm l Vmn C In10 Jc for H gt 1012 w and equal to the friction formulation for Chezy otherwise The gradient for the bottom friction also changes according to the adaptations resulting from the friction formulations For a detailed derivation the reader is again referred to Verlaan 1994 and Brouwer 1995c Table 5 1 shows how the gradient changes as a result of other friction formulations Version 1 0 October 25 1996 159 Technical documentation of WAQAD friction absolute percentual formulation White Colebrook table 6 3 Adaptations of the bottom friction gradient as a result of other friction formulations where grcis the bottom friction gradient gre v 2H fruwl in u points and gre w 2H frvwl in v points e The percentual gradient is calculated by percentual grad absolute grad friction e The adaptations according to White Colebrook are only calculated if 1 0129 H 12 W 5 6 3 Implementation According to other friction formulations the coefficients fruu fruv frvu and frvv do not change However the coefficients fruwl and frvwl do change they acquire an additional factor Manning 4 3 White Colebrook 1 36 if H gt 10129 w CIn10 12 and 1 otherwise The changes in the gradient resulting from the friction f
6. Flow models used to be calibrated fully manually by experts This trial and error approach was very time consuming and laborious In recent years advanced mathematical techniques in the field of optimal control theory have been successfully applied to model calibration The adjoint model has become particularly popular In the early nineties RIKZ developed WAQAD the automatic calibration program based on the adjoint method The theory underlying WAQAD was derived by Ten Brummelhuis 1992 and implemented by J R Brouwer This first version of WAQAD was designed in order to calibrate the two dimensional shallow water flow model WAQUA WAQUA models were calibrated automatically by estimation of depth parameters bottom friction parameters or the wind stress coefficient Cy WAQAD was subsequently extended to calibrate WAQUA models by estimating the boundary conditions Later a second version of WAQAD was implemented in a SIMONA environment resulting from WAQUA in SIMONA This version calibrates WAQUA in SIMONA models by estimating the depth bottom friction and boundary parameters However estimation of the Cd coefficient was left out during this conversion This second version is the current WAQAD version except for some extensions described by Brouwer 1995a and Brouwer 19950 This chapter describes the 2D adjoint model as it is implemented in the current version of WAQAD version 2 1 see the User s Guide WAQAD Brouwer 1996b
7. The u contributions in the continuity equation are eliminated by means of the u momentum equation with a substitution There results for each row a tridiagonal matrix equation which consists solely of G contributions 6 C Cm 6 3 28 6 m G 6 Sinus This matrix equation is solved with Gaussian elimination double sweep method Stelling 1984 After here the u momentum equation can be solved by simple substitution v momentum in the first stage The v momentum equation in the first stage of the ADI method is solved on the basis of the solution of the model equations at time t The term is discretized as Vin Vna 3 29 2A Y Here the difference v is approximated implicitly it is expressed in v velocities at the new time level Therefore the structure of the matrix equation will alter Without the cross advective term uv and the viscosity terms all the v momentum equations for a column in the domain form a tridiagonal matrix equation As mentioned earlier the well known double sweep algorithm can be applied to perform this task 31 Technical documentation of WAQAD 32 Difficulties arise with the implicit approximation of the cross advective term uv Um n Vna 3 30 in which TIT i eb Lus cba PU Pls deal 4Vm t n Vm 2 n if Umn gt 0 2AX SA Vmn SV 4 Vm 2 n if Umn lt 0 2A X This is the so called second order upwind approximation
8. V V Adjoint Advective Terms As mentioned earlier advective terms only occur in both momentum equations For each term at each stage in the computation process the adjoint counterpart will be deduced VVy The v momentum equation is solved in the first stage of the ADI method Here v is approximated implicitly The discretisation of vv is as follows T Vin ntt n 3 34 2AY Applying the variational formalism Acad A V m n 1 2A Y V m n 1 Bi V m n 1 2A Y m n As 3 35 In constructing the adjoint v momentum equations the unknown terms with A y and A v are important The derivation will be done step by step for a fixed location m n far away from the boundaries and impeded points At the surrounding grid points the normal equations are solved 33 Technical documentation of WAQAD In the conservation of momentum in the y direction at point m n 1 the term vv occurs as V m n 27 Vion Vm n 3 36 M 2 VER Introduce Lagrange multipliers and apply variational formalism Vm n 1 v mna MEL RES AVin n T vv na 2A Y 3 37 Vm n 27 V m n Vm n 1 uh A mmt 2A Y Javan Saa Javna Remark the dots in equation 3 37 and the following equations denote the other non advective adjoint terms Inthe conservation of momentum in the y direction at point m n 1 the term vv occurs as nn Ynn 3 38 Vm n 1 Introdu
9. amit eat R 68 3 3 THE ADJOINT SOLVER cccccccccseseeceeeeeeeeseeeeeeeseaeeeeeees 72 3 3 1 Method for solving the adjoint model equations 72 3 3 2 AdVeCctlon oo a o D e d ie 80 3 3 3 Dryfall Sore de Tot uon oso tata 91 9 4 93 3 4 1 Bottom friction nen ire e e beee eE a a 93 34 2 Depth irni ash ep 94 3 4 3 Boundary conditions eeee 98 3 5 REMARKS ON WAQAD WITH RESPECT TO WAQUA 100 4 THE 3D 101 4 1 INTRODUCTION 0ccccccseseceeeeeeeeeeseeeeeceeseeeseeseaeeeeeees 102 4 2 THE ADJOINT MODEL 103 4 3 THE 109 4 4 THE GRADIENT edoan aaan eaaa apania eaa anaie nnn nnn 118 AA Bottorm friction uiri a aG 118 44 2 Depth nenii ni ei Gita i o Wie aas iea 119 4 4 3 Vertical VISCOSIty tede 122 4 5 REMARKS ON WAQAD WITH RESPECT TO TRIWAQ 123 5 POSSIBILITIES OPTIONS AND SHORTCOMINGS 125 5 1 IDENTIFIABILITY AND THE USE OF THE PENALTY 126 5 1 1 Introductions coire oe ecce eee 126 5 1 2 Some examples incre 126 5 1 3 Splitting the calculation
10. Tfimona Technical documentation WAQAD Werkdocument RIKZ OS 96 151X Rijkswaterstaat Ministerie van Infrastructuur en Milieu Technical documentation WAQAD The automatic calibration program Version 1 0 October 25 1996 Maintenance see www helpdeskwater nl waqua Copyright Rijkswaterstaat Preface Version 1 0 October 25 1996 iii Technical documentation of WAQAD Preface The technical documentation of WAQAD was written as part of the project Extensions of WAQAD contract RKZ 232 commissioned by the National Institute for Coastal and Marine Management RIKZ and carried out by SIMTECH an independent engineering consultant specializing in civil and software engineering It is a comprehensive compilation of the mathematical and technical aspects of WAQAD the automatic calibration system The following persons contributed text J R Brouwer SIMTECH J v d Linden RIKZ E E A Mouthaan TUD I D M Rozendaal SIMTECH M Verlaan RIKZ E V L Kuijper RIKZ supervised and coordinated the editing and rewriting The English text was revised by J Burrough Boenisch Science Editing amp Translation This technical documentation is intended to meet the needs of persons wanting to know about e the mathematical and theoretical background of WAQAD e the WAQAD system itself and how the mathematical model is implemented Persons interested in using the WAQAD system in order to be able to calibra
11. The extensions described in this section can be used in several ways Probably the most important use is to depict the spatial gradient for depth and bottom friction to select rectangles or triangles for the optimization process Using these depictions the situation in which a non zero gradient cancels over a selected rectangle can be avoided This is important since a zero gradient leads to no adaptation while the criterion can still be lowered using other triangles Furthermore the magnitude of the gradient in space gives an indication of the sensitivity of the criterion to changes and it is probably possible to make a detailed estimation of the parameter In areas where the gradient is small a detailed estimation is probably not possible and larger rectangles should be selected The application of these rules is shown in an example in Verlaan 1994 Another application is to use the gradient as a function of time for the selection of the simulation period The actual gradient used in the minimization procedure is the average over time of the function This average corresponds to the mean of the function Usually the function will show deviations from this average in time The random fluctuations show the propagation of the noise and indicate if the time interval is long enough to obtain a sufficiently accurate estimate of the gradient The systematic fluctuations indicate model imperfections If these systematic fluctuations change in tim
12. 1 where g can be shown to be 1 0 1 sino 140 This is for infinitely long series For finite series the truncated convolution is divided by the sum of the corresponding g I In two dimensions the problem can be solved analogously but it becomes slightly more difficult to handle of all exceptions becomes Implementation Several programs are used at RIKZ for more common variables like water levels None of these programs is currently able to cope with the new variables introduced in this section After comparing various alternatives a sequence of two programs was chosen First the data file is read with a FORTRAN program and the filtered data are written to an ASCII file In the second step a script written in MATLAB language is used to plot the pictures to the screen or the printer When this program is started a menu appears on the screen and the user can use the mouse to select plots and analyse the data More than one window for plots can be used at the same time and several plots can be displayed in one window Also a number of options for the presentation of the data are offered At the moment the FORTRAN program ADPIC is fully implemented There is also a working test version VIZ of the MATLAB script See the WAQAD User s Guide for more information on the ADPIC and VIZ programs Version 1 0 October 25 1996 153 Technical documentation of WAQAD 5 4 4 Practical use and interpretation of results
13. 2 AAK Version 1 0 October 25 1996 The 3D adjoint model After eliminating the adjoint u velocities from their adjoint momentum equations 4 12 in the adjoint conservation of mass 4 4 a tridiagonal system in adjoint water levels is formed which can be resolved by the well known double sweep method Stelling 1984 Am ve m 1 Bn vc Cm vc Dn 4 13 in which Am 1 3 E 9 AA Kona Gea l 28 Kn 1 AG ECC A d g CC ha g AA i Cae UE AE K 1 we 2 K AE 118 Technical documentation of WAQAD 114 and vu v EAE Dm zat T uel F if 0g v ut u EN elected iy HE E v ut u h 2h m E m Vy k 1 elo Ox Vuk hy hi h K 1 S n n Vy u u 4 2h h eec dis k d W 5 il x 9 h thi K vo s ua K xt 1 Uo o0 wk Al k 2 kalv ynie m w a v J vy V he 21 lor zonk e P E eyes en K 1 v v vi 2n h Zlo o a EE L a n ds n ehe IS vy v v i 2 k 1 LP E 9 uou los ox Kv Kc hi o Veal ti lox ok Kvv X he Now can be solved and thereafter v m can be computed per layer with the adjoint u momentum equation Note that the bottom roughness formulation equals
14. AU1 At i If GEOUUX m n 1 AU1 Umer Unio 7 10 2AX AU1 D Uhn 7 Unin 2AX If GEOUUX m n 2 7 11 AUT un Un iin if Un n gt 0 AX 0 otherwise AU1D Uinn Unan 20 0 otherwise If GEOUUX m n 3 7 12 AU1 msnzum ity AX 0 otherwise u u m 1 n m n if lt 0 AX 0 otherwise AU1D 177 Technical documentation of WAQAD So far we handled the degenerations in the adjoint advective terms caused by impeded points in the adjoint u impulse equation In the next part we focus on the deformations of the second order difference operator S and S ee Degenerations in the adjoint v impulse equation associated with caused by impeded geometry In the adjoint v impulse equation associated with Avi the second order difference operators S central and S upwind appear In this appendix the boundary treatment of these operators in the adjoint computation process is indicated These approximations depend on the dry wet status of the surrounding points Therefore the time dependent local flow situation has to be determined first to figure out which approximation has been used in the computation process of WAQUA Sj Um n and Sus We define for each u velocity point the logical variables True or False KBB KB KT and KTT with which the arrays GEBVUY and GETVUY are determined KBB m n KHV m n 2 gt 0 a KHU m n
15. These v velocities at the new time level are perpendicular to the ones which set up the tridiagonal structure Therefore these equations are solved column by column in a Gauss Seidel iterative way If a v velocity at the new time level from a nearby column is needed in the computation process and has not yet been computed the most recent value should be taken instead This enables the equations to be resolved again by solving tridiagonal matrices over grid lines This so called predictor corrector procedure will now be explained in more detail First of all the direction has to be determined whether the tridiagonal systems over columns should be solved from left to right or from right to left This depends on the dominant flow direction of u defined by q and the iteration number defined by q 0 if Y 120 1 4 if Yuso a q q 4 1 4 83 3 2 2 3 Version 1 0 October 25 1996 The 2D adjoint model The second order upwind approximation q 2 is 3459 4v 3 4 S9 70 2AX l 3 32 g u veh vig 8 9 5 0 2 The maximum number of iterations is two These v momentum equations are solved column by column in the dominant flow direction of the u velocity If the sign of u is constant then these equations are solved in one iteration otherwise a second step is necessary This step proceeds in the opposite direction V Vna V VV V 3 33 V V
16. Verlaan 1994 In this section the derivation of the friction formulations is described briefly since part of the derivation in Verlaan 1994 has already been described in chapters 3 and 4 of this documentation The friction terms for the Chezy formulation are given by equations 3 11 and 3 14 Z una Y vmn m n Quy aussi ven n frUV mn g Um n Vmn E TaD n fruwl mn eee t Waal frvu ae M C Was Y vmn fiw 9 Unn a r HAR nn n 2 2 frvwl a ER Vm n 7 Vmn In these formulations C is the Chezy coefficient where Cu is defined at a u point and Cv at a v point For the other friction formulations coefficient C has a different dependence on H e Chezy formulation C constant e Manning formulation 6 ont M where M is Manning s parameter Possibilities options and shortcomings e White Colebrook formulation C 18 log max 1 0129 where W is the White Colebrook parameter The coefficients fruu fruv frvu and frvv do not depend on the friction formulation chosen The coefficients fruwl and frvwl however do depend on the friction formula since they depend on the total water depth H For the Manning and White Colebrook formulation they become e Manning 7 4 g Umn ee y vi y 30 fruwl in 4 g Vmn y Vn y 3 frvwl mn e White Colebrook
17. Version 1 0 October 25 1996 35 Technical documentation of WAQAD 36 e Inthe conservation of momentum in the y direction at point m n 1 in the second stage the term vv occurs as f Vian E Vmn 2 V 3 47 cu 2 Y a Introduce Lagrange multipliers and apply variational formalism v 2 a o M AME cde v j 2A Y n 2 3 48 Vin n Vs n 2 Vin n 1 AV aci MEL LN AV inn at 2AY i 2 In 3 44 3 46 and 3 48 we focus on the terms with AVE and Av The latter can be used for constructing the adjoint equation for v s In 3 48 v zi ME s 3 49 m In 3 44 v AVI RU LL UBL Lo E PM 3 50 n 3 44 v A A 2AY 3 50 In 3 46 vv MM Javan 2 3 51 Finally some parts of the adjoint v momentum equations can be determined To compute the gradient the cost function has been augmented with additional terms each of which is the product of an undetermined multiplier and a model equation After differentiating this Lagrange function the summation is rearranged some additional terms are now the product of an unknown derivative and an equation in Lagrange multipliers In the other additional terms no unknown derivatives of discrete model variables occur This procedure enables the derivative of the cost function to be evaluated providing the Lagrange multipliers satisfy the adjoint equations Rearranging the summation in 3 40 3 41 3 42 3
18. for instance the model is going to be used to compute sediment or pollution transport the current velocity is very sensitive to the tidal average of the water level that is prescribed at the boundary If only measurements of the water level are available and the average water level at the boundary is to be calibrated the problem can be very ill conditioned To illustrate this assume no tide is present and the effect is dominated by the following part of the shallow water equations in one dimension eG gu _ 5 g ox CH Two observations can be made from this equation The velocity is highly dependent on the bottom friction value which is not known very accurately and also the velocity is very sensitive to the difference in water level between two points Substituting the values C 50 s m H 10 m then a 10 change in the average velocity u 0 1 to u20 11 ms 4 corresponds to a difference of of about one centimetre per 100 kilometres This kind of accuracy is very difficult to obtain with water level measurements and current velocity measurements can greatly improve the performance of the calibration In general automatic calibration performs best when as many reliable measurements as possible are included since all other variables in the model will have to be guessed from the available measurements and the model which sometimes means that variables for which no measurements are available will not be accu
19. if man SO 3 123b or Version 1 0 October 25 1996 The 2D adjoint model 3 1230 au vy Yun T mi 0 Saran ee x vas if Umer n 0 3 123d tas ER Jif Onion 0 3 123e dm eee es 2 il on 2 6 u 2 In the adjoint v momentum equation associated with A Vm n Ve oia Foy 3 124 En 2AN nano Enon 2 Js In the adjoint u momentum equation associated with Ann b us vu Kin 2AE Km 1 n 2 Um 1 n TRE n 1 Wu mani T TR Az Mu mt nt 3 125 i capis am vu Jans Vin v nn A AV n En n 2 1 2An Enni 1 2An M 3 126 vu Jansi vu Janis V m n 2 En n1 12An En n 2 12An TE E v 19 v in a SI l 3 127 vy Imin 1 T Wied 1 Vinci 59 Technical documentation of WAQAD vv mna Sel Vm n vy ment Vmetn i i 3 128 v 154 Vm 1 n 1 v Jic ASK Vm n 1 vn 1 Vs n 2 tua e N if gt 0 3 129 vu nme En n2 2 Vm n 2 a une Aino i Feo 3 129b MATS Enna 2 f AN 1 3v n es GSO j En n 2Nn M or 3 129c Td 3Vn 2 t v Jnn mn if 0 2 van vob Amna i T 3 129d nai Enn 2An j E T 1 n2 y gat if lt 0 3 129 vu 1 Bs 2 2 e Remark The boundary treatment and the degenerations caused by impeded points are ignored in this s
20. ram m n m n f E j a j 2 n6 2 E g edi nin g Una Urn ma 0 Kiin Ag Go o DhstC in which ur u velocity at time t At m n e The coupled equations in the second stage This is the second step of the ADI method in v direction The WAQUA variables at time t 2At and the u velocities at time t At are known All v velocities and water levels at time t At are calculated such that the continuity equation is 1 cat em Gin 3 8 Kin En AY Yin Vn Km n AY inating Kmn 11 0 1 uy uy Ad H Jm n Um Enn B H Ent 1 T An the v momentum equation is 1 v vV X At m n m n 2 5n g omisi Sms 9g Vrau m n y Vin y 0 Ei An G Dina TC m n PET n RT 3 9 22 Version 1 0 October 25 1996 The 2D adjoint model Step 2 The differentiated WAQUA equations The WAQUA equations are differentiated to C u and v applying the variational formalism e The differentiated v momentum equation at the first stage from t gt t 72At for vv mn is 1 AVinn 7 A Vm T nA m n A A g Ginea Cmn 3 10 Enn An frVUmn Ki frVV m n Ann frvwl m n in which frvu contribution of bottom friction at v velocity point to the adjoint u momentum equation at time t frvv contribution of bottom friction at v velocity point to the adjoint v momentum equation at time t frvwl contribution of bottom friction
21. 2 Kmr Krnn 1 Ee nt An Contribution of the v momentum equation in the second stage in m n 1 M m n 1 Enn 1 An 2 vy Jana Contribution of the v momentum equation in the first stage from t gt t ZAt in m n 1 9I Per m n 1 Enn 1 An 2 Contribution of the continuity equation in the first stage in m n 1 1 vem 2 Vm Knn KO n 4 ES o An Contribution of the WAQUA boundary condition 65 Technical documentation of WAQAD 66 Jmn By adding the contributions and rearranging the variables the adjoint boundary condition is described by 1 9 frvwl m n 1 Jana 4 Ce ane 2 1 F 2 Vma Knn Jes ve mn 1 An Ki n 1Eh n 1 Lower open boundary As for the upper open boundary the u momentum equation does not contribute to the boundary The continuity equation and the v momentum equation do Consider point m n to be on the boundary The contributions to the boundary are as follows Adjoint equation for AG 0 Adjoint equation for AC Contribution of the continuity equation in the second stage from t 2At gt t in m n 1 1 vc T g VmnKmn Kl net Ea An Contribution of the v momentum equation in the second stage in m n g M EAn 2frvwl Contribution of the v momentum equation in the first stage from t gt t 72At in m n g 1
22. 2 gt 0 KHV m 1 n 2 gt 0 KB m n KHV m n 1 gt 0 a KHU m n 1 gt 0 KHV m 1 n 1 gt 0 KT m n KHV m n gt 0 a KHU m n 1 0 KHV m 1 n 0 KTT m n KHV m n 1 gt 0 a KHU m n 2 gt 04 KHV m 1 n 1 0 Hereafter the arrays GEBVUY and GETVUY can be filled E ee 1 2 2 2 1 t E NN 1 a mil cael m Hi tr Table Il Definition of the arrays GEBVUY and GETVUY 178 Version 1 0 October 25 1996 Appendices The time dependent geometrical situation can now be distinguished and consequently the proper discretization in WAQUA first second order upwind or central or zero Stelling 1984 The boundary treatment of S Um n Jand S Uhn Jis as follows Determine the dry wet status of the u velocity points by the arrays GEBVUY and GETVUY Table at time t 72At Unne t E me for situation i S Up for situation ii iv 7 13 i 2AY 0 for situation v ix The upwind difference operator depends on the flow pattern of the v velocity field LEVINE O n 4 n n Bumn dus for situationGEBVUY 2 2AY S u 7 Cunn Unna for situation GEBVUY 1 WM i AY i 0 for situation GEBVUY 0 Ue t 4Ur n J Urine 2AY for situation GETVUY 2 pu Uninet 7 Urn 7 1 5 S Unn AY for situation GETVUY 1 0 for situation GETVUY 0 Now we focus on another part of the righthand side of the adjoint v
23. 5 n C z frvwl minl KE a 7 Ka net Bnet Enn An 2 An 1 g 1 EL Dn At 2frvwl ma ill DD g 1 IE An 2frvwl mallk Vv m n DD 1 esee ve mn vint An Kin Esa Enna 1 enfe ve mn1 E ve mn An Ke net Enna If at whole time levels water level measurements have been recorded the residuals of the observed and computed water levels are totalled to coefficient Dn Version 1 0 October 25 1996 77 Technical documentation of WAQAD The obtained system of matrix equations is transformed into a tridiagonal system 2 D B C Ve m1 B 2 B C Ve m2 E A B C Ve m3 E As 4 By a Cg a Dis A By T This system of equations is written as the following equivalent systems of equations De Ji _ Di 1 C R Vein D 1 C2 Ve m2 Ds 4 G Velma 1 Cn 1 D 1 jo Xx E Dar or Jmn eo Jmn Di vn These coefficients are determined recursively C CG E E Ba AnCy 1 A ne Dn An Dia Di Ba AnCn 1 a _ C a _ D and p Ci B Di B 78 Version 1 0 October 25 1996 The 2D adjoint model Solving backwards using the same formula gives ve mnt Dar vc m n D a C ere The adjoint v velocities are then calculated by substituting the adjoint water levels into the adjoint v momentum equation vc Jr Jin nd Eau
24. En vy Jun DDna CC AD UXD The adjoint u momentum equation at the second stage from t t gt t is an explicit equation XN vu Jn T vi fren frvu mnll vy Jos T vv mn 0 The adjoint v velocities at time t are calculated straightforwardly by substituting the adjoint u velocities at time t ZAt and the adjoint v velocities at time t Remarks e The solution method applied to solve the adjoint equations is anumerically stable method The structure of this method is similar to the numerical method used in WAQUA For more information on this subject see Stelling 1984 Some tests regarding the numerical stability of WAQAD are described in Mouthaan 1992 lf uncertain parameters of the WAQUA boundary conditions are estimated the corresponding adjoint variables at the boundary must be calculated The adjoint boundary conditions are solved after solving the adjoint model equations In WAQAD the adjoint boundary conditions are solved in module AD GRKB see subsection 3 4 3 In order to reduce the computational effort the adjoint water level in WAQAD is divided by the WAQUA transformation factors K and E in advance The adjoint water level variable is implemented in WAQAD as 79 Technical documentation of WAQAD 80 3 3 2 vc Jmn 7 This is important when the adjoint water levels are printed or plotted as output of WAQAD Advection This subsection describes the m
25. TRIWAQ by means of a constant value the algebraic model or the k eps model see the WAQUA User s Guide WAQAD can estimate the vertical viscosity parameters for all the options mentioned 4 3 Version 1 0 October 25 1996 The 3D adjoint model The adjoint solver This section deals with the method for solving the adjoint equations deduced from the Chavent method for the simplified three dimensional shallow water flow model This solution method also has a two stage time splitting structure However only the first stage is considered since the other stage is virtually identical due to the underlying ADI method Stelling 1984 Similar to the approach in the forward model the approximation of the vertical viscosity terms in the adjoint momentum equations is fully implicit For detailed information on the numerical treatment of the discretized equations in the forward model see Lander et al 1996 The inverse model equations are solved in reverse order compared with the three dimensional shallow water flow model In the first stage in WAQAD the adjoint continuity equation and the adjoint u momentum equation are solved in module AD SUV3 first and then the adjoint v momentum equation is solved in module AD VYD3 In the second stage in WAQAD the adjoint continuity equation and the adjoint v momentum equation are solved in module AD SVU3 first and then the adjoint u momentum equation is solved in module AD UXD3 Below the modu
26. The coupled equations in the second stage from t 2At gt t The WAQAD variables at time t 2At are known All adjoint v velocities and adjoint water levels at time t are calculated such that the adjoint continuity equation is Version 1 0 October 25 1996 27 Technical documentation of WAQAD 1 f XA Jan i voa ral vy Jos T E vy Jui Wi rs Edid Enn 1 j T vema t va vc nnt venns 3 21 Am Kin Eia Kl net Bane 1 2 mine Jui T E vc Jos vena Ki n 1 Ehn Kan En frvwl mn vv Judd vost 0 the adjoint v momentum equation is 1 p At vy Jans Gaia n nKmn vc Jut vela vc n n 1 ve inne 3 22 An Kis Esa Ki net Ent frvv mnl vy Ja aaa 0 e The adjoint u momentum equation in the second stage The WAQAD variables at time t At and the adjoint v velocities at time t 72At are known All adjoint u velocities at time t are calculated such that zal dau vu mn fen frvUnoll w nn idee J 3 23 Remarks e The adjoint model equations are solved backwards in time The initial condition for solving the equations is equal to zero e The adjoint model is activated by the misfits between the observed data and the corresponding model values Residuals of observed and computed water levels affect the adjoint continuity equation If a water level measurement has been recorded at a full ti
27. Time F F External forces in the x and y directions A two stage time splitting method is used to integrate the model equations for the two dimensional shallow water flow model The method is based on one used by Stelling 1984 Because the model is used in a wide range of flow conditions the discretization method used must be unconditionally stable and have at least second order accuracy The discretization in Version 1 0 October 25 1996 The 2D adjoint model time is analogous to an ADI Alternating Direction Implicit method Each stage describes the propagation of the model state in half a time step In both stages the terms of the model equations are solved in a consistent way with at least second order accuracy in space For example if a term in the model equation is approximated implicitly in time in one stage this term will be approximated explicitly in time in the other stage The spatially and temporally discretized equations for the two dimensional hydrodynamic model are the starting point for deducing the corresponding adjoint equations which are derived by following the next three steps e Determine the constraints in order to transform the constrained optimisation problem into an unconstrained one These constraints are the non linear discretisations of the shallow water equations described by 3 1 3 3 e Derive the linearised equations This is done in the following way regard all WAQUA variables in the co
28. Version 1 0 October 25 1996 Possibilities options and shortcomings Coupling harmonic constituents Introduction The boundary conditions for many model schematizations are described by harmonic constituents During the calibration process these boundary conditions are usually modified The calibration can be done by hand or by using WAQAD The angular velocities of certain constituents are very similar for example S2 30 00 hour K2 30 08 hour A minimum simulation of 360 30 08 30 00 hours 4500 hours half a year is required to separate these constituents Using WAQUA models a simulation period of approximately one month is a practical limit for the computational effort If the simulation period is too short another way must be found to determine the amplitude ratio and the phase difference between constituents S2 and K2 One possible partial solution to this problem is to couple the two constituents and then consider them as a single constituent The HATYAN program is available for coupling constituents when calibration is done by hand This section discusses how constituents can be coupled by using WAQAD for automatic calibration Theoretical background In this subsection the coupling of two constituents is described However the results described still hold if more than two constituents are coupled 145 Technical documentation of WAQAD 146 Assumptions Let en be two constituents
29. an area can be found by summing the derivatives of all the points in the area If we select an area with average derivative zero the gradient over this area will be zero even if the derivatives themselves are not Clearly this is not what we want since the parameter appears to be optimal because the derivative is zero But other areas can be selected that improve the model further The contributions to the gradient can be split in time as well as in space The resulting function of time can be interpreted as the effect on the gradient if the parameter were changed during a short period around the time in question For some parameters this interpretation may seem hypothetical but the information obtained can be useful even in these cases The programs described in this section enable graphs to be made of the gradient of parameters split as a function of time One direct application of the gradient split as a function of time is when selecting the period for the estimation First notice that only the average of the function is important for the gradient itself To obtain an accurate estimate of the 150 5 4 2 Version 1 0 October 25 1996 Possibilities options and shortcomings gradient the average to variance ratio should be sufficiently large When the signal to noise ratio is large a short period is enough but when the signal to noise ratio is small a long period is needed Of course this ratio is not the only aspect for th
30. are in which Um 1 n Am n l 2AX Bm n Um 1 n GA E Judi E Jaai 2AX Uis n i AV m n1 vu Jina 12A Y Janet 12A Y 12A Y AV v vu Jana Vu hus LLL y vy Jmn 9 Mia vy Jodi EE Vinin i v Jain yeaa v mn 1 16 Vinn mnt Sx Vm n vv Juss Gal Vinet n j vy Justi 154 Vm 1 n 1 vy 1s Vm n 1 12A Y 3 159 3 160 3 161 3 162 3 163 3 164 3 165 3 166 Remark Note that the coefficients A B C and D that are calculated in module ADSVU are not the same as those calculated in module ADUXD 89 Technical documentation of WAQAD Adding equation 3 111 to equation 3 158 yields Bin Brn 3 167 SVna ify 50 2 i M p 5 if v lt 0 2A Y me Adding equation 3 111 to equations 3 163 to 3 166 yields Dae Dios 3 168 ex if Vo geo 0 vu und E if Vm 0 vu m n if vi lt 0 hz SAY Vm n 2 The system is solved iteratively using a Jacobi method The computed values from the previous loop are used to compute the marked adjoint variables in equation 3 168 90 3 3 3 Version 1 0 October 25 1996 The 2D adjoint model Dryfall In this subsection an overview of the drying procedures in the WAQAD modules AD SUV AD VYD AD SVU and AD UXD is given AD SUV If dry
31. at v velocity point to the adjoint continuity equation at time t According to Ten Brummelhuis 1992 and Brouwer 1995b the bottom friction contribution is given by M g Um n Vm n frvu in 2 n n 42 2 C Dm T Cmn Vn us y 2 Vm n E few S m gan Dian ed Unn Vmn imt S Vra Formula 3 11 describes the Chezy formulation For further information on the friction formulation see section 5 6 of this documentation Remark The overlined variables in equations 3 10 and 3 11 indicate that they are interpolated to v velocity points 23 Technical documentation of WAQAD e The differentiated continuity equation in the first stage for voa 1 XA AG an A Cmn 1 1 T fe der A Ki Bs nel H mnEnn Au mat ou mn Enn AG mn AG m 1 n 3 12 1 1 B AR H Enn Au sant U m tin En 1n A Cray A Ch 1 1 H Jmn Kmn A Vmnt Vmn Kmn A Cmn A Gmn An 2 1 1 H Jmn Kn n A Vmn A tz Vm n 1 Kennet A Cari A Gmn An 2 e The differentiated u momentum equation in the first stage for Vu mn iS on AW wn AUmn z fmn AV mn x g E 3 13 Kin AS 1 At n 6 frUU mn AU mn frUV mn Av ma fruwl ma AG inn in which fruu contribution of bottom friction at u velocity point to the adjoint u momentum equation at time t 72At fruv contribution of bottom friction at u velocity
32. eigenvectors scaled by the square roots of the eigenvalues It is clear that the occurrence of an I problem can be characterized by the presence of large eigenvalues of the inverse In practice large is probably best described relative to the other eigenvalues During the iteration process of the BFGS algorithm an increasing part of the inverse of the Hessian is estimated The occurrence of large eigenvalues in this estimation seem a reasonable indication of the occurrence of an l problem figure 5 7 Cost function for measurements to be calibrated and control measurements o figure 5 8 Eigenvectors and eigenvalues For more information on the use of the Hessian see subsection 5 4 6 5 1 6 1 Version 1 0 October 25 1996 Possibilities options and shortcomings Regularization Methods for giving the cost function a more regular progress without causing I problems are often called regularization methods One such method will be described briefly in this subsection The other method called the penalty will be described in detail in the next subsection Reduction of the number of parameters If an I problem is caused by choosing too many parameter boxes for depth or friction it can be solved by reducing the number of boxes for example by fusing them Suppose that the calculation from l problem 2 has a cost function shown by figure 5 2 Fusing the two boxes results in both parameters being adapted b
33. equation in the second stage from t 72At gt t At at point m 1 n E 2 2 v y g MEUS Vista 1 C in ACen y Dinan t dire Contribution of the WAQUA continuity equation in the second stage from t 72At gt t At at point m n vee 5 Ulin Erin 2 Vin Knn Ki Ehn AE An Contribution of the WAQUA continuity equation in the second stage from t 72At gt t At at point m 1 n 1 veNnetn B 2 Umn Enn 2 Km 1 n Kmetn Eta AG An Version 1 0 October 25 1996 The 2D adjoint model e Contribution of the WAQUA continuity equation in the second stage from t 72At gt t At at point m 1 n 1 n 1 vena D 2 Umn i Enn 73 Vinetin Km 1 n Kivi Eine AG An e Contribution of the WAQUA continuity equation in the second stage from t 72At gt t At at point m n 1 n 1 vena 2 Umn Enn 7 2 Vmn Kmn Kinn 1 Ent AG An The expression of the gradient for depth in the first stage becomes 2 Y fruwa voa frUW met m n k Vu mon frUWl ma Vu net fruwl mae von Unn Enn vemen Uhn Enn J Kan Een AG Kaen Eun AG v zs Enna y vena Unit Enna y Kl met Bet AG Ki net Eh n1 ve mn Vm n Kmn ve inetn Vm 1 n Km 1 n Kis Ehn An Kmetn Ens An Vm 1 n Km 1 n J vom 7 Vm nKm n Kher n1 Erne An Kan Ehn An 97 Technical documentation of WAQAD And the expression of the grad
34. equation yields the same structure with other matrix coefficients Jian DDin 1 n7 CCm 1 n vc Jon E vc Jin 3 132 in which H Jmn 2 ax Um 1 n 7 Um 1 n At 2AX 2 fru is H uon m Um 4n vu ys _ At i 2AX DDm 1 n m 2 F 2 fruu Um 1 n Um 1 n At 2AX v CB GE Rn AX z frUUmn ums Uns At i 2AX For the sake of completeness the adjoint continuity equation advection included associated with AG will also be given 2 troy ire Ve Jmn D vc Jon g UY 1 XU Qna eit Vu Dan nant PR E Eve dnt ve Dan Insta ve Jnn 3 135 zU 1 1 1 1 AX C ve Dan cs vonga ras fruwl an E vu Vaan vu n 70 m 1 n Version 1 0 October 25 1996 81 Technical documentation of WAQAD 82 In this equation three adjoint water levels and two adjoint u velocities at the new time level occur They form a pentadiagonal matrix equation together with the adjoint u momentum equation if the corresponding equations for all grid points in a row of the domain are written down ve doce Ve v s ma Vu Vu Ve v a Ve Yum Ve ve osa Ve sa Ve 3 134 The adjoint u contributions in the adjoint continuity equation are eliminated by means of the adjoint u momentum equation with a substitut
35. m n 32g Vinay in Winn Qs y D e 20 The expression of the gradient for bottom friction in the first stage becomes ZEN A2 2 2 Vv mn g Vm n JT Vmn m n k ian y Dina Grn g Um in SUL js yn Vu m n F 3 6 mn Dm 2 C nih And the expression of the gradient for bottom friction in the second stage becomes 2 ydus 8 m mal u Vu E Uc 093 m n k Cink y Dim zh g Vm n us 1 4 mn y Date vy Thich m n The gradient for bottom friction is calculated in WAQAD module AD GRCZ 3 4 2 Depth In this subsection the expression of the gradient for depth is derived The parameter to be estimated is a correction factor for the depth in point m n An overview of the contributions of the WAQUA momentum equations and the WAQUA continuity equation is given first Then the expression of the gradient is deduced 94 Version 1 0 October 25 1996 The 2D adjoint model Contribution of the WAQUA v momentum equation in the first stage from t gt t 72At at point m n vv 9 Vm n urs y H Vm y Vv m n 3 pam 2 0n5 Dinan Gu Contribution of the WAQUA v momentum equation in the first stage from t gt at point m 1 n y E g Vm 1 n y y V 1 ACen y Dn MR j Contribution of the WAQUA u momentum equation in the first stage from t gt t 72At at point m n Ua aa
36. momentum equation Only the important formulas are presented below In the first stage the difference v is approximated by a second order upwind difference implicit This results in the following approximation of uv Um n 3 87 in which Oh 1 Um n Umn 4 Um n ae ee F Um 1 n l PAT ph 3Vmn 4 Vito Vm2 n if Un 0 2A X Qh Vm n 3Vm n 4 Vmin if Um n lt 0 2A X Version 1 0 October 25 1996 The 2D adjoint model Differentiating this expression with respect to an arbitrary model parameter yields A Um n 4S Vmn Aunt 1S Vs KS MS 3 88 A Um 1 n 1 i amp Vn A Um n1 i amp Vian added to A ae Hnn A Vm 1 n PIER 2 2 u E m n gt 0 m 2 n 2AX Um n or added to 3 Un n 4 Um n AVin n AV n 2A X 2A X 3 Un n Av if lt 0 m 2 n 2AX Um n Equal to amp Vina tus S AVinn First we concentrate on the counterpart of this term in the adjoint umomentum equation Therefore we examine the v momentum equations solved in the first stage of the ADI method at the locations m n m 1 n m 1 n 1 m n 1 these equations contain factors with the differential AU mn after introducing Lagrange multipliers and applying variational formalism Assembling yields Aus Winn 4S Vin Gps vv ere 49 Model Juni 49 Vinni 2 3 89 The expression between braces contributes to t
37. momentum equation at the first stage from t At gt t t is an explicit equation Am XN Lov T vv mn i fmn Vu m n P 0 The adjoint v velocities at time t 2At are calculated straightforwardly by substituting the adjoint v velocities at time t At and the adjoint u velocities at time t 72At AD SVU The adjoint v momentum equation and the adjoint continuity equation at the second stage from t At t are implicitly coupled equations The adjoint v velocity at time t is expressed as a function of the unknown adjoint water levels at time t vc Jn Dp m n DDr i i Ko net net Kin Ehn Coefficients CC and DD are determined by H Jmn Km n A PU ES m n At frvv mol CC 7 1 ou frvv PARU em DDr A 5 fr m n a N H Kn ve mn ven 1 Kh nEh A age n n n n t1 ni Gat frvv mnl The adjoint v velocity at time t is eliminated from the adjoint continuity equation This results for each row in a tridiagonal system of unknown adjoint water levels at time t An ve mn 1 ii C D 76 The 2D adjoint model Coefficients A B C and D are calculated by 1 1 g 1 p el Km n 1 An _frvwl m n ICG Kn n 1 ES n4 Enn 1 An 2 An zea 1 g Bi n E 7frvwi m n Bes An CC jp y IE m atti rele e 1 alvin Kin Vm n 1 Kw nil 1 An 45 At l Km 1 g 1
38. n 2 8 ANE ET 2 Vm 1 n Vm 2 n if Um n gt 0 AS 1 3 117 amp vna 7 2d SS if Um n SO AG At source code level there is no difference between WAQUA in curvilinear or spherical coordinates Everything is fixed by the choice of the variables Kj and Emn the distances between depth points in the and directions see subsection 3 2 1 From this together with the mathematical derivation in the previous chapters it follows that the adjoint advective terms for curvilinear models and models with spherical coordinates are compare with equations 3 100 to 3 111 57 Technical documentation of WAQAD 58 In the adjoint u momentum equation associated with Aur n TM MM y Um 1 n 7 Um 1 n dee y f Umi4n 7 Um 4n 3 118 LER e Ka 2 Km n 2 In the adjoint v momentum equation associated with Vna vv Jana oe vv Jn Bios Eso 2AN 777 7177 il B Mn 3 119 En n 1 2 End 2 NUS Js 2n Um 2 n 4 bw Jiz 4n Km 2 n 1 AE Km 1 n 12 j f B 3 120 vv 4u m 1 n _ w jm U m 2 n Km 1 n 12 A Km 2 n 12 vu Jan a Uhn a oy Unin t Ia vu epee Jmn 2 S 1 mf e Um n uetus Um 1 n v Jn a S unn vu m 1 aS in 3 122 Joss is Um 1 n 1 vu je is Um n 1 nN i aon ido 50 3 1232 e w ant
39. there are unknown parameters To determine the gradient of the cost function expression 2 2 is differentiated with respect to one correction factor for an uncertain parameter All discrete model variables are dependent on these uncertain parameters The model state is a difficult implicit function of the correction factors Therefore the chain rule must be applied J _ k gk 06r 2 m ul OF 06 F u OF Ov OF Ee E 2 3 mak Op Op op aq F Ou OF ane TL p OU Ep OVdp op Y XE OF OF OF Ou OF OF ar OuoOp OvOp Op This expression still contains unknown terms for the gradient namely all differentials of the water levels and the velocities with respect to an incremental change in the control parameter For example a local change in the bottom roughness also affects the tidal flow in other areas In principle changes in one parameter influence the water level Technical documentation of WAQAD and the velocity current throughout the entire domain of the flow model These influences are a priori unknown Looking more closely at the terms in this expression of the gradient 2 3 it can be seen that The differentials of the model state at each location for every time step with respect to the parameter to be estimated are unknown These terms are u ov Op Op Op The differentials of the ter
40. value instead of solving the momentum equation or instead of solving the continuity equation in the case of a water level boundary condition Terms with state variables at interior points do not appear in the linearizations of the WAQUA boundary conditions In principle adjoint boundary variables do not appear in the adjoint model equations for the interior state variables either This is equivalent to solving the adjoint model equations with homogeneous WAQAD boundary conditions Only when parameters in the WAQUA boundary conditions are estimated is it necessary to determine the contribution to the gradient for the adjoint state variables at the boundary This can be interpreted as a WAQAD boundary condition The WAQAD boundary conditions are solved after solving the adjoint model equations for the interior state variables Drying and flooding of tidal flats in WAQUA can be interpreted as the movement of the closed boundary conditions This is expressed in WAQAD by the movement of the corresponding homogeneous boundary conditions Finally it must be remembered that the adjoint equations are solved backwards in time when implementing the adjoint drying and flooding procedures The values of the WAQUA state variables are therefore known at every time step To discover the times of drying in WAQUA it is not necessary to check the local water depths with the drying criteria described A simple test to check if the velocity at a certain point of t
41. velocities in either the x or y direction Due to the drying and flooding procedure this geometric information varies with time We are going to examine the influence of impeded geometry on the numerical treatment of the advective term uu in the computational procedures of WAQUA Firstly we consider the alterations in the adjoint u impulse equation associated with Aur Degenerations in the adjoint u impulse equation associated with AU caused by impeded geometry An array GEOUUX is constructed from which the time dependent local flow situation can be determined In WAQUA impeded u velocity points are denoted by KHU lt 0 0 for permanent dry and 1 for temporary dry These variables are part of the solution flow and may change from time to time KHU m 1 n KHU m n KHU m 1 n GEOUX m n Table Definition of array GEOUUX The situation GEOUUX 1 is already handled in section 3 2 2 This is the case were the surrounding points are wet Than a central difference is used to approximate u otherwise a first order upwind scheme or zero This information is important for the construction of the adjoint 174 Appendices advective terms since there is a one to one relationship between WAQUA and its inverse model GEOUUX m n 2 In this case the term uu is discretized at the first stage as Um n 7 Um 1 n Umn AX JA 7 1 Introducing Lagrange multipliers and applying variational formalism y
42. we examine the u momentum equations solved in the second stage of the ADI method at the locations m n m 1 n m 1 n 1 m n 1 these equations contain factors with the differential Av after introducing Lagrange multipliers and applying variational formalism Assembling yields l 1 AV mn e vu Jr 4S Umn vu dr m 1n ig u m 1n 3 3 68 1 vu gS Unin intl 4 ig u mnn The expression between braces contributes to the adjoint v momentum equation which is solved at the first stage It only affects the right hand side of the equation not the matrix structure The appearance of the differentials Au in equation 3 67 also implies contributions to the adjoint u momentum equation solved in the second stage Since the discretisation of uy is dependent on the flow direction of the v velocity two situations have to be distinguished at each location After differentiating we focus on the equations at the locations which contain factors with the differential Aur At location ur n 2 in this stage vuy is discretized as V m n 2 Sy U m2 Consequently we find Vs n 2 E Au aes if yv gt 0 nn 2 Y 0 7 3 69 The 2D adjoint model At location m n 1 in this stage vuy is discretized as Vm net S Uns Consequently we find aut oS if FO 2AY man 0 if go 20 3 70 e At location m n in this stage vuy is discretized as V m n S Umn Consequently w
43. 1 0 October 25 1996 85 Technical documentation of WAQAD 86 Without the advective terms the existing adjoint v momentum equation at this stage looks like lache quu quis Ee sa Gidea load eso U bat ve rmt vy Yan 170 AY V6 Jm n VV6 Jm n mnl jAVv m n Vv Jmn Adding 3 106 to this existing adjoint v momentum equation results in the same matrix structure in which H m n 0 2 l y AY a at m n 1 Vm n 1 E iiis 2 Y 2 a J Vial Ap tn 2A Y Vv Jing DDm n 1 2 frvv Vinni At ian 2 Y At ZA 2 Vin wet Vin n AY refe m n 1 m n 1 The 2D adjoint model For the sake of completeness here is the adjoint continuity equation associated with AG n ve ee E Jo Semn T Jl 7 vv Jai 7 Vv Jmn Van Vv m n 1 1 2Vmn i Y Y a ASTU mat ve Jmn 7 eme ve 3 152 4 2 Vm 1 1 1 A de AY t vc Jona 7 Jmn vc vc pm mmu 00 ES In 3 152 three adjoint water levels and two adjoint v velocities at the new time level occur Together with this adjoint v momentum equation they form a pentadiagonal matrix equation if for all grid points in a column of the domain the corresponding equations are written down v2 Ve v Vent Vy v Wnt Ve ve v Ven v 3 153 ve vh ve vy Vener Ve v
44. 135 Technical documentation of WAQAD 136 More precisely If p is the parameter vector and y t p are the model results computed with these parameters at time t and y t are the measured values then the quadratic cost criterion was J p 0 5 X y t p y t 5 4 teT On the other hand if we assume that the measurement errors of y t are Gaussian and identically distributed with mean 0 and standard deviation o the Maximum Likelihood estimator can be written as 1 y t p yo t L p e x 5 5 tet 2 Taking minus the logarithm of this expression reduces it to the old criterion for all c Bayesian estimation is a logical extension of Maximum Likelihood estimation that is based on Bayes s rule which reads P A B P B A P A P B 5 6 This can be seen easily from the definition of conditional probability P A B PAOB P B 5 7 if we apply this to our estimation problem we get Pyt p y t P p p te D TT PEP y P b p dp 5 8 Version 1 0 October 25 1996 Possibilities options and shortcomings if for the prior distribution P p p we assume that it is Gaussian distributed with mean 0 and covariance C and note that the denominator does not depend on p By also taking the minus logarithm we can obtain J p 5 y t p y t p C p 2 teT o 5 9 This looks very much like the Maximum Likelihood estimator The only difference is the penalt
45. 29 5 3 The separate transit times and with them the depths are not to be estimated Each combination of depths with the same total transit time has the same value of the cost function which means that there is an I problem If the phase at the boundary does not have to be estimated it is not necessary for Xm to be at the boundary This means that this problem can also occur at another place in the model I problem 3 harmonic constituents that are difficult to separate It is a known fact that harmonic constituents with frequencies which are close together can only be separated using long series of measurements This is still a problem when WAQAD is used It is regarded here as an l problem See section 5 3 for more information on solving this problem I problem 4 more boundary points than stations It is possible to prove that estimation of the constituents at the boundary leads to an l problem if there are more boundary points than measurement points However the proof is very technical mathematically and has therefore been omitted Splitting the calculation The number of parameters to be estimated is sometimes so very large that is is tempting to break down the calibration process into a number of steps The resulting separate problems are smaller and therefore require less computational time For example one could first calibrate the constituents of the open boundary and next the bottom depth and friction Note however tha
46. 2AX v Jinin 2AX Remark the marked adjoint variables in these equations are the ones at the new time level e Inthe adjoint v momentum equation associated with 4 M vc y E MERO VN 3 101 2AY l 4u TU mn vv m 2n E Jn 4n 12AX _ 3 102 4U m n i U m 2 n vv Join 12AX Jaian e vu Jnn on Una 2 roca ey CUL in 3 103 vu Jr m 1 n 1 15 Ur yu m n 1 15 Unn Um n bond Um 1 n i eh v i de 1 8 104 vu 194 Um 1 n 1 T vu aoe aS Um n 1 t vu 2a mee if Um 2 n lt 0 3 105a 54 The 2D adjoint model vv inant E 3 105b PE Ea uo aud 2AX or 3 105c HY AK 3Un n vv mni if gt 0 v Jnnt 2AX Um n Yan ELT if 1 70 3 105d wee 2AX m 1 n T n i vv nan if Um 2 n 0 3 105e Inthe adjoint v momentum equation associated with Ning Inthe adjoint u momentum equation associated with A Um n Jmn E Uc oral vu Jost 2AX i Um 1 Umi Dr A N usa Vu Jmm sou Imetm 3 107 2AX dn chm T m 1 ni Vin n AV m n t vu Jans DIY vu ee EY AV m net Vm m2 v Juin J sio 12A Y 3 108 Version 1 0 October 25 1996 55 Technical documentation of WAQAD 56 v nin 49 Vinn T vv Jos S Mi t v 164 sint v s iol
47. 49 3 50 and 3 51 yields Version 1 0 October 25 1996 The 2D adjoint model T a ounce 2AY 3 52 2AY Vinn 7 Misi T Jmn 2AY 3 53 v Vn n 1 ead E ne 2 Y where denotes the state variable at time t 4At The expressions between braces are part of the adjoint v momentum equations which are solved in the computing procedures in the first and second stage 37 Technical documentation of WAQAD VUy This cross advective term occurs in the conservation of momentum in x direction In the first stage of the ADI method uy is approximated by a weighted central difference explicit In the second stage another difference scheme is used namely a second order upwind difference implicit The phase and amplitude errors of the combined approximation are discussed in Stelling 1984 At this first stage in the u momentum equation at location m n the term is discretized as follows Vno S Un 3 54 in which l1 Vm n Vmn i Vm n Vm 1 n Vim 1 n 1 Vim n 1 Um n 2 4 Um n7 4Un n17 Um n 2 Sh Um n 12AY Applying variational formalism yields V mn AV mn A Um n2 t Um n 1 77 3T 12A Y 12 Y 8 55 A nin V m n A m n 12AY 4 42A Y A Um n 1 iS Um t AM natn aS Um n tac 195 Um n AV mnn iS Um n Equal to Van S AUm n AV mn S Cnm In this express
48. 5 16 This equation can also be written in matrix notation AQ GA 5 17 Where Q is a column vector that contains the Q x y the columns of G contain the weighting functions and A is a column vector with all the parameters Several types of weighting functions have been suggested in the literature such as point functions rectangles triangles splines and other functions fitted to the problem A description of all the weighting functions mentioned is given in Verlaan 1994 Rectangles and triangles are implemented in the WAQAD system These weighting functions are described below figure 5 11 Rectangles e Rectangles Every weighting function is one of a set that is composed of rectangles and is zero outside this set These are the weighting functions that were initially used in the WAQAD system Their great advantage is that they are easy to compute and many algorithms that use the weighting functions can take advantage of the special structure Their disadvantage is that at certain edges the 141 Technical documentation of WAQAD 5 2 2 142 adaptations jump from one value to another These discontinuities may cause problems since waves are partially reflected at these edges These reflections are an artefact of the model They do not exist in reality The discontinuities are also very different from the physical situation where there are no jumps along straight lines JLRS LON ELLIO SOOO LISS SSN LI
49. 5 4 1 Introduction The gradient of the criterion with respect to the parameters is an important element of the procedure for estimating parameters using a cost criterion and the adjoint model In fact the main purpose of the adjoint model is to compute the gradient of this criterion This gradient gives information on the direction in which the parameters should be adapted to decrease the value of the criterion which in turn brings about an improvement of the model This is in fact how the gradient is used for parameter estimation Another interpretation of the gradient of the criterion is the local sensitivity of the quality of the model to the corresponding parameters If a large element of the gradient is associated with a certain parameter this indicates that a small change in this parameter has a large effect on the model If this parameter represents an area of the depth parameter in the model the model is sensitive to a change of the depth in that area This enables different areas to be compared and a simple map to be made of the sensitivity over the spatial domain One of the purposes of the programs described in this section is to provide an automatic procedure for making these maps As a result of the chain rule of derivatives the effect of combining two areas parameters into one is the weighted using the area average of the two elements of the gradient If we have a map in which the areas consist of points the gradient for
50. IMONA environment to these standards in combination with the complexity of the model is quite a job The best option is to employ the Chavent method on the discretized model equations to obtain the discrete gradient i e the exact gradient of the discretized criterion To understand the mathematics of the adjoint model it is necessary to understand the hydrodynamics of the forward model The basic concepts of the latter reappear in the mathematical description of the inverse model The depth averaged shallow water equations are described in Cartesian coordinates by Version 1 0 October 25 1996 15 Technical documentation of WAQAD The u momentum equation which describes the conservation of momentum in x direction ha2 2 fag erg xl Ot Ox oy OX C H 3 1 _ pw Ox oy The v momentum equation which describes the conservation of momentum in y direction 2 2 ON ON gro V esu 05 aN ot Oy Ox Oy C H 3 2 ox oy The continuity equation which describes the conservation of mass 954 2 Hu 2 Hv 0 3 3 Ot Ox oy in which u v Depth averaged velocity components in the x and y directions Water elevation above plane of reference d Water depth below plane of reference H Total water depth d C g Acceleration due to gravity C Chezy coefficient for model bottom roughness Viscosity coefficient f Coriolis force t
51. KHV at time t At Table V GEOVVY V 0 f vss v Jani TARNA Table VII The degenerations in equation 3 101 182 Appendices For convenience we give again the adjoint v impulse equation associated with DN cn including the advective terms neglecting the boundary treatment and the dryfall procedure and the iterative aspect w Jansi Dmn 7 16 with Am n DX 7 17 Bnn 7 18 Oc 7 19 2 Ds lt Ww Yin finn UV inn vu Jan u Imn T 7 20 v5 Jus Vrasakvv Jet 2AY i AU mtn vy Vicon ed vy Drea 12AX AU ms n Us n vy Join NONE vy ee ae 4 g vu Jn a ur 4 Ur tn D Us iicet aie z Uhn 7 Jos 194 Um n vu dus 194 Um 1 n 7 eae 1s Um 1 n 1 v ASK Um n 1 Version 1 0 October 25 1996 183 Technical documentation of WAQAD 184 Incorporating the boundary treatment and the dryfall procedure yields still without the iterative aspect Amn vy EM Bain vv Js Cun vv Junii 7 Dm 7 25 p The coefficients Amn and Cmn can be determined from Table VI e One part of the righthand side D is prescribed by Table VII These elements have to be subtracted e The boundary treatment of the second operators S and S can be read from Table They appear in the formula for Dr e One part of t
52. LLIES ASSIS EE op IIL LI LISS S Lge SLL hives ESS SSS go a TESS SSS SSS SOON rr tots LLL i 7 2 SSS 50 EES NS EEE EES ESS 50 20 30 figure 5 12 Triangles Triangles where the weighting functions are pyramids The resulting adaptation is a linear interpolation of the values on the edges of a set of triangles in space This type of representation is well known for its use in finite element methods for partial differential equations With some caution the adaptations at the boundary are always continuous An additional advantage is that the geometry can be followed more easily using triangles than using rectangles In the next part of this section the extension of the WAQAD system for triangles is described Triangles have been chosen because the geometry is easy to follow and the algorithms are relatively simple Theoretical background in the case of triangles Linear interpolation between the three sides of a triangle is not difficult when baricentric coordinates are used If the three sides of a triangle are given by the Cartesian coordinates X1 1 X2 Y2 and X3 y3 then the baricentric coordinates for this triangle are given by the equations aX4 2 2 3X3 X MY F A2yo FY 5 18 o tis 71 5 2 3 5 2 4 Version 1 0 October 25 1996 Possibilities options and shortcomings where A1 A2 A3 are the baricentric coordinates In baricentric coor
53. Modelling R 3712 NETH RC Rand Corporation Santa Monica 1989 Mooij M N A Calibration of shallow water flow models using current velocities RIKZ werkdocument RIKZ OS 96 140X 1996 Mouthaan E E A Geadjungeerde droogval en onderstroom procedure voor WAQUA in Dutch RIKZ werkdocument RIKZ OS 94 160X 31 december 1992 Mouthaan E E A Geadjungeerd systeem voor rechtlijnig WAQUA in Dutch RIKZ werkdocument RIKZ OS 94 113X 1994 Mouthaan E E A Faserapport fase 2 ECAWOM in Dutch TU Delft faculteit TWI juli 1995 RIKZ TU Delft WL de Voorst Data assimilation with altimetry techniques used in the Continental Shelf Model BCRS Report 94 08 november 1994 Rozendaal D M J Karbaat and J R Brouwer Faserapport fase 2 Uitbreidingen van WAQAD Triangularisatie Componentensplitsing in Dutch SIMTECH oktober 1995 Rozendaal D M Implementatie en testen 3D adjoint in Dutch RIKZ werkdocument RIKZ OS 96 125X 30 mei 1996 Stelling G S On the construction of computational methods for shallow water flow problems Rijkswaterstaat Communications 35 Rijkswaterstaat The Hague The Netherlands 1984 Stelling G S A K Wiersma and J B T M Willemse Practical aspects of accurate tidal computations Journal of Hydraulic Engineering p 802 817 september 1986 Verlaan M Some extensions to the calibration program WAQAD TU Delft faculteit TWI Report 94 111 1994 Verlaan M Identificeerbaarhe
54. Mund vy inna S Vm n vy mrna Sel Vida t ere is Vm n 1 v Jana isi Vs n 1 a 3 L 5 l l gt J AET T T I AR zu 3 A I 3 5 o ces Ee TESI or sm micum ee pec ES Note that these two expressions are symmetrical with adjoint advective terms in the first stage of the adjoint 3 109 3 110 3 1112 3 111b 3 1116 3 111d 3 111e the computation procedures This is due to the ADI structure in WAQUA 3 2 2 4 Version 1 0 October 25 1996 The 2D adjoint model Curvilinear and spherical coordinates In this subsection the advective terms for curvilinear and spherical coordinates are indicated If the cross terms are ignored the advective terms can be derived from the advective terms for WAQUA in rectangular coordinates In this case it is only a matter of multiplying by coefficients The structure of the equations remains the same For curvilinear models Um n Um n Um 4 sii n7 Unt 112 UR KA AE PS Vm n 3 113 EL on 2An The second order difference operators 1 2 Us n7 Um n1 1 Um n 2 Um n 2 an m a 3 114 SUR zl 2An 4An a 3 Um n 2 ond i Um n 2 if Vm n gt 0 An 1 En 3 115 Sih Um n n 3 Um n 2 su seo if y nO An 2 1 S Vmin 1 2 Venet Vim 1 n Ju 1 Vio n Vm 2 n 3 116 Km
55. This allows the equations to still be resolved by solving tridiagonal matrices over grid lines because in the adjoint cross advective part the latest computed adjoint u velocities Gauss Seidel or the ones originating from the previous iteration level Jacobi are used The 2D adjoint model UUx The u momentum equation and the conservation of mass are simultaneously solved in the first stage of the ADI method Here the difference u is approximated explicitly The discretization of uu in this equation is as follows us C 3 75 After differentiating AU ae Mene 3 76 In the second stage of the ADI method the difference ux is approximated implicitly Here the discretization of uu is ETAN 7 Um 1n 3 77 2AX After differentiating u S u Au 2 Au Aut m 1 n m 1 n Jig ul m 1 n m 1 n 3 78 2AX 2AX Since the adjoint treatment is completely analogous to the term vvy only the most important formulas will be given We focus on those u momentum equations at different locations in which terms with Au and are found after applying the variational formalism and introducing the proper Lagrange multipliers for equations 3 76 and 3 78 In the first stage in the u momentum equation at location m 1 n we find us n BE FN poe 3 79 Version 1 0 October 25 1996 45 Technical documentation of WAQAD 46 e Inthe first stage in the u mome
56. a jt 2h h x nlo od w wn nne Y F amp x ox vu 106 Version 1 0 October 25 1996 The 3D adjoint model Adjoint conservation of momentum in x direction per layer in the first stage vu va h ve y 08 At K NES E Vv vu k 1 vu 1 1 h hk 1 h ry pn Suk ds Jag hth L h hc Adjoint conservation of momentum in y direction per layer in the first stage l r y fy g w A ml hi h PETN t _ tusks Jeg Z h hi h in which adjoint water level vu adjoint velocity in x direction for a layer adjoint velocity in y direction for a layer Remarks e Note that the vertical viscosity coefficient and the adjoint velocity in y direction are denoted by the same symbol However the adjoint velocity in y direction is marked by braces v e The adjoint equations in the second stage are virtually identical and can easily be deduced from the ones in the first stage by Stelling 1984 the interchange of amp and u and v the alteration of sign of the Coriolis coefficient the transformation of the time level with half a timestep e The adjoint equations for the simplified three dimensional shallow water equations are deduced with the Chavent method Since uncertain boundaries are not prescribed
57. and the System Documentation of WAQAD Brouwer 1996a The 2D adjoint model 3 2 The adjoint model equations This section focuses on the mathematical deduction of the adjoint model equations for depth averaged shallow water flow models The method used here to derive the adjoint equations is the Chavent method outlined in chapter 2 The adjoint model equations are derived in section 3 2 1 Sections 3 2 2 to 3 2 4 describe advection boundary conditions and dryfall respectively 3 2 1 Derivation of the adjoint model equations The mathematical description of water flow consists of a system of coupled differential equations which are physically based on the conservation laws for mass and momentum These shallow water equations are deduced from the well known Navier Stokes equations Accordingly the adjoint model is also mathematically described by differential equations Adjoint codes for mathematical models can be constructed in various ways One way is to deduce the adjoint equations in continuous form from the continuous model equations and apply a finite difference scheme to obtain the discretized adjoint equations and the discretized gradient Chavent advises against using this approach Another way is to apply adjoint model compilers which generate the adjoint code automatically Giering et al 1995 The great disadvantage of this is that the source code must be remoulded to the standards of the compiler Rewriting the S
58. assumed This property also holds for uncertain boundary conditions Finally the adjoint boundaries are determined by the adjoint solution for the inner points Using the adjoint technique has several advantages e The adjoint equations are always linear even though there may be nonlinearities in the model equations This is because the differential operator is linear e Computation is efficient since the computational effort for the determination of the gradient of the criterion is not influenced by the number of control parameters Integrating the adjoint equations backwards in time is more or less comparable with one forward model simulation e The inverse model is flexible in the sense that the number of parameters is easily changed It is even possible to estimate different spatially varying correction factors for bottom roughness and harmonic constituents simultaneously The only difference is the way the adjoint solution is combined with the original model output e With respect to the parameters to be estimated the Chavent method has two advantages over the conventional finite differences method for determining the gradient of the cost function it computes the exact gradient and it is efficient There is only one disadvantage of using the adjoint technique e The computed water levels and velocities of each grid point at each time step during the calibration period have to be stored for use in the adjoint system A shortage of d
59. ated such that the continuity equation is 1 C ees t 3 5 1 1 Ki E AB Hin U m nEmn B Hin U m tn En 1 n 1 H Vm n Km n H ni Vm n 1KKm n1 0 An The 2D adjoint model and the u momentum equation is Crain Cin a 1 E At N Urn 7 Umin fmnVmn 3 6 0 9 Km n in which CH CH Jas n g ust Y Gs AE C Este total water depth at u velocity point Ve Om nt Gmet n tdm nt Dm n1 total water depth at v velocity point Ve m n t Gm ntt dmt dm transformation coefficient K at water level point Ya Kmn itKm n transformation coefficient E at water level point Ya Em 1 n Emn transformation coefficient K at u velocity point VA Km neat Kms n 1t Km t n Kmn water depth in u velocity point Y2 dm ntdm n 1 chezy coefficient in u velocity point and where the following terms are interpolated to u velocity points n m n Gih mn Vm n n Version 1 0 October 25 1996 VA Vm nt Um n 1tUm 1 n 1 Um 1 n Ve m nt Gm Yalfm nV mnt fm n1 Vm n t hme nt Vm n 1 fm 1 nVm 1 nl 21 Technical documentation of WAQAD The u momentum equation at the second stage from t 72At gt t At This is the second step of the ADI method in u direction The WAQUA variables at time t 2At are known All u velocities at time t At are calculated such that 1 n u u fmn Van
60. ator depends on the flow pattern of the u velocity field gt p If 0 3Vin n n Wea f gt for situation GELUVX 2 van Waan T 7 31 S v oe for situation GELUVX 1 i AX 0 for situation GELUVX 0 Version 1 0 October 25 1996 187 Technical documentation of WAQAD M 3vi t 4v A Vir Mii eae 2AX Vinn Mt 7 32 Q amp Yhn semi AX 0 for situation GERUVX 0 Another part of the righthand side of the adjoint u impulse equation under consideration forms equation 3 108 Table IX shows the degenerations in case of impeded geometry The necessary arrays GEBVUY and GETVUY Table II have to be deduced from the flow field at time t 72At GEBVUY GETVUY situation Situation iiy iv Table IX The degenerations in equation 3 108 188 Appendices Then the part of the right side originating from the second order upwind difference in WAQUA Table X shows the deformations caused by impeded geometry The arrays GEBVUY and GETVUY Table have to be deduced from the flow situation at time t 7 At SEBY u point m n 2 V mn 2 gt 0 u point m n 1 Table X The degenerations in equation 3 108 The elements corresponding to u point m n contribute to the main diagonal of the matrix equation The coefficients have to be added to Bmn The other elements influence the right side since an iterative
61. be done is to save all the required information The situation when one 151 Technical documentation of WAQAD 152 wants to split the gradient in time is a little more intricate but the solution is almost the same and again all that needs to be done is to store the required information Filtering the data The raw data obtained by storing the information during the computations of the adjoint model have a regular part and some high frequency irregularities or noise A discrete time filter is used to obtain a smoother picture The only parameter needed for the filter is the cut off frequency Most of the noise has a frequency close to the one corresponding to the time between the data The input should be given as a fraction of the Nyquist frequency this means that if one is given as an input all frequencies are passed and output data is the same as input data A good indication of this fraction is the inverse of the number of simulation steps between two time steps with measurement data Let a k be a series of data then the Fourier transform is given by a o de a k and the inverse transform is given by T a k ze alodo x 21 If we cut off all the frequencies above then the Fourier transform of this filtered data b becomes o fis Op 0 o gt Possibilities options and shortcomings 5 4 3 In the time domain this transformation can be described by the convolution b k 2 g l a k
62. ccuracy may occur In order to be able to calculate a reliable gradient using the finite difference method some research must be performed on the shape of the parabola WAQAD offers the possibility of calculating the gradient by using either the adjoint method or the finite difference method The latter method has one advantage It can be used to check the gradient calculation of the adjoint method by checking the correctness of the implementation when new types of parameters are introduced such as vertical viscosity or salinity The finite difference method will mainly be used as control tool or to investigate if it is useful to calibrate a certain type of parameter Its use as a calibration tool is not advised for the following reasons The gradient cannot be calculated exactly During every iteration the shape of the cost function for each parameter must be searched in order to get an indication of the step size to be used The computational time required is dependent on the number of parameters In addition to an initial WAQUA run a WAQUA run for each parameter is required in order to calculate the gradient This in contrast with the adjoint method which only requires one single adjoint run in order to calculate the complete gradient vector For more information on using the finite difference method see Brouwer 1995b 163 Technical documentation of WAQAD 164 5 8 Velocity and discharge boundaries Until recen
63. ce Lagrange multipliers and apply variational formalism Vm n 1 m a d AV vv Joint 2AY m n 1 3 39 Vis V m n m n 2 A Vm n Vm n 1 Avi 2AY 2AY In 3 35 3 37 and 3 39 we focus on Av and since we are interested in the adjoint conservation of momentum in y direction at location m n Vm n 1 In 3 39 vy c7 3 40 Edid 3 40 34 The 2D adjoint model In 3 35 vy J A Vent 7 Vino 3 41 Vv Jmn Vm AY In 3 37 vv Vans zl en Jv z 3 42 A fragment of the adjoint v momentum equations can be determined using this as will be demonstrated at the end of this section In the second stage the v momentum equation is solved simultaneously with the conservation of mass In this second stage of the ADI method v in the conservation of momentum in y direction for location m n is approximated explicitly The discretization of vvy is Vm n 1 7 Vinn 1 Vinn 3 43 i 2AY Applying the variational formalism again Vina 7 V AV sar AN Av m n 1 m n 1 v m n 1 m n 1 3 44 oT 2A Y i 2AY e Inthe v momentum equation for location m n 1 in this second stage vv occurs as Vias E Vinn V 3 45 m n 1 2AY Introduce Lagrange multipliers and apply variational formalism Vinn d weg o AW oe wy Jost 2A Y m n 3 46 Vii 7 Vm n Vin n 1 AVE a AV ott 2AY 2
64. characterized by a set of weighting functions 5 2 gi x y The adaptations made to the model can always be written as a linear combination of these functions AQ x y ug X y 5 11 If the gi form a linear set of basic functions over all x y then there is a linear one to one correspondence between A and the AQ We can rewrite the above operation as a matrix multiplication with AQ denoting a vector that stores all AQ x y and containing all Ai AQ Gi 5 12 The covariance of the can now be calculated by E AX GE AQAQ G 5 13 Possibilities options and shortcomings 5 1 6 3 Since we usually only want to use a few functions g we do not have a linear basis and thus we cannot invert G But we can estimate from AQ using the least squares method GG G AQ 5 14 Using this formula we estimate the variance of X from the variance in AQ E 44 GG GE AQAQ G G G 5 15 This is the method used in the program to compute the covariance matrix needed for the penalty The symmetry of covariance matrices can be exploited to reduce the time needed for computation The algorithm becomes slightly more complicated due to the indirect storage of the WAQUA in SIMONA and WAQAD programs The indirect specification of the covariance matrix for the penalty has some important advantages over a direct specification The parameters that have to be given by the user have a clear interpretation and the pena
65. cs of the model must not be violated The constrained optimization problem is transformed into an unconstrained one by introducing Lagrange multipliers Bryson and Ho 1975 and Chavent 1980 J6 2 amp op amp P Xe 2 2 t F C u V p b R 6 u v p mn mn k The functions state the conservation of mass and Fy state the conservation of momentum eastward and northward These functions describe the physical processes of the Version 1 0 October 25 1996 The adjoint method model the propagation in time of the model state They consists of finite difference methods for the continuous model equations Stelling 1984 Additional terms are appended to the cost function each of which is the product of an undetermined multiplier and a model equation Each model equation for each discrete model variable at each grid point for every time step must be associated with its own multiplier the adjoint variable and contribute one term To simulate the flow the implicitly coupled model equations have to be solved for each grid point for each time level in the simulation period Any value for an adjoint variable is allowed in 2 2 since the model state satisfies the model equations More precisely after the model simulation the corresponding model equation holds at each grid point for every time step F 70 F 0 9 0 Note that there are exactly the same number of equations as
66. ctive term uu is discretized in this first stage as Usu Um 4n 3 25 Um n i 2AX The difference u is approximated explicitly This means using the u velocities at the old time level which were determined in an earlier stage of the computation process Consequently the structure of equation 3 24 augmented with 3 25 remains the same only the coefficients will alter The cross advective term vu is discretized in this stage as V m n S Um n 3 26 in which Ie 6 m Vm n V mn Vm 1 n Vm 1 n 14 Vm ni S Cun m Um n 2 7 4 Um n 17 Um n 2 12AY The difference uy is approximated explicitly with a weighted central difference For a discussion of the phase and amplitude errors of this approximation see Stelling 1984 In the preceding procedure these v velocities were computed and due to the explicit approximation of uy the structure of equation 3 24 remains unaltered when equation 3 26 is added With the exception of the boundary treatment incorporating the advective terms at this stage of the computation process leaves the structure of the equations unimpaired Version 1 0 October 25 1996 The 2D adjoint model For each row in the domain of the model all the u momentum equations and the continuity equations form a pentadiagonal matrix equation together with the boundary conditions 6 Um 2 6 u Caer U Um 6 u Cm U 6 6 Um 6G 6 u Cm U 6 6 Umer 3 27
67. dent defined for eddy viscosity in the second stage is described by A va U V u s Zunch n con irs hy h h h kat v n v Vit 4 5 Version 1 0 October 25 1996 The 3D adjoint model Remarks on WAQAD with respect to TRIWAQ The remarks on WAQAD with respect to TRIWAQ are The 3D WAQAD bottom friction terms fruu3D fruv3D frvu3D and frvv3D described at the end of section 4 3 contain a shift of half a time level with respect to the bottom friction terms in TRIWAQ In the term frwl3D the shift of half a timelevel with respect to TRIWAQ is corrected Horizontal and vertical advection is not included in WAQAD but it is in TRIWAQ Horizontal viscosity is not included in WAQAD but it is in TRIWAQ Just as in the case of WAQUA additional features of TRIWAQ are not included in WAQAD This may lead to the same problems described in section 3 5 The features of concern are External forcings like wind discharges and density gradient Transports STRESS2d barriers and weirs Smoothing boundary values Depth multiplier In TRIWAQ there are two ways to prescribe layers fixed layers or sigma layers WAQAD assumes that no fixed layers are prescribed in the forward model Just as WAQUA TRIWAQ writes map data only at every whole time step and not every half time step In WAQAD the map data at whole time step
68. dinates linear interpolation can be written as Q x y Q x Y1 2 2 As Q X4 y 5 19 Using the set of linear equations 5 18 and formula 5 19 Q can be evaluated at a certain point these evaluations can then be used to construct the weighting functions When the appropriate routines are rewritten in terms of weighting functions the transfer to triangles is straightforward and is achieved by changing the weighting functions Computation of the gradient in the case of triangles When a different set of weighting functions is used then the computation of the gradient needs to be changed Let Q x y be a parameter that depends on spatial coordinates x y and gi x y i 1 m a set of weighting functions for the adaptations in the parameters The relation between the adaptation of Q x y and the coefficients A of the weighting functions gi x y is given by AQ x y g y 5 20 If the gradient of the criterion with respect to Q x y is known then the gradient with respect to the coefficients of the weighting functions can be computed by ad _ 8 a 2 9 aK y 5 21 So in contrast to the use of rectangles the summation of the gradient J 0Q x y which was over the rectangle is now over the whole computational grid and there is an additional weighting function However in both cases the computation can be done using weighting functions Implementation in the case of triangles For the computatio
69. e a v The adjoint v contributions in the adjoint continuity equation are eliminated by means of the adjoint v momentum equation with a substitution The result is a tridiagonal matrix equation for each column consisting solely of adjoint G contributions ve ves Ve Vo Us 3 154 This matrix equation is solved with Gaussian elimination double sweep method Stelling 1984 The adjoint v Version 1 0 October 25 1996 87 Technical documentation of WAQAD 88 momentum equation can then be solved by simple substitution AD UXD In procedure AD UXD the adjoint u momentum equation is solved to determine v m The existing adjoint u momentum equation at this second stage of the adjoint computation process reads 1At 2 3 155 Without the advective terms diagonal matrix equations over a grid line in x direction are obtained Bs Vu Jmn Din 3 156 in which Bn 2 At 2 T m Dm n At vu Jan fmn 2 frvu m n vy Jah vy ae Adding 3 107 3 110 omitting 3 111 to this existing adjoint u momentum equation yields a tridiagonal matrix equation for a row in the domain of the model vu 1 vs v 3 157 vu v l Version 1 0 October 25 1996 The 2D adjoint model This banded system could be solved over the gridlines in x direction by standard procedures The corresponding coefficients of this matrix equation
70. e case of 3D models e The horizontal advection terms are solved in WAQUA using the Gauss Seidel method In WAQAD these horizontal advection terms are solved using a Jacobi method e Additional WAQUA and TRIWAQ features are not included in WAQAD The features of concern are External forcings like wind discharges and density gradient Transports STRESS2d barriers and weirs Smoothing boundary values Depth multiplier In TRIWAQ there are two ways to prescribe layers fixed layers or sigma layers WAQAD assumes that no fixed layers are prescribed in the forward model In the past more differences existed between WAQAD and WAQUA For example horizontal advection was not a feature of WAQAD But when the horizontal advection was also included in WAQAD the calibration process really improved It is expected that WAQAD will perform better if other missing features are assimilated too Therefore it is recommended that WAQAD is fully consistent with WAQUA and TRIWAQ Furthermore it is recommended that WAQAD is extended to the use of water velocity measurements This has already been tested by means of a special test version of WAQAD see Mooij 1996 The results are promising 165 Technical documentation of WAQAD 166 WAQAD is extended to estimate the standard deviations of the measurements In theory it is possible to use the Maximum Likelihood method to determine the optimal standard deviations of th
71. e it is important for the interval used for calibration to be a good reflection of the interval that will be used with the model for prediction If the model is to be used to predict tidal maxima and minima for example and there is in spring tide and neap tide fluctuation in the gradient then the interval for calibration should be a month or a multiple of a month These effects are illustrated with some examples in Verlaan 1994 5 4 5 Use of the Hessian Together with the normal output of the program the Hessian and its eigenvalues are given see subsection 5 1 4 The Hessian the matrix of all second derivatives gives us information on the shape of the surface of the criterion at the local minimum From maximum likelihood theory it can be shown that the inverse of the Hessian converges to the a posteriori covariance of the parameters and thus the inverse Hessian can be used to give an error bound to the estimate The estimate of the inverse Hessian should be used with care Only when the estimate converges to the true parameter is it certain that the covariance is correct If this is the case the error of the parameters will be asymptotically normal and all kinds of tests may then be used 154 5 4 6 Version 1 0 October 25 1996 Possibilities options and shortcomings The estimate of the inverse Hessian is generated by the second order algorithm that is used in the minimization This algorithm only gives an accurat
72. e Navier Stokes equations The simplified i e advection and horizontal viscosity are ignored three dimensional hydrodynamic equations for continuity and motion may be written in Cartesian coordinates as De Goede 1992 B So yo at OX OZ OV _ OG O OV MT g VW ot Oy OZ OZ uez voz ot OR et in which u v Velocity components in x and y directions G Water elevation above plane of reference f Coriolis coefficient g Acceleration of gravity Ny Vertical viscosity coefficient The two stage time splitting method ADI as used in the 2D model is used to integrate the model equations of the three dimensional shallow water flow model However a slightly different approach is used for the vertical viscosity terms they are always taken implicitly in time Because the coefficients vary in time in the vertical terms the unconditional stability cannot be reached with central time derivatives Lander et al 1996 The discretized equations in space and time for the simplified three dimensional hydrodynamic model are the starting point for the deduction of the corresponding adjoint equations Since the two stages of a time step are virtually identical only 103 Technical documentation of WAQAD 104 one stage will be considered For an interior point in the field not enclosed by impeded geometry the flow is described by Conservation of momentum per layer in y d
73. e estimate in the limit In practice the result is often close to the limiting result after a number of iterations equal to the number of parameters plus one This result implies that in order to use the inverse Hessian estimate for error bounds on the estimate of the parameters the number of iterations should be large enough for the convergence of this estimate Another and probably more important use of the inverse Hessian is to detect problems in the estimation problem When the problem specified by the user has almost the same value of the criterion for some parameter in neighborhood it is not possible to give accurate estimates for the parameters see also subsection 5 1 4 Conclusions The gradient of the criterion used in WAQAD for estimating parameters can be split either in time or in space The resulting functions can be used to analyse the model They can also be useful in selecting the parameters for the input of WAQAD but additional experience is needed to make this a practical and reliable procedure Programs ADPIC and VIZ were designed and implemented to produce plots of these functions It is currently possible to make plots of the gradients of bottom friction depth vertical viscosity and boundary components for every model that runs under WAQAD for SIMONA 155 Technical documentation of WAQAD 5 5 The use of current velocity measurements Velocity measurements can be of great use for the calibration of models If
74. e find Au SVmn if F 0 DAY m n 3V mn ig Aur lt 0 3 71 Um n 2AY I Vm n Atlocation m n 1 in this stage vuy is discretized as Vs ni S Us Consequently we find 0 if Van 0 AV if Vin nt lt 0 3 72 2AY iu At location m n 2 in this stage vuy is discretized as V m n2 Sy Uinn 2 Consequently we find 0 if V mn27 0 n v m n 2 MD Au if v lt 0 3 73 m n 2AY if V m n 2 Since the computation procedure in the second stage performs the computation of to make Au vanish the expressions above have to be translated in time this means that they can be derived from the same expressions in the previous stage Version 1 0 October 25 1996 43 Technical documentation of WAQAD 44 It follows from equations 3 69 to 3 73 that a part of the adjoint equation associated with A um n looks like ae d eden Day f Yona 3 74 n T E Jann E if Vm ntl d 0 TA n y i e Vu m n JAY f vs or n Ms 2 if Vm n EO if Vain O Vm n DT ms if Vm n2 0 These adjoint u variables at time t are perpendicular to the ones in in the computation procedure in the second stage which set up the tridiagonal structure if the advective terms were ignored Therefore in imitation of WAQUA an iterative solution scheme is proposed a Gauss Seidel or Jacobi method
75. e for the initial cost function e Calculate the gradient of the cost function e Use a Quasi Newtonian method in combination with a line minimization to modify the correction factors which decrease the cost function e Stop if you are satisfied or start again if you are not satisfied with the result It is well known that the minimum of a function is characterized by a zero gradient in all directions The adjoint technique can be applied to determine the gradient of the criterion efficiently Unlike other assimilation methods it has the advantage that the optimal estimate is consistent with the underlying physics of the model The gradient can also be determined by another technique called finite difference see section 5 7 However using this method the gradient cannot be calculated exactly Therefore the use of the finite difference method is restricted and not described in this section In the next section the adjoint method is illustrated by means of the mathematical description Technical documentation of WAQAD 2 2 The mathematical description When calibrating a shallow water flow model one has to find the values for the uncertain parameters in the hydrodynamic model which minimize the cost function As mentioned above the lower the cost function the better the agreement between computed and observed water levels The water levels computed by the flow model are implicit functions of the uncertain parameters The grad
76. e in space receive a penalty The result of this last type is a more regular distribution of the parameters in space A combination of these two types of penalties usually creates a better conditioned problem and solves the problem of unrealistic adaptations In the case of the first example the cost hardly changes from the adaptations at the Danish coast while the penalty becomes large which results in a large total cost cost plus penalty The minimum total cost will now be attained with smaller adaptations In the second example large adaptations that give the same cost will now be given a large total cost and thus these adaptations will not be made by the algorithm This shows on 5 1 6 2 Version 1 0 October 25 1996 Possibilities options and shortcomings one hand that there will be no unrealistic adaptations but on the other that there still is not enough data to make the correct adaptations The two examples above were of course unrealistic But problems like this often occur in less obvious forms Often but not always they occur when we want to estimate too many parameters with too few data The second motivation for using a penalty can be that this enables the user to select which parameters are known well or which parameters are known with less accuracy Often one knows from the way in which the data are obtained which parameters in the model are known accurately and which are not For example depths can be measu
77. e measurements The ability to estimate these standard deviations automatically would also reduce the number of parameters that the user has to specify However this estimation is computationally somewhat different from the other parameters but not trivial the derivation is much more complex Experiments will be needed to determine the usefulness of this procedure since it can become instable To implement the possibility of producing warnings if the estimation procedure does not work properly e g if the problem is not identifiable This information can be obtained from the Hessian the standard deviations and the value of the criterion The automatic detection of these kinds of problems could save many hours of waiting for a result that then shows that the input was wrong Appendices 7 Appendices 7 1 Documenthistory June 3 1996 October 25 1996 First edition Version 1 0 October 25 1996 167 Technical documentation of WAQAD 7 2 References Boogaard H F P van den Inregelen van wiskundige modellen op basis van besturingstheorie in Dutch WL rapport 62 10 Z107 104 augustus 1988 Brouwer J R In Dutch Faserapport fase 1 Uitbreidingen in WAQAD in SIMONA Nieuw iteratief proces Problemen met kromlijnige modellen Satelliet data RSDS uitvoer in Matlab formaat in Dutch SIMTECH 1995a Brouwer J R Faserapport fase 1a Uitbreidingen van WAQAD Halve grid verschuivingen frictie termen Snelheids e
78. e selection of the period for instance the characteristic times of the model are important too It is important to note that the gradient only gives local information When the model is changed eg due to WAQAD the gradient can be very different The computation of the gradient split in time or space is explained below However it should be stressed that this procedure is largely intended as a research tool because not all implications have been worked out and neither can a clear cut algorithm for its use be given Theoretical Background Chain rule The computation of the gradient split in space or time is based on the chain rule for derivatives This rule is given by GOD Mn V 29 09 5 23 Ox O y OX evaluated for y gi x Let Q x y be a parameter that depends on the spatial variable and parameters 2 that represent the combined effect over an area This effect can be represented by the weighting functions gi x y The relation between these functions is AQ x y Xig X y 5 24 The gradient of the criterion J A 1 m can now be written as OJ OX OQ X y g Xy 5 25 For the rectangular or triangular method used in WAQAD gi x y is O or 1 and the gradient is just the summation over the rectangle If one looks at the way the gradient is computed with the adjoint method the summation 5 25 already forms a natural part of the computation All that needs to
79. e variables Km n and Emn are determined by Kmn gVV 74 g4 the WAQUA transformation coefficient in amp direction Enn guu Guan the WAQUA transformation coefficient in direction Time dependent variables In the discretized equations state variables at different time levels are denoted as follows the state variables at time level t 2At are denoted by a single apostrophe The state variables at time level ttAt are denoted by two apostrophes If there is no apostrophe the state at time level t is meant Advection and viscosity The advection and viscosity terms have been ignored in the discretized equations given in this subsection The horizontal advection terms of the discretized equations which are taken into account in the 2D adjoint model are described in subsection 3 2 2 The viscosity terms are not taken into account in the 2D adjoint model Changes arising from the implementation in WAQAD The discretized equations in this section were described in Mouthaan 1994 However some modifications regarding the formulation of bottom friction were made during the implementation of the equations in WAQAD and these have been incorporated in this documentation the equations described in this documentation are according to the WAQAD implementations Let us now perform the three steps The 2D adjoint model Step 1 The WAQUA equations n axis H water level point velocity point
80. ecified which causes the depth values to be modified When calibrating depth parameters this option will definitely lead to erroneous results WAQAD has no knowledge of this multiplier and will deliver WAQUA adapted depth values as if the depth multiplier is 1 WAQUA writes map data only at every whole time step and not at every half time step In WAQAD the map data at whole time steps is interpolated to half time steps Version 1 0 October 25 1996 The 3D adjoint model The 3D adjoint model This chapter goes into the adjoint model for three dimensional shallow water flow in detail The adjoint model equations are derived in section 4 4 The adjoint solver as it is implemented in WAQAD is discussed in section 4 3 The derivation of the gradient is given in section 4 4 An overview of the remarks on WAQAD with respect to TRIWAQ is shown in section 4 5 101 Technical documentation of WAQAD 4 1 Introduction One major shortcoming of the two dimensional shallow water flow model WAQUA is the fundamental loss of information about vertical structures of velocities This property is inherent in the 2D depth averaged shallow water equations The first three dimensional model used at Rijkswaterstaat was proposed by Leendertse 1989 It used a rectangular grid in both horizontal and vertical planes Stelling also proposed a three dimensional model based on so called sigma coordinates In 1992 Rijkswaterstaat RIKZ built a three dime
81. edures is described below The procedures are described for the U part of the ADI method applied to the interior points of the model At the beginning of the first stage of the ADI method in the u direction first the local water depth is calculated at every velocity point Um n The local waterdepth in a u velocity point is defined by pA Che oi Qna dmn The flooding procedure If the local waterdepth exceeds the threshold value VAR which is given in the WAQUA user input and if velocity point Umn was not included in the calculation for the 68 Version 1 0 October 25 1996 The 2D adjoint model previous time step because of dryfall that velocity point will now be involved in the calculation After this check the continuity equation and the u momentum equation are solved for this row After this row has been completed the drying procedure at velocity points is activated In this first stage of the ADI method in the u direction the drying procedure only checks the u velocity points In addition it is possible to active drying procedures that check water level points The drying procedure in velocity point Um n If the local water depth at the velocity point is less than 72 VAR the point Um is no longer involved in the calculation The dryfall threshold is half the flooding threshold This prevents points from becoming active for one time step and inactive again for the next e Optional drying procedures at wat
82. en the adjoint u momentum equation is solved in module AD UXD An overview of the WAQUA and WAQAD computation procedures is given below The WAQUA procedures VYD computes v SUV computes UXD computes u SVU computes The WAQAD procedures AD SUV computes AD VYD computes AD SVU computes AD UXD computes to make AC Au vanish vy to make AC Av vanish to make Au vanish vv to make Av vanish vu Vu 72 The 2D adjoint model Version 1 0 October 25 1996 In the next part of this subsection the solution method in the WAQAD modules AD SUV AD VYD AD SVU and AD UXD is described AD SUV The adjoint u momentum equation and the adjoint continuity equation in the first stage from t At gt t AAt are implicitly coupled equations The adjoint u velocity at time t is expressed as a function of the unknown adjoint water levels at time t 72At Comer E vna vu nm DDm CCm Kin t n En Ean The coefficients CCm 1 and DDm are determined by FEE UA AE fruur DEN 1 fruu v 1 j m n JA Vu m n DDm 1 AAt o o 1 At EE Enn vena ven Ag 1 fru Uhn Kan Eis Khiri n En 45 At The adjoint u velocity at time t 72At is eliminated from the adjoint continuity equation This results in a tridiagonal system of unknown adjoint water levels at time t 72 t for each row Am ve
83. ent is ignored throughout this section See appendix A for information on boundary conditions and dryfall with advection Subsection 3 2 2 3 gives a derivation of the adjoint advection terms In that subsection and subsection 3 2 2 2 the advective terms are described using rectangular coordinates This is followed in subsection 3 2 2 4 by an overview of the differences between rectangular curvilinear and spherical coordinates Advective terms in WAQUA In this subsection the computation procedures for solving the u and v momentum equations in the first stage are described The computation procedures for solving the u and v momentum equations in the second stage consist of similar symmetrical expressions and have therefore been omitted from this subsection u momentum in the first stage 29 Technical documentation of WAQAD 30 In this stage the discretized continuity equation and the discretized u momentum equation are solved simultaneously based on the model state at time t and the results of the preceding procedure which yields the v velocities at time t AAt If the advective and the viscosity terms are omitted the resulting discretized u momentum equation has the following structure see formula 3 6 1 X At U m n Um n fmn V m n 3 24 g Eng E m n g TEN TUR i T v j 0 AG Dec The conservation of momentum in x direction contains three unknown variables The adve
84. er level point Cnn In WAQUA the user has three options for defining local water depths at water level points The user also has the option to skip this procedure The definitions of the local water depths in water level points are 1 maximum criterion Cmn tJ MaX dm n dm 1 n 1 max dai dia 2 averaged criterion Grae A dm 1 n dmn 1 3 minimum criterion Ginn 22 MIN dm ns dm 1 n 1 min daa s dm n 1 4 no check at water level points If the user chooses one of the first three options and if the corresponding local water depth at the water level point is less than 72 VAR the two adjacent velocity points uj en ug are no longer involved in the calculation If one or more points fall dry the continuity equation and the u momentum equation are solved again for this row with these velocity points as stationary points This means that during this step these velocity points are considered to be closed boundary points Normally the new velocities are calculated by means of the momentum equations In the case of dryfall a velocity blockade is placed in order to give the new velocity a zero value 69 Technical documentation of WAQAD 70 Before calculating the second stage of the ADI method in u direction the u velocity points are checked for dryfall This drying procedure has already been checked during the calculation process though Next the new u velocities are calculated according to the previou
85. est model performance can be achieved when the cost function reaches its minimum To evaluate the cost function a model simulation has to be run for a given value of the parameters to provide the computed water levels The value of the cost function varies with the values chosen for the mathematical parameters The cost function may be controlled by these uncertain parameters The gradient i e the direction in which the parameters have to be adjusted so that the cost function will decrease provides information crucial for searching efficiently for a more accurate estimate of the uncertain parameters The minimum of the cost function can now be found iteratively Ultimately Version 1 0 October 25 1996 The adjoint method the optimal estimate of the uncertain parameters is determined In other words the mathematical model is calibrated The following steps are important in off line data assimilation e Define the parameters in the mathematical model which have to be calibrated In the case of a 2D simulation these may be correction factors for the model values of the bottom friction harmonic constituents in the boundary conditions or depth In the case of a 3D simulation they may be correction factors for the model values of the bottom friction the depth or the vertical viscosity coefficient e Runa simulation with the original parameter values in the model and the correction factors initially set at zero This gives a valu
86. fall occurs in WAQUA in the first stage from t gt t 7AAt at a u velocity point then the corresponding u velocity point in WAQAD in the first stage from t At gt t 72At must fall dry This means that at this stage the system of equations must be solved considering this point as a closed WAQAD boundary condition Instead of solving the adjoint u momentum equation the adjoint u velocity is set to equal zero AD VYD If dryfall occurs in WAQUA in the first stage from t gt t 72At at a v velocity point then the corresponding v velocity point in WAQAD in the first stage from t At gt t 72At must fall dry In short the adjoint v velocity is set to equal zero 91 Technical documentation of WAQAD 92 AD SVU If dryfall occurs in WAQUA in the second stage from t 7 2At gt t at a v velocity point then the corresponding v velocity point in WAQAD in the second stage from 2 t must fall dry In short the adjoint v velocity is set to equal zero AD UXD If dryfall occurs in WAQUA in the second stage from t 7 2At gt t at a u velocity point then the corresponding u velocity point in WAQAD in the second stage from t 72At t must fall dry In short the adjoint u velocity is set to equal zero The 2D adjoint model 3 4 3 4 1 Version 1 0 October 25 1996 The gradient When the adjoint model equations are solved the gradient can be calculated with this solution In the following subsectio
87. formation coefficient Ve Em 1 t Gm E VA Ema tEmuitEmmdtEmdqadt Coriolis force bottom friction contribution in u velocity point to the adjoint u momentum equation g 2 Umn Van 2 u T5 2 n C mnt Chal Um Vn y bottom friction contribution in u velocity point to the adjoint v momentum equation g Um n 2 u Cu Dan Gr n y uma Van bottom friction contribution in u velocity point to the adjoint continuity equation g Um n y Van C Dinn t aap ig bottom friction contribution in v velocity point to the adjoint u momentum equation g Tu Vm n 2 C Dm ae 1 3 Via y bottom friction contribution in v velocity point to the adjoint v momentum equation g Gas j 2 Vm n y C Dian Gundy ann Y Venn bottom friction contribution in v velocity point to the adjoint continuity equation g Vna uns vn C Apes ecl 2 r Ci hk ux vn g Uk Vic Ci hk af uc v 171 Technical documentation of WAQAD 172 frvu frvv A u m n eleg Uog n g mea C hk uz y vi g uc 2 wy C hi uz y v d onm nf ots e n ni g V 2 iz C hi external force acceleration of gravity layer thickness total water depth d C HE Ya Cmn Cm 1 n Omn Om n 1 H Ya Cm nt Gm ne1 dm 1 n dm n layer index distance between depth poi
88. g the variational formalism we find m moli p monic 3 59 12A Y Um Finally by rearranging the terms in the summation a part of the adjoint u momentum equation can be determined Pe AV n n 1 A caw aie us 2 DE S 2 el vu nns 12 Y va nns EJ 9 60 Version 1 0 October 25 1996 39 Technical documentation of WAQAD The expression between braces is part of the adjoint conservation of momentum in x direction which is solved in the computation procedure of the second stage Now we focus on the u momentum equations in those locations in which terms with Av __ are found descended m n from this cross advective term solved in the first stage of the ADI method e In the first stage of the conservation of momentum in x direction at location m n after introducing Lagrange multipliers and applying the variational formalism we find vu Yin 4 a S unna AVan 3 61 e Inthe first stage of the conservation of momentum in x direction at location m 1 n after introducing Lagrange multipliers and applying the variational formalism we find vu nein ond S aA 3 62 e Inthe first stage of the conservation of momentum in x direction at location m 1 n 1 after introducing Lagrange multipliers and applying the variational formalism we find voa cd S ass AVR 3 63 e Inthe first stage of the conservation of momentum in x direction at locat
89. has been used in the computation process of WAQUA S Vinn and S Vv Jin 3 109 and 3 110 We define for each v velocity point the logical arrays KLL KL KR and KRR with which the arrays GELUVX and GERUVX are determined KLL m n KHU m 2 n gt 0 a KHV m 2 n gt 0 a KHU m 2 n 1 0 KL m n KHU m 1 n gt 0 A KHV m 1 n gt 0 A KHU m 1 n 1 20 KR m n KHU m n 0AKHV m 1 n 0 KHU m n 1 0 KRR m n KHU m 1 n 0AKHV m 2n gt 04 KHU m 1 n 1 gt 0 Hereafter the arrays GELUVX and GERUVX can be filled Appendices GELUVX KR m n KRR m n GERUVX m n m n m 2 2 KLL m n KL m n GE oats te SY oe pep Table VIII Definition of the arrays GELUVX and GERUVX With this information the local flow situation can be determined and consequently which approximation for the v term has been used in WAQUA first second order upwind or central or zero Stelling 1984 This is necessary for the proper form of the adjoint advective terms in the adjoint computation process The boundary treatment of S Vn n and S vi is as follows e Determine the dry wet status of the v velocity points by the arrays GELUVX and GERUVX Table VIII at time t ran 4 Vinetn 4Vina n Vm 2 n for situation i v 12AX 7 30 S vVhn E for situation ii iv 0 for situation v ix The upwind difference oper
90. he adjoint u momentum equation which is solved in the computation procedure of the second stage It only affects the right hand side The appearance of the differentials Av in equation 3 88 also implies contributions to the adjoint v momentum equation 49 Technical documentation of WAQAD solved in the computation procedure of the first stage Since the discretization of v is dependent on the flow direction of the u velocity two situations have to be distinguished at each location After differentiating we focus on the equations at the locations which contain factors with the differential Avi At location m 2 n in this stage uv is discretized as Um 2 n amp Vin 2n Consequently we find 0 if 7 0 u 2 n f AV mEn if lt 0 3 90 m n 2 A X Um 2 n At location Ane 1 n in this stage uv is discretized as Um 1 n amp Consequently we find 0 if Umea Os 4u 1 n if A if lt 0 3 91 Vian 2AX Um 1 n At location m n in this stage uv is discretized as Um n amp Yinn Consequently we find Avi Suma if 0 2AX mn 3u m Av if lt 0 3 92 mn 2AX Um n ve e Atlocation imet n in this stage uv is discretized as Um t n amp Vin Consequently we find 4Um 1 n 0 2A X meee 0 3 93 50 Version 1 0 October 25 1996 The 2D adjoint model At location m 2 n in this stage uv is discreti
91. he righthand side Dmn is prescribed by Table Ill These elements have to be subtracted Now we deal with the iterative aspect caused by the addition of equation 3 105 Table IV shows the degenerations in case of impeded geometry The arrays GERUVX and GELUVX Table VIII have to be calculated for the flow field at time t The elements corresponding to v point m n contribute to the main diagonal of the tridiagonal matrix equation The coefficients have to be added to Bmn dependent on the geometrical situation The other elements in Table IV contain variables which are at right angles with the ones which set up a tridiagonal structure The iterative aspect use values from the previous loop at time t At in the first iteration This way these elements have to be subtracted from the righthand side Dm n Applying the double sweep algorithm for solving the tridiagonal matrix equation yields a first guess for all adjoint v variables at time t t The same strategy can be applied in the second iteration All coefficients remain the same except the other elements in Table IV This first guess is used for the adjoint v variable at time t 72At Applying the double sweep algorithm results in a second guess It is expected that two iterations are enough for convergence Appendices Degenerations in the adjoint v impulse equation associated with A v caused by impeded geometry In procedure AD SVU e Determine for all v velocity poin
92. he second operators has to be used 3 114 3 117 In Table IIl and Table IV the adjoint v velocity has to be 1 multiplied by the corresponding K In Table VI and Table VII the v velocity has to be multiplied by the corresponding in The variables AV and AVD has to be multiplied by n 7 26 In Table IX and Table X the adjoint u velocity has to be 1 multiplied by the corresponding E In Table XI and Table XII the u velocity has to be multiplied by the corresponding z 193
93. homogeneous boundary conditions for the adjoint model 107 Technical documentation of WAQAD 108 are assumed Therefore calibration of boundary condition parameters is not implemented in WAQAD in the case of three dimensional models WAQAD takes account of dryfall in the case of three dimensional models in the same way as in the case of two dimensional models The adjoint drying and flooding procedure described in subsection 3 2 4 is also applied to the 3D adjoint model equations The adjoint model is activated by the misfits between the observed data and the corresponding model values Residuals of observed and computed water levels affect the adjoint continuity equation If a water level measurement has been recorded at full time levels the residual has to be added to the right hand side of the corresponding adjoint continuity equation in the second stage Differentials of layer thicknesses were transcribed into differentials of water levels Since hk z 1 z and z o H 6 it follows that Ahk 5 AC AG It is assumed that no fixed layers are prescribed in the forward model A central time integration of the vertical terms in the momentum equations of the forward model is assumed Since it is possible for the approximation of the vertical viscosity to be fully implicit the Euler implicit variant has in fact become redundant There are several options to define the vertical viscosity coefficient in
94. id en gebruik van penalty voor afregelen met WAQAD in Dutch RIKZ werkdocument RIKZ OS 96 108X februari 1996 Version 1 0 October 25 1996 169 Technical documentation of WAQAD 170 7 3 List of symbols I gt aay nou M amp E ur variable at time t At variable at time t At variable at time t 2At water elevation above plane of reference averaged sea surface V2 m nt Gm V2 m nt Gm Sud E S An end horizontal curvilinear coordinates grid distance in direction grid distance in amp direction adjoint water level J17 vc An u i7 vc AG vc Jus ve mat adjoint velocity in x direction va vu m 1 AS gn Vu vu De vu Jiasi t vu adjoint velocity in y direction vy vy 1 An 3L GO Wht tO nein nat vertical viscosity coefficient UO cU U Ue X ha h XZ h h Verc Y VV X h h X h h h phase of harmonic constituent angular velocity of harmonic constituent Version 1 0 October 25 1996 fruu fruv fruwl frvu frvv frvwl fruu fruv Appendices amplitude of harmonic constituent chezy coefficient for model bottom roughness water depth below plane of reference D Ve dm n dm n 1 D Y2 dm 1 ntdmn distance between depth points in direction guu Jg the WAQUA trans
95. ields vu Yon AX 7 2 y A Um n 7 Um 1 n E AX At the second stage as Unn CIs gt 0 7 3 AX Introducing Lagrange multipliers and applying variational formalism yields vu y ths Aut n Unn E Uia bd AX u E AX e GEOUUX m n 3 In this case the term uux is discretized at the first stage as Um 1 n 7 Um n Umn AX if Unn lt 0 7 5 Version 1 0 October 25 1996 175 Technical documentation of WAQAD 176 Introducing Lagrange multipliers and applying variational formalism yields Um 1 n 7 Um n va Yon TN AX 7 6 ju A Um 1 n 7 7 n AX At the second stage as m n Unn u unes i Lg lt 0 7 7 m n AX m n Introducing Lagrange multipliers and applying variational formalism yields u u s Allg A 7 8 AMT dy 7 AU n des s AX It follows that the adjoint u impulse equation associated with AU n solved in procedure AD SUV boundary effects and dryfall and advective terms included looks like Determine for all u velocity points the array GEOUUX Table at time t 72At Version 1 0 October 25 1996 Appendices vu s DDnm a n7 GGsdal 7 7 9 with 2 CCn 1 1 2 HS AX fruu AUT At 2 n ae fruu n AU1D 7 Jmn DDm 1 n tru AU At HE egt sto AX zx fruu
96. ient for depth in the second stage becomes 1 Y Uus frvwl ma vs FVW metin 2 m n k vy frvwl m n vy Juris frvwl m 1 n Jin Um n Ena vc Jub Umn Ena Kin Kien En AG dcin Sp ed Erat Khin Ehr A ve ond Un nct Emn KO net ned AG Jii Vm n Kmn ve Vm 1 n Km 1 n Kin Ean An Km t n Erin An 7 Vm 1 n Kntn m n 1 Vm nKm n Kien et net Ko nsi Enn An The gradient for depth is calculated in WAQAD module AD GRDP 3 4 3 Boundary conditions In this subsection the expression of the gradient for the boundary condition is derived The parameters to be estimated are correction factors for the amplitude and phase of the boundary condition in point m n An overview of the contributions of the WAQUA continuity and momentum equations was given in section 3 2 3 The only terms containing the boundary condition parameters are Ve Jinn TA Gmn AA Comat Asir ot vc AG AA Coat 2 1 Asir oXt 720 6 Ad 98 Version 1 0 October 25 1996 The 2D adjoint model The expression for the gradient of the amplitude and phase at the first stage becomes vy coolt J At OA 1 Venn A sir o t At And the expression for the gradient of the amplitude and phase at the second stage becomes I y cos at OA aU Jan A sint 6 In WAQAD mod
97. ient of the criterion gives information about the direction positive or negative and the size small or big of the adjustments for each individual uncertain parameter An adjoint method can be used to determine this gradient The time needed to compute the gradient with the adjoint technique is more or less the same regardless of the number of parameters and is comparable with the computation time needed for one simulation run of the flow model The water level elevation and the velocity currents u and v flowing respectively east and north are the unknown variables in the mathematical model equations and are described by the laws of conservation of mass and by the conservation of momentum together with the model s boundary conditions These descriptors are implicitly coupled in nonlinear partial differential equations These shallow water equations contain some uncertain parameters for example bottom friction If an additional correction factor is defined for the bottom roughness at a given location the adjoint method can compute the derivative of the cost function with respect to this parameter only This is comparable with a sensitivity analysis for this parameter If more than one control parameter has to be estimated the gradient becomes a vector The mathematical derivation of the adjoint equations and concomitantly the determination of the gradient will be outlined below While searching for the optimal values the physi
98. ime is exactly equal to zero will be sufficient 71 Technical documentation of WAQAD 3 3 3 3 1 The adjoint solver This section deals with the method for solving the deduced adjoint equations for the two dimensional shallow water flow model The method for solving the adjoint equations also has a two stage time splitting structure In subsection 3 3 1 the method for solving the adjoint equations without advective terms is described Subsection 3 3 2 gives the modifications to solving the adjoint equations when advective terms are included Subsection 3 3 3 gives an overview of solving the adjoint model equations in the case of dryfall Method for solving the adjoint model equations The adjoint model equations are solved in the reversed order compared with WAQUA In the first stage in WAQUA the v momentum equation is solved in module VYD first and then the continuity equation and the u momentum equation are solved in module SUV In the first stage in WAQAD the adjoint continuity equation and the adjoint u momentum equation are solved in module AD SUV first and then the adjoint v momentum equation is solved in module AD VYD In the second stage in WAQUA the u momentum equation is solved in module UXD first and then the continuity equation and the v momentum equation are solved in module SVU In the second stage in WAQAD the adjoint continuity equation and the adjoint v momentum equation are solved in module AD SUV first and th
99. impulse equation associated with AV Equation 3 102 originates from the discretization of the term uv at the second stage in the v impulse equation The mathematical 179 Technical documentation of WAQAD derivation has been performed under the assumption that all surrounding points were always wet Boundary effects and dryfall cause degenerations in the discretization of the term uv in WAQUA In this case Table Ill has to be used The necessary arrays GELUVX and GERUVX Table VIII have to be deduced from the WAQUA velocities at time t At For example if GELUVX m 1 n 2 and GERUVX m 1 n 1 follows from the flow field at time t At situation ii in Table VIII than the indicated degeneration takes place in equation 3 102 Situation i in Table IIl corresponds to the case were all surrounding points are wet GELUVX GERUVX Situation iiy iv point point 180 Table IIl The degenerations in equation 3 102 Appendices point m 1 n lt 0 point m n Umn gt O m 1 n Um 1n gt 0 v point m 2 n Um 2 n gt 0 A part of the righthand side of the adjoint v impulse equation consists of equation 3 105 a sort of adjoint upwind difference Table IV shows the alterations in case of impeded geometry Therefore the arrays GERUVX and GELUVX Table VIII have to be calculated for the flow field at time t _ Um 2 n 2AX Table IV The degenerations in equatio
100. in time as follows C t A sir o t o Co t A sin o t in which A42 amplitude C12 angular velocity 2 12 phase 6512 Since the angular velocities of both constituents are not identical the phase difference of both constituents can vary during the simulation The longer the simulation period lasts the bigger the variation can become In this section it is assumed that if the variation of the phase difference becomes more than 45 the constituents are too far apart to be separated Let S be the simulation period The next simple test must then hold 6 0 S in 5 22 Such a test does not have to be implemented for the amplitude ratio between the two constituents for this will stay the same during the period that is considered the amplitude is not a function of S Coupling the constituents WAQAD calculates the gradient For the coupled constituents the gradient can be determined in two ways 1 By deriving the gradient for the coupled constituents all over again 2 Byusing the gradients already derived in a considered way The very laborious first method is not applied since using existing knowledge an elegant derivation of the gradient can be found Gradient of the phase If 2 and constituents and amp are coupled together the phase can be calibrated in two ways Version 1 0 October 25 1996 Possibilities options and shortcomings 1 both constituents are
101. int equation at each grid point for each time step Similarly to the forward model we can distinguish an inverse continuity equation and inverse momentum equations If we select the Lagrange multipliers such that they satisfy the adjoint equations the unknown terms in this expression for the gradient will be multiplied by zero More precisely the adjoint variables vu and v should be chosen such that at each grid point at every time step the corresponding adjoint equation holds ak ko 4 a ps Ce if measurement is available Fel VesVusVv O if measurement is not available AF ul Ve Vue v 70 AF v vesVusVv 9 The unknown terms vanish leaving a simple expression for the gradient of the criterion The gradient for this selection of the adjoint variables can be computed os op m n k OF Op OF SE i k ee vc jm Vu yen The adjoint method is a so called inverse model The adjoint equations are solved going backwards in time starting with a homogeneous initial condition Actually they describe the propagation reverse in time of the co state The adjoint system is activated by the disparity between the measurements and the computed model counterparts inside the system Fortunately the co states in the interior are independent of the adjoint variables located at the boundaries therefore homogeneous adjoint boundary Technical documentation of WAQAD conditions can be
102. ion The result is a tridiagonal matrix equation for each row consisting of adjoint contributions Ve m1 Yelm v Welma Ve Ve Velm Ve 3 135 This matrix equation is solved with Gaussian elimination double sweep method Stelling 1984 The adjoint u momentum equation can then subsequently be solved by simple substitution AD VYD In procedure AD VYD the adjoint v momentum equation is solved to determine v mn for all inner points of the domain of the model So far the adjoint advective terms have been omitted The existing adjoint v momentum equation at this first stage of the adjoint computation process reads A Jan Man 1 fs st frav T Gs Cuan 70 3 136 Version 1 0 October 25 1996 The 2D adjoint model Without the advective terms a diagonal matrix equation over a grid line in y direction is obtained Brin vv Dis 3 137 in which Bn 2 At 2 n En Dm n At vy nn ji fmn 1 fru Veal Jo Jon Adding 3 101 3 104 omitting 3 105 to this existing adjoint v momentum equation yields a tridiagonal matrix equation for a column in the domain of the model 3 138 83 Technical documentation of WAQAD This banded system could be solved over the grid lines in y direction by standard procedures The corresponding coefficients of this matrix equation are in which Vm n 1 mn Ac 3 140 Am 2A Y B B
103. ion m n 1 after introducing Lagrange multipliers and applying the variational formalism we find vu ewe VT 15 Um n 1 eee 3 64 Finally by rearranging the terms in the summation we obtain AV n 4 Sy us s vu Jrisiin 2 S Um 1 n v 1S uni na v Jana is Um n 1 3 65 40 The 2D adjoint model The expression between braces is part of the adjoint conservation of momentum in y direction which is solved in the computation procedure of the first stage In the second stage the difference uy is approximated by a second order upwind difference implicit This results in the following approximation of vuy Vinny Uma 3 66 in which me NG tee Vm n V mn 2 Vm n Vm 1 n Vm 1 n 1 Vim n 1 ES if V ma 0 2AY Stn us s dez if V ors 0 2 Y Differentiating this expression with respect to an arbitrary model parameter yields AV an 49 Um n Tt AV ntn iS Urn AM e aat 49 DUM AM ie 45 Uhn added to 35 45 V m n AU V m n i 2 Y 2AY TT mn if m n 2 Zee V m n or added to 35 AS AUT AUS AU a4 T t 2 i 2 AUT sco if vn 0 Equal to 3 67 Sy Unt vi SC AUS Version 1 0 October 25 1996 41 Technical documentation of WAQAD 42 First we concentrate on the counterpart of this term in the adjoint v momentum equation Therefore
104. ion several derivatives of discrete model variables are found located at different grid points The entire formula has to be multiplied by the corresponding adjoint variable To build up the adjoint equation this summation of products has to be rearranged one derivative of a discrete model variable located at one grid point has to be multiplied by a formula in adjoint variables First of all we concentrate on the u momentum equations at different locations in which terms with A are found descended from this cross advective term solved in the first stage of the ADI method 38 The 2D adjoint model In the first stage in the conservation of momentum in x direction at location m n 2 after introducing Lagrange multipliers and applying the variational formalism we find V mn 2 m i A 5 3 56 Jm 12A Y Um n In the first stage in the conservation of momentum in x direction at location m n 1 after introducing Lagrange multipliers and applying the variational formalism we find AV mnt f sd X mie 3 57 12 Y u 2 In the first stage in the conservation of momentum in x direction at location m n 1 after introducing Lagrange multipliers and applying the variational formalism we find LAS ee EE Ae ee 3 58 om 12A Y Um e Inthe first stage in the conservation of momentum in x direction at location m n 2 after introducing Lagrange multipliers and applyin
105. ions the equations resulting from the process above are averaged over the vertical coordinate from bottom to surface This integration introduces a number of new terms one of which is the stress turbulent and viscous on the bottom A new assumption has to be made here to obtain a closed set of equations The bottom stress is often assumed to be directed opposite to the depth averaged flow and proportional to the square of the magnitude In the past measurements have indicated that the friction depends on the local water depth and several empirical or mainly empirical formulas for this dependence are widely used especially the Chezy Manning and White Colebrook formulations In WAQUA in SIMONA one of these three formulations for the dependence of friction on the local water depth has to be chosen Since all three formulations contain a parameter usely known the calibration procedure should be able to deal with all three types of parameters The Chezy formulation had already been implemented in WAQAD The friction terms in this documentation in chapters 3 and 4 were therefore described for the Chezy formulations Recently the WAQAD system was extended to use Manning and White Colebrook formulations too In this section the extension of WAQAD to the other formulations is described 157 Technical documentation of WAQAD 5 6 2 158 Theoretical background For a detailed theoretical derivation of the friction formulations see
106. irection in the first stage v v 5D Hf S gsal wk 4 1 The vertically integrated continuity equation in the first stage HUE HVK hy At K 4 2 Conservation of momentum per layer in x direction in the first stage u g ets 8 U 4 3 X A fv X Go 4 3 in which H Total water depth U Depth averaged horizontal velocity component in x direction V Depth averaged horizontal velocity component in y direction Horizontal curvilinear coordinates K gvv the WAQUA transformation coefficient in direction see subsection 3 2 1 E guu the WAQUA transformation coefficient in direction see subsection 3 2 1 h Layer thickness at u point also denoted as hu h Layer thickness at v point k Layer index K Number of layers Zk Position of layer interface At Integration time step Remarks The state in the new time level t 72 t is denoted by a single apostrophe The absence of an apostrophe indicates the state in time level t For convenience only indices which deviate from m n k are denoted The terms considered in these forward model equations do not The 3D adjoint model require a special boundary treatment Therefore the corresponding adjoint equations do not change near boundaries e Inthe TRIWAQ system a grid definition in the vertical plane needing the prescription of K 1 layer interfaces is used The layer thic
107. isk space can be a serious problem especially in the case of a 3D simulation However the technique of check pointing can be used to overcome this problem Using this technique the calibration period is split into smaller periods that are calibrated separately and subsequently when the first period is calibrated the data estimated in this period is used as restart data in the next period which is then calibrated and so on For more information on this subject see the WAQAD User s guide Once the gradient has been computed the realm of minimising functions can be entered As mentioned before the gradient is a directional vector in the parameter space Quasi Newtonian formulae or conjugated gradient methods alter this direction to search for a better estimate Nowadays experts generally agree that the BFGS Broyden Fletcher Goldfarb Shanno updating formula is superior to others for medium sized lt 50 problems of parameter estimation A line minimisation finds the values of the parameters along the search direction that have the lowest value of the cost function With exact line minimisation a Quasi Newtonian Version 1 0 October 25 1996 The adjoint method method has second order convergence For a thorough treatment of this topic the interested reader is referred to Fletcher 1987 A short description of line minimisation is given below line search When a new search direction is found the stepsize that minimises the cos
108. its two dimensional version The 3D adjoint model AD VYD3 The adjoint conservation of momentum per layer in y direction in the first stage 4 6 can be written in the following form ak vv 1 Dk vv dk vv deed d 4 14 in which _ n Vv eer 12 hia hk 1 h 1 n Vv Vv wA EN na m hh hth n Vv s 1 h hk acero RT In more compact form A vw d 4 15 in which Guys d d jy 3 dk and An n n D 12 hk 1 hk 1 h C 8 Version 1 0 October 25 1996 115 Technical documentation of WAQAD 116 This tridiagonal system in adjoint v velocities can be solved by the well known double sweep method For this the layer thicknesses are determined from the computed total water depth and the sigma coordinates and at the exterior layers a term associated with wind stress or bottom roughness is also defined The two dimensional approach is used here Part of the vertical viscosity term still remains on the main diagonal associated with these layers Remarks Since no uncertain boundary conditions are assumed it is not necessary to solve the adjoint boundary conditions Therefore calibrating boundary condition by using 3D adjoint is not implemented in WAQAD The drying and flooding procedure in the 3D case is implemented in WAQAD in the same way as in the 2D case This procedure i
109. kness hx is defined as z 4 Zk Notice that zo equals the water level and equals the bottom with opposite sign The sigma coordinates o are obtained using the following transformation figure 4 1 Sigma transformation The layer thickness is defined in a relative way as a constant part of the total water depth However it is also possible to use layer thicknesses defined in an absolute way The Chavent method as described in chapter 2 is applied to these forward model equations 4 1 4 3 to obtain the corresponding adjoint equations Since the derivation has already been described in detail for the two dimensional case in chapter 3 the derivation for the three dimensional case is omitted here The results of the derivation are summarized Version 1 0 October 25 1996 105 Technical documentation of WAQAD Adjoint conservation of mass in the first stage y T S vu y Y 1 gt At K qunm m oxi Gut Cuts k 2 wh Ch yet 4 4 vc y K E m it gt x Lu N n n hi Lh dS eec k 1 u k x k 10 h h hy HS sueco sane an e 6 ER u u 7 ok vu x l E 2 x N x N er 2p k 1 A 7 vik A n t nv ng Es we A Wea 1 ID v Vi
110. l the minimum with the required accuracy is determined There are different methods for calculating the gradient Using the adjoint model it is possible to determine the complete gradient vector for an undetermined number of parameters in one simulation run This method on which WAQAD is based is fully described in this documentation Finite difference is the classical method for calibrating models The minimization of this method shows the same iterative progress as WAQAD except that the gradient is calculated differently WAQAD calculates the gradient for each selected parameter in one single adjoint simulation run Using the finite difference method a WAQUA simulation run must be started for each selected parameter in order to determine the gradient Using the finite difference method the gradient of a selected parameter is calculated as follows e Determine value JO of the cost function of the initial WAQUA run see the WAQAD User s Guide for calculation of the cost function e Adapt the values of the selected parameters with a prescribed value S and determine the value J1 of the cost function of the WAQUA run with adapted input Calculate the gradient G JO J1 S Step size S must be chosen small enough so that gradient G approximates the direction coefficient of the parabola However if the step size is too small problems of machine 162 Version 1 0 October 25 1996 Possibilities options and shortcomings a
111. les of the first stage are discussed AD SUV3 and AD VYD3 109 Technical documentation of WAQAD 110 AD SUV3 The adjoint conservation of mass appears to be a pentadiagonal matrix equation Combined with the adjoint u momentum equation a tridiagonal system in adjoint water levels is obtained which can easily be solved and afterwards the adjoint u velocities can be determined The starting point of the description of the solution method is that the adjoint conservation of momentum per layer in x direction 4 5 can be written in the following form Aele bn vu je t 4 7 V6 m k RIA Em k vado fm k Vu uoti dm k in which ge hE m k AE 1 ae E Dm k Vy Vv X yh hath Zh h h WE Cm k AE n k 7 7 i nei 1 h 5 fm k Wy 1 he heh j2 At The layer thicknesses can be calculated from the total water depth and the sigma coordinates Version 1 0 October 25 1996 The 3D adjoint model In more compact form A J d vc n a vc 4 8 K E ES in which T 6 ion d d doe dul a a 852 e and Dans LPS n2 bo LN A L Dink ts fm dT Ww 25 he hi he Ww m K 7 16 hk 4 hk 1 h In the uppermost layer a term associated with wind stress is defined however the wind stress in WAQAD is supposed to be eq
112. lita Bn Cm Vb udin Dm Coefficients Am Bm Cm and Dm are calculated by 1 a Urin En n L 35 166 2 Am m 1 n Kinin Khana 2 4 AE 73 Technical documentation of WAQAD Bn Km 1 n AG 1 i2 1 1 1 fruw KB 2 Unn Ena Cn Kiedin En 1 Un at g 1 i fruwl Ue ae a TUW lh oll ve Ta Enn _2 von T AS _2 Kin Ean Kaen En 1 3 Un tn En 1 n AS This system of matrix equations is transformed into a tridiagonal system 74 Kn t n Ein Kn Ean Tus Umn En E truw lnn CCm P 2 C ou T Dara mf Version 1 0 October 25 1996 The 2D adjoint model which is written as the following equivalent system of equations 16 a E 22 di Ve Jan Ds 1 Gs a 1 Cu n Duci D Dnt or von Cin vean Dm vm These coefficients are determined recursively Cn Cn Bm AmCmn 1 zDs Aubas 3 Bn An Cs C Ci and D4 D Solving backwards using the same formula gives v Jimin Dmr Dm T 6 The adjoint u velocities are then calculated by substituting the adjoint water levels into the adjoint u momentum equation von v n DDm x i i Eun Ean 75 Technical documentation of WAQAD AD VYD The adjoint v
113. litude of harmonic constituent Qj angular velocity of harmonic constituent i i phase of harmonic constituent i The discretized equation in space and time for the boundary condition described by 3 130 is the starting point for the deduction of the corresponding adjoint boundary conditions Similar to the derivation of the adjoint model equations the adjoint boundary conditions can be derived by following the next three steps e The boundary conditions described above are constraints just like the WAQUA model equations Remark Since the boundary conditions are not defined by state variables at interior points of the model and since itis solely a function of model parameters the deduction of the adjoint model equations in the previous section does not change However some terms concerning the adjoint variables at the boundary must be added e The terms to be added are derived from the linearized WAQUA equations by looking to see which terms contribute to the boundary e The adjoint boundary conditions are then deduced by rearranging the contributing terms and assigning them a value in such a way that all AG do not contribute to the expression of the gradient 61 Technical documentation of WAQAD 62 In the next part of this subsection the three steps are performed for each open boundary in order to derive the corresponding adjoint boundary conditions Left hand open boundary The velocity parallel to the bounda
114. ls are frequently used to calculate detailed flow patterns and to estimate water levels so that water movement can be simulated These computational flow models often contain empirical parameters which are not known very accurately The most crucial phase in the development of a hydrodynamic model is the calibration Generally a model can be calibrated on the following parameters bottom roughness bottom topography and harmonic constituents in the boundary condition All of these parameters are known only to a certain degree of accuracy and need to be optimally adjusted in the calibration process The model builder is free to adapt these parameters within certain bounds Observations are used for comparison with the model outcome and empirical parameters in the model are adjusted in order to improve the model s performance Flow models used to be calibrated fully manually by experts This trial and error approach was very time consuming and laborious In recent years advanced mathematical techniques in the field of optimal control theory have been successfully applied to model calibration In a process generally called data assimilation information derived from observations is combined with the model outcome in order to improve the model s performance The adjoint model developed by this process has become particularly popular It is an inverse modelling technique in which misfits between computed and observed values are retranslated into adj
115. lty is independent of the parameterization The latter means that when for instance a rectangle is divided into two smaller rectangles and the adaptations of the two rectangles are equal the penalty equals the penalty of the original situation In this case also any difference between the adaptation of the two rectangles receives an additional penalty Practical use Although it is very important to understand the mathematics of the penalty in order to derive the formulas to use the penalty it is sufficient to understand what it does rather than how it does it In general the penalty will not allow adaptations that are too large What is considered too large can be controlled by the standard deviations and characteristic lengths that are required in the input Usually it will suffice to take reasonable values known from experience or intuition The result of the adaptation process should mostly not be very sensitive to the exact values However if the adaptations after a first run are too large or too small for the user this can be changed How the input of WAQAD affects the adaptations is explained in detail in the WAQAD User s Guide and Verlaan 1994 In Verlaan 1996 the use of the penalty is also illustrated by means of an experiment Version 1 0 October 25 1996 139 Technical documentation of WAQAD 5 1 7 Conclusion The WAQAD program always gives output to the input of a calibration problem Also the calculated solutio
116. me level the residual has to be added to the right hand side of the 28 3 2 2 3 2 2 1 3 2 2 2 Version 1 0 October 25 1996 The 2D adjoint model corresponding adjoint continuity equation 3 21 at the second stage e Homogeneous boundary conditions are assumed for calculation of the adjoint state variables in the interior model field Only if uncertain parameters of the WAQUA boundary conditions are estimated do the corresponding adjoint variables at the boundary have to be calculated afterwards This is described in subsection 3 2 3 Advection Introduction Advective terms can play an important role in modelling the flow in estuaries Stelling 1984 Such terms occur only in the conservation of momentum These terms are uu and vuy in the u momentum equation and vv and uv in the v momentum equation In the first stage in WAQUA the v momentum equation is solved first and then the continuity equation and the u momentum equation are solved together In the second stage in WAQUA the u momentum equation is solved first and then the continuity equation and the v momentum equation are solved together In subsection 3 2 2 2 the structure of the solution method of the discretized WAQUA equations with advection terms is described For simplicity the external forces are set equal to zero In order to obtain a unique solution to the discretized equations a set of boundary conditions has to be prescribed The boundary treatm
117. meters to be estimated can be changed For example it is impossible to estimate the bottom depth using only a few measurement stations bu it is possible to estimate the average bottom depth instead Reducing the number of parameters will generally lead to a better estimation problem And using a method called the penalty can also solve some l problems 5 1 2 Some examples The following examples using a simple linear one dimensional problem illustrate a number of l problems that can occur during the WAQAD calibration process Figure 5 1 shows a one dimensional problem with an open boundary on the left hand side and a solution that converges to zero at infinity on the right hand side This problem illustrates for example a Kelvin wave that progresses counter clockwise through the North Sea with the x axis following the coast line The equations are described by a simple linear version of the shallow water flow equations o6 pou 0 ot OX ou OE f g ot ox Du C 0 Hosir ot lim G X t 0 5 1 126 Version 1 0 October 25 1996 Possibilities options and shortcomings The solution to this problem at constant depth in first order is given by ae C x t Ho sn i e t x 5 2 with a linear friction with coefficient f f is small In this solution a harmonic wave is recognized which progresses from the left to the right at velocity V9D en which is slightly atten
118. ms in the discretised model equations with respect to each individual model state at a certain point of time can be determined after simulating the tidal flow for the calibration period in the case of a non linear underlying forward model These are efc eR OG du v OF OF OF v OF OF OF du v The differentials of the terms in the discretised model equations with respect to the control parameter can be determined after the simulation These are oF OF Op dp Op The difference between the measurements and the computed water levels can be determined after the simulation This term is So known and unknown terms occur in this expression for the differential of the cost function with respect to the parameter The adjoint variables and v are unrestricted Any value is allowed This gradient can be computed by the Chavent method Rearranging the summation in 2 3 yields Version 1 0 October 25 1996 The adjoint method 0J 06 T A uz 4F Vu Vv i m n Op m n k Op vc Tunt c NT Oui i 2 Ve Vav m n k Op B QN 2 4 E X AE VesVusVv m n k Op OF k OF k F 2 v jus anis T vy Jan M g Op Op gt p The AF functions describe the so called adjoint equations Every differential of a discrete model variable with respect to a control parameter is now associated with an adjo
119. n 3 105 The terms corresponding to v point m n contribute to the main diagonal The other contribute to the righthand side sign changes due to the Jacobi method The last part left in these adjoint advective terms of this adjoint v impulse equation is equation 3 101 It consists of adjoint variables at the new time level which set up the matrix structure and adjoint variables at the old time level To determine the degenerations caused by impeded geometry an array GEOVVY has to be constructed from which the time dependent local flow situation can be determined Table V In WAQUA impeded v velocity points which inhibit the flow of water are denoted by KHV lt 0 Version 1 0 October 25 1996 181 Technical documentation of WAQAD KHV m n 1 KHV m n KHV m n 1 GEOVVY m n 20 gt 0 50 0 lt 0 0 1 gt other Table V Definition of array GEOVVY The situation GEOVVY m n 1 corresponds to the case were the surrounding v points are wet Now we focus on the other three situations The coefficients Amn and Cmn from the corresponding tridiagonal matrix equation can be read from Table VI The array GEOVVY has to be constructed from the variables KHV at time t Table V jp a edi Vmn Amn V point m n 1 Cmn V point m 1 2AY Table VI The degenerations in equation 3 101 For another contribution to the right side the array GEOVVY has to be constructed from the variables
120. n debietranden Eindige differentie in Dutch SIMTECH 1995b Brouwer J R Faserapport fase Ill Onderhoud WAQAD Visualisatie satelliet data en aanpassingen frictietermen in Dutch SIMTECH 1995c Brouwer J R e a WAQAD System s documentation SIMTECH 1996a Brouwer J R e a User s guide WAQAD SIMTECH 1996b Brummelhuis P G J ten Parameter estimation in tidal models with uncertain boundary conditions TU Twente Ph D Thesis 1992 Bryson A E and Y C Ho Applied Optimal Control Hemisphere New York 1975 Chavent G Identification of distributed parameter systems about the output least square method its implementation and identifiability Proc 5th IFAC Symposium on Identification and System Parameter Estimation Vol I New York Pergamon Press 1980 Fletcher R Practical methods of optimization New York John Wiley and Sons 1987 Giering R and T Kamanski Recipes for adjoint code construction Max Planck Institute f r Meteorologie 1995 Goede E D de Numerical Methods for the Three Dimensional Shallow Water Equations on Supercomputers Ph D Thesis University of Amsterdam The Netherlands 1992 Lander J W M P A Blokland and J M de Kok The three dimensional shallow water model TRIWAQ with a flexible vertical grid definition RIKZ werkdocument RIKZ OS 96 104X Rijkswaterstaat 1996 168 Appendices Leendertse J J A New Approach to Three Dimensional Free Surface Flow
121. n is always better if a comparison is made using the measurements that are assimilated by WAQAD However the calibrated model is not always really better It is therefore very important to be able to analyse in which cases WAQAD can be used in a useful way or how the problem can be changed in such a way that WAQAD is useful Using a penalty enables a number of problems to be avoided Large unrealistic adaptations are no longer possible and furthermore large differences in adaptations for points that are close to each other will be prevented The penalty gives us the freedom to indicate which parameters are known accurately or not If desired this can be used to stress certain aspects of the model Although in many cases the results will now be acceptable to the user it remains important to have many well spread data available Only parameters that affect the results of the simulation at the points and times of the measurements can be estimated 140 5 2 5 2 1 Version 1 0 October 25 1996 Possibilities options and shortcomings The use of weighting functions for parametrization Introduction When trying to identify a parameter here denoted by Q that depends on spatial coordinates a set of weighting functions needs to be selected to specify the types of adaptations that are allowed for changing the initial estimate The adaptation becomes a linear combination of the weighting functions gi x y AQ x y X219 X y
122. n n 3 141 7 Vm n 1 Pe s IM 3 142 Cn 2 Dm n Dm n 3 143 T MeV v p n Misi Vy Jani 3 144 2AY AU mtn vy Je 12AX 12AX Z _ 3 145 Au n n U m 2 n vv Jis 12AX Denm 12A X j Urs 7 n m Uds j S i v 1 iach K 3 146 D Jose 16 Um 4n 19 Uninet a vu Jes 194 Um n Jaan 194 Um 1 n 3 147 cine 1 Um 1 n 1 vu are 1 Um n 1 Remark Note that the coefficients A B C and D that are calculated in module ADSUV are not the same as those calculated in module ADVYD 84 The 2D adjoint model Adding equation 3 105 to equation 3 139 yields Ban Bn 3 148 SUm n if 3 lt 0 2AX ne E JUm n if 7 gt 0 2AX Adding equation 3 105 to equations 3 143 to 3 147 yields Dm n Das 3 149 Um n S x vy ira if Um 2n 0 Um 1 n S vna E if qo 0 Um 2 n vy revo nen if 5 0 The system is solved iteratively using a Jacobi method The computed values from the previous loop are used to compute the marked adjoint variables in equation 3 149 AD SVU In procedure AD SVU the adjoint conservation of mass and the adjoint v momentum equation are simultaneously solved to determine and v Since this computation procedure is analogous to the procedure AD SUV only the important formulas will be given below Version
123. n of the weighting functions it is important to identify triangles with common sides are identified Therefore the storage of the triangles is indirect First the coordinates of the points that are sides are given together with a unique number The triangles are now described by enumerating the numbers of the sides 143 Technical documentation of WAQAD To ensure continuity on the boundaries of the region where adaptations are made it is possible to make sides of triangles inactive as a parameter This option can also be used to deactivate other sides when needed for some reason The implementation of triangles in WAQAD is described in more detail in Rozendaal 1995 5 2 5 Conclusion The use of other basic functions can improve the performance of the WAQAD system by dealing with the problem of discontinuities that are introduced when using indicator functions on rectangles as weighting functions These new weighting functions also make it easier to follow the model geometry When fully implemented it will probably be an improvement to the WAQAD system It will be necessary to evaluate the performance on a real life model since it cannot be guaranteed that these new weighting functions will not introduce unwanted effects At present rectangles and triangles are implemented in WAQAD For a detailed description of the input requirements for these weighting functions see the WAQAD User s Guide 144 5 3 5 3 1 5 3 2
124. ns the expression of the gradient is derived for three cases e the parameter to be estimated is a correction factor for the bottom friction coefficient subsection 3 4 1 e the parameter to be estimated is a correction factor for the depth subsection 3 4 2 e the parameters to be estimated are correction factors for the amplitude and phase of the boundary condition subsection 3 4 3 Bottom friction In this subsection the expression of the gradient for bottom friction is derived The parameter to be estimated is a correction factor for the Chezy coefficient at u velocity point m n and v velocity point m n For bottom friction in the case of Manning and White Colebrook see section 5 6 An overview of the contributions of the WAQUA u and v momentum equations is given first Then the expression of the gradient is deduced e Contribution of the WAQUA v momentum equation in the first stage from t gt t 72At at point m n 2 Vm n re y T Vmn y g 94867 Wm n m n e Contribution of the WAQUA u momentum equation in the first stage from t gt t 72At at point m n 2g uus Y Chay Date van m n e Contribution of the WAQUA u momentum equation in the second stage from t 72At t At at point m n 2g uus Y 15 hay Dic van m n 93 Technical documentation of WAQAD e Contribution of the WAQUA v momentum equation in the second stage from t 72At t At at point
125. nsional system TRIWAQ with curvilinear coordinates in the horizontal plane and a general description of layer interfaces in the vertical plane This prescription can be made by defining the layer thickness either relatively or absolutely in relation to the total water depth For an extended description of this three dimensional shallow water flow model see the technical documentation of TRIWAQ Lander et al 1996 The adjoint model for the simplified three dimensional shallow water equations was derived by Mouthaan 1995 in order to calibrate TRIWAQ This chapter describes the 3D adjoint model as it is implemented in the current version of WAQAD version 2 1 see the WAQAD User s Guide Brouwer 1996b and the System Documentation of WAQAD Brouwer 1996a For more information on the implementation see Rozendaal 1996 102 4 2 Version 1 0 October 25 1996 The 3D adjoint model The adjoint model equations In this section the adjoint model for a simplified three dimensional hydrodynamic sea model will be presented The adjoint equations were derived using the Chavent method outlined in chapter 2 But before presenting these the hydrodynamics of the forward model are described As in the case of the 2D model the mathematical description of three dimensional water flow consists of a system of coupled differential equations which are physically based on the laws for the conservation of mass and momentum deduced from th
126. nstraints as functions of the unknown model parameters Then differentiate the constraints to the model parameters e Rearrange the linearised equations by ordering the contributions of each multiplication factor These contributions are the adjoint equations which must be solved The three steps are performed in the next part of this subsection First some remarks about the discretized equations e Rectangular spherical and curvilinear coordinates The discretisations of the model equations is dependent on the model grid used In WAQAD there are three possibilities the grid points can be given in rectangular co ordinates spherical co ordinates or curvilinear co ordinates Rectangular and spherical co ordinates are considered to be special cases of curvilinear co ordinates The following grid dependent variables are used in the discretized equations AE grid distance in amp direction An grid distance in direction Kmn distance between depth points in amp direction Ema distance between depth points in direction In the case of rectangular co ordinates A amp equals Ax and An equals Ay whereas Km and Emn are both equal to 1 In the case of spherical co ordinates the variables Km n and Em are determined by Technical documentation of WAQAD Kmn RA Emn Rcos AX in which R radius of the earth grid angle in the northward and eastward directions In the case of curvilinear coordinates th
127. nts in amp direction JVV 9 the WAQUA transformation coefficient Km n 1tKm n K n Ya Km n 1 Km 1 n 1 Km 1 ntKmn horizontal grid index in m direction horizontal grid index in y direction time integration timestep horizontal velocity component in x direction Un u U Unga FUR 454 Ug 4 n 1 depth averaged horizontal velocity component in x direction K 1 Yuh H U U AG horizontal velocity component in y direction Version 1 0 October 25 1996 I Von n lt II Von Zk Appendices als I V Vna Vine tin 1 Vine n I depth averaged horizontal velocity component in y direction K 1 vh H a V Va An position of layer interface 173 Technical documentation of WAQAD 7 4 Advection in the case of bounderies and dryfall Up till now we assumed that no degenerations take place in the discretization of the advective terms in WAQUA although this discretization depends on the dry wet status of surrounding points In fact it was assumed that the surrounding points were all wet The presence of boundaries and the dryfall procedure are responsible for the deformations of the advective terms These alterations have consequences for the adjoint advective terms in the inverse model That is the subject of this appendix 7 4 1 Closed Boundaries Closed boundaries are land water boundaries These are represented by zero
128. ntum equation at location m n we find Cae E 3 80 e Inthe first stage in the u momentum equation at location m 1 n we find Ut ai JY n esu ET 3 81 n the second stage in the u momentum equation at location m 1 n we find Ur n vu ar TUN SC Jus um 3 82 Inthe second stage in the u momentum equation at location m n we find Urin Ta ju au som 3 83 Inthe second stage in the u momentum equation at location m 1 n we find Un n vu Jovis p E x Jos EE 3 84 Rearranging the summation in equations 3 79 to 3 84 and thereby constructing parts of the adjoint equations yields Um 1 n 7 Um 1 n A on B yg Umn Joa 2AX 3 85 U m n 7 Um 4n 2AX Version 1 0 October 25 1996 The 2D adjoint model no 2AX 3 86 Um 1 n Vu n Umnet n Vu Ye 2AX The expressions between braces are part of the adjoint u momentum equations which are solved in the computation procedures in the first stage 47 Technical documentation of WAQAD 48 UVx This cross advective term occurs in the v momentum equation In the first stage of the ADI method a second order upwind difference is used to approximate v implicitly In the second stage a weighted central difference approximates vx explicitly The treatment is analogous to the other cross advective term vuy in the u
129. odifications that must be made as a result of the advective terms when solving the existing adjoint computation procedures see subsection 3 3 1 The description is given per WAQAD module AD SUV AD VYD AD SVU and AD UXD The degenerations caused by the boundary treatment and impeded points which may inhibit the flow of water have been left out These are discussed in appendix A In this subsection the advective terms are described using rectangular coordinates For an overview of the differences between rectangular curvilinear and spherical coordinates see subsection 3 2 2 4 AD SUV In procedure AD SUV the adjoint conservation of mass and the adjoint u momentum equation are simultaneously solved to determine Jmn and vu mn for all inner points in the domain of the model In imitation of module SUV in WAQUA the adjoint u momentum equation is used to eliminate the adjoint u velocity from the adjoint conservation of mass This leaves a tridiagonal system of equations in adjoint water levels for a row which can be solved by the well known double sweep algorithm Stelling 1984 After this the adjoint u momentum equations can be resolved Without the advective terms the adjoint u momentum equation in this stage looks like necs l HY elke un Feo S Ost 16 131 an Fly ve Xan The 2D adjoint model Adding equation 3 100 as described in subsection 3 3 1 to the existing adjoint umomentum
130. off line data assimilation It is used to make best estimates of uncertain parameters in a mathematical model The mathematical theory for identifying parameters in distributed systems is based on Chavent 1980 and the methods for computing the shallow water equations are based on Stelling 1984 The National Institute for Coastal and Marine Management RIKZ has developed the technique of applying Chavent s method to models based on Stelling s method for conventional in situ data thereby combining advanced theory with important practical applications The technique is outlined below The cost function J p provides an objective way of assigning preference to certain values of the parameters This function is a measure of the model error Usually the J p is chosen according to Xp db AN 2 1 m n k with E i Observed water level at location m n at time k p Vector of control variables ac Model counterpart of observation The summation in this quadratic criterion runs over all locations in the domain and over all time levels in the simulation period for which measurements are available The control variables can be defined as additive or percentual correction factors for model parameters The set of parameters p is optimized in order to minimize the cost function J The more accurate the estimates of the parameters the better the agreement between computed and observed water levels and the smaller the cost function The b
131. one can be used for calibration The cost function of the other set could be reproduced in a graph in order to check the first set The second data set can also be obtained by splitting one data set in time or into stations Figure 5 7 shows an example of both graphs Using the steepest descent method for BFGS too the steepest direction of the descent will be followed in the first iteration In these directions iterations the model will really improve and the cost function of the control data will also decrease After a number of iterations no steepest directions remain and then the direction of the line of the descent will also be searched In this direction the cost function will decrease slightly but the improvement is only cosmetic because of the observation errors In fact the model will get worse and the cost function of the control data will increase Using estimation of the Hessian The convexity of the cost function at its minimum can be searched by means of the second order derivative of the cost function Hessian The contour lines of the cost function in the neighbourhood of the minimum are approximately elliptical The directions of the axes of these ellipses are formed by the eigenvectors of the inverse of the Hessian and the ratios of the lengths of the axes are formed by the ratios 131 Technical documentation of WAQAD 132 of the square roots of the eigenvalues of the inverse of the Hessian Figure 5 8 shows the
132. ore complicated if a combination of parameters cannot be estimated For example trying to estimate depth friction and the phase of the boundary conditions simultaneously from data from only waterlevelstation can lead to the following problem If the tidal wave arrives too early at the measurement station this can be due to the wave velocity or to the phase at the boundary Now the adaptation routine does not know which of the two parameters to adapt In an attempt to solve this problem the method will change both slightly This part is inherent to the problem and more data is needed to solve it The biggest problem however is that it is possible to adapt the depth and the phase at the boundary such that the model predictions at the measurement point do not change This kind of problem can give very large and unrealistic values for the parameters Various ways of dealing with these problems have been suggested in the literature Verlaan 1994 and Verlaan 1996 Most can be written as the addition of an extra term to the cost function This additional term will be called penalty from now on This penalty often takes one of two forms One way is to penalize large adaptations for example by adding a penalty that is proportional to the square of the change in the parameters The second is meant for parameters that describe the spatial distribution of a certain property like depth Here usually differences between adaptations that are clos
133. ormulations are described by the factors shown in table 6 1 160 Possibilities options and shortcomings These extensions for Manning and White Colebrook friction formulations are implemented in WAQAD For more information on this implementation see Brouwer 1995c Version 1 0 October 25 1996 161 Technical documentation of WAQAD 5 7 Finite difference calculations Models can be calibrated by selecting input parameters which may be adapted in such a way that the prediction results are the best possible The square mean of the error error prediction observation is usually chosen as the critirion for the prediction value A certain input parameter can be adapted using steps of different sizes For each step a simulation is run and the corresponding prediction value is calculated In this way a number of points are calculated with X step size and Y prediction value These points form the cost function of the corresponding parameter to be calibrated The cost function is parabolic A model is calibrated by searching for the minimum of the cost function of the selected input parameter This minimization is carried out by calculating the gradient direction coefficient of the tangent of the parabola for X 0 for each parameter after which a search process is started During this search a number of simulations are run with different step sizes for the adaptions of the parameter in the direction of the gradient This is done unti
134. point to the adjoint v momentum equation at time t 72At fruwl contribution of bottom friction at u velocity point to the adjoint continuity equation at time t 72At 24 The 2D adjoint model The bottom friction contribution is given by 2 U mn 24 ve 2 fruu mn a vin a E Ci na Cae U mn vin A n fru mn 9 pris 9114 C ane U m n ye Vim ae 7 g Dice fruwl pn Formula 3 14 describes the Chezy formulation For further information on the friction formulation see section 5 6 of this documentation Remark The overlined variables in equations 3 13 and 3 14 indicate that they are interpolated to u velocity points This notation will be used in this documentation from now on e The differentiated u momentum equation in the second stage from t 72At for is Auna fmn A Vma 2 g AEG 3 15 Kenn AB n fru Eich AU mn fru V in A V m n fruw loin A in with fruu fruv and fruwl given by 3 14 Version 1 0 October 25 1996 25 Technical documentation of WAQAD e The differentiated continuity equation in the second stage for mn is 1 45 At mn 7 AGinn 1 i 1 ES Jmn Enn 2 Una Enn AGinn A nas 3 16 1 xn 15 Eno AUS in Um 1n Exi ACh 1n A rn 1 SUIS AVin 2 m in Kmn Aban Ab mn 1 X 7 H
135. r n 1 Km n 1 AV m m1 An j j 2 Vinn Kmn 1 AG nt AGinn The differentiated v momentum equation in the second stage for v Jm n i e 1 ave EAM C E Oa ins fm fo AUS IE x 3 17 En An n an n n Ac frvur Aur FIV Vinn AVmn frvwil AC with frvu frvv and frvwl given by 3 11 but at time t At instead of time t Step 3 The adjoint equations The coupled adjoint equations in the first stage from t At gt t ZAt The WAQAD variables at time t At are known All adjoint u velocities and adjoint water levels at time t 2At are calculated such that 26 The 2D adjoint model the adjoint continuity equation is 1 F n XA ven 8 vena 4 g vu m in T viria E vu id vila AC Kin t n Kn 1 20m vas vean van 3 18 AG Kin En i 1 vsus rn _ va ve fan Kun Ein Kin En fruw as vu m in vu the adjoint u momentum equation is 1 XA vu E vun Hn vena vena vein Ve metn 9 19 AE Kin En Kin Bein vi 5 170 e The adjoint v momentum equation in the first stage The WAQAD variables at time t At and the adjoint u velocities and adjoint water levels at time t ZAt are known All adjoint v velocities at time t 72At are calculated such that Tl so 1 9 20 5
136. rate Verlaan 1994 describes in detail extending WAQAD to use velocity measurements However this extension has not yet been implemented in WAQAD Research is currently in progress on using velocity measurements in a separate test version of WAQAD see Mooij 1996 156 5 6 5 6 1 Version 1 0 October 25 1996 Possibilities options and shortcomings Friction formulations Introduction The friction term is an important part of the two and three dimensional shallow water equations It is very indirectly related to physical concepts Very briefly the origin of this term can be described as follows Basically the shallow water equations as approximated by WAQUA in SIMONA originate from Newton s law F MA and the conservation of mass In the presence of viscosity the flow becomes turbulent for higher Reynold numbers Since finite difference methods as used in WAQUA in SIMONA can only compute smooth solutions a direct solution of the turbulent flow is impossible Often one is only interested in the average flow and not in the turbulent fluctuations In this case it is possible to average the equation over a time interval to obtain time averaged flow In the process of averaging new terms are created that represent the turbulent energy To obtain a closed set of equations an assumption about this energy has to be made Often this is in the form of diffusion approximation To obtain two dimensional depth averaged equat
137. reases explosively when the number of parameters increases Because of the measurement noise the cost function will differ slightly from its theoretical minimum value see figure 5 3 It is obvious that even a small measurement noise will lead to a large variation in the minimum point in the case of an problem The solution estimation will vary especially in the direction of the descent since a slight vertical change will have a large effect on the minimum because of the measurement noise This increases the noise in this direction figure 5 3 Change in the cost function for an I problem due to noise I problem 2 two depth boxes between two measurement stations Version 1 0 October 25 1996 Possibilities options and shortcomings Consider the situation of I problem 1 In this case there is a second water level station at x x We assume that x20 measurements are made at the boundary which means that the phase at the boundary is supposed to be known The depth is divided into two boxes one box from Xm tO Xbetween NO measurements are made at Xpetween and one box from Xpetween tO Xn Xm lt Xbetween lt Xn The phase of the signal at x is now determined by the phase at the boundary known to a greater or lesser extent and the transit time from the boundary to xn The transit time is the sum of the transit times in the two boxes The phase in x is thus o O Xbetween 7 Xm OX K etwesri n 0 D 9 J4D
138. red quite well while bottom friction measurements are difficult to obtain This also implies that larger adaptations are acceptable for bottom friction values than for depth values By adjusting the penalty all this information can be used The disadvantage of using a penalty is that one or more parameters for this penalty need to be specified Often these settings can only be obtained by experimenting with different values and selecting a value that gives reasonable results The method described here however requires parameters to be in a form that can easily be interpreted and the values can often be guessed instantly from experience or information about the accuracy of the data used The arguments in this introduction show that the penalty can be useful for solving a number of problems However not all problems are solved automatically Some can only be solved by using more data or data at better chosen points in space or at better times Theoretical background The cost criterion used in WAQAD is quadratic and can be viewed as a Maximum Likelihood estimator In this framework the addition of a penalty is simply the transition from a Maximum Likelihood esitmator to a Bayesian estimator The prior probability distribution needed for the Bayesian estimator can under the assumption that it is Gaussian be characterized by its mean and covariance This covariance is computed with a model to assist the user in building this matrix
139. ry equals zero Therefore the v momentum equation does not contribute to the boundary The continuity equation and the u momentum equation do contribute Consider point m n to be on the boundary The contributions to the boundary are as follows e Adjoint equation for AC Contribution of the continuity equation in the first stage from t gt t ZAt in m 1 n E venin 2 Enn Kaen AG Contribution of the u momentum equation in the first stage in m n 8 QI nuwt vu mn Ki AE 2 m n Contribution of the u momentum equation in the second stage from t 72At t At in m n 9 vss Ki AE 2 m n Contribution of the continuity equation at the second stage in m 1 n UM Umn Ena m 1 n 2 Khen En AG Version 1 0 October 25 1996 The 2D adjoint model Contribution of the WAQUA boundary condition By adding these contributions including the contribution of the WAQUA boundary condition and rearranging the variables the adjoint boundary condition is derived 1 f vna Ga efr IE Gu I 1 2 Uma ven AE Kien Bisin e Adjoint equation for A Cmn ve mn 0 Right hand open boundary As for the left hand open boundary the v momentum equation does not contribute to the boundary The continuity equation and the u momentum equation do Consider point m n to be on the boundary The contributions to the boundar
140. s v NE 164 Vini 191 Technical documentation of WAQAD e The coefficients Am n and can be determined from Table XI e One part of the righthand side D is prescribed by Table XII These elements have to be subtracted e The boundary treatment of the second operators S and S can be read from Table VIII They appear in the formula for Ds e One part of the righthand side Dmn is prescribed by Table IX These elements have to be subtracted e Now we deal with the iterative aspect caused by the addition of equation 3 111 Table X shows the degenerations in case of impeded geometry The arrays GEBVUY and GETVUY Table VIII have to be calculated for the flow field at time t The elements corresponding to u point m n contribute to the main diagonal of the tridiagonal matrix equation The coefficients have to be added to Bm n dependent on the geometrical situation The other elements in Table X contain variables which are at right angles with the ones which set up a tridiagonal structure The iterative aspect use values from the previous loop at time t 72At in the first iteration This way these elements have to be subtracted from the righthand side Dm n Applying the double sweep algorithm for solving the tridiagonal matrix equation yields a first guess for all adjoint u variables at time t The same strategy can be applied in the second iteration All coefficients remain the same except the other elemen
141. s implemented in the 3D modules AD SUV3 AD SVU3 AD VYD3 and AD UXD3 according to subsection 3 3 3 In the lowermost layer of the 3D model bottom roughness is assumed The same approach is applied to the 3D case as in the 2D case This means that the bottom friction terms for the 3D model are calculated in the same way as for the 2D model The bottom friction contributions to the u and v momentum equations are calculated in a similar way by equations 3 11 and 3 14 uy tery 3d g Uk vie Ci hk ux Gi ad _ 1 Uk Vk M hr Mae Oa quy 8 2 vK hk U P v However the bottom friction contribution to the adjoint continuity equation is calculated by the last two terms of coefficient Dm of equation 4 13 Version 1 0 October 25 1996 The 3D adjoint model g uk luk lg frwl C h 5 n JG uz ni The term frwl is calculated at u points as well as v points every half time step In the term frwl the shift of half a time step with respect to TRIWAQ is corrected 117 Technical documentation of WAQAD 118 4 4 4 4 1 The gradient Once the adjoint state has been determined the most difficult task of the inverse modeling technique has been accomplished The gradient of the error function can now be computed easily with the adjoint solution In the following subsections the expression of the gradient is derived for the following three ca
142. s is interpolated to half time steps 123 Version 1 0 October 25 1996 Possibilities options and shortcomings Possibilities options and shortcomings In this chapter the possibilities options and shortcomings of WAQAD regarding the following subjects are discussed Identifiability and the use of the penalty section 5 1 The use of weighting functions for parametrization section 5 2 Coupling harmonic constituents section 5 3 Analysis of the gradient in time and space section 5 4 The use of current velocity measurements section 5 5 Friction formulations section 5 6 Finite difference calculations section 5 7 Velocity and discharge boundaries section 5 8 125 Technical documentation of WAQAD 5 1 Identifiability and the use of the penalty 5 1 1 Introduction When a model is calibrated automatically it is possible for several combinations of parameters to result in the same minimal value of the cost function This is called an identifiability problem I problem In such cases it is impossible to find the right values of the parameters because there are not enough observations or the observations are not appropriate to estimate the parameters In general one can indicate locations at which observations could be obtained to solve the I problem However it is not practical to try to solve every l problem with extra data Instead of solving the I problem the input that is the description of the para
143. s solely a function of model parameters 3 2 4 Dryfall In this subsection the drying and flooding procedures in WAQUA are described first Then the adjoint drying and flooding procedures are derived For more information on drying and flooding during tidal computations see Mouthaan 1992 and Stelling 1986 Drying and flooding procedures in WAQUA WAQUA probably does not calculate the water flow very accurately in parts of estuaries with a very low water level Drying procedures have been developed to prevent calculations with physically unrealistic values Several procedures for drying and flooding of tidal flats are applied in WAQUA in order to solve the shallow water equations correctly These procedures decide which grid points are included in the next calculation The decisions are made on the basis of the calculated water level elevations The criterion is the local water depth Several definitions of the local water depth are possible in the case of a water level point In the case of a velocity point there is only one definition Checks of local water depths at velocity points are always performed in the calculation Th local waterdepth must be positive as otherwise the continuity equation and the momentum equations will not be solved correctly In addition itis also possible to check local water depths at water level points This is important especially in the case of transport simulations The drying and flooding proc
144. s state variables After this the drying procedures do not check the u velocity points any more since the new water level elevations have not yet been calculated All these procedures have one thing in common the decision about which points should be involved in the next calculation step is made on the basis of the water depths If points are excluded the velocity points are made stagnant These velocity points are regarded by WAQUA as closed boundary points with zero values If they are re incorporated in the calculation they are no longer considered to be closed boundary points The alternating drying and flooding of the tidal flats means that moving closed boundaries can be created in WAQUA The adjoint drying and flooding procedures The dynamic WAQUA equations are important for deriving the adjoint equations The starting point for deducing the adjoint variables is the way that WAQUA calculates new state variables at a certain point of time If drying occurs in WAQUA the adjoint equations at the corresponding time steps must be modified The whole process of drying and flooding in WAQUA can be reduced to the following essential property either the velocities are calculated by the normal momentum equations or they are determined by imposing zero values As a consequence of this the exact formulations of the different drying and flooding procedures are not important to WAQAD Only the time steps during which velocities points are trea
145. ses e the parameter to be estimated is a correction factor for the bottom friction coefficient subsection 4 4 1 e the parameter to be estimated is a correction factor for the depth subsection 4 4 2 e the parameter to be estimated is a correction factor for the vertical viscosity coefficient subsection 4 4 3 Bottom friction The differential of an additive control parameter defined for Chezy in the first stage is given by as LS 2g viet FEGE E v 2g m Jj v vk m Tu hk The differential of an additive control parameter defined for Chezy in the second stage is described by 2 vef Pt v 22 ue V ux PHOS V k C hb C m For bottom friction in case of Manning and White Colebrook see section 5 6 The 3D adjoint model 4 4 2 Depth The differential of an additive control parameter defined for bathymetry in the first stage is given by n Jue KE oe K E Jo K a k 1 u i m Vy u Pp ox i vu X h hL h y ES f 5 7 n vy 447 U Ay 2h pos ox SOF a 15 l Deca o Sethe NU Pea 5 h h thi 1 K vi v tn oid vy 5 k 2 E hi4 h K n y ny x m Vy k 1 2h pa oxl wv hh h Y m Vy vy Via Y 2h hea tor vy X h y h hives Site oi 4l v v
146. shifted over the same phase Ad that is Gt A Sin a t gt Ad G t A sin ost Ad 2 both constituents are shifted over the same period T that is C t Arsin a t T d Co t A sin o t T gt The second method is chosen for WAQAD for in this way the calculated signal will stay the same with respect to its shape The gradients of the phase are then calculated by 0J gt oJ OQ ti ab alti 06 2 mo COS ox ti T 0J Gf md rr Bett A COS cs t T 2 and the gradient to T UN este bolt oT m OC ti OT oT D 5 A sin o t T 4 Z AcSintos t T 6 mem Aim COL t T A29 cogdozlti T From this it follows that 0J 0J 0J t o aT 0 06 Both constituents can thus be summarized with weights en for and C2 respectively 147 Technical documentation of WAQAD 148 The adaptations of the phases are applied to Ag and Ag in the following way AQ AT and Ads 2AT It is useful and numerically desirable to normalize T with such that d aT wis described by the mean of o and which gives aJ _ vein I a Ob PUES t TA the weights are therefore 5 and 2 they will be e approximately 1 The adaptations in this case are AQ 7 2 Ad Ad 2 Ad with the new parameter Ad Gradient of the amplitude The amplitude is calibra
147. solution method is proposed Version 1 0 October 25 1996 189 Technical documentation of WAQAD 190 To determine the deformations caused by impeded geometry for equation 3 107 an array GEOUUX has to be determined Table I from the variables KHU at time t 72At The coefficients Am n and Cmn from the corresponding tridiagonal matrix equation can be read from Table XI GEOUUX Table XI The degenerations in equation 3 107 For another contribution to the right side the array GEOUUX has to be constructed from the variables KHU at time t 7 At Table GEOUUX tomal Table XII The degenerations in equation 3 107 Version 1 0 October 25 1996 Appendices For convenience we give again the adjoint u impulse equation associated with Au including the advective terms neglecting the boundary treatment and the dryfall procedure and the iterative aspect Ann vu Jaen Bs Jmn Crh vu metn Dm n 7 33 with Am n Un 7 34 Bn n I 7 35 Ds Uns 7 36 2 i Dm n At vu Jui E fin z frvUm nl v Jon vy m n 7 37 REA GU ieee e Jain 2AX 7 38 AU su vu Jig 12A Y vu n n 12A Y AV m n 1 Vn n 2 ae ee e 4 vu Jama 12 Y 12 Y 7 39 ee 18 Vinn vy Vrain 7 M KP 7 40 w MET Vmin S Vu macs a Vinna p vy Vian zs TA vy Vai 7 41 B v m 1 n 1 i
148. t function J has to be determined This search is performed in a forward or backward direction depending on the value of the cost function J computed after the first estimation When the value of the new computed J is larger than the original J the optimum value is in between and the search is performed backwards Otherwise the search is performed forwards The forward search is performed with equal stepsizes defined by the parameter parcor that is the initial parameter correction step size The backward search is performed by reducing the previous stepsize to half its size Suppose that cost functions J4 Jo Jn are computed during one iterative search process the optimal step is then computed as follows e if J2 J4 the search is performed forwards The minimum value can be computed as soon as Jn gt Jn 1 with n22 eR step parcor dha h 2h The 2D adjoint model 3 The 2D adjoint model This chapter describes the adjoint model for two dimensional shallow water flow in detail The adjoint model equations are derived in section 3 2 The adjoint solver as it is implemented in WAQAD is discussed in section 3 3 The derivation of the gradient is given in section 3 4 An overview of remarks on WAQAD with respect to WAQUA is shown in section 3 5 A description of the three dimensional adjoint model will be given in chapter 4 Version 1 0 October 25 1996 13 Technical documentation of WAQAD 3 1 Introduction
149. t such a procedure could conceal the I problem 129 Technical documentation of WAQAD 130 For example look at the cost function in figure 5 2 and consider the problem of calibrating parameter p only Figure 5 4 shows the corresponding cost function with a clear minimum Minimizing of the second parameter figure 5 5 using the minimum for the first parameter results in a clear minimum again and the whole I problem seems to have disappeared However this is not so the problem is still there The advantage is that the observation errors are less amplified resulting in adaptations which will not be unrealistically large However it is not possible to recognize an l problem right away unless the calibration process is iterative after calibrating parameter p parameter p is calibrated again and so on Figure 5 6 shows the progress of the iterative procedure 0 01 002 003 004 005 006 007 0 08 0 09 0 1 p2 figure 5 5 Cost for parameter 2 only Version 1 0 October 25 1996 Possibilities options and shortcomings figure 5 6 Progress of the iterative procedure Validation of the calculation process Since it is impossible to determine beforehand if an I problem will occur it would be practical to ascertain if all is still going well during the optimization process the execution of the calculation Two methods seem feasible at this moment Using two data sets If two data sets are available
150. te models automatically are referred to the WAQAD User s Guide Brouwer 1996b Persons starting to learn about the calibration process are advised to read some other books on this subject before reading this document For example the following documents are recommended Inregelen van wiskundige modellen op basis van besturingstheorie Van den Boogaard 1988 e Data assimilation with altimetry techniques used in the Continental Shelf Model RIKZ 1994 e User s guide WAQAD Brouwer 1996b Diana Rozendaal SIMTECH AUTOMATISERING BV Contents Contents 1 INTRODUCTION pn 1 2 THE ADJOINT 3 2 1 INTRODUCTION itn eren d a aaa paani 4 2 2 THE MATHEMATICAL 6 3 THE 2D ADJOINT MODEL eee 13 8 1 INTRODUCTION 3 cuss opes eps ede ete teg 14 3 2 THE ADJOINT MODEL EQUATIONS cere 15 3 2 1 Derivation of the adjoint model equations 15 9 2 2 AGVECTOMN sey 29 3 2 2 1 Introduction eines 29 3 2 2 2 Advective terms in 29 3 2 2 3 Adjoint Advective 33 3 2 2 4 Curvilinear and spherical coordinates 57 3 2 3 Boundary conditions eese 60 3 2 4 Dryfall
151. ted as closed boundary conditions are important A universal adjoint drying and flooding procedure can be based on this characteristic The adjoint equations derived without the drying and flooding procedure are the starting point of the next analysis Suppose that at the first stage from t t 72 t velocity point Um is no longer involved in the WAQUA calculation due to an drying procedure which may be arbitrary The following can be derived e Unn 0 instead of being calculated by the u momentum equation The 2D adjoint model Version 1 0 October 25 1996 The equation u 0 does not explicitly contain any m n unknown model parameter The contribution to the calculation of the gradient is therefore equal to zero irrespective of the fact that this constraint is modelled as a known or unknown boundary condition in the reverse calculation process vuln disappears from all the adjoint equations therefore Vu mn has a zero value in all these equations The universal adjoint drying procedure is therefore similar to the WAQUA drying procedure Vu mn 0 instead of being calculated by the adjoint u momentum equation The preceding reasoning is in fact the adjoint approach for velocity boundaries in WAQUA that is boundary conditions in WAQUA correspond with homogeneous boundary conditions in WAQAD In this situation WAQUA state variables at boundary points are also determined by assigning the variable a certain
152. ted by calibrating variables and amp in the following way C t C t A1 1 amp sin ost gt A 1 amp sin ost 5 5 3 3 Version 1 0 October 25 1996 Possibilities options and shortcomings The gradients of the amplitude are LEE NC 0 amp TT OC ti 08 S 0b 64 sinat Man A sin ct 0J WN a E Haein Pene be If 5 and constituents C en C are coupled amp 6 holds The gradient to is equal to WS x J oy er A i A Sir oxt 6 A2sir ast from which it follows that 8 a 06 b 0 amp The gradients for the percentual adaptation of the amplitude can thus be summarized without using weights The adaptations are AG AE Implementation The constituents are coupled in WAQAD as follows e If the test defined by 5 22 fails for the constituents to be coupled the preprocessor of the program gives a warning e Regarding the existing calculation of the gradient an extra step is implemented in order to be able to summarize some gradients by using weights Just before adaptation the coupled variables amp and are decoupled to amp i amp 2 1 62 Then the adaptation is applied to these decoupled variables For more information on the implementation of coupling constituents see Rozendaal 1995 149 Technical documentation of WAQAD 5 4 Analysis of the gradient in time and space
153. the characteristic distance of the error Let Q be a parameter that depends on the spatial variables x and y that take only discrete values Assuming that all these stochastic variables are jointly Gaussian with expectation 0 then the probability density is fully described by the variances of Q x y the correlations between any two points and the zero mean We model the correlation of the error in Q denoted by AQ between the points x1 y1 and X2 y2 as xrx2Y y4 y2 cor AQ x y Q Q y e a 5 10 137 Technical documentation of WAQAD 138 The parameter d serves as a Characteristic distance see figure 5 10 This distance indicates over what order of length the correlation is significant The method used for measurements for example of depths may indicate reasonable values for this In our experience one tenth of the model diameter is a good default when no accurate information is available When other coordinates are used the distance in the numerator of the exponent of equation 5 5 may be defined differently 2 distance characteristic distance figure 5 10 Correlation and characteristic distance The correlation in the a priori errors is important for two reasons The errors are not considered to cancel each other out when averaged over an area and the resulting penalty is large when two points close together are adapted in different directions The parameters i that are used by the program can be
154. tion standards to make maintenance easier This documentation contains a complete description of the mathematical and technical aspects of WAQAD Chapter 2 contains a general description of the adjoint method on which WAQAD is based In chapter 3 the adjoint model for two dimensional shallow water flow is described in detail the adjoint model equations are derived the adjoint solver as it is implemented in WAQAD is discussed and the derivation of the gradient is given At the end of chapter 3 an overview of the remarks on WAQAD with respect to WAQUA is given Chapter 4 focuses on the adjoint method for three dimensional shallow water flow Chapter 5 comprises a discussion of the potential options and constraints of WAQAD in relation to identifiability and the use of the penalty the use of weighting functions for parametrization coupling harmonic constituents analysis of the gradient in time and space the use of current velocity measurements friction formulations finite difference calculations and velocity and discharge boundaries Finally an overview of recommendations on WAQAD is given in chapter 6 Version 1 0 October 25 1996 The adjoint method The adjoint method The adjoint method is described in general in this chapter It is outlined in section 2 1 and illustrated by means of the mathematical description in section 2 2 Technical documentation of WAQAD 2 1 Introduction The adjoint method is very important in
155. tly WAQAD was only able to recognize closed and open water level boundaries Open velocity and discharge boundaries were treated as open water level boundaries Because of this the velocities prescribed by WAQUA at open boundaries were treated in the adjoint model equations as if they were velocities in the interior field This wrong approach has been corrected in WAQAD in order to gain completeness Thus open boundaries are now treated by WAQAD in the following way e lf an open boundary is a water level boundary then the water level points at this boundary do not take part in the adjoint calculation The nearest neighbour velocity points do take part in the adjoint calculation fan open boundary is not a water level boundary but a velocity or discharge boundary then the water level points and the velocity points at the boundary do not take part in the adjoint calculation For more information on the implementation of velocity and discharge boundaries see Brouwer 1995b Version 1 0 October 25 1996 Recommendations Recommendations Remarks on WAQAD with respect to WAQUA and TRIWAQ were made in chapters 3 and 4 respectively An overview of these remarks is listed below The WAQAD bottom friction terms contain a shift of half a time step compared with the bottom friction terms in WAQUA e Horizontal viscosity is not included in WAQAD both for 2D and 3D models e Horizontal advection is not included in WAQAD in th
156. ts in Table X This first guess is used for the adjoint u variable at time t Applying the double sweep algorithm results in a second guess It is expected that two iterations are enough for convergence 7 4 2 Open Boundaries Open boundaries are introduced to restrict the size of the domain of the problem The numerical boundary procedures for an open water level boundary are described by Stelling 1984 The solution at the inner points is not greatly influenced by the order of the advection discretization near the open boundary Stelling 1984 At inflow uu is approximated by a zero value otherwise first order upwind The influence on the inner points of this approximation is negligible as was found by practical experiments Stelling 1984 Therefore it is 192 Appendices proposed to neglect the boundary treatment of the terms uu and vvy in the adjoint computation process Example if m n is a left water level boundary point then GEOUUX m n and GEOUUX m 41 n equal zero The treatment of GEOVVY near a water level boundary is similar 7 4 3 Spherical and curvilinear coordinates The boundary treatment and the degenerations caused by impeded points is similar to the case for rectangular WAQUA For curvilinear models they are Version 1 0 October 25 1996 In all equations AX must be replaced by A amp and AY by Am The variables AU1 and AU1D have to be multiplied by d qa 7 9 Kin n The definition of t
157. ts the array GEOVVY Table V at time t vy ar DDm n1 CG neal Jn E Dinan 7 26 with m n CCh net 2 AY m n t AV x frVV m 2 2 fwn AVD ic DDm n 1 2 mnt AV E frvv m Hina ve Jmn ve 2 RAV x frvv m f GEOVVY m n 1 AV Vm n 17 Vm n 1 7 27 2AY Nb aa g Vinn AVD 2AY If GEOVVY m n 2 e iac Via if van gt 0 7 28 AY 0 otherwise Vis V AVD Vian Vana if vmn gt 0 AY i 0 otherwise Version 1 0 October 25 1996 185 Technical documentation of WAQAD 186 If GEOVVY m n 3 Vin mei Vink T Nimo Vma tz lt Q 7 29 AY QO otherwise ssi Vna if Vm n lt O AY 0 otherwise AVD So far we handled the degenerations in the adjoint advective terms caused by impeded points in the adjoint v impulse equation Degenerations in the adjoint u impulse equation associated with A um caused by impeded geometry In the adjoint u impulse equation associated with Aum the second order difference operators S central and S upwind appear Here the boundary treatment of these operators in the adjoint computation process is indicated These approximations depend on the dry wet status of the surrounding points Therefore the time dependent local flow situation has to be determined first to figure out which approximation
158. u velocity point v o depth point d m axis figure 3 1 The shaded points in the figure all have indices m n For an interior point in the field see figure 3 1 not enclosed by any impeded geometry the discretized two dimensional shallow water flow is described by The v momentum equation in the first stage from t gt t AAt This is the first step of the ADI method in v direction The WAQUA variables at time t are known All v velocities at time t At are calculated such that n V mn 7 Vmin fmn Um n 3 4 g Gnin g Vmn us J Vmn r 0 Ean An in which Vm v velocity at time t M mn v velocity at time t 72At D is water depth at v velocity point Ve dm 1 n dmn En transformation coefficient E at v velocity point Ya Em 1ntEmntEmn1tEm 1 n 1 ds space dependent Coriolis force C Chezy coefficient at v velocity point Version 1 0 October 25 1996 19 Technical documentation of WAQAD 20 and where the following terms are interpolated to v velocity points i z VA Um nt Um n 1FUm 1 n 1 Um 1 n n E es n Ve m nt mn n i fm n Um n x V4 fm nUm n t fm n 4Um net fm n 1U m1 n1 m1 nU m 1 0 The coupled equations in the first stage This is the first step of the ADI method in u direction The WAQUA variables at time t and the v velocities at time t 2 t are known All u velocities and water levels at time t At are calcul
159. ual to zero In the lowermost layer a term associated with bottom roughnouss is defined the two dimensional approach is used for the adjoint momentum equations in the bottom layer Here on the main diagonal part of the vertical viscosity term still remains The right hand side of this equation 4 8 is independent of vu so the tridiagonal matrix A can be inverted which formally leads to Qm at at yy Umi at e atd 9 A 3 m 111 Technical documentation of WAQAD 112 In practice it is not necesary to compute the whole matrix A the three tridiagonal systems can be solved by a double sweep This strategy decouples the adjoint u velocity points in the vertical Using equation 4 9 again the momentum equation per layer can be derived m 4 1 AAK K pit nct CCK Ine kgs a 0 in which AAK k th element of A a k th element of A c DDK k th element of A d These equations have to be added up per layer from bottom to top in order to eliminate the adjoint u velocities at the implicit time level in the adjoint continuity equation The vertically summed adjoint u momentum equation is then gt AAK v mk K 4 11 ane CCK k7 DDK gt ns gt ha which will be noted as vna AA m CO DD 4 12 j K E v Jmk js K E n in which AA
160. uated bottom at D figure 5 1 A simple problem I problem 1 phase at the boundary and depth Suppose there is only one water level station in X Xm At this station a sinusoidal signal obscured by noise is measured Equation 5 2 shows that the phase of this signal is determined solely by the phase at the boundary and the transit velocity of the wave that is the depth The amplitude of this signal is also dependent on the depth but not on the phase at the boundary This means that in principle the depth can be estimated in this way when the friction is known If the friction is unknown or not accurately known then every combination of the depth and the phase at the boundary that leads to the correct phase at the measurement point will minimize the cost function This means that there is an l problem Figure 5 2 shows the progress of the cost function schematically In general an l problem can be recognized from a graph if it shows a long drawn out decline with at least one direction in which the cost hardly varies Unfortunately in practice it is not possible to produce such a picture for several reasons For example evalution of the cost function at one point requires much computational time Performing the dozen evaluations necessary in the case of two dimensional 127 Technical documentation of WAQAD 128 pictures is already stretching the limit of computational possibilities And the number of evalutions required inc
161. ubsection See appendix A for more information on this subject 3 2 3 Boundary conditions In this subsection the adjoint boundary conditions are deduced It is only necessary to solve these adjoint boundary conditions is only necessary when estimating uncertain parameters in the WAQUA boundary conditions We will begin by describing the WAQUA boundary conditions In order to obtain a unique solution to the WAQUA model equations 3 1 3 3 a set of boundary conditions is prescribed Two types of boundaries are considered open and closed At closed boundaries such as coastlines and dams the boundary condition yields 60 Version 1 0 October 25 1996 The 2D adjoint model v 0 This means that no inflow or outflow can occur through these boundaries Open boundaries are non physical mathematical boundaries that are to delimit the domain of the model In the case of open boundaries two boundary conditions are required In WAQUA the first requirement is fulfilled by placing velocity blockades at the open boundary This guarantees that the velocity parallel to the boundary has a zero value The second requirement is satisfied by imposing a water level boundary condition in an explicit way This boundary condition consists of a series of cosines of harmonic constituents C t YA COS wit i 8 130 in which water level elevation at time t C t Co averaged sea surface n number of harmonic constituents amp
162. ule AD GRKB the boundary conditions described by section 3 2 2 are solved first Then the gradient for the amplitude and phase is calculated 99 Technical documentation of WAQAD 100 3 5 Remarks on WAQAD with respect to WAQUA The remarks on WAQAD with respect to WAQUA are The WAQAD bottom friction terms fruu fruv fruwl described by equation 3 14 and the terms frvu frvv frvwl described by equation 3 11 contain a shift of half a time step compared with the bottom friction terms in WAQUA The horizontal advection terms are solved in WAQUA using the Gauss Seidel method In WAQAD these horizontal advection terms are solved using a Jacobi method Horizontal viscosity is not included in WAQAD but it is in WAQUA Additional features of WAQUA which are not mentioned in Stelling 1984 are not included in WAQAD This may lead to problems For example External forcings apart from open boundary values Forcings like wind discharges and density gradient influence the results of the WAQUA computations WAQAD does not take these influences into account which may lead to an unreliable gradient Transports STRESS2d barriers and weirs are taken into account by WAQUA but not by WAQAD Smoothing boundary values WAQUA has the option to smooth the boundary values during the play in period This option is not supported by WAQAD Depth multiplier In the simulation input file a depth multiplier can be sp
163. ustments of the model input The more accurate the estimates of the parameters the better the agreement between computed and observed water levels The adjoint method enables models to be calibrated automatically and efficiently Different parameters can be estimated simultaneously in an optimal way Using the adjoint method an error function is introduced to quantify the misfit of the model outcome This cost function is minimized by using an iterative algorithm in each step of the iteration an improved vector of control variables is searched for in a direction in the parameter space computed from the gradient of the cost function The adjoint model computes this gradient vector efficiently The adjoint method guarantees complete consistency with the dynamics RIKZ developed WAQAD the automatic calibration program based on the adjoint method The program calibrates two dimensional shallow water flow models WAQUA by estimating depth bottom friction or boundary condition Technical documentation of WAQAD parameters WAQAD can also calibrate three dimensional shallow water flow models TRIWAQ by estimating depth bottom friction or vertical viscosity parameters WAQAD has been implemented in a SIMONA environment One of the goals of the SIMONA project is to develop a generic framework in which it is easy to make connections between various mathematical modelling tools The project also produces programming standards and documenta
164. w v Viet 1 VR Version 1 0 October 25 1996 119 Technical documentation of WAQAD 120 The 3D adjoint model The differential of an additive control parameter defined for bathymetry in the second stage is described by n 6 fuf q K E b K E v H Meh h k 1 T hy K ay i S 17 oxl vu y Ua d Mr z 2l vn ox 2 Ok vi Lu N 1 h h qe ce Sew n JS y vv used Ok Ok 1 vu Xh h hy Y an K w Va V T V PS oi idC vv oh Sous ont es ede k 2 y hi h y JS lv y vy V v 21 hi k V RM A h Y h hia 6 K 1 VEY v vy i e Yoh h th ay y Bleed 1 VU C2 Lox 1 vk 6 9 y C hk Y Version 1 0 October 25 1996 121 Technical documentation of WAQAD 4 4 3 122 Vertical viscosity The differential of an additive control parameter which is time independent defined for eddy viscosity in the first stage is given by F K 1 5 yo u uj 1 zn e TE iM es vv Via V V PT hie thy n w v Via hy h a k 2 The differential of an additive control parameter which is time indepen
165. y Ea m n m n Gili XC Dnt e Contribution of the WAQUA u momentum equation in the first stage from t gt t 72At at point m n 1 v 3 g CMT y T os u m n 1 d 2 Cna y 8 y m n 1 Contribution of the WAQUA continuity equation in the first stage from t gt t ZAt at point m n ve mn 2s Enn uT Kmn Kis Ehn AG An Contribution of the WAQUA continuity equation in the first stage from t gt t 72At at point m 1 n venta T 205 Enn Km 1 n Ban AG An Contribution of the WAQUA continuity equation in the first stage from t gt t ZAt at point m 1 n 1 1 vem 2 Unni Enna 72 Vm 1 n Km 1 n Benet AG An 95 Technical documentation of WAQAD 96 Contribution of the WAQUA continuity equation in the first stage from t gt t At at point m n 1 1 vent 2 Um n4 Enn 735 Vmn Km n Kinet ant AG An Contribution of the WAQUA u momentum equation in the second stage from t 72At gt t At at point m n ge ete Gal Gd Fi ums u X2 u PrE 2 Cmn Dinn t Conn Contribution of the WAQUA u momentum equation in the second stage from t 72At gt t At at point m n 1 D g TAN Came y uis i Zl y Dm T y vila Contribution of the WAQUA v momentum equation in the second stage from t 72At gt t At at point m n 9 Nirn ms inn ME 2 Cma Dhat Cmn Contribution of the WAQUA v momentum
166. y are as follows 63 Technical documentation of WAQAD 64 Adjoint equation for AG Contribution of the continuity equation in the first stage from t gt t 2At in m 1 n 1 vean m 1 n En 1 n En AG Contribution of the u momentum equation in the first stage in m 1 n i a E I of Kiss AE 2 m 1 n Contribution of the u momentum equation in the second stage from t 7 At t At in m 1 n 4 A I on n Ki 2 m 1 n Contribution of the continuity equation in the second stage in m 1 n d U vein gum Ent Kin Ein Ag Contribution of the WAQUA boundary condition va By adding the contributions and rearranging the variables the adjoint boundary condition yields vo 9 fruwll AE 2 Vu in tn T m 1 n Jide piae cae AG Exin The 2D adjoint model Version 1 0 October 25 1996 e Adjoint equation for A Cmn Jui 0 Upper open boundary The velocity parallel to the boundary equals zero Therefore the u momentum equation does not contribute to the boundary The continuity equation and the v momentum equation do contribute Consider point m n to be on the boundary The contributions to the boundary are as follows e Adjoint equation for AG vn 0 e Adjoint equation for A Cmn Contribution of the continuity equation in the second stage from t 2At t in m n 1 vc pom
167. y p C p It can be shown that this penalty is always positive So instead of minimizing the cost we now try to minimize the total cost which means cost plus penalty If the prior is taken uniform over p independent of p the Bayesian estimator reduces to the Maximum Likelihood estimator This can be interpreted as there being no prior information on parameters or all parameters are equally probable before the results of the simulation are known The routines for using penalty functions also provide some help in generating the covariance matrix C The model used for parameters that in principle could vary for all x and y coordinates is based on the spatial distribution of the parameters These parameters for example the values for depth and bottom friction are parametrized using rectangles or triangles But this is just for the program For the original data we can characterize the error structure by the variance of the parameters at a point and the correlation between two points This correlation to what extent the points depend on each other decreases with increasing distance between the points Thus if we measure the depth somewhere near Scheveningen with an error of one metre too deep we cannot say anything about the error of a measurement obtained somewhere near Bergen Norway but a depth measurement obtained near Hoek van Holland will probably be to deep too The range within which two measurements influence each other is called
168. y the same value In this way only the line P1 P2 is searched Figure 5 9 shows the resulting cost function The l problem is solved in this way The problem is solved in a mathematical sense but it is of course extremely important that the constant adaptation over the large fused box is a reasonable assumption 0 5 0 6 0 7 0 8 0 9 1 pizp2 figure 5 9 Cost function for the joined parameters Penalty for deviation from prior estimates Introduction In some cases the standard quadratic cost does not lead to acceptable results In such cases the adaptations of some parameters are unreasonably large while others are quite acceptable This is not a feature of the algorithms used but it is inherent to the problem The adaptations of the parameters are based solely on the values of the waterlevels measured at certain points in the model Sometimes it is attempted to adapt a parameter that 133 Technical documentation of WAQAD 134 cannot be estimated well from the data For example it is clearly not possible to estimate the depth of the bottom near Denmark from data on the English and Dutch coasts since changing the depth at the Danish coast has hardly any influence on water levels at the Dutch coast Data obtained near the area in question are needed to solve this problem Alternatively the estimation of the depth near Denmark can be omitted This will make the other estimations behave better The situation becomes m
169. zed as Um 2 n amp Vm 2n Consequently we find u 2 n T s AN if gt 0 m n 2 X Um 2 n 0 ieee 3 94 It follows from equations 3 91 to 3 94 that a part of the adjoint equation associated with looks like vv m 2n Ge if Um2n lt 0 3 95 AUm n ME vv ds i if Um 1 n s 0 F vv mn EJ if i5 05 or vw Jinn E f GeO H vy Jin EJ if p vv ier ps if 0 These adjoint v variables at time t 2At are perpendicular to the ones in the computation procedure of the first stage which set up the tridiagonal structure if the advective terms were ignored Therefore in imitation of WAQUA an iterative solution scheme is proposed a Gauss Seidel or Jacobi method This allows the equations to still be resolved by solving tridiagonal matrices over grid lines because in the adjoint cross advective part the latest computed adjoint v velocities Gauss Seidel or the ones originating from the previous iteration level Jacobi are used 51 Technical documentation of WAQAD 52 The weighted central difference of uv at the second stage looks like U mn amp x 3 96 in which 4 Umn Umn UC RE Um 4n Us 4n44 E 4 vmetn7 4 Vm4n7 Vm 2 n Sd una 12AX Applying variational formalism yields AV E a AV nn au a 12A X 12A X AU n u m n 2 m 1 n 1 2A X m 2 n 1 2A 3 97 AUmn FX

User's Guide - Applicaties Helpdesk Water

Contents

Download Pdf Manuals

Related Search

Related Contents