(1)
A. Atzeni1, A. Balzano1, R. A. Falconer2, B. Lin2 and Y. Wu2
1 Dipartimento di Ingegneria del Territorio, Sezione Idraulica, Università di Cagliari
Piazza d’Armi, 16, I09123 Cagliari, Italy
Phone: +39 070 675 5304, Fax: +39 070 675 5310, E-mail: balzano@unica.it
2 School of Civil Engineering, Cardiff University
PO Box 686, Newport Road, Cardiff CF24 3TB, Wales, UK
Phone: +44 (0) 29 2087 4280, Fax: +44 (0)
29 2087 4597
E-mail: {falconerra,linbl,wuy2}@cardiff.ac.uk
Abstract: Three two-dimensional shallow water flow models and four methods for simulating flooding and drying processes are compared and evaluated. The comparison has been carried out for the significant case study represented by Poole Harbour, Dorset, UK, a natural harbour whose wet plain surface is nearly halved between high and low tide. After a brief description of the hydrodynamic models, preliminary simulations are carried out, with each model being implemented with its own wetting and drying algorithm. On this base, one model is chosen to implement all four methods for flooding and drying in it for testing them. The four methods are then compared pointing out both rough and detailed differences of results. It is shown that two of the methods tested, as expected, actually yield what appears to be an effectively inaccurate evaluation of the storage volume of the water body, resulting in significant modification of the peak values of the computed results. Further simulations carried out with the other two methods yield significant improvement of the results compared to the original couplings between hydrodynamic model and flooding and drying method. In particular, it is shown that the depth tolerance value, often concerned with in several criteria for the checking of flooding and drying conditions, can be actually chosen to be one order of magnitude smaller than previously used in the DIVAST model, with use of the WWU method proposed by Balzano (1998) (10-2 m or less). This result is achieved without wiggles or other irregularities of merely numerical origin being introduced in the computed results. Also, bottom elevations are accounted for which are consistent with actual bathymetry, implying a correct evaluation of the storage volume of the cell in the context of the linear bottom elevation model assumed.
Keywords: moving boundaries, numerical models, tidal flows, wetting and drying
Accurate simulation of free surface flows on domains with moving boundaries is needed in several engineering applications, including flood propagation, storm surge, wave swash and tidal flows predictions. However, while in the first two only the advancement of a wave front needs be accurately represented for practical purposes, its alternate recession, with consequent drying of the bottom, also needs to be properly modelled in wave swash and tide- or wind-generated flows, through several wetting and drying cycles. In this paper only modelling of tidal and wind flows will be considered. However, relevant results can, in principle, be applicable to modelling of the advancing front of a storm surge or river flood, these phenomena being susceptible of modelling, as tidal flows, in the framework of the shallow water equations for hydrostatic flow conditions, at least at a first approximation level.
In the case of tidal or wind currents, further difficulties are associated with the computation of advection-diffusion of scalar tracers in domains subject to flooding and drying, even small mass unbalances or computational inaccuracies potentially being able to induce large relative errors in the tracer mass balance in the wetting and drying cells. As a consequence, concentration oscillations caused by numerical methods for advection computation could be amplified, and possibly lead to instabilities.
At the moment, after several research contributions (Reid and Bodine, 1968, Leendertse 1970, Benqué et al., 1982, Stelling, 1984, Stelling et al., 1986, Falconer, 1986, Falconer and Owens, 1987, Leendertse, 1987, Falconer and Chen, 1991, Casulli and Cheng, 1993, Molinaro et al., 1994, Defina et al., 1994, Balzano, 1998), only methods based on the Eulerian discretization approach seem to have been used in industrial numerical models for simulation of such processes in real world applications. Although more attractive from a theoretical point of view, methods based on co-ordinate transformations involving the time variable seem to have been tested for schematic cases only.
A systematic comparison of ten methods suitable for practical applications, seven of which taken from literature, was carried out by Balzano (1998). The analysis was restricted to methods ensuring strict mass conservation, in view of applications to water quality models. The methods were tested on two schematic test cases for which some important features of the solution is known, and on a third test case for which the analytical solution is known. To different extents, all seven methods taken from literature showed inconsistencies of various kinds. Models in which the water elevation at the common side of two adjacent cells was always computed as the average value of the values at the cell centres showed a major physical inconsistency in the wetting and drying process of idealised 1D basins with a bottom irregularity, suggesting the need for a more refined evaluation of the total water depth at the side between two cells.
It was also pointed out that the criteria for the declaration of the wet state of a cell resulted in different estimations of effective retention volumes, ranging from maximum underestimation with the Leendertse (1987) (LEE) method to maximum overestimation with the Casulli and Cheng (1993) (MLU) method.
Of the three proposed new methods, all of which behaved consistently in the above test cases, the most accurate (VRS, Variable Retention Surface) was based on an iterative procedure to achieve full mass conservation and was identified as presently unpractical for industrial applications in long-term simulations. The other two proposed methods (WWU, Weighted Water level Upstreaming, and WUT, Water level Upstreaming with Threshold) yielded similar results and proved efficient for use in practical applications.
Indeed, the question remained open if the above evident inaccuracies of the tested methods (and associated improvements achievable with the newly proposed methods) in a detailed analysis on the schematic test cases would have their counterparts in equally evident inaccuracies and related improvements in field applications, in which, for instance, the actual knowledge of the bottom topography itself is not free from uncertainties. Also, it is of interest ascertaining possible improvements related to the use of different kinds of the whole hydrodynamic model rather than of specific algorithms for the simulation of wetting and drying only.
A comparison has thus been carried out between numerical models and relevant methods for handling flooding and drying, on a significant case study located in the south coast of England, namely Poole Harbour, Dorset, where extensive flooding and drying problems exist.
The numerical models used in this study are based on the depth integrated shallow water wave equations, including the x and y direction momentum equations and the continuity equation respectively, one possible form of which is:
(1)
(2)
(3)
where u, v are the depth averaged velocity components in x and y directions; U = hu and V = hv are the corresponding unit discharge components; h is the total depth, t is the time, b is the momentum correction factor for the non-uniform vertical velocity profile; f is the Coriolis parameter, g is the gravitational acceleration, h is the water surface elevation above datum; r is the fluid density; tbx and tby are the bed shear stress components in x and y directions; tsx and tsy are the surface (wind) shear stress components in x and y directions, and nt is the depth averaged turbulent viscosity.
The three numerical models tested actually use different equations from those written above as to some terms and/or the computational variables used, as briefly shown in the following:
(1) DIVAST_FDM (Depth Integrated Velocities And Solute Transport) is a finite difference models of the ADI type, using the above equations but in the form containing the u, v, h variables alone. Time centring of the advection terms is achieved via iterations, ensuring second order time accuracy to the numerical scheme, which is also second order accurate in space and unconditionally stable. Full details of the governing finite difference scheme can be found in Falconer (1986);
(2) DIVAST_FVM is a finite volume version of the latter, of the ADI type as well. In order to obtain equations in full conservation form, the water elevation gradient term is arranged as follows in the x momentum equation and similarly in the y momentum equation:
(4)
where d is the undisturbed water depth. The solution techniques are similar as those used in the finite difference version of the model. The model is also second order accurate both in time and space and unconditionally stable. Full details of the numerical scheme can be found in Wu and Falconer (1998);
(3) THELMA_MOS is a finite difference model of the operator splitting, semi-implicit type based on exactly equations (1-3), but with a unit momentum correction factor, the effects of vertical non-uniformity being included in the diffusion term as is often customarily. The model is a modification of the former THELMA model (Two-dimensional Hydrodynamic Eulerian-Lagrangean Model), in which the original Lagrangean advection step has been replaced with the Eulerian MOSQUITO scheme (Balzano, 1999), which is free from numerical diffusion. The model is second order accurate in space and unconditionally stable. However, due to the operator splitting technique used it is only first order accurate in time. Further details of the operator splitting, semi implicit scheme are given in Balzano (1996).
All three models are implemented on the usual Arakawa-C staggered grid and are strictly mass conserving.
Poole Harbour is a large tidally flushing natural coastal harbour, where extensive flooding and drying occurs, with the plan-form area being reduced by approximately 50% during spring tide. This natural harbour is located on the south-east coast of England and connected to the English channel by a narrow entrance with a width of only 300 m. A plan view of Poole Harbour area is shown in Fig. 1.
Fig.1 Map of Poole Harbour
Inflows occur into Poole Harbour and Holes Bay from the Rivers Frome and Piddle, and from sewage treatment works at Poole, Lytchett Minster and Keysworth, with the latter of these sewage works discharging into the harbour via the River Piddle. The model area was represented horizontally using a mesh of 65 ´ 57 uniform grid squares, with a length of 150 m. The bathymetry data for Poole Harbour were obtained from the Admiralty Chart (No. 2611), together with additional data provided by Poole Harbour Commissioners.
All the models were always started from rest, with the initial water levels across the basin being set everywhere to the open seaward boundary level at the start of simulations. The tidal boundary data were obtained from Poole Harbour Commissioners and Hydraulics Research Ltd., with the river flows and additional calibration data being provided by Wessex Water plc.
As a first check the results given by the three models with the same parameter values were compared in order to point out relevant characteristics of accuracy. Figures 2 and 3 show computed water elevations and velocities at Poole bridge. DIVAST_FDM and DIVAST_FVM are implemented with the FAC (Falconer and Chen, 1991) method for wetting and drying, THELMA_MOS with the WWU (Balzano, 1998) method.
Fig.2 Water elevations at Poole Bridge computed with DIVAST_FDM, DIVAST_ FVM and THELMA_MOS.
Fig.3 Velocities at Poole Bridge computed with DIVAST_FDM, DIVAST_FVM and THELMA_MOS.
Whereas elevations computed by DIVAST_FDM and THELMA_MOS are nearly identical, DIVAST_FVM shows a pronounced overestimation of minimum elevations. Further discrepancies are shown by DIVAST_FVM results compared to those obtained by THELMA_MOS and DIVAST_FDM as to velocities, with the latter model yielding sharper velocity peaks compared to THELMA_MOS. Although definitive conclusions of the above results are left for the future in the absence of accurate calibration of the models, the differences between DIVAST_FDM and THELMA_MOS results could be well ascribed to the different order of the time truncation error as mentioned before.
From this point on the methods for flooding and drying simulation were tested as implemented in DIVAST_FDM.
The flooding and drying methods compared are those reported in Falconer and Chen (1991)(FAC), Balzano (1998)(WWU), Leendertse (1987)(LEE), and Cheng and Casulli (1993)(MLU). The last three methods were all tested by Balzano (1998), while the first one is a further development of method (FAO) tested in the latter paper. Methods FAC and WWU are quite similar in the analysis of the wet/dry state of the cell; however, with WWU water depth at the common side of two cells are computed based on a weighted average of the relevant cell centres’ elevations instead of a simple average as in FAC.
LEE and MLU have been chosen among several existing methods tested in Balzano (1998) because, as shown in the paper above, they represent “limit cases” in the sense that the criteria used for declaring the wet/dry state of a grid imply maximum effective over-estimation of bottom elevations by LEE and under-estimation by MLU. Correspondingly, effective under-estimation and over-estimation of the storage volume of the water body is expected.
This is clearly shown in Figure 4, where elevation time series at the downstream reach of Wareham Channel computed with the four methods and the same parameter values are plotted. Whereas FAC and WWU give quite the same results, LEE underestimates low tide elevations since upstream areas cease first to contribute to the total discharge during the ebb tide. At the same time, upstream cells are still contributing to the total discharge with MLU, since they get dry only when the water elevation at the cell centres falls equal or below the lowest bottom side elevation (in practice, this occurs only asymptotically).
Even more dramatic is the effect on predicted
velocities at Poole Bridge shown in Figure 5, where MLU yields significantly
higher peak values compared to FAC and WWU, while LEE appears never to allow the
flooding of Holes Bay. This last is a consequence of making reference to the
highest side bottom elevation to compute a control depth with the cell water
elevation. Actually, in a cell located in the narrow channel connecting Holes
bay to the central zone of the water body, an average cell centre depth under
the chart datum of about 2 m is obtained from averaging four side depths, of
which the smallest is equal –2.0 m. Since this latter is the value used to
compute the control total depth of the cell with the LEE method, the cell is
never allowed to be flooded.
Fig.4 Water elevations at downstream reach of Wareham Channel computed with DIVAST with methods FAC, WWU, LEE, and MLU for handling of flooding and drying.
All the simulation results shown so far were computed imposing in FAC, WWU and LEE a threshold depth value equal to 0.1 m, which is also of the order of magnitude of typical values used in previous field studies with FAO and FAC (Falconer 1987, 1991), namely 0.1¸0.2 m. No threshold value needs be given for MLU. In figure 6 it is shown the effect of imposing a threshold value one order of magnitude smaller, that is 0.01 m, with WWU. This results in the sharpening of the velocity peaks, which may prove effective in a detailed calibration of the model, which is left for future work, however. Also, in principle reducing the depth tolerance implies a better evaluation of the actual storage volume of the water body, which is referred to the average cell bottom elevation instead of the maximum or minimum one as in LLE and MLU. Moreover, even at this very small depth threshold value the resulting time series of velocities appears to be still smooth and free from discontinuities of merely numerical origin.
Fig.5 Velocities at Poole Bridge computed with DIVAST with methods FAC, WWU, LEE, and MLU for handling of wetting and drying. Depth threshold values equal to 0.1 m for FAC, WWU, and LEE.
Fig.6 Velocities at Poole Bridge computed with DIVAST_WWU with values of the depth tolerance equal to 0.1 m and 0.01 m.
Although a detailed calibration of models for the case study of Poole Harbour is left for future work, the comparison of the four methods for wetting and drying seem to show qualitative results which are consistent with what was expected. In particular, two of the methods for wetting and drying tested actually yield what appears to be an effectively inaccurate evaluation of the storage volume of the water body, resulting in significant modification of the peak values of the computed results. On the other hand, using the WWU algorithm (Balzano, 1998) in the DIVAST model, the depth tolerance value, often concerned with in several criteria for the checking of flooding and drying conditions, can be chosen to be one order of magnitude smaller than previously used in model applications (10-2 m or less). This result is achieved without wiggles or other irregularities of merely numerical origin being introduced in the computed results. Also, bottom elevations are accounted for which are consistent with actual bathymetry, implying a correct evaluation of the storage volume of the cell in the context of the linear bottom elevation model assumed, as opposed to LEE and MLU algorithms.
Balzano, A. (1996). A semi-implicit model for 2D tidal flow currents. Internal Report, CRS4-TECH-REP-96/14, Cagliari, Italy.
Balzano, A. (1998). Evaluation of methods for numerical simulation of wetting and drying in shallow water flow models. Coastal Engineering, 34, 83-107, Elsevier.
Balzano, A. (1999). MOSQUITO: an efficient finite difference scheme for numerical simulation of 2D advection. Int. J. for Numerical Methods in Fluids, 31, 481-496.
Benqué, J.P., Cunge, J.A., Feuillet, J., Hauguel, A., and Holly, F.M. Jr (1982). New method for tidal current computations. J. Waterways, Port, Coastal Ocean Division, ASCE 108, 396-417.
Casulli, V., and Cheng, R. (1993). Semi-implicit finite difference methods for three-dimensional shallow water flow. Int. J. Numerical Methods Fluids, 15, 629-648. Wiley, New York.
Defina, A., and D’Alpaos, L. (1994). A new set of equations for very shallow water and partially dry areas suitable to 2D numerical models. Proc. Specialty Conference on modelling of flood propagation over initially dry areas, ASCE, 72-81.
Falconer, R.A. (1986). A water quality simulation study of a natural harbor. J. Waterways, Port, Coastal Ocean Eng. Division, ASCE 112 (1), 15-34.
Falconer, R.A., and Chen, Y. (1991). An improved representation of flooding and drying and wind stress effects in a two-dimensional tidal numerical model. Proc. Instn. Civ. Engrs., Part 2, 91, Dec., 659-678.
Falconer, R.A., and Owens, P. H. (1987). Numerical simulation of flooding and drying in a depth-averaged tidal flow model. Proc. Instn. Civ. Engrs., Part 2, Mar, 161-180.
Leendertse, J.J. (1970). A water quality simulation model for well-mixed estuary and coastal seas, Vol. 1, principles of computation. RAND publication, S. Monica, February, RM-6230-RC.
Leendertse, J.J. (1987). Aspects of SIMSYS, a system for two-dimensional flow computations. RAND publication R-3572-USGS, February.
Molinaro,
P., Di Filippo, A., and Ferrari, F. (1994). Modelling of flood wave propagation
over flat dry areas of complex topography in presence of different
infrastructures. In: Proc. Specialty Conference on Modelling of Flood
Propagation Over Initially Dry Areas, ASCE, 209-228.
Reid, R.O., and Bodine, B.R. (1968). Numerical model for storm surges in Galveston Bay. J. Waterways, Harbors Division, ASCE, 94 (WW1), 33-57.
Stelling, G.S. (1984). On the construction of computational models for shallow water equations. Rijkswaterstaat communication No. 35/1984.
Stelling, G.S:, Wiersma, A.K., Willemse, J.B.T.M. (1986). Practical aspects of accurate tidal computations. J. Hydraulic Eng., ASCE, 112 (9), 802-817.