MATHEMATICAL TECHNIQUES & COMPUTER IMPLEMENTATION OF A NEW APPROACH FOR MODELLING OF ESTUARIES

1

NAGHMEH KESHAVARZI-ROONIZI

RAHIM GHORASHI

VAHID NASSEHI

Department of Chemical Engineering

LoughboroughUniversity

Loughborough, LE11 3TU, UK

1

Abstract: -Water quality monitoring and management in tidal waterways need accurate predictive hydrodynamic and pollutant transport models. However, the development and application of realistic three dimensional numerical models to natural water systems require prohibitive amounts of hydrographic data and their use in practice is restricted by very high computational costs. To avoid these problems simplified one or two dimensional models, which are derived using depth or cross sectional averaging, are commonly used to simulate flow and mass transport in estuaries. The approximate nature of these models means that the quality of the results obtained by them can be varied and in cases where significant three dimensional effects are present they cannot be reliably used. In this paper we outline the development and computer implementation of a new approach for modelling solute transport in estuaries based on the representation of tidal hydrodynamics and solute mixing in terms of multi-mode velocity and concentration fields. The vertical components of the filed variables in the flow and transport equations corresponding to multi-mode system are solved separately from the lateral components and thus such a model can yield information in all three space dimensions while economising computational costs and hydrographic data requirements.

Key-Words: -Computer Implementation,Modelling, Estuaries, Multi-mode.

1 INTRODUCTION

In natural water systems water quality is affected by the interactions of a wide range of hydrodynamic and physico-chemical phenomena. Therefore computer models developed to simulate hydro-environmental conditions in these systems use sophisticated techniques and architectures. In particular, in systems such as estuaries and coastal zones, where crude, large scale averaging is not acceptable, the modelling is often based on sophisticated mathematical procedures which can preserve accuracy of the simulations. However, the necessity of maintaining computingeconomy as well as accuracy requires the use ofmodern IT in conjunction with reliable mathematical techniques.Over the past decade there had been a number of publication which deal with both aspects of mathematical modelling and computer implementation of complex phenomena encountered in non-liner problems of the type which may be seen in estuaries. For example see Soleimani fard and Vahidian kamyad (2003) and Wei Sun et al, (2003). The present work is an attempt to create a system for complex non-liner flow problems which improves mathematical modelling and computer applications.

In this paper we outline a novel mathematical method which can generate reliable simulations for the transport of solute pollutants in tidal waterways. This model avoids using full three dimensional computations and hence is very cost effective. We explain that the most optimum way of implementing the present model is via the use of an IT system which provides a flexible and user friendlyenvironment.We further describe the basic architecture of the ITsystem developed and provide an example of the application of the model.

2 MATHMATHICAL MODEL

The derivation of multi-modal models for tidal water ways is based on the spectral (i.e. eigen function) expansion of general three dimensional hydrodynamic and solute transport equations. In general , the velocity or concentration modes satisfy Sturm-Liouville type equations. For example, the velocity modes relate to the following equation

µ(m)Φ(m)=0, (1)

withΦ(m) = 0 on = 0 (estuary bed) and = 0 on = 1 (water surface), where Φ(m)are the velocity modes and are the associated eigenvalues, as indicated is a normalised vertical co-ordinate direction ,is a dimensionlessshape function characterising eddy viscosity profile of flow in the tidal water way. The expansion of the velocity components, for example, in the longitudinal (i.e. x) direction are hence posed as

u = (x ,y, t) Φ(m)(σ) (2)

etc. Similarly, replacing the eddy viscosity term with an appropriate eddy diffusivity, the concentration modes of solutes can be represented as

C = (x ,y, t)(m)(σ) (3)

Substitution of the above expansions into the general hydrodynamic and transport equations yields the modal equations used in the present work. In this model each mode represents a cross section of the water way at a certain depth and therefore solution of a number of two dimensional sets of governing equations provides the required vertical profile of the velocity and concentration fields. Detailed derivation of the multi-modal equations for tidal water ways is available in the published literature,(Smith, 1995). The governing equation of salt distribution which form the basis of the spectral expansions used in this work is as follows

= 0 on σ = 0,1. (4)

Here , H is water depth, u, υ and W are the longitudinal, lateral and vertical components of velocity, q and Q fresh and salt water discharges and K are the vertical and horizontal eddy diffusivity, respectively. In equation denotes a normalised vertical direction in the three dimensional space. After the derivation of the concentration modes ψ(m)(σ) and associated eigenvalues λ(m) via the Sturm-Liouville equation the modal expansion of the concentration variable is obtained asshown in equation 3.

Insertion of the modes into equation gives rise to the following modal equations which are solved to obtain the required results.

=H

-Hdσ

(5)

However, in the present work instead of exclusively relying on situations where an analytical solution for equation (1) is available, the general case where this equation is solved numerically is considered. This solution is based on the utilisation of a mesh independent diffuse approximation technique (Mokhtarzadeh, 1998). The expanded cross sectional equations, on the other hand, are solved using a Lagrange-Galerkin scheme whose details have been reported previously (Nassehi and Kafai, 1999).

3 IMPLEMENTATION

To achieve the goals of flexibility, practicality, speed, multi-user capability and computing economy the use of the described model is based on its implementation via an Information Processing Tool. The details of this system was presented in a recent international conference on Industrial Simulation (N.Keshavarzi-Roonizi et al, 2004). The following is a brief description of the system within the context of the present application.

The IT tool used here is constructed as a combined system consisting of three distinct modules.

  • Front-end, provides a network node for communication between the users and the software. Quantitative results generated by the model are treated in a back-end module and the final result is provided in clearly tabulated or graphical formats. A log is kept by the system which records the particulars of the application for future comparisons and evaluation of the efficiency of the solution generated for a given problem. The user is guided by the front end towards inserting appropriate input data. In view of the complexity of the present application the importance of this stage cannot be over emphasized. In this application the user should be guided, firstly, to insert hydrographic data consisting of water way geometry and bed profile, and time dependent fresh water and tidal flows into the flow channel as well as bed friction, wind stress and other hydraulic parameters. These data change continuously in a moving boundary domain according to flood and ebb tides and hence the system should track and up-date the flow domain. To avoid extra complications, as mentioned earlier, we have adopted a mesh independent technique to solve the mode generator equations for velocity and concentration and the solution of the expanded equations is based on a Lagrange-Galerkin approach. The front end facilitates data handling whilst monitoring the suitability and relevance of the data inserted by the users in using the data generator (simulation) software.
  • Data-generator, is a systematic library of various modelling software that use FORTRAN and symbolic programming as well as genetic algorithms to solve the model equations. Specific applications such as optimization of the bed friction values (i.e. model calibration) is carried out via genetic algorithms (Passone et al, 2004). The front end interface feeds the input data to the modelling sections initiating numerical computations. At the end of the computations the simulation results are passed to the back-end module for further processing required before presentation of the answers.
  • Back-end, provides a network node for the import and processing of the data (i.e. results, answers, evaluations) generated by the software in the data generator module. The processed results are logged and returned to a window in the front-end for the utilization by the users. It should be noted that the simulations results obtained via the data-generator mainly consist of large amounts of numerical data and hence the back–end module re-organizes the data as preformatted tables or graphical out put. The back-end module can also communicate directly with the front-end at the preliminary stages of an application (i.e. problem simulation). For example, it is designed to provide answers regarding the validity of a given input on the basis of its self-consistency or logical status. This feature is mainly utilized to guide the untrained user to obtain required information from the system with ease and efficiency. In figure 1 the general architecture and basic procedures of the present IT system is shown.

4 APPLICATION AND SAMPLE RESULTS

The applicability of the developed model to realistic situations is demonstrate by the modelling of salt intrusion in the Tay estuary.

This is as a typical solute transport case study in a large three dimensional tidal water way where distinct vertical variations in the solute level is expected to exist casting doubt in the suitability of using a simple depth averaged model.

The Tay estuary in Scotland, UK consists of the rivers Tay and Earn, Figure 2. It has the highest fresh water discharge of any British estuary, which consequently affects the mixing between fresh and salt water over a significant length. The estuary is 54 km and the front of salt intrusion extended for 45 km while the limit of tidal motion in near Perth, 50 km from the estuary mouth.

The estuary is mesotidal-macrotidal with variation of the tidal range between 2.8 and 4.7 m. Hydro-dynamically influenced by its geomorphology, the Tay estuary is classified as “complex”(Buck and Davidson, 1997).

The flow channel geometry in Tay has resulted from the combined effects of glaciations, river erosion and sea-level changes. Hence, the Tay Estuary is quite meandering and characterised by 165 km of shoreline. The dominant flow channel on its southern side is wide and shallow and includes a number of sand over-banks. The width and depth of estuary are not regular. From the entrance of the estuary at Buddon Ness, which is 1.9 km wide at low water, the width of the estuary increases exponentially inland to a maximum of 3 km at Dundee. At high water, the width of the estuary reaches a maximum of 5 km at Birkhill (22.5 km from the estuary mouth) with an exponential decrease from this point towards the estuary head.

Fig. 1 Schematic diagram of the IT system

Fig. 2 The Tay Estuary

Because of its typical flow behaviour and morphology, the Tay estuary has been the subject of many hydrological and statistical analyses in the last four decades (e.g. McManus, 1968) .

Considerable variation in the quantity and the velocity of water have been observed in the flow channel along the estuary. During the propagation of the tidal wave the phase lag can reach a maximum of 3 and 1.5 hours at low and high water, respectively. The propagation of the tidal wave is accompanied by the deformation in its shape, which increases with the phase lag. The irregular geometry of the estuary also affects the mixing and salt transport along the channel. While the differences in surface and bottom salinity in upper parts of the estuary are relatively small (3g/kg), salinity distributions in the vertical as well as longitudinal directions, measured at different locations in lower reaches, show significant asymmetries ( Buller et al, 1972).

As the brief description of the physical environment of the Tay estuary indicates the geometrical and hydro-dynamical complexity of this tidal water way is such that any accurate simulation of physical phenomenon in its environment require extensive data handling as well as elaborate and sophisticated mathematical modelling. Therefore the simulation of salt intrusion , in particular in lower reaches of the estuary where significant longitudinal and vertical variations in salinity levels are expected to exist, provides a very appropriate case study to demonstrate the applicability of the present model and IT system. An additional factor influencing this selection is the possibility of using accurately measured salinity diffusivity data for spring and neap tides based on extensive 1972 field survey of Tay estuary (Nassehi and Williams, 1987).

As mentioned earlier after the verification of the suitability of the input data inserted by the user the front end passes the data to the simulator module and computations begin. The simulation stage starts with the solution of mode generator equations (i.e. equation 1 for the velocity and a similar equation in terms of diffusivity for salt concentration). Generated modes are used to obtain modal expansions of the field unknowns via equations 2 and 3. In practice, at this stage the values of coefficients in the multi-modal forms of the hydro-dynamic and salt transport models which have been written in a FORTRAN code are determined. Therefore after this stage the multi-mode equations (i.e. equation 5) can be solved by the finite element method to generate the field variables.

The hydro-dynamic equations for velocity components are initially solved assuming an average salinity throughout the estuary and inserting a pre-determined set of bed friction values (i.e. Manning’s coefficients). The initial determination of the bed friction coefficient is based on the types of sediments in the flow channel bottom and bed morphology. After the initial solution of the flow equations the salt transport equation is solved. The cycle is then iterated by updating the salinity and friction coefficients so that flow equations yield results which are close to a selected set of observed data for a given tidal event. This is step is commonly known as the calibration step in the simulation of tidal phenomena (Nassehi and Bikangaga, 1993). The adjustment of bed friction is carried out automatically via a genetic algorithm (Passone, 2002) which significantly speeds up the model calibration. After the calibration the model is verified by running it for a different tidal event whilst all of the physical parameters are kept similar to the values determined by the calibration and not altered.

The simulation vertical salinity profiles for the calibration event ( Spring tide of 12 June, 1972) at stations B and C (Figure 2, corresponding to Broughty Ferry and Newport along the flow channel) are shown in figure 3. As this figure shows significant variations in the salt concentration throughout the tidal cycle in the water column at these locations are predicted. These profiles are compared with depth average salinities to provide a quantitative measure for the deviation from reality that will result if the salinity changes in the water column is ignored.

The results obtained for the verification event (i.e. neap tide of 20 June, 1972) at stations B and C are shown in figure 4. Despite the significant reduction in water level changes during the neap tide noticeable variations in the vertical salinity profile throughout the tidal cycle is again predicted.

The predicted asymmetry in the salinity profiles between stations B and C and also with respect to time are as expected for a complex estuary such as the Tay.

Fig. 3 Salinity Profiles, Spring tide of 12 June 1972

Fig. 4 Salinity Profiles, Neap tide of 20 June 1972

Our numerical experiments showed that to obtain results comparable with those shown in figure 3 and 4 via the use of full 3-D computations, computer CPU time of two order of magnitudes larger than the time used in the present computations is required. It should also be noted that in most cases 3-D computations require parallel processing systems.Therefore another advantage of our method is that useful quantitative results about 3-D dispersion of solutes in tidal water systems can be obtained withoutemploying costly super computers.

5 CONCLUSIONS

The novel modelling approach presented in this paper is shown to yield accurate solute concentration profiles in a water column (i.e. bed to surface direction) in tidal water ways. In this approach full three dimensional computations are avoided and hence the model provides a very cost effective scheme for obtaining results for complex tidal water ways under realistic conditions. Practical problems involving high cost of 3-D computations can be solved by this method. To overcome the complexity of the implementation of the developed multi-stage software and to make it accessible to users who may have a minimal knowledge of mathematical modelling the scheme should be used in conjunction with an information processing tool. The basic structure of the proposed IT system is explained and its application to a realistic case study is presented. The results obtained by the model for the salinity profiles in lower reaches of the Tay estuary reveal detailed information regarding solute concentration variations during a tidal cycle in the water column at selected locations along the estuary. Critical evaluation of the generated results show their theoretical consistency and close comparison with realistically expected data. The developed model is a general tool and can be applied to simulate a variety of hydro-environmental phenomena in complex estuaries and coastal zones.

6 ACKNOWLEDGMENT

The authors wish to thank Engineering and Physical Research Council of UK (EPSRC) and LoughboroughUniversity for PhD studentships, which enabled the development of the model and information processing system described in this work.