Home

SHYFEM Finite Element Model for Coastal Seas User Manual

1. take the value 1 and the hollow circles where they are 0 The contributions anm to the system matrix are therefore different from 0 only in elements containing node m and the evaluation of the matrix elements can be performed on an ele ment basis where all coefficients and unknowns are linear functions of x and y This approach is straightforward but not very satisfying with the semi implicit time step ping scheme for reasons explained below Therefore an other way has been followed in the present formulation The fluid domain is still divided in triangles and the water levels are still defined at the nodes of the grid and represented by piecewise linear interpolating functions in the internal of each element i e E Cmm m 1 K However the transports are now expanded over each triangle with piecewise constant non continuous form functions y over the whole domain We therefore write U UWy n 1 J where n is now running over all triangles and J is the total number of triangles An example of y is given in the lower right part of Fig la Note that the form function is constant 1 over the whole element but outside the element identically 0 Thus it is discontinuous at the element borders Since we may identify the center of gravity of the triangle with the point where the trans ports U are defined contrary to the water levels Gm which are defined on the vertices of the triangles the resulting grid may be seen as a staggered grid
2. Nitrogen and phosphorous are then returned to the organic compartment ON OP via phy toplankton and zooplankton respiration and death After mineralization the organic form is again converted into the dissolved inorganic form available for phytoplankton growth The DO mass balance is influenced by almost all of the processes going on in the system The reaeration process acts to restore the thermodynamic equilibrium level the saturation value while respirations activities and mineralization of particulated and dissolved organic matter consume DO and of course photosynthetic activity produces it Other terms in cluded in the DO mass balance are the ones referring to redox reactions such as nitrification and denitrification The reaeration rate is computed from the model in agreement with ei ther the flow induced rate or the wind induced rate whichever is larger The wind induced reaeration rate is determined as a function of wind speed water and air temperature in agreement with O Connor 1983 while the flow induced reaeration is based on the Co var method Covar 1976 i e it is calculated as a function of current velocity depth and temperature The dynamic of a generic herbivorous zooplankton compartment ZOO meant to be repre sentative of the pool of all the herbivorous zooplankton species is followed and accordingly the subroutines relative to phytoplankton organic matter nutrients and dissolved oxygen which were influenc
3. is describing linear friction or a Chezy or Strickler form of friction default 0 iczv Normally R is evaluated at every time step iczv 1 If for some reason this behavior is not desirable iczv 0 evaluates the value of R only before the first time step keeping it constant for the rest of the simulation default 1 The value of A may be specified for the whole basin through the value of czdef For more control over the friction parameter it can be also specified in section area where the friction parameter depending on the type of the element may be varied Please see the paragraph on section area for more information 15 Physical parameters The next parameters describe physical values that can be adjusted if needed rowass Average density of sea water Default 1025 kg m roluft Average density of air Default 1 225 kg m grav Average gravitational acceleration Default 9 81 ms 7 Wind parameters The next two parameters deal with the wind stress to be prescribed at the surface of the basin The wind data can either be specified in an external file ASCII or binary or directly in the parameter file in section wind The ASCII file or the wind section contain three columns the first giving the time in seconds and the others the components of the wind speed Please see below how the last two columns are interpreted depending on the value of iwt ype For the format of the binary file please see the relative section If b
4. itvd Type of advection scheme used for the transport and diffusion equation Normally an upwind scheme is used 0 but setting the parameter itvd to 1 choses a TVD scheme This feature is still experimental so use with care Default 0 16 dhpar Horizontal diffusion parameter general Default 0 chpar Horizontal diffusion parameter for the tracer Default 0 diftur Vertical turbulent diffusion parameter for the tracer Default 0 difmol Vertical molecular diffusion parameter for the tracer Default 1 0e 06 Temperature and salinity The next parameters deal with the transport and diffusion of temperature and salinity itemp Flag if the computation on the temperature is done A value different from 0 computes the transport and diffusion of the temperature De fault 0 isalt Flag if the computation on the salinity is done A value different from 0 computes the transport and diffusion of the salinity Default 0 temref Reference ambient temperature of the water in centigrade Default 0 salref Reference ambient salinity of the water in psu practical salinity units per mille Default 0 thpar Horizontal diffusion parameter for temperature Default 0 shpar Horizontal diffusion parameter for salinity Default 0 Output for concentration temperature and salinity The next parameters define the output frequency of the computed concentration temperature salinity to file They also define the internal time step to
5. 2 2 The Model The model uses the semi implicit time discretization to accomplish the time integration In the space the finite element method has been used not in its standard formulation but using staggered finite elements In the following a description of the method is given 2 2 1 Discretization in Time The Semi Implicit Method Looking for an efficient time integration method a semi implicit scheme has been chosen The semi implicit scheme combines the advantages of the explicit and the implicit scheme It is unconditionally stable for any time step At chosen and allows the two momentum equations to be solved explicitly without solving a linear system The only equation that has to be solved implicitly is the continuity equation Compared to a fully implicit solution of the shallow water equations the dimensions of the matrix are reduced to one third Since the solution of a linear system is roughly proportional to the cube of the dimension of the system the saving in computing time is approximately a factor of 30 It has to be pointed out that it is important not to be limited with the time step by the CFL criterion for the speed of the external gravity waves At lt A Vs where Ax is the minimum distance between the nodes in an element With the discretization described below in most parts of the lagoon we have Ax 500m and H 1m so At 200 sec But the limitation of the time step is determined by the worst case For example for
6. 73 74 76 353 350 349 1374 1154 1160 1161 408 409 786 795 Send paved nas old chezy values 0 Io 25e 2 74 75 2 21 23 350 346 Sh 20 lt 28 LLISA 1153 4 5 6 Send Figure 4 1 Example of a parameter input file 12 title free one line description of simulation name_of_simulation name_of_basin end Figure 4 2 Example of section title file extension RST No restart file is written with idtrst equal to 0 itrst Time to use for the restart If a restart is performed then the file name containing the restart data has to be specified in restrt and the time record corresponding to itrst is used in this file idtres Time step and start time for writing to file RES the file containing resid itmres ual hydrodynamic data idtrms Time step and start time for writing to file RMS the file containing itmrms hydrodynamic data of root mean square velocities idtflx Time step and start time for writing to file FLX the file containing dis itmflx charge data through defined sections The transects for which the dis charges are computed are given in section flux of the parameter file Model parameters The next parameters define the inclusion or exclusion of certain terms of the primitive equations ilin Linearization of the momentum equations If ilin is different from 0 the advective terms are not included in the computation Default 1 itlin This parameter decides how the advective non linear
7. terms are com puted The value of 0 default uses the usual finite element discretiza tion over a single element The value of 1 choses a semi lagrangian approach that is theoretically stable also for Courant numbers higher than 1 It is however recommended that the time step is limited using itsplt and coumax described below Default 0 iclin Linearization of the continuity equation If iclin is different from 0 the depth term in the continuity equation is taken to be constant Default 0 nrand Compute system matrix every nrand iterations This parameter can be used to speed up computations if the system matrix does not depend crucially on the varying water depth e g in deep waters It is recom mended to leave it to the default value of 1 compute at every iteration Default 1 The next parameters allow for a variable time step in the hydrodynamic computations This is especially important for the non linear model ilin 0 because in this case the criterion for stability cannot be determined a priori and in any case the time integration will not be unconditionally stable The variable time steps allows for longer basic time steps here called macro time steps which have to be set in idt It then computes the optimal time step here micro time step in order to not exceed the given Courant number However the value for the macro time step will never be exceeded 13 itsplt Type of variable time step computation If this value
8. Ax 100 m and H 40 m the time step criterion would be Ar lt 5 sec a prohibitive small value The equations 1 3 are discretized as follows ae _ Gh 4 au U i a vnt V I o 2 6 n 1__yyn d n l fn gH FEY aye 4 X 0 2 7 yrtl y 0 n l arn a te E Ryg 0 2 8 With this time discretization the friction term has been formulated fully implicit X Y fully explicit and all the other terms have been centered in time The reason for the implicit treatment of the friction term is to avoid a sign inversion in the term when the friction parameter gets too high An example of this behavior is given in Backhaus 1 If the two momentum equations are solved for the unknowns U and V we have ua 1 G SE GA n 1 n 1 urea v AtgH MX 2 9 1 a gnt Tv n 1 n 1 vert v AtgH 4 5 MY 2 10 If were known the solution for U and V could directly be given To find we insert 2 9 and 2 10 in 2 6 After some transformations 2 6 reads Mie n Dinh OC e 5 1 2 gp ZU gt 4 a gt 2 11 2 AtR oU ov Ana eR dx x oy Ar ax Y 2 1 AtR Ox x The terms on the left hand side contain the unknown C the right hand contains only known values of the old time level If the spatial derivatives are now expressed by the finite element method a linear system with the unknown is obtained and can be solved by stan
9. C ratio C P C ratio O C ratio grazing efficiency fraction of ON from algal death fraction of OP from algal death fraction of dissolved inorganic phosphorous mineralization of dissolved OP rate constant mineralization of dissolved OP rate temperature constant CBOD half saturation constant for oxidation algal respiration rate constant algal respiration rate temperature constant optimal value of light intensity for phytoplankton growth light extinction coefficient sediment oxygen demand rate constant sediment oxygen demand temperature constant optimal temperature value Table 6 3 Parameters 31 T C water temperature Toir C air temperature Osat mg L DO concentration value at saturation Io lux day incident light intensity at the surface H m depth Vol m volume Vel m sec current speed Wind m sec wind speed Table 6 4 Variables for a site They have to be set therefore only once at the beginning of the simulation Once set these parameters are available to the EUTRO module as global parameters For every box in the discretized domain horizontal and vertical and for every time step the main program calls the subroutine EUTROOD Inside EUTROOD the differential equations that describe the bio chemical reactions are solved with a simple Euler scheme The values passed into EUTROOD can be roughly divided into 4 groups The first group is made out of the aforementioned constants that re
10. Chapter 2 Equations and resolution techniques 2 1 Equations and Boundary Conditions The equations used in the model are the well known vertically integrated shallow water equations in their formulation with water levels and transports dU ac _ gt tela RU X 0 2 1 av ac G 2 Ha RV Y 0 2 2 ae Ww XW va a Von i where C is the water level u v the velocities in x and y direction U V the vertical integrated velocities total or barotropic transports 6 E v udz v vdz h h g the gravitational acceleration H A the total water depth h the undisturbed water depth t the time and R the friction coefficient The terms X Y contain all other terms that may be added to the equations like the wind stress or the nonlinear terms and that need not be treated implicitly in the time discretization following treatment The friction coefficient has been expressed as gv u2 y2 R gt C H 2 4 with C the Chezy coefficient The Chezy term is itself not retained constant but varies with the water depth as C k H 6 2 5 where ks is the Strickler coefficient In this version of the model the Coriolis term the turbulent friction term and the nonlinear advective terms have not been implemented At open boundaries the water levels are prescribed At closed boundaries the normal ve locity component is set to zero whereas the tangential velocity is a free parameter This corresponds to a full slip condition
11. all values is 0 i e if no file with boundary values is supplied and no constant is set the value of 0 is imposed on the open boundary conz Constant boundary values for concentration temperature and salinity temp respectively If these values are set no boundary file has to be supplied salt Default 0 The next two values are used for constant momentum input This feature can be used to describe local acceleration of the water column The values give the input of momentum in x and y direction The unit is force density m s 7 In other words it is the rate of volume m s7 times the velocity m s to which the water is accelerated These values are used if boundary condition ibtyp 4 has been chosen and no boundary input file has been given If the momentum input is varying then it may be specified with the file boundn In this case the file boundn must contain three columns the first for the time and the other two for the momentum input in x y direction umom Constant values for momentum input Default 0 vmom Section wind In this section the wind data can be given directly without the creation of an external file Note however that a wind file specified in the name section takes precedence over this section E g if both a section wind and a wind file in name is given the wind data from the file is used The format of the wind data in this section is the same as the format in the ASCII wind file 1 e three columns wi
12. concentration values psu and Temperature values temp deg C for the initialization conz Files with tracer concentration values for the initialization 37 Bibliography 1 Jan O Backhaus A semi implicit scheme for the shallow water equations for aplica tion to shelf sea modelling Continental Shelf Research 2 4 243 254 1983 2 Kurt C Duwe and Regina R Hewer Ein semi implizites gezeitenmodell fiir wattge biete Deutsche Hydrographische Zeitschrift 35 6 223 238 1982 3 G Grotkop Finite element analysis of long period water waves Computer Methods in Applied Mechanics and Engineering 2 2 147 157 1973 4 B Herrling Computation of shallow water waves with hybrid finite elements Ad vances in Water Resources 1 313 320 1978 5 Bruno Herrling Ein finite element modell zur berechnung von Tidestr mungen in astuarien mit Wattflachen Die K ste 31 102 113 1977 6 K P Holz and G Nitsche Tidal wave analysis for estuaries with intertidal flats Advances in Water Resources 5 142 148 1982 7 Michael Kwizak and Andr J Robert A semi implicit scheme for grid point atmo spheric models of the primitive equations Monthly Weather Review 99 1 32 36 1971 8 A Schoenstadt A transfer function analysis of numerical schemes used to simulate geostrophic adjustment Monthly Weather Review 108 1248 1980 9 C Taylor and J Davis Tidal and long wave propagation a finite element approach Com
13. is given in figure 4 2 The first line of this section is a free one line description of the simulation that is to be car ried out The next line contains the name of the simulation in this case name_of_simulation All created files will use this name in the main part of the file name with different exten sions Therefore the hydrodynamic output file extension out will be named name_of_simulation out The last line gives the name of the basin file to be used This is the pre processed file of the basin with extension bas In our example the basin file name_of_basin bas is used The directory where this files are read from or written to depends on the settings in section name Using the default the program will read from and write to the current directory Section para This section defines the general behavior of the simulation gives various constants of pa rameters and determines what output files are written In the following the meaning of all possible parameters is given Note that the only compulsory parameters in this section are the ones that chose the duration of the simulation and the integration time step All other parameters are optional Compulsory time parameters This parameters are compulsory parameters that define the period of the simulation They must be present in all cases itanf Start of simulation Default 0 itend End of simulation idt Time step of integration Output parameters The following parameters
14. organic P F 1 Velocity coefficient ws 10 m day Sinking velocity To 20 C Optimal temperature value Table 6 8 Sediment parameters Vol m volume Volsea m sediment volume H m total depth of water column sec time step Table 6 9 Sediment variables 6 6 1 Parameters for the Water Quality Module Compulsory bio parameters These parameters are compulsory parameters that define if the water quality module is run and what kind of output is written ibio Flag if the computation on the temperature is done The model writes at each time step the state variable values in the the bio output file it smed Flag if the average minimum maximum file of variables bio salinity temperature is done if itsmed 1 the model writes sav tav output files of the corresponding variables Boundary conditions Boundary conditions have to be given in a file in every section Sbound bio2dn File name that contains boundary conditions for concentration of the water quality state variables The format is the same as for the file boundn The unit of the values given in the second and following col umn 9 data columns for EUTRO must the ones of the variable Initial conditions Initialization of variables are done by file The files can be created by the progam laplap They have to be given in section name bio File with concentration values of water quality variable to be used for the initialization salt Files with salinity
15. the text to be written and the last entry on the line is the text for the annotation A small example of an annotation would be legend 3000 2000 12 Adriatic Sea 3200 2500 12 Lido Inlet Send where the two annotations are written with a font size of 12 points The text Adriatic Sea is written at 3000 2000 so that the lower left corner of the A of Adriatic is at the specified position Section name Directory specification This parameters define directories for various input and output files basdir Directory where basin file BAS resides Default datdir Directory where output files are written Default tmpdir Directory for temporary files Default defdir Default directory for other files Default 25 Chapter 6 The Water Quality Module by Donata Melaku Canu Georg Umgiesser Cosimo Solidoro The coupling between EUTRO and FEM constitute a structure which is meant to be a generic water quality for full eutrophication dynamics The Water Quality model is de scribed fully in Umgiesser et al 2003 6 1 General Description The water quality model has been derived from the EUTRO module of WASP released by the U S Environmental Protection Agency EPA Ambrose et al 1993 and modified It simulates the evolution of nine state variables in the water column and sediment bed including dissolved oxygen DO carbonaceous biochemical oxygen demand CBOD phytoplankton carbon and chlor
16. transects As can be seen the nodes of the transects can be given on one line alone first transect two transects on one line transect 2 and 3 spread over more lines transect 4 and a last transect 21 Chapter 5 Post Processing There are several routines that do a post processing of the results of the main routine The most important are described in this chapter Note that in the model framework no program is supplied to do time series plots However there are utility routines that will extract data from the output files These data will be written in a way that it can be imported into a spreadsheet or any other plotting program that does the nice plotting 5 1 Plotting of maps with plotmap 5 1 1 The parameter input file for plotmap The format of the parameter input file is the same as the one for the main routine Please see this section for more information on the format of the parameter input file Some sections of the parameter input file are identical to the sections used in the main routine For easier reference we will repeat the possible parameters of these section here Section title This section must be always the first section in the parameter input file It contains only three lines An example has been given already in figure 4 2 The only difference respect to the para section of the main routine is the first line Here any description of the output can be used It is just a way to label the parameter file
17. where the unknowns are defined on different locations This kind of grid is usually used with the finite difference method With the form functions used here the grid of the finite element model resembles very much an Arakawa B grid that defines the water levels on the center and the velocities on the four vertices of a square Staggered finite elements have been first introduced into fluid mechanics by Schoenstadt 8 He showed that the un staggered finite element formulation of the shallow water equa tions has very poor geostrophic adjustment properties Williams 12 13 proposed a similar algorithm the one actually used in this paper introducing constant form functions for the velocities He showed the excellent propagation and geostrophic adjustment properties of this scheme The Practical Realization The integration of the partial differential equation is now performed by using the subdivi sion of the domain in elements triangles The water levels are expanded in piecewise linear functions m m 1 K and the transports are expanded in piecewise constant func tions Yn n 1 J where K and J are the total number of nodes and elements respectively As weighting functions we use y for the momentum equations and m for the continuity equation In this way there will be K equations for the unknowns one for each node and J equations for the transports one for each element In all cases the consistent mass matrix has been substi
18. H la P PE t af T e de eh Jen Seni J all dt 7 Tae ra le Se Equation 6 7 represents the limiting function given by Di Toro and used in the EUTRO program of WASP Therefore equation 6 7 is just the Steele limiting function but using la the average incident light intensity over daylight hours instead the instantaneous incident light intensity used in the original Steele formula 6 6 6 3 2 Light attenuation formula by Smith EUTRO uses also a light limitation formula by Smith In the manual this is given as I i li i PO jet 3 et 6 8 33 where now Is is the optimal light intensity which is not a constant but is continuously adjourned by the program Again I is the instantaneous light intensity at the surface I is given as TU Tt sin 6 9 2f f between 0 and f daylight and J 0 otherwise Averaging 6 9 over one whole day 0 1 gives I frl t day Ti T lr sin dt I and averaging only over daylight hours gives a dava 1 sf al Tt I To ghi _ sin dt f fso 2f of f This shows that in equation 6 9 is the average light intensity at the surface over one whole day ITOT in Eutro and that J Z f is the average light intensity at the surface over daylight hours This must be taken into account when the Di Toro formulation is used Note that in EUTRO the actual formula used is J 0 91 f where the parameter 0 9 prob ably accounts for some losses dur
19. K 2 12 where um is the coefficient of the function m and K is the order of the approximation In case of linear finite elements it will just be the number of nodes of the grid used to discretize the domain To find the values um we try to minimize the residual that arises when u is introduced into L multiplying the equation L by some weighting functions and integrating over the whole domain leading to f ynL u dQ f Wal Ont dO tn f YnL 0m dQ 2 13 Q Q Q If the integral is identified with the elements of a matrix anm we can write 2 13 also as a linear system AnmUm 0 n 1 K m 1 K 2 14 Once the basis and weighting functions have been specified the system may be set up and 2 14 may be solved for the unknowns un Wn Figure 2 1 a form functions in domain b domain of influence of node i Staggered Finite Elements For decades finite elements have been used in fluid mechanics in a standardized manner The form functions m were chosen as continuous piecewise linear functions allowing a subdivision of the whole area of interest into small triangular elements specifying the coef ficients um at the vertices called nodes of the triangles The functions m are 1 at node m and 0 at all other nodes and thus different from 0 only in the triangles containing the node m An example is given in the upper left part of Fig 1a where the form function for node i is shown The full circle indicates the node where the function
20. Model In the following an overview is given on running the model SHYFEM The model needs a parameter input file that is read on standard input Moreover it needs some external files that are specified in this parameter input file The model produces several external files with the results of the simulation Again the name of this files can be influenced by the parameter input file 4 1 The Parameter Input File The model reads one input file that determines the behavior of the simulation All possible parameters can and must be set in this file If other data files are to be read here is the place where to specify them The model reads this parameter file from standard input Thus if the model binary is called hp and the parameter file param str then the following line starts the simulation hp lt param str and runs the model 4 1 1 The General Structure of the Parameter Input File The input parameter file is the file that guides program performance It contains all nec essary information for the main routine to execute the model Nearly all parameters that can be given have a default value which is used when the parameter is not listed in the file Only some time parameters are compulsory and must be present in the file The format of the file looks very like a namelist format but is not dependent on the compiler used Values of parameters are given in the form name value or name text If name is an array the following
21. S The routines and their usage are the following nosmaniav f It generates a file containing for each node of the spatial domani aver age minimum and maximum values of the specified variable of the whole simulation nosmaniqual f It generates a water quality file from the elaboration of the state vari able It computes for each node and at each time step a water quality index that can be chosen between two suggested indexes Vismara QualityV and TRIX a well known quality index applied to the water quality definition at coastal seas and estuaries 34 Class Var 1 2 3 4 5 O2sat 90 110 70 90 or 110 120 50 70 or 120 130 30 50 lt 30 or gt 130 BODS lt 3 3 6 6 9 9 15 gt 15 NH3 lt 0 4 0 4 1 1 2 2 5 lt 5 Table 6 5 Classification of Water Quality Indices These indices can be computed using the definitions in Table 6 5 and the following equa tions QualityV class 02sat class BOD5 class NH3 auxtl phyto 30 1000 auxt2 nh3 nox 1000 auxt3 opo4 op 1000 TRIX log10 auxt1 o2satp auxt2 auxt3 1 5 1 2 nosmanintot f generates a file of total inorganic Nitrogen as sum of NH3 and NOx nosdif f computer for the chosen state variable the difference between the values at two times step nosdiff f computer the difference between the variable outputs of two simulations nosmanibod5 f computes the BODS values from the CBOD outputs as bod5 chod 1 exp 5 parl 64 14 nh3 1 exp 5
22. SHYFEM Finite Element Model for Coastal Seas User Manual Georg Umgiesser ISDGM CNR S Polo 1364 30125 Venezia Italy Version 4 93 December 2 2005 Contents Disclaimer 1 Introduction 2 Equations and resolution techniques 2 1 Equations and Boundary Conditions 22 THE Model 100 3h Pe ee AS 2 2 1 Discretization in Time The Semi Implicit Method 2 2 2 Discretization in Space The Finite Element Method 2 2 3 Mass Conservation 2 2 4 Inter tidal Flats 2 000 3 Pre Processing 3 1 The pre processing routine vp 3 2 Optimization of the bandwidth 3 3 Internal and external node numbering 4 The Model 4 1 The Parameter Input File 4 1 1 The General Structure of the Parameter Input File 4 1 2 The Single Sections of the Parameter Input File 5 Post Processing 5 1 Plotting of maps with plotmap 5 1 1 The parameter input file for plotmap 6 The Water Quality Module 6 1 General Description LL 6 2 The coupling 6 3 Light limitation 00 6 3 1 Light attenuation formula by Steele and Di Toro 6 3 2 Light attenuation formula by Smith 6 4 Initialization LL 6 5 Post processing 6 6 The Sediment Module 6 6 1 Parameters for the Water Quality Module Bibliography ii 10 10 10 11 22 22 22 26 26 27 33 33 34 34 35 37 38 Disc
23. So the computed trans ports are actually different from the transports inserted formally in the continuity equation The continuity equation is therefore not satisfied These contributions of nodes lying further apart could in principle be accounted for In this case not only the triangles Q around node i but also all the triangles that have nodes in common with the triangles Q would give contributions to node i namely all nodes and elements shown in Fig 1b The result would be an increase of the bandwidth of the matrix for the computation disadvantageous in terms of memory and time requirements Using instead the approach of the staggered finite elements actually only the water levels of elements around node i are needed for the computation of the transports in the triangles Q In this case the model satisfies the continuity equation and is perfectly mass conserving 2 2 4 Inter tidal Flats Part of a basin may consist of areas that are flooded during high tides and emerge as islands at ebb tide These inter tidal flats are quite difficult to handle numerically because the elements that represent these areas are neither islands nor water elements The boundary line defining their contours is wandering during the evolution of time and a mathematical model must reproduce this features For reasons of computer time savings a simplified algorithm has been chosen to represent the inter tidal flats When the water level in at least one of the thre
24. The other two line with the name of simulation and the basin are used to open the files needed for plotting Section para These parameters set generic values for the plot Note that the only compulsory parameter in this section is iwhat a parameter that deter mines what to plot All other parameters are optional 22 iwhat x0 y0 x1 yl x0leg yOleg xlleg ylleg dxygrd velfac trafac typls typlsf velref traref velmin tramin bndlin Flag that determines what to plot default 0 0 Nothing plotted ci Plot basin grid and isolines of depth Plot velocities Plot transports Plot water levels Plot concentration Plot temperature Plot salinity o Nl A aA A Q N Plot rms velocity Lower left corner of the plotting area Default is whole area Upper right corner of the plotting area Default is whole area Lower left corner of the area where the legend is plotted Upper right corner of the area where the legend is plotted Grid size if the results are interpolated on a regular grid A value of 0 does not use a regular grid but the original finite element grid for plotting Default 0 Not used anymore Typical length scale to be used when scaling velocity or transport ar rows If dxygrd is given this length is used and typls is not used If not given it is computed from the basin parameters Default 0 Additional factor to be used with typls to determine the length of the
25. a by Steele and Di Toro The well known light limitation function proposed by Steele is given as ae e lo where is the light intensity and the optimal light intensity P is the limiting function and takes values between 0 and 1 In this form P is a function of depth z and of time t P P z t since the light intensity depends on depth and time J I z t The depth dependence of J can be written as P i 6 3 I z lje 6 4 where J is the incident light intensity on the surface still dependent on time and k is an extinction coefficient Inserting 6 4 into 6 3 gives the equation I _ li ok I li kz P z t pene be e le Ken 6 5 o o We now compute the average of P over the water column This gives ss Yt ef H g liek P t af P z t dz af e e D dz With the substitution x e and therefore dx dz ke kx the integral can be transformed into Solving this integral gives finally e aie e dik e li KH i P t kH e lo kH e lo kH e lo e bj 6 6 0 0 This is the depth integrated form of the Steele limiting function for instantaneous light If we want to work with the average light over one day then equation 6 6 can be easily averaged over one day If T is the averaging period one day f the fraction of the day with daylight and J is approximated with a step function 0 at night and J during daytime then we can write dl 1 e la KH la e Ia k
26. be used with the time integration idtcon Time step and start time for writing to file CON concentration TEM itmcon temperature and SAL salinity istot Frequency of internal time step of the solution of the transport and dif fusion equation Normally at every external time step of the hydro dynamic equations one transport diffusion internal time step is exe cuted If the external time step is too long the solution of the transport diffusion equations with the same time step may lead to instabilities These instabilities can be avoided if more internal time steps are ex ecuted in one external step istot gives the number of internal time steps to be executed in one external step Default 1 Section name In this sections names of directories or input files can be given All directories default to the current directory whereas all file names are empty i e no input files are given Directory specification This parameters define directories for various input and output files basdir Directory where basin file BAS resides Default datdir Directory where output files are written Default 17 tmpdir Directory for temporary files Default defdir Default directory for other files Default File names The following strings enable the specification of files that account for initial conditions or forcing bound File with initial water level distribution This file must be constructed by the utility routin
27. d by the main routine The reason why the last two parameter groups are handled differently from each other has also to do with the fact that the third group is highly variable in time and space Variables like current velocity change with every time step and are normally different from element to element The fourth group is very often only slowly variable in time light wind and can very often be set constant in space Therefore these values can be set at larger intervals and do not have to be changed when looping over all the elements in the domain The overall flow of information during one time step is the following First the hydrody namic model resolves the momentum and continuity equation to update the current veloc ities and water levels After that the physical temperature and salinity and bio chemical scalars are advected and diffused Once this advection step has been handled the new load ings and forcing terms are set up and then EUTROOD is called for the bio chemical reactions Note that the operator splitting technique which decouples the advective and diffusive transport from the source term allows for different time steps of the two processes for a more efficient use of the computer resources Default values of the water quality parameters are already set in the code Owner specific parameters for the water quality model should be written in the subroutine param_user 32 6 3 Light limitation 6 3 1 Light attenuation formul
28. dard methods Once the solution for Catt is obtained it can be substituted into 2 9 and 2 10 and these two equations can be solved explicitly In this way all unknowns of the new time step have been found Note that the variable H also contains the water level through H h C In order to avoid the equations to become nonlinear is evaluated at the old time level so H h and H is a known quantity 2 2 2 Discretization in Space The Finite Element Method While the time discretization has been explained above the discretization in space has still to be carried out This is done using staggered finite elements With the semi implicit method described above it is shown below that using linear triangular elements for all unknowns will not be mass conserving Furthermore the resulting model will have propa gation properties that introduce high numeric damping in the solution of the equations For these reasons a quite new approach has been adopted here The water levels and the velocities transports are described by using form functions of different order being the standard linear form functions for the water levels but stepwise constant form functions for the transports This will result in a grid that resembles more a staggered grid in finite difference discretizations Formalism Let u be an approximate solution of a linear differential equation L We expand u with the help of basis functions m as U mlm m 1
29. deal with the output frequency and start time to external files The content of the various output files should be looked up in the appropriate section The default for the time step of output of the files is 0 which means that no output file is written If the time step of the output files is equal to the time step of the simulation then at every time step the output file is written The default start time of the output is 0 idtout Time step and start time for writing to file OUT the file containing the itmout general hydrodynamic results idtext Time step and start time for writing to file EXT the file containing hy itmext drodynamic data of extra points The extra points for which the data is written to this file are given in section extra of the parameter file idtrst Time step and start time for writing the restart itmrst 11 Stitle benchmark test for test lagoon bench venlag Send Spara itanf 0 itend 86400 idt 300 ireib 2 ilin 0 href 0 23 iczv 1 Send boundl kbound 73 74 76 ampli 0 50 period 43200 phase 10800 zref 0 end Sbound2 kbound 353 350 349 iqual 1 Send Sbound3 kbound 1374 1154 1160 1161 iqual 1 Send Sname wind winl8sep win Send MAREOGRAFI PER TARATURA 13 133 99 259 328 772 419 1141 1195 1070 1064 942 468 1154 Sextra MAREOGRAFI SEZIONI PER TARATURA 13 133 99 259 328 772 419 1141 1195 1070 1064 942 468 1154
30. e nodes of an element falls below a minimum value 5 cm the element is considered an island and is taken out of the system It will be reintroduced only when in all three nodes the water level is again higher then the minimum value Because in dry nodes no water level is computed anymore an estimate of the water level has to be given with some sort of extrapolation mechanism using the water nodes nearby This algorithm has the advantage that it is very easy to implement and very fast The dynamical features close to the inter tidal flats are of course not well reproduced but the behavior of the method for the rest of the lagoon gave satisfactory results In any case since the method stores the water levels of the last time step before the ele ment is switched off introducing the element in a later moment with the same water levels conserves the mass balance This method showed a much better performance than the one where the new elements were introduced with the water levels taken from the extrapolation of the surrounding nodes Chapter 3 Pre Processing The pre processing routine vp is used to generate an optimized version of the file that de scribes the basin where the main program is to be run In the following a short introduction in using this program is given 3 1 The pre processing routine vp The main routine hp reads the basin file generated by the pre processing routine vp and uses it as the description of the domain where the h
31. e zinit wind File with wind data The file may be either formatted or unformatted For the format of the unformatted file please see the section where the WIN file is discussed The format of formatted ASCII file is in standard time series format with the first column containing the time in seconds and the next two columns containing the wind data The meaning of the two values depend on the value of the parameter iwtype in the para section rain File with rain data This file is a standard time series with the time in seconds and the rain values in mm day qflux File with heat flux data This file must be in a special format to account for the various parameters that are needed by the heat flux module to run Please refer to the information on the file q lux restrt Name of the file if a restart is to be performed The file has to be produced by a previous run with the parameter idtrst different from 0 The data records used in the file for the restart must be given by time itrst Lagrangian trajectories LAGR Section bound These parameters determine the open boundary nodes and the type of the boundary level or flux boundary At the first the water levels are imposed on the second the fluxes are prescribed There may be multiple sections bound in one parameter input file describing all open boundary conditions necessary Every section must therefore be supplied with a boundary number The numbering of the open boundaries must be i
32. ed by such a modification The grazing has been described by means of a type II functional relationship as it is usu ally done for aquatic ecosystems However the possibility to select a type III relationship as well as to maintain the original parameterisation of constant zooplankton has been in cluded The zooplankton assimilates the ingested phytoplankton with an efficiency EFF and the fraction not assimilated ecologically representative of faecal pellets and sloppy feeding is transferred to the organic matter compartments dotted lines Fig Finally zooplankton mortality is described by a first order kinetics The code has been written by adopting the standard WASP nomenclature system and the choice between the different available functional forms is performed by setting the index IGRAZ A choice of 0 the default value corresponds to the original EUTRO version giving the user the ability to chose easily between the extended version or revert to the original one 6 2 The coupling Mathematical models usually describe the coupling between ecological and physical pro cess by suitable implementation of an advection diffusion equation for a generic tracer reads d O f O 2 0 0 O K V7 0 K dt dz sc de dz where U is the average components of the velocity the are the tracers which compose the entire vector of the biological state variable and F is a source term T and J indi cate respectively water tempe
33. en from phytoplankton and zooplankton death nitrification denitrification mineralization of ON mineralization of OP source of inorganic phosphorous from algal death sink of inorganic phosphorous for algal growth source of organic phosphorous from phytoplankton and zooplank ton death oxidation of CBOD Table 6 2 Functional Expression Description 29 C1 OC K1D PHY Zsink 31 NIT2 x 32 NIT1 32 DO1 KA Osat DO 33 DO2 PN GP1 x PHY x OC 34 DO3 1 PN GPI PHY ui x324 15 1 5 NC DO4 OC x RES PHY 36 N2 4 N1 37 NH3 NOX PN KN NH3 KN NOX iy 38 NH3 KN NH3 NOX KN NOX RES K1RC KIRT 39 SOD SODI x SODT 40 Lnu min X1 X2 mult X1 X2 41 __NH3 NOX da EN TNH3INOK 42 OPO4 Ao Fopos tOPO4 a Liignt lo x e KE H peli pad 44 KA F Wind Vel T Tair H 45 source of CBOD from phytoplank ton and zooplankton death sink of CBOD due to denitrification reareation term dissolved oxygen produced by phy toplankton using NH3 growth of phytoplankton using NOX respiration term oxygen consumption due to nitrification ammonia preference algal respiration sediment oxygen demand minimum or multiplicative nutrient limitation for phytoplankton growth nitrogen limitation for phytoplank ton growth phosphorous limitation for phyto plankton growth light limitation for phytoplankton growth re areation coeffici
34. ent Table 6 2 continued Functional Expression Description 30 K1D 0 12 day KGRZ 1 2 day KPZ 0 5 mg C L KDZ 0 168 day K1C 2 88 day KIT 1 068 KN 0 05 mg N L KP 0 01 mg P L KCnir 0 05 day KT 1 08 Knit 2 0 mg O2 L KCaenit 0 09 day KTygenit 1 045 Kaenit 0 1 mg O2 L KNCnin 0 075 day KNT nin 1 08 KDC 0 18 day KDT 1 047 NC 0 115 mg N mg C PC 0 025 mg P mg OC 32 12 mg 02 mg C EFF 0 5 FON 0 5 FOP 0 5 FOPO4 0 9 KPC nin 0 0004 day KPT min 1 08 KBOD 0 5 mg O2 L K1RC 0 096 day K1RT 1 068 Is 1200000 lux day KE 1 0m SODI 2 0 mg O2 L day m SODT 1 08 To 20 C phytoplankton death rate constant grazing rate constant half saturation constant for phytoplankton in graz ing zooplankton death rate phytoplankton growth rate constant phytoplankton growth rate temperature constant nitrogen half saturation constant for phytoplank ton growth phosphorous half saturation constant for phyto plankton growth nitrification rate constant nitrification rate temperature constant half saturation constant for nitrification denitrification rate constant denitrification rate temperature constant half saturation constant for denitrification mineralization of dissolved ON rate constant mineralization of dissolved ON rate temperature constant oxidation of CBOD rate constant oxidation of CBOD rate temperature constant N
35. for the application of very shallow basins or well mixed estu aries Storm surge phenomena can be investigated also This two dimensional version of the program is not suited for the application to baroclinic driven flows or large scale flows where the the Coriolis acceleration is important Finite elements are superior to finite differences when dealing with complex bathymetric situations and geometries Finite differences are limited to a regular outlay of their grids This will be a problem if only parts of a basin need high resolution The finite element method has an advantage in this case allowing more flexibility with its subdivision of the system in triangles varying in form and size This model is especially adapted to run in very shallow basins It is possible to simulate shallow water flats i e tidal marshes that in a tidal cycle may be covered with water during high tide and then fall dry during ebb tide This phenomenon is handled by the model in a mass conserving way Finite element methods have been introduced into hydrodynamics since 1973 and have been extensively applied to shallow water equations by numerous authors 3 9 5 4 6 The model presented here 10 11 uses the mathematical formulation of the semi implicit algorithm that decouples the solution of the water levels and velocity components from each other leading to smaller systems to solve Models of this type have been presented from 1971 on by many authors 7 2 1
36. format is used name valuel value2 valueN The list can continue on the following lines Blanks before and after the equal sign are ignored More then one parameter can be present on one line As separator blank tab and comma can be used Parameters arrays and data must be given in between certain sections A section starts with the character followed by a keyword and ends with end The keyword and end must not contain any blank characters and must be the first non blank characters in the line Other characters following the keyword on the same line separated by a valid separator are ignored Several sections of data may be present in the input parameter file Further ahead all sec tions are presented and the possible parameters that can be specified are explained The 10 sequence in which the sections appear is of no importance However the first section must always be section title the section that determines the name of simulation and the basin file to use and gives a one line description of the simulation Lines outside of the sections are ignored This gives the possibility to comment the param eter input file Figure 4 1 shows an example of a typical input parameter file and the use of the sections and definition of parameters 4 1 2 The Single Sections of the Parameter Input File Section title This section must be always the first section in the parameter input file It contains only three lines An example
37. gen ON mg N L Li Inorganic Phosphorous Q OPOA seq OPO4 yes 6 OPO4 mg P L Q OP sed OPyes OPyink 7 Organic Phosphorous OP mg P L Q ONsead ONsink ONyes id Sediment Organic Nitrogen NH3 res ONsea mg N L Q OPsea OPsink OPres sed Sediment Organic Phosphorous OP Ores OP sea mg N L Table 6 6 Sediment Mass balances NH3res Volsea KNCsed Mineralization of organic nitrogen in x KNT 2 1 x ONsed sediment ON VANONI 2 Resuspention of organic nitrogen in sediment ONgink Vol ee 3 Sink of organic nitrogen from the wa ON FPON ter column OPO4 gg Volsea KPCsed a KPTT T x OP 4 Mineralization of organic phospho rous in sediment OP reg Volsj KP s Eg OPj 5 Resuspention of organic phospho rous in sediment OP ing Vol eee 6 Sink of organic phosphorous from OP x FPOP the water column Time scale for sinking processes of ty H ws d organic N HM 8 Time scale for sinking processes of organic P Table 6 7 Sediment functional expressions 36 KNC seq 0 075 Mineralization of sediment ON rate constant KNT 1 08 Mineralization of sediment ON rate temperature constant KPC ea 0 22 Mineralization of sediment OP rate constant KPT 1 08 Mineralization of sediment OP rate temperature constant KNyes 0 1 Fraction of sediment depth resuspended day KPres 0 1 Fraction of sediment depth resuspended day FPON 0 5 Fraction of particulate organic N FPOP 0 5 Fraction of particulate
38. h This choice will lead to a nearly optimum num bering of the nodes and will be by all means good results If however you decide to do a manual optimization please follow the online instructions in the program 3 3 Internal and external node numbering As explained above the elements and nodes of the basin are re numbered in order to opti mize the bandwidth of the system matrix and so the execution speed of the program However this re numbering of the node and elements is transparent to the user The pro gram keeps pointers from the original numbering external numbers to the optimized num bering internal numbers The user has to deal only with external numbers even if the program uses internally the other number system Moreover the internal numbers are generated consecutively Therefore if there are a total of 4000 nodes in the system the internal nodes run from 1 to 4000 The external node numbers on the other side can be anything the user likes They just must be unique This allows for insertion and deletion of nodes without having to re number over and over again the basin The nodes that have to be specified in the input parameter file use again external numbers In this way changing the structure of the basin does not at all change the node and element numbers in the input parameter file Except in the case where modifications actually touch nodes and elements that are specified in the parameter file Chapter 4 The
39. he actual boundary Note that the value of iqual must be smaller than the number of the actual boundary i e boundary i must have been defined before The next parameters give a possibility to specify the file name of the various input files that are to be read by the model Values for the boundary condition can be given at any time step The model interpolates in between given time steps if needed The grade of interpolation can be given by intpol All files are in ASCII and share a common format The file must contain two columns the first giving the time of simulation in seconds that refers to the value given in the second column The value in the second column must be in the unit of the variable that is given The time values must be in increasing order There must be values for the whole simulation i e the time value of the first line must be smaller or equal than the start of the simulation and the time value of the last line must be greater or equal than the end of the simulation boundn File name that contains values for the boundary condition The value of the variable given in the second column must be in the unit deter mined by ibtyp i e in meters for a level boundary in m s for a flux boundary and in m4 s for a momentum input zfact Factor with which the values from boundn are multiplied to form the final value of the boundary condition E g this value can be used to set up a quick sensitivity run by multiplying all discharge
40. hting of the new time level of the transport terms in the continuity equation Default 0 5 ampar Weighting of the new time level of the pressure term in the momentum equations Default 0 5 Coriolis parameters The next parameters define the parameters to be used in the Coriolis terms icor If this parameter is 0 the Coriolis terms are not included in the com putation A value of 1 uses a beta plane approximation with a variable Coriolis parameter f whereas a value of 2 uses an f plane approxi mation where the Coriolis parameter f is kept constant over the whole domain Default 0 dlat Average latitude of the basin This is used to compute the Coriolis pa rameter f If not given the latitude in the basin file is used If given the value of dlat in the input parameter file effectively substitues the value given in the basin file Depth parameters The next parameters deal with various depth values of the basin href Reference depth If the depth values of the basin and the water levels are referred to mean sea level href should be 0 default value Else this value is subtracted from the given depth values For example if href 0 20thena depth value in the basin of 1 meter will be reduced to 80 centimeters 14 hzmin Minimum total water depth that will remain in a node if the element becomes dry Default 0 01 m hzoff Total water depth at which an element will be taken out of the compu tation because it becomes dry Defaul
41. ing the integration The corresponding variables in EU TRO are FDAY for f and IAV for Ig 6 4 Initialization This section describes the interpolation of data for the initialisation of the model To create a file with initial conditions the program laplap can be used The program is called as laplap lt namefile dat This makes a laplacian interpolation of specified data contained in the namefile dat This data file should have the first line empty and shold contain two colums containing respectively node number and data values for the node It generates two files laplap nos and laplap dat The first one can be used to check if the procedure has been conducted well creating a map with the plotting procedure see Postprocessing section The dat file name should be given in the section name of the str file to initialize the model You can create initialisation files for temperature salinity wind field and biological vari ables If you want to initialise the biological model with biological data you should create a single data file merging the 9 data files one for each variable using the inputmerge f routine 6 5 Post processing This section shows how to generate derivate variables The post processing routines elaborate the water quality outputs to generate derivate vari ables They allow to generate variables such as averages both in time and in space sum differences and water quality variables such us Vismara TRIX and BOD
42. is 0 the time step will kept constant at its initial value A value of 1 devides the initial time step into possibly equal parts but makes sure that at the end of the micro time steps one complete macro time step has been executed The last mode itsplt 2 does not care about the macro time step but always uses the biggest time step possible In this case it is not assured that after some micro time steps a macro time step will be recovered Please note that the initial macro time step will never be exceeded Default 0 coumax Normally the time step is computed in order to not exceed the Courant number of 1 However in some cases the non linear terms are sta ble even for a value higher than 1 or there is a need to achieve a lower Courant number Setting coumax to the desired Courant number achieves exactly this effect Default 1 idtsyn In case of itsplt 2 this parameter makes sure that after a time of idt syn the time step will be syncronized to this time Therefore setting idtsyn 3600 means that there will be a time stamp every hour even if the model has to take one very small time step in order to reach that time This parameter is useful only for itsplt 2 and its default value of 0 does not make any syncronization These parameters define the weighting of time steps in the semi implicit algorithm With these parameters the damping of gravity waves can be controlled Only modify them if you know what you are doing azpar Weig
43. l direction and the application of the boundary conditions The typical use of the new EUTRO module is as follows the main program first sets all parameters needed in EUTRO through the call to EUTRO_INI These parameters are the kinetic constants of the reactions that are described in EUTRO and are considered constant 28 GPP GP1 PHY DPP DP1 x PHY GRZ KGRZ x PHY PHY KPZ 400 GP1 Lnut Liign K1C KIT T7 DP1 RES K1D GZ EFF x GRZ DZ KDZ x ZOO Zineff 1 EFF GRZ Zsink Zineff DZ Naig NC DPP x 1 FON Natg2 PN NC GPP NOaig 1 PN NC GPP ONaig NC DPP FON Zein T Tp N1 KC KTE x NH3 Ki DO T T NITI KCaenitK i 0 x NOX x Kido ON1 KNCpin KNT x ON OP1 KPCpin KPT x OP OPaigi PC DPP x 1 FOP OPaig2 PC GPP OP 1g3 PC DPP x FOP Zyink KDT 1 20 DO KBOD DO OX KDC x CBOD x 10 11 12 13 14 15 17 18 19 20 21 22 23 24 25 26 27 28 29 30 phytoplankton growth phytoplankton death grazing rate coefficient phytoplankton growth rate with nu trient and light limitation phytoplankton and death rate respiration zooplankton growth rate zooplankton death rate grazing inefficiency on phytoplankton sink of zooplankton source of ammonia from algal death sink of ammonia for algal growth sink of nitrate for algal growth source of organic nitrog
44. laimer Copyright c 1992 1998 by Georg Umgiesser Permission to use copy modify and distribute this software and its docu mentation for any purpose and without fee is hereby granted provided that the above copyright notice appear in all copies and that both that copyright notice and this permission notice appear in supporting documentation This file is provided AS IS with no warranties of any kind The author shall have no liability with respect to the infringement of copyrights trade secrets or any patents by this file or any part thereof In no event will the author be liable for any lost revenue or profits or other special indirect and consequential damages Comments and additions should be sent to the author Georg Umgiesser ISDGM CNR S Polo 1364 30125 Venezia Italy Tel gt 39 41 5216875 Fax 1 39 41 2602340 E Mail georg lagoon isdgm ve cnr it ii Chapter 1 Introduction The finite element program SHYFEM is a program package that can be used to resolve the hy drodynamic equations in lagoons coastal seas estuaries and lakes The program uses finite elements for the resolution of the hydrodynamic equations These finite elements together with an effective semi implicit time resolution algorithm makes this program especially suitable for application to a complicated geometry and bathymetry This version of the program SHYFEM resolves the depth integrated shallow water equations It is therefore recommended
45. maximum or reference vector This is the easiest way to scale the ve locitiy arrows with an overall factor Default 1 Reference value to be used when scaling arrows If given a vector with this value will have a length of typls typlsf on the map or in case dxygrd is given dxygrd typlsf If not set the maximum value of the velocity transport will be used as velref traref Default 0 Minimum value for which an arrow will be plotted With this value you can eliminate small arrows in low dynamic areas Default 0 Name of file that gives the boundary line that is not part of the finite element domain The file must be in BND format You can use the program grd2bnd pl to create the file from a GRD file Default is no file 23 ioverl Create overlay of velocity vectors on scalar value With the value of 0 no overlay is created 1 creates an overlay with the velocity speed The value of 2 overlays vertical velocities 3 water levels and 4 overlays bathymetry Default 0 inorm Normally the horizontal velocities are plotted in scale The value of inorm can change this behavior A value of 1 normalizes velocity vec tors all vectors are the same length whereas 2 scales from a given minimum velocity velmin Finally the value of 3 uses a logarithmic scale Default 0 Section color The next parameters deal with the definition of the colors to be used in the plot A color bar is plotted too icolor Flag that determines the type of col
46. ncreasing The number of the boundary must be specified directly after the keyword bound such as bound1 or bound 1 kbound Array containing the node numbers that are part of the open boundary The node numbers must form one contiguous line with the domain ele ments to the left This corresponds to an anti clockwise sense At least two nodes must be given 18 ibtyp Type of open boundary 0 No boundary values specified 1 Level boundary At this open boundary the water level is imposed and the prescribed values are interpreted as water levels in meters 2 Flux boundary Here the discharge in m s has to be prescribed 3 Internal flux boundary As with ibtyp 2 a discharge has to be imposed but the node where discharge is imposed can be an in ternal node and need not be on the outer boundary of the do main For every node in kbound the volume rate specified will be added to the existing water volume This behavior is differ ent from the ibtyp 2 where the whole boundary received the discharge specified 4 Momentum input The node or nodes may be internal This feature can be used to describe local acceleration of the water column The unit is force density m s 2 In other words it is the rate of volume m3 s times the velocity m s to which the water is accelerated iqual If the boundary conditions for this open boundary are equal to the ones of boundary i then setting iqual i copies all the values of boundary i to t
47. nkton growth and described by the same one step kinetic law More specifically the influence of inorganic phosphorous and nitrogen availability on phytoplankton growth nutrients uptake is simulated by means of Michealis Menten Monod kinetics Reactions 42 and 43 Table 6 2 Phytoplankton uptakes nitrogen both in the forms of ammonia and nitrate but ammo nia is assimilated preferentially as indicated in the ammonia preference relation Reaction 38 Table 6 2 The influence of temperature is given by an exponential relation Reaction 13 Table 6 2 while the functional forms for the limitation due to sub optimal light con 26 dition can be chosen between three alternative options namely the formulation proposed by Di Toro et al 1971 and the one proposed by Smith 1980 Di Toro and Smith sub routines Reaction 44 Table 6 2 and the Steele formulation Steele 1962 that can use hourly light input values The choice between different available functional forms Ditoro Smith and Steele is made by setting the index LGHTSW equal to 1 2 or 3 The new version is therefore able to simulate diurnal variations depending on light intensity such as night anoxia due to phytoplankton respiration during nighttime Finally the two frequently used models for combining maximum growth and limiting fac tors the multiplicative and the minimum or Liebig s model are both implemented and the user can choose which one to adopt Reaction 41 Table 6 2
48. of two independent contributions 00 ot 90 d0 a dt 6 2 phys biol and it might be convenient in writing a computer code to devote independent modules to computation of each of them Indeed most of the modern water quality programs do have at least conceptually a modular structure In this way the same code can be used for simulating different situations by switching off the module referring to the reactor term the transport of a purely passive tracer is reproduced while a 0D close and uniformly stirred biological system is simulated if the module referring to the physical term is not included Finally the inclusion of both modules gives the evolution of tracers subjected to both physical and biogeochemical transformation in a representation that depending upon the parameterisation of the physical module can be 1 2 or 3 dimensional The whole water quality module is contained in a file weut ro f and the call to EUTRO is made through a subroutine call that is done from the main program through an appropriate interface There is a clean division between the hydrodynamic motor parameters used by the model and the resolution of the differential equations and the ecological model as evidenced by the overall structure of the modules It is the responsibility of the main module to implement the time loop administration the advective and diffusive transport of the state variables both in the horizontal and vertica
49. ophyll a PHY ammonia NH3 nitrate NOX organic nitrogen ON organic phosphorus OP orthophosphate OPO4 and zooplankton ZOO The interacting nine state variables can be considered as four interacting systems the car bon cycle the phosphorous cycle the nitrogen cycle and the dissolved oxygen balance Fig Different levels of complexity can be selected by switching the eight variables on and off in order to address the specific topics The evolution of phytoplankton concentration Reaction 1 Table 6 1 is described by the anabolic and the catabolic terms plus a grazing term related to zooplankton concentration Reaction 10 11 and 12 Table 6 2 which however is treated as a constant in the original version The anabolic term Reaction 10 Table 6 2 is related to light intensity temperature and concentration of nutrients in water while the catabolic term Reaction 11 Table 6 2 depends on temperature Phytoplankton growth is described by combining a maximum growth rate under optimal conditions and a number of dimensionless factors each ranging from 0 to 1 and each one referring to a specific environmental factor nutrient light availability which reduces the phytoplanktonic growth insofar as environmental conditions are at sub optimal levels Phy toplankton stochiometry is fixed at the user specified ratio so that no luxury uptake mech anisms are considered and the uptake of nutrients is directly linked to the phytopla
50. or table to be used 0 stands for gray scale 1 for HSB color table Default 0 isoval Array that defines the values for the isolines and colors that are to be plotted Values given must be in the unit of the variable that will be plotted i e meters for water levels etc color Array that gives the color indices for the plotting color to be used Ranges are from 0 to 1 The type of the color depends on the vari able icolor For the gray scale table 0 represents black and 1 white Values in between correspond to tones of gray For the HSB color table going from 0 to 1 gives the color of the rainbow There must be one more value in color than in isoval The first color in color refers to values less than isoval 1 the second color in color to values be tween isoval 1 and isoval 2 The last color in color refers to values greater than the last value in isoval x0col Lower left corner of the area where the color bar is plotted yOcol xlcol Upper right corner of the area where the color bar is plotted ylcol dval Difference of values between isolines If this value is greater then 0 the values for isolines and the respective colors are created automatically without need to specify the single values in arrays isoval and color Default 0 faccol Factor for the values that are written to the color bar legend This en ables you e g to give water level results in mm faccol 1000 Default 1 ndccol Decimals after the decimal poin
51. oth a wind file and section wind are given data from the file is used The wind stress is normally computed with the following formula T Pacplulu pacp u w 4 1 where pa Po is the density of air and water respectively u the modulus of wind speed and u u the components of wind speed in x y direction In this formulation cp is a dimen sionless drag coefficient that varies between 1 5 107 and 3 2 10 3 The wind speed is normally the wind speed measured at a height of 10 m iwtype The type of wind data given default 1 0 No wind data is processed 1 The components of the wind is given in m s 2 The stress t 7 is directly specified 3 The wind is given in speed m s and direction degrees A direction of 0 specifies a wind from the north 90 a wind from the east etc 4 As in 3 but the speed is given in knots dragco Drag coefficient used in the above formula The default value is 0 so it must be specified Please note also that in case of iwt ype 2 this value is of no interest since the stress is specified directly Concentrations The next parameters deal with the transport and diffusion of a conserva tive substance The substance is dissolved in the water and acts like a tracer iconz Flag if the computation on the tracer is done A value different from 0 computes the transport and diffusion of the substance Default 0 conref Reference ambient concentration of the tracer in any unit Default 0
52. par2 To run one of the postprocessing routine write the name of the routine and enter 6 6 The Sediment Module The sediment buffer action on the biogeochemical cycles could be very important espe cially in the shallow water basins and during the storm surge events The routine wsedim introduced in April 2004 aims to address the resuspention sinking dynamics of nitrogen and phosphorous This routine can be switched on and off as needed by the user setting the bsedim parameter true or false in the bio3d routine It is called after the the eutrophication subroutine It allows to follow the dynamics of two additional variables OPsed and ONsed that simu late the evolution of Nitrogen and Phosphorous detritus in sediment that are not subjected to advection diffusion processes These two variables interact with the Nitrogen and Phosphorous cycle as decribed by the equations in Table 6 6 When the wsedim subroutine is switched on OP ON NH3 and OPO4 are updated at each time step in agreement with those equations The resuspension is a linear function of the water velocity calculated by the hydrodynamic model at each box as written in Table 6 7 The amount of the sinking nutrients depends on specific prossess parameters as given in Table 6 7 and on the depth of the underlying column 35 oS _ ot O S sea General Reactor Equation O NH3 sea NH3 yes 3 Ammonia NH3 mg N L Q ON sea ONres ONsink 5 Organic Nitro
53. plete or are missing altogether the program terminates with an error 3 2 Optimization of the bandwidth The main task of routine vp is the optimization of the internal numbering of the nodes and elements Re numbering the elements is just a mere convenience When assembling the system matrix the contribution of one element after the other has to be added to the system matrix If the elements are numbered in terms of lowest node numbers then the access of the nodal pointers is more regular in computer memory and paging is more likely to be inhibited However re numbering the nodes is absolutely necessary The system matrix to be solved is of band matrix type Le non zero entries are all concentrated along the main diagonal in a more or less narrow band The larger this band is the larger the amount of cpu time spent to solve the system The time to solve a band matrix is of order n m where n is the size of the matrix and m is the bandwidth Note that m is normally much smaller than n If the nodes are left with the original numbering it is very likely that the bandwidth is very high unless the nodes in the file GRD are by chance already optimized Since the bandwidth m is entering the above formula quadratically the amount of time spent solving the matrix will be prohibitive E g halving the bandwidth will speed up computations by a factor of 4 The bandwidth is equal to the maximum difference of node numbers in one element It is the
54. present the kinetic constants and other pa rameters that do not vary in time and space The second group represents the state variables that are actually modified by EUTROOD through the bio chemical reactions These variables are transported and diffused by the main routine and are just passed into EUTROOD for the description of the processes After the call no memory remains in EUTROOD of these state variables They must therefore be stored away by the main routine to be used in the next time step again The third and fourth groups of values have to do with the forcing terms They have been divided in order to account for the different nature of the forcing terms The third group consists of the hydrodynamic forcing terms that are directly computed by the hydrodynamic model and parameters related to the box discretization They consist of water temperature salinity current velocity and depth and type of the box Here the type identifies the position of the box surface water column sediment which is needed for some of the forcings to be applied These variables are passed directly into EUTRO through a parameter list The last group contains other forcing terms that are not directly related to the hydrodynamic model These consist of the meteorological forcings wind speed air temperature ice cover light climate surface light intensity day length and sediment fluxes These parameters are set through a number of commodity functions that are calle
55. puters amp Fluids 3 125 148 1975 10 Georg Umgiesser A model for the Venice Lagoon Master s thesis University of Hamburg 1986 11 Georg Umgiesser and Andrea Bergamasco A staggered grid finite element model of the Venice Lagoon In J Periaux K Morgan E Ofiate and O C Zienkiewicz editors Finite Elements in Fluids Pineridge Press 1993 12 R T Williams On the formulation of finite element prediction models Monthly Weather Review 109 463 1981 13 R T Williams and O C Zienkiewicz Improved finite element forms for the shallow water wave equations International Journal for Numerical Methods in Fluids 1 81 1981 38
56. rature and Irradiance level while w represent the downward U VO w r err 6 1 27 Q PHY GPP DPP GRZ Q ZOO GZ DZ Q NH3 Nagi ON1 Naig2 N1 Q NOX N1 NOgig NIT 1 Q ON ONaig ONI Q OPO4 OPag1 OP1 OPaig2 Q OP OPaig3 OP1 Q CBOD C1 OX NIT2 Q DO DO1 DO2 DO3 General Reactor Equation 1 Phytoplankton PHY mg C L 2 Zooplankton ZOO mg C L 3 Ammonia NH3 mg N L 4 Nitrate NOX mg N L 5 Organic Nitrogen ON mg N L 6 Inorganic Phosphorous OPO4 mg P L 7 Organic Phosphorous OP mg P L Carbonaceous Biological Oxygen Demand CBOD mg O2 L 9 Dissolved Oxygen DO mg O2 L D04 N2 OX SOD Table 6 1 Mass balances flux rates sinking velocity for the tracer and K and K are the eddy coefficients for horizontal and vertical turbulent diffusion The term F includes the contributions of the biological biogeochemical activities and the whole biological state vector is explicitly considered in the last term of equation 6 1 without a spatial operator As far as the biologically induced variations are regarded the fate of each tracer in every location x y z is tightly coupled to other tracers in the same location but is not directly influenced by processes going on elsewhere Therefore in this approximation the global temporal variation of any tracer state variable conservative or not can be split into the sum
57. refore important to re number the nodes in order to minimize this number However there exist only heuristic algorithms for the minimization of this number The two main algorithms used in the routine vp are the Cuthill McGee algorithm and the algorithm of Rosen The first one starting from one node tries to number all neighbors in a greedy way optimizing only this step From the points numbered in this step the next neighbors are numbered This procedure is tried from more than one node possibly from all boundary nodes The numbering resulting from this algorithm is normally very good and needs only slight en hancements to be optimum Once all nodes are numbered the Rosen algorithm tries to exchange these node numbers where the highest difference can be found This normally gives only a slight improvement of the bandwidth It has been seen however that if the node numbers coming out from the Cuthill McGee algorithm are reversed before feeding them into the Rosen algorithm the results tend to be slightly better This step is also performed by the program All these steps are performed by the program without intervention by the operator if the automatic optimization of bandwidth is chosen in the program vp The choices are to not perform the bandwidth optimization at all GRD file has already optimized node number ing perform it automatically or perform it manually It is suggested to always perform automatic optimization of the bandwidt
58. s by a factor without generating a new file Default 1 conzn File name that contains values for the respective boundary condition tempn i e for concentration temperature and salinity The format is the same saltn as for file boundn The unit of the values given in the second column must the ones of the variable i e arbitrary unit for concentration centi grade for temperature and psu per mille for salinity 19 intpol Order of interpolation for the boundary values read through files Use for for stepwise no interpolation 2 for linear interpolation The default is cubic interpolation 4 The next parameters can be used to impose a sinusoidal water level tide or flux at the open boundary These values are used if no boundary file boundn has been given The values must be in the unit of the intended variable determined by ibtyp ampli Amplitude of the sinus function imposed Default 0 period Period of the sinus function Default 43200 12 hours phase Phase shift of the sinus function imposed A positive value of one quar ter of the period reproduces a cosine function Default 0 zref Reference level of the sinus function imposed If only zref is specified ampli 0 aconstant value of zref is imposed on the open boundary With the next parameters a constant value can be imposed for the variables of concentration temperature and salinity In this case no file with boundary values has to be supplied The default for
59. t 0 05 m hzon Total water depth at which a dry element will be re inserted into the computation Default 0 10 m hmin Minimum water depth most shallow for the whole basin All depth values of the basin will be adjusted so that no water depth is shallower than hmin Default is no adjustment hmax Maximum water depth deepest for the whole basin All depth values of the basin will be adjusted so that no water depth is deeper than hmax Default is no adjustment Bottom friction The friction term in the momentum equations can be written as Ru and Rv where R is the variable friction coefficient and u v are the velocities in x y direction respectively The form of R can be specified in various ways The value of ireib is choosing between the formulations In the parameter input file a value A is specified that is used in the formulas below ireib Type of friction used default 0 0 No friction used 1 R is constant 2 is the Strickler coefficient In this formulation R is written as R 44 with C ksH 6 and A ky is the Strickler coefficient In the above formula g is the gravitational acceleration u the modulus of the current velocity and H the total water depth 3 A is the Chezy coefficient In this formulation R is written as R u and A C is the Chezy coefficient 4 R A H with H the total water depth qll 5 R Xkl czdef The default value for the friction parameter A Depending on the value of ireib the coefficient
60. t for the values written to the color bar legend Use the value 1 to not write the decimal point Default 1 legcol Text for the description of the color bar This text is written above the color bar 24 Section Sarrow The next parameters deal with the reference arrow that is plotted in a legend The arrow regards the plots where the velocity or the transport is plotted x0arr Lower left corner of the area where the reference arrow is plotted yOarr xlarr Upper right corner of the area where the reference arrow is plotted ylarr facvel Factor for the value that are written to the arrow legend for the velocity factra and transport This enables you e g to give velocities in mm s facvel 1000 Default 1 ndcvel Decimals after the decimal point for the values written to the arrow leg ndctra end velocity and transport Use the value 1 to not write the decimal point Default 2 legvel Text for the description of the arrow legend velocity and transport legtra This text is written above the arrow legend arrvel Length of arrow in legend in velocity or transport units If not given arrtra the arrow length will be computed automatically Default 0 Section legend In this section annotations in the plots can be given The section consists of a series of lines that must contain the following information The first two values x y give the position of the annotation The third value gives the font size in points for
61. tegrals is zero If now the expressions for U V are introduced we obtain a system with again only the water levels as unknowns db ocr dP det n 1 2 N N fotto nta A a by dt aby AC n 2 N N ewe dQ At 2 ag f H Fata dd oP de ar 2 1 0 f U N yn qQ 2 15 Qo dx oy dDy dDy Ar 2 f X Y dQ ar 2ja Xt SP Dd Here we have introduced the symbol as a shortcut for n 1 1 ArR The variables and unknowns may now be expanded with their basis functions and the com plete system may be set up 2 2 3 Mass Conservation It should be pointed out that only through the use of this staggered grid the semi implicit time discretization may be implemented in a feasible manner If the Galerkin method is applied in a naive way to the resulting equation 2 11 introducing the linear form func tions for transports and water levels and setting up the system matrix the model is not mass conserving This may be seen in the following way see Fig 1b for reference In the computation of the water level at node i only and transport values belonging to triangles that contain node i enter the computation full circles in Fig 1b But when in a second step the barotropic transports of node j are computed water levels of nodes that lie further apart from the original node i are used hollow circles in Fig 1b These water levels have not been included in the computation of the water level at node i
62. th the first specifying the time in seconds and the other two columns giving the wind data The interpretation of the wind data depends on the value of iwtype For more information please see the description of iwtype in section para Section Sextra In this section the node numbers of so called extra points are given These are points where water level and velocities are written to create a time series that can be elaborated later The output for these extra points consumes little memory and can be therefore written with a much higher frequency typically the same as the integration time step than the complete hydrodynamical output The output is written to file EXT 20 The node numbers are specified in a free format on one ore more lines An example can be seen in figure 4 1 No keywords are expected in this section Section flux In this section transects are specified through which the discharge of water is computed by the program and written to file FLX The transects are defined by their nodes through which they run All nodes in one transect must be adjacent i e they must form a continuous line in the FEM network The nodes of the transects are specified in free format Between two transects one or more 0 s must be inserted An example is given in figure 4 3 Sflux 1001 1002 1004 0 35 37 46 0 0 56 57 58 0 407 301 435 0 89 87 Send Figure 4 3 Example of section flux The example shows the definition of 5
63. tuted with the their lumped equiva lent This was mainly done to avoid solving a linear system in the case of the momentum equations But it was of use also in the solution of the continuity equation because the amount of mass relative to one node does not depend on the surrounding nodes This was important especially for the flood and dry mechanism in order to conserve mass Finite Element Equations If equations 2 9 2 10 2 11 are multiplied with their weighting functions and integrated over an element we can write down the finite element equations But the solution of the water levels does actually not use the continuity equation in the form 2 11 but a slightly different formulation Starting from equation 2 6 multiplied by the weighting function y and integrated over one element yields n 1 n n 1 n enero ada f o tu oy 2 dQ 0 If we integrate by parts the last two integrals we obtain f Oy Crt E dQ 4 Seon U IPN Cyn dQ 0 Q a dx oy plus two line integrals not shown over the boundary of each element that specify the normal flux over the three element sides In the interior of the domain once all contribu tions of all elements have been summed these terms cancel at every node leaving only the contribution of the line integral on the boundary of the domain There however the boundary condition to impose is exactly no normal flux over material boundaries Thus the contribution of these line in
64. ydrodynamic equations have to be solved The program vp is started by typing vp on the command line From this point on the program is interactive asking you about the basin file name and other options Please follow the online instructions The routine vp reads a file of type GRD This type of file can be generated and manipulated by the program grid which is not described here In short the file GRD consists of nodes and elements that describe the geometrical layout of the basin Moreover the elements have a type and a depth The depth is needed by the main program hp to run the model The type of the element is used by hp to determine the friction parameter on the bottom since this parameter may be assigned differently depending on the various situations of the bottom roughness This file GRD is read by vp and transformed into an unformatted file BAS It is this file that is then read by the main routine hp Therefore if the name of the basin is lagoon then the file GRD is called lagoon grd and the output of the pre processing routine vp is called lagoon bas The program vp normally uses the depths assigned to the elements in the file GRD to determine the depth of the finite elements to use in the program hp In the case that these depth values are not complete and that all nodes have depths assigned in the GRD file the nodal values of the depths are used and interpolated onto the elements However if also these nodal depth values are incom

SHYFEM Finite Element Model for Coastal Seas User Manual

Contents

Download Pdf Manuals

Related Search

Related Contents