- Email: [email protected]
The simulation of a continuously deforming lateral boundary in problems involving the shallow water equations
Corsairs and ~quids Vol. 10, No. 2, pp. 105-116, 1962 Printed in Great Britain.
0045.-79301821020105-12503.0010 Pergamon Press Ltd.
THE SIMULATION OF A CONTINUOUSLY DEFORMING LATERAL BOUNDARY IN PROBLEMS INVOLVING THE SHALLOW WATER EQUATIONS B. JOHNS Department of Meteorology, University of Reading, England and S. K. DUBE, P. C. SINHA, U. C. MOHANTYand A. D. RAO Centre for Atmospheric Sciences, Indian Institute of Technology, Delhi, India
(Received 5 August 1981) kbstract--A finite-differencemethod is given for the treatment of a continuously deforming lateral fluid boundary in computer simulations of systems involving the shallow water equations. The technique is described in the context of a numerical storm surge model for the east coast of India. Using a forcing wind-stressdistributionrepresentative of the 1977Andhra cyclone, a comparisonis made between simulations using models both with and without inland intrusion of water along the coast of Andhra Pradesh.
I. INTRODUCTION The treatment of deforming fluid boundaries is an essential part of computer simulations in many dynamical problems. A classical area of note is in the context of surface wave propagation where the position of the free surface is unknown a prior/. One of the difficulties here is that the free surface boundary conditions have to be applied at a level which may not coincide with any of the computational levels in a scheme of numerical solution. In the case of shallow water wave propagation, these problems may be overcome by firstmaking a coordinate transformation to ensure that the free surface always corresponds to a fixed computational level in the numerical scheme as described in[I,2]. A similar difficultyarises in problems involving the deformation of lateralfluid boundaries. An important example of this occurs in numerical storm surge prediction models. Frequently, the lateralboundaries in these are taken to be vertical side-walls through which there is no flux of water and many such studies have been undertaken. In actuality, however, the water will usually move continuously inland and the use of idealized vertical side-walls may lead to a misrepresentation of the surge development. In particular, a simplistic horizontal inland extrapolation of the predicted surge height at the fixed side-wall in order to estimate the distance of inland penetration does not in itselftake account of the dynamics of the intrusion process. The recognition of this has led several authors to develop numerical models that allow the flow of water across the initialposition of the model coastline. This has been accomplished in [3,4] by arranging for the advancing water front to move discontinuously from one grid-point to another according to a prescribed criterion. However, to obtain acceptable accuracy, this procedure requires an extremely fine grid-spacing which may not be otherwise necessary. More recently, a finiteelement simulation of shallow water flow in continuously deforming regions has been proposed by Lynch and Gray[5]. The attractions of their method are considerable because it accommodates a continuous motion of the shoreline. However, the computational overheads characteristic of finite element procedures may make the method expensive. In the present work, we describe a finite-difference method which also models a continuously moving lateral boundary. Although admitting of application in various contexts, we apply it to the numerical simulation of the surge generated by the 1977 Andhra cyclone which struck the east coast of India and led to disastrous inland flooding. An attempt to simulate this event numerically has previously been made by Johns et al. [6]. Their treatment incorporated fixed vertical side-wails through which there is no flux of water. However, they did employ a curvilinear representation of the lateral boundaries which allowed these to be modelled in a fairly realisticway. Using this framework, finite-differencecalculations were undertaken with wind-stress forcing data appropriate to the 1977 Andhra cyclone. 105
106
B. JOHNS
The method we employ here is a generalization of the curvilinear boundary technique described in[6] and is conceptually similar to the free surface treatment given in[I]. Using an upgraded representation of the data for the wind-stress forcing field, we have illustrated the method by comparing simulations of the surge with and without a deforming fluid boundary along the coast of Andhra Pradesh. Our conclusion is that the incorporation of a deforming coastline in storm surge models may lead to an important change in the prediction detail compared with that from a fixed coastline model. 2. FORMULATION All conditions are referred to a set of rectangular Cartesian coordinates in which the origin, O is within the equilibrium level of the free surface. Ox points towards the east, Oy towards the north and Oz is directed vertically upwards. The displaced position of the free surface is given by z = ~(x, y, t) and the position of the sea-floor by z = - h ( x , y). A western coastal boundary (the east coast of India) is situated at x = bl(y, t) and an eastern open-sea boundary is situated at x = b2(y). Southern and northern open-sea boundaries are at y = 0 and y = I, respectively. This configuration is shown in Fig. 1. The depth-averaged components of velocity, u and v, then satisfy
__aU+ u au + au 3~ + l_~ {,r~ _ kpu (u2 + v2),12} at -~x V ~ - y - [V = - g ax Hp
(2.1)
u2}.
(2.2)
and
_av+uav+ 3t
av
-~x
_
V ~y + [U =
a_~+ 1 {~.~_kpv(u2+v2)
gay
np
In (2.1) and (2.2), [ denotes the Coriolis parameter, the pressure is taken as hydrostatic and H is the total depth ~ + h. (~, t~) denotes the applied surface wind-stress and the bottom stress is parameterized in terms of a quadratic law. p denotes the water density and the friction coefficient, k, is taken as 2.6 x 10-3. The kinematical condition at the western coastal boundary requires either that the coastline consist of a vertical side-wall through which there is no flux of water or, if there is inland intrusion of water, that the shoreline move with the same velocity as that of its constituent water particles. In the first case, we require Obi u-v-~:--O at x = b t ( y , 0 )
(2.3)
_
oy
25*
y=t
BANGLADESH
-~ .-:...:f""~
INDIA/~//~/'-~" ' "v~!,. x : b,(y, t ) ~ /
". B U R M A
//
/i
f" iI
- - 20"
i"
!"~
c'"
-- 15"
X= bz(y)
I
,,
i'/,o %... '"
I
80°
....
°
d
"~ - - - ~ x
I
85" F i g . 1. T h e
I
[
90" analysis
95" area.
~uu" 5"
Problemsinvolvingthe shallowwaterequations
107
and, in the second Oy at x=b~(y,t).
(2.4)
Thus, in (2.3) the value of b~ is invariant in time and is always equal to its initial value. In (2.4), the position of the coastline changes and the kinematical condition must be applied at its a priori unknown location. At the eastern open-sea boundary, a radiation condition is applied that is identical to that used by Johns et al. [6] in their coastal zone model. This leads to
062_ ( ~ ) 1/2 u-v~y ~=0
at x=b2(y).
(2.5)
Conceptually similar radiation conditions are applied at the southern and northern open-sea boundaries to yield v+(~)
112 ,=0
at y=O
(2.6)
at y = l .
(2.7)
1/2 v-(~)
~=0
To accompany (2.1) and (2.2), the equation of continuity is applied in the vertically integrated form
o/-/+o__ Ot
Ox (Hu)+
(Hv) = O.
(2.8)
When bt is temporally invariant, the above formulation is the same as that of the coastal zone model described in [6] except that the northern boundary is now an open-sea boundary.
3. COORDINATE TRANSFORMATION A coordinate transformation is introduced which is a generalization of that used in[6] to facilitate the numerical treatment of an irregular boundary configuration. We now write _ x - b~(y, t)
-
b(y,t)
(3.1)
where b(y, t) = b2(y)- b~(y, t).
(3.2)
Thus, the western and eastern boundaries correspond, respectively, to ~ = 0 and ~ = 1. Equations (2.1) and (2.2) may firstly be written in flux form and then, with ~, y and t as new independent variables, may be transformed to Oti 0 . . ti)+ 0 ( v ~ ) _ f ~ = _ O~ b'r~_kfi . 2 ± . 2 , , / 2 ~-+-~(U 0), gH o!~+ P H ( u -,v j
(3.3)
and
p - -~ (u 2 + v2) 112,
(3.4)
106
B. JOHNS
where
(3.5) Hbu} = Hbv " =
(3.6)
Equation (2.8) is readily found to yield
o~(Hb)+
(Hb U)+
05)=0.
(3.7)
From (2.3), (2.4) and (3.5) we see that the kinematical condition at the curvilinear boundary is always satisfied provided that U=0
at
~=0.
(3.8)
The radiation condition (2.5) at the eastern boundary is satisfied if 1/2
As in [6], the numerical procedure will be to solve (3.3)-(3.7) in a rectangular computational area 0 ~<~ ~< 1, 0 ~< y ~< I. On satisfying (3.8) and (3.9), the conditions (2.3)-(2.5) on the western and eastern coastlines are automatically fulfilled. The conditions (2.6) and (2.7) are satisfied directly when implementing the numerical scheme. In the case of inland intrusion of water, it remains to establish a procedure for determining the changing position of the coastline. The kinematical condition (2.4) has already been satisfied by taking U = 0 at ~ = 0 and a further requirement is that the depth of the water be zero at the coastline. This leads to H=0
at
x=b~(y,t),
(3.10)
or, equivalently,
(,(~, = O, y, t) + h[x = b~(y, t), y] = O.
(3.11)
Depending on whether b,(y, t),~ ba(y,0), we may then either interpolate or use new inland orographical data to fix the value of h[x = b,(y, t), y]. For present purposes, however, the simplest procedure is to differentiate (3.11) with respect to t. This leads to O--~-(~= O, y, t) +
s=O
(3.12)
where S=
• =b,(y.,)"
(3.13)
If s is prescribed, (3.12) yields a prognostic equation for b,. The simplest case corresponds to a constant value of s when (3.12) immediately integrates to ba(y, t) = bl(y, 0) _ 1 ~(~ = 0, y, t).
(3.14)
Then, if ~(~ = 0, y, t ) < 0, the sea-surface at the shoreline is depressed and the shoreline has
Problemsinvolvingthe shallowwater equations
109
consequently receded from its initial position. If ~(~ = 0, y, t ) > 0, there is a positive surge; the elevation at the shoreline is raised above its equilibrium level and there is a corresponding inland penetration of water. The numerical procedure to be described in Section 4 leads to an updating scheme for ~(~ = 0, y, t). Our formulation may, therefore, be used to monitor the continuous movement of the shoreline. This procedure is based rigorously on the fulfilment of the physically correct coastline conditions (2.4) and (3.10) and allows full account to be taken of the dynamics of the intrusion (or recession) process. 4. NUMERICAL PROCEDURE We define discrete coordinate points in the ~-y plane by =
¢i
=
(i -
1 ) A s ¢,
i
=
1, 2 ....
m
;A¢ = ll(m ll(n
y = y~ = (j - l)Ay, j = I, 2 .... n ;Ay =
-
-
1)'~ I) J"
(4.1)
A sequence of time instants is defined by t = tp = p A t ,
p = O, 1 . . . . .
(4.2)
X~.
(4.3)
For any variable, X, we write
X(~i,
Y~,
tp)
=
In the ~-y plane we use a staggered grid where there are three distinct types of computational point. With i even and j odd, the point is a ~-point at which ~ is computed. If i is odd and j is odd the point is a u-point at which both u and U are computed. If i is even and j is even, the point is a v-point at which v is computed. Along s¢ = 0 we, therefore, have only u-points and at each of these we have U = 0 identically. We choose m to be even so that ~ = 1 consists of both ~-points and v-points. We also choose n to be odd thus ensuring that there are only ~-points and u-points along y = 0 and y = i. This arrangement of grid-points is shown in Fig. 2. With this scheme, it will be noted that computed elevations are not carried along the coastline ~ -- 0. In fact, an updating procedure is never invoked along ~ -- 0 where (in the case of inland intrusion) the depth is zero. This feature in our scheme avoids problems relating to the singular behavior of (2.1) and (2.2) when H = 0 to which reference has been made in[5]. In order to overcome this difficulty, it was suggested there that a physical re-formulation of the bottom friction may be required near the shoreline--a view to which we do not necessarily subscribe. In order to describe the finite-difference equations, we define difference operators by
A'X = (X~+' - X~)IAt
1 (4.4)
80~ = (X,P+~,j - X~-~,j)/(2AO t. 8yX = (X~+~ - X£~-d/(2,~y)J
.T:5
•
[
3=4 J=3
()
s-,
3=2
.T=I
Ay • I:1 Fig.
~p
I:2 1:3 1:4 1:5 I=6 2. Grid-pointarrangement.
110
B. JOHNS
Averaging operations are defined by
~ = ~(xr+,.j + xr-,.~) 1 ~?' =
'(xb+,
+
x~J-,) t
(4.5)
and a shift operator is defined by
EtX = vP'+l ,~,j .
(4.6)
The discretization of (3.7) is then based upon A,(Hb) + a~ (/4~b U) + ar 03) = 0.
(4.7)
Equation (4.7) is used in the updating procedure at points corresponding to i--2(2)m- 2 and j = 3(2)n - 2. The elevation at j = 1 is determined from (2.6) which, in practice, is applied at j = 2 and replaced by l(g~
'/2
~/
(W"+W')+vf2=o,
(4.8)
thus leading to an updating procedure for the elevation along the southern open-sea boundary. Likewise, the elevation at j = n is determined from (2.7) which is then replaced by 1( g ],,2 p+l X \z----.,,,.._,/ ( W ' + v~,,.-2,v~._, = O.
_
2
(4.9)
The elevation along ~ = 1 determined from (3.9) which is replaced by 1
"
g
xl12
-2 [-~\h,.-1,J
(~mjP+l+ ~m-2,j)p+* b Pj IJ,,-I,j p = O,
(4.10)
thus yielding updated elevations for i = m. In contrast with[6], the time-dependent nature of b should be noted in (4.10). Additionally, and most importantly, the application of (4.7)-(4.10) implies that the position of the grid points is changing with time in the physical x-y plane although they are fixed in the computational 6-y plane. This means that the grid-point values of h in the computational plane will be changing with time whereas they are independent of time at fixed points in the physical plane. Accordingly, after each updating of the prognostic variables it is necessary to re-compute new values of hit in the computational plane. This can only be done when the change has been determined in b~ corresponding to the current advance of the solution through a time At. This change may be computed from (3.12) by applying
(b,),~÷' = (b,),~ - ~ (¢V I - ¢~j).
(4.11)
The value of ~ 1 in (4.11) is not carried in our computational scheme. It is, therefore, determined by linear extrapolation from the adjacent ~-points by using ~¢j+1= ~(3~j+1_ ~j+l).
(4.12)
Values of ~,~+~ at the ~-points are determined from the prognostic variable Hb in (4.7) by application of
~+1 = _ hit + (Hb )~+llb~.
(4.13)
With the values of (bl)~ +1 and b~ +~ known, the x-coordinates of points in the physical plane at
Problems involving the shallow water equations
Ill
the advanced time level corresponding to the fixed computational points ~ = ~ are deduced from (3.1) in the form x
=(bDf+1+~bf +l.
(4.14)
Values of hij at the computational points in the ~-y plane corresponding to (4.14) may then be found by linear interpolation between values of h at the nearest adjacent points in the physical plane at which the bathymetry is known. After this procedure, the values of hij so determined are available for the next complete advance of the solution through a time At. However, before the next advance in ~, ti must be updated by applying a discretization of (3.3) based upon A,ti + Be(0e t1~) + 8y(~~t~r) - .f~'6Y= - gE,(H ~~ ) -~ br~x P _
k[u2 + (~y)2]1/2
Et(a)
E,(/-~ ~)
(4.15)
Equation (4.15) is used to update fi at the u-points. Updated values of u may then be deduced from fi and ~ by using u
fi b/~.
(4.16)
Equation (4.15) is applied for i - - 3 ( 2 ) m - 1 and j = 3(2)n-2. By virtue of (3.8), the boundary value of O at i = 1 (referenced by the averaging operator in (4.15)) is identically zero. When applied at j = 3 and j = n - 2, the averaging operator references values of u at j = 1 and j = n and to overcome this difficulty, an appropriate one-sided definition of 8y is then used. Similarly, when i = m - 1, the averaging operator in (4.15) references a value of O outside the analysis area and, at this position, we also use a one-sided definition of 8~. With this procedure, values of t~ (and hence u) may be updated at all the interior u-points. The discretization of (3.4) is based upon
A,~ + ,~(u' ~) + 8 ~ ( ~ ~) + E,q ~i~)
+ bTg _ k[(ueY)2 + v2] I/2 Et (1~)
p
(4.17)
E, (/-I y)
Equation (4.17) is used for updating ~ at the v-points after which values of v are deduced by applying
v
b/~y.
(4.18)
Equation (4.17) is applied for i = 2 ( 2 ) m - 2 and j = 2 ( 2 ) n - 1. One-sided definitions of ~y are used when j = 2 and j = n - 1. After the updating of u and v, the updated value of U may be obtained by applying (3.5) in the discretized form 1 O = ~ - [ u - ( 8 , bl+~8, b ) - ( S y b l + ~ S y b ) ~eY].
(4.19)
The ensuing values of U at the u-points are then available for the next complete advance of the solution through a time At. Details relating to the stability characteristics of this computational scheme are given in [6]. In fact, stability is only conditional upon the time step being limited by the space increment and CAFVol.10,No. 2-.-B
112
B.JOHNS
the gravity wave speed. However, with excessive shearing of the coastal boundary, that condition may become severely restrictive. It will, of course, be noted that the updated values of the elevation in the computational ~-y plane correspond to the elevation at moving points in the physical x-y plane. Accordingly, after a prescribed period of integration, the programmed procedure is designed to interpolate from the grid-point values of ~ in the computational plane to the initial positions of the grid-points in the physical plane. This then enables a comparison to be made of the time history of the sea-surface elevation between simulations with and without the incorporation of continuous inland intrusion. 5. NUMERICAL EXPERIMENTS Numerical experiments are performed using the analysis area depicted in Fig. 1. Thus, we consider that part of the Indian Ocean and Bay of Bengal lying between 6 and 22°N. On average, the model has its fixed eastern boundary at about 300 km from the east coast of India. The surgegenerating capacity of the cyclone is recorded for about 3 days before landfall at the coast of Andhra Pradesh. During this period, the cyclone is moving along an approximately north-westerly track (shown in Fig. 1) with a translation speed of between 13 and 14 km hr -I. In our simulations, the earliest recorded position of the centre of the cyclone (corresponding to t -- 0 when the system is at rest) is about 1000 km south-east of the point of landfall. The time of landfall then corresponds to t = 70 hr. An idealized cyclonic wind field is simulated by applying an empirically-based representation of the azimuthal wind speed in the form
i Vo(~-, ./" )312
for r ~
IV0 exp ( ~ - ~ )
for r > R
V =-
] •
(5.1)
In (5.1), V0 is the maximum wind speed, R the radius of maximum wind, r is the distance from the centre of the cyclone and ~ is a length scale fixing the areal extent of the far wind-field. Using reports from the India Meterological Department, the available wind field data may best be represented by choosing Vo = 70ms -l, R = 80km and a = 240km. The associated windstress is calculated by a quadratic law with a uniform friction coefficient of 2.8 Equation (5.1) should be contrasted with the corresponding representation used in [6]. It will be observed that (5.1) yields a reduced value for the wind speed at distances from the cyclone centre for which r ~>R. This is an upgraded representation compared with that in [6] and is in closer accord with estimated wind speeds at distances from the cylone centre in excess of 250 km. Dube et al. [7] have shown that by calculating the applied wind-stress forcing from (5.1), the bulk of the surge response is confined to the Andhra coast, its quantitative value there being little changed in comparison with the simulations given in[6]. The finite-difference grid in the computational plane is chosen with m = 10 and n = 49. With our choice of open-sea boundaries, the east-west space increment varies between 18 and 40 km and Ay--36.5 km. Tests reported in[6] for a non-deforming grid confirm the adequacy of this resolution. A generally representative bathymetry has been used with the depth increasing from between 0 and 10 m at the coastline to in excess of 500 m in the southerly regions near the open-sea boundary ~ = 1. In actuality, the maximum depth of water is greater than that used in the model. However, from the point of view of wind-generated Surges, 500 m is deep and an increase in the depth has little effect on the surge response, the only consequence being to make the computation more expensive by further limiting the maximum admissible value of At. We restrict the region of inland intrusion to the coastline of Andhra Pradesh. This implies that variations in the coastline position are confined to values of ] between 21 and 35 inclusive. If j is outside this range, inland intrusion does not occur and the coastline consists of a perfectly reflecting vertical side-wall at which we have to satisfy (2.3) or, equivalently, (2.4) with s = ~. Thus, with the initial configuration of the coastline we take hL~= 0 for 21 ~ j ~ 35. During the updating procedure this means that re-computed values of hij in the computational plane are only required for 21 ~<] ~ 35.
× 10-3.
Problemsinvolvingthe shallowwater equations
113
In the numerical experiments, we have used constant values of the inland gradient, s. Our experience has shown that these must be compatible with the sea-floor gradient in the off-shore region. We have successfully performed simulations with s = 2 x 10-4, 3 x lO-4 and 5 x 104, the equilibrium depth at the first off-shore grid-point being, respectively, 2.5, 5 and 10 m. Off-shore depths in excess of these values have been found to lead to numerical problems that are probably associated with an unacceptably large discontinuity in the bottom gradient between the inland and sea-floor regions. With s = 2 x 10-4, Fig. 3(a) gives the variation of the coastline displacement at two stations (Kavali and Divi) that are, respectively, south and north of the point of landfall in Andhra Pradesh. This is calculated from the deforming coastline model (MD). Negative values of the displacement correspond to inland intrusion and positive values correspond to a seaward recession of the coastline. We also give an estimate of the coastline displacement based upon a computation made with the fixed coastline model (MF). This estimate is derived from a horizontal inland (or seaward) extrapolation given by the value of g(~ = 0, y, t)/s. We note that MD predicts a maximum inland intrusion at Kavali of 36km at t = 60hr whereas MF yields a significantly greater maximum intrusion of 42 km at t = 55 hr. Consequently, calculations based on MF result in the predicted maximum inland inundation occurring some 5 hr earlier than in MD. The nature of the resurgence phase is also noteworthy. This commences for t > 70 hr when the centre of the forcing cyclone leaves the analysis area. At Kavali, we see that MD predicts that the coastline has returned to its initial position at t -- 74 hr after which it recedes seaward by about 17 km during the next 7 hr. Subsequently, the position of the coastline again moves landward and is within about 5 km of its initial position at t = 100 hr. This resurgence phase should be compared with the extrapolative calculation based on the use of MF which reveals that the maximum seaward recession is 37 km at t = 80 hr. During the 100 hr of model integration, there is only a slight recession of the coastline at Divi and comparable inundation is predicted by both MD and MF for t > 40 hr. The maximum intrusion in MD again occurs about 5 hr later than in MF. It is informative to compare the intrusion curves in Fig. 3(a) with the variation in time of the
_
I'~
8 --
,'X
/ I
I I I
6--
I
I
I
,,~ 2/,
IIi ~ I~. i
4 --
o
ii
"}
.....
_z 12
~o ~
I
-"-2:.'..:, /I _./if_..-'" -..'X
x"::, ' [ /'
~m 8
-
",
~\-~ If..."! l l: :
'\:"111 .F I
-24
--
i \
',J
iV I
-36
20
40
"./
60
~0
100
20
TIME IN HOURS
~ 60 TiME IN HOURS
80
100
Fig. 3. (a) Variationof coastlinedisplacementwith s = 2 × 10'4 calculatedfromMD and MF. Kavali(MD) ; Kavali(kiF). . . . . ; Divi (kiD)..... ; DiviOdF) .... Co)Variationof sea surfaceelevationat initial position of coastline with s = 2 × 10" calculated from MD and kiF. Kavali (MD) ; Kavali (MF) ..... ; Divi (kiD) ..... ; Divi (kiF) ....
114
B. JOHNS
sea-surface elevation calculated from MD and projected onto the initial position of the coastline. This information is given for s = 2 × 10-4 in Fig. 3(b) together with the variation of the coastal elevation computed from MF. At Kavali, we note that MD predicts a maximum sea-surface elevation of almost 7 m at t = 60 hr which coincides with the time of maximum inland intrusion. On the other hand, MF predicts a maximum elevation of approximately 8.5 m at t = 55 hr. Consequently, MF overpredicts in comparison with MD presumably because of the unrealistic piling up of water at the vertical side-wall boundary. In MD this piling up does not occur since the water may flow freely inland thus reducing the maximum elevation in MD at the initial position of the coastline. The resurgence phase at Kavali must be interpreted with due caution when using the derived value of the sea-surface elevation at the initial position of the coastline. This is because the initial position of the coastline will ultimately become dry during recession in consequence of which both MF and MD will lead to a fictitious negative depth of water. This feature is apparent for MF when t >/68 hr and for MD when t/> 75 hr. In the case of MD, we have already noted that the coastline has returned to its initial position at t-~ 74 hr which is consistent with our interpretation of the extrapolated, and therefore fictitious, negative depth of water at the initial position of the coastline. Of the two models, however, MD is clearly the more realistic during this phase since its formulation is not based upon the assumption of there always being water at the initial position of the coastline. A comparison of these predictions with observational estimates of the flooding in Andhra Pradesh documented in[6] is interesting. Undoubtedly, MD predicts an excessive distance of inland surge penetration (reckoned, in actuality, to be a maximum of 25 km) together with an overprediction of the maximum surge height at the initial position of the coastline. This is possibly a consequence of a misrepresentation of the off-shore bathymetry and a value of s that is too small. On the other hand, it is encouraging to note that MD predicts the maximum inland flooding to occur some 5 hr later than MF. This latter prediction is, indeed, in closer accord with estimates during the flooding in Andhra Pradesh than is the prediction based upon MF. In the next experiment we have taken s = 3 x 10-4. Figure 4(a) gives the intrusion and recession curves calculated from both MD and MF. The maximum inland penetration at Kavali computed from MD is about 23 km compared with 25 km in MF. These both occur at t -~ 55 hr. An examination of the entry into the resurgence phase reveals that MF yields a maximum seaward recession of 26 km compared with 16 km in MD. This overestimate in MF might be
8I
30
-!
2O
i
o z z
~o a. m
! i
N -tO
-20
:'. -'~'"
I
m
1 2O
"1
40 60 TIME IN HOURS
]
I
80
100
-8
I
20
I
40
60
~'-'1
TIME IN HOURS
Fig. 4. (a) s = 3 × 10-4: (b) s = 3 × l0 -4
80
J l
100
115
Problems involving the shallow water equations
r~
8 --
!1 z
"x
./".\
0
0 ¸
i
- - ~ \ .
,..J
/
-4 I
-
\~. :;: .i
.5
I l I l l i vl I
2-8 -12
I
20
I
40 60 TIME IN HOURS
I
80
I
100
I 20
l
40
IL}l
60
80
1 100
TIME IN HOURS
Fig. 5. (a) s = 5 × 10-4: (b) s = 5 x 10-4
attributed to the artificial reflection of the flood wave from a vertical side-wall compared with the more realistic simulation of the run-down of water in MD. It is noteworthy that the maximum inland penetration predicted by MD and MF at Divi is, respectively, 16 and 14 km. However, the time of maximum inundation is attained more than 5 hr later in MD compared with MF. The resurgence phase at Divi is not characterized by any significant seaward recession of the coastline from its initial position. A consequence of this is that the peak flood wave is longerqasting at Divi and inland inundation is predicted to persist for 55 hr compared with about 45 hr at Kavali. In Fig. 4(b), we give the corresponding variation of the sea-surface elevation projected onto the initial position of the coastline. It is seen that MD and MF yield respective peak elevations of about 6.7 and 7.4 m at Kavali. At Divi, MD produces a peak elevation of 4.8 m whilst MF leads to 4.3 m. The peak in MF occurs some 5 hr earlier than in MD. Finally, we consider the case of s = 5 × 10-4 with an off-shore depth of 10 m. The intrusion and recession curves are given in Fig. 5(a) from which we see that the differences between MD and MF are far less marked. With this parameter setting, it is clear that MF (relative to MD) will yield an acceptable estimate of the extent of inland flooding. At Kavali, for example, both models predict an inland penetration of 11 km at t = 55 hr. The corresponding variation of the sea-surface elevation at the initial position of the coastline is given in Fig. 5(b). As with the intrusion curves, the differences between MD and MF are not significant. As is to be expected the two models are tending to produce coincident results as the gradient, s, is increased. Acknowledgement--This work was supported by the U.K. Overseas Development Administration under the auspices of the Imperial College Delhi Committee.
REFERENCES 1. B. Johns, The modeling of tidal flow in a channel using a turbulence energy closure scheme. J. Phys. Oceanogr. 8, 1042-1049 (1978). 2. B. Johns and R. J. Jefferson, The numerical modeling of surface wave propagation in the surf zone. J. Phys. Oceanogr. 10, 1061-1069 (1980). 3. R. O. Reid and B. R. Bodine, Numerical model for storm surges in Galveston Bay. Proc. Am. S. Civil Engng, J. Waterways Harbors Div. 94, 33-5"}. (1968).
116
B. JOHNS
4. R. A. Flather and N. S. Heaps, Tidal computations in Morecambe Bay. Geophys. J.R. Astr. Soc. 42, 489-517 (1975). 5. D. R. Lynch and W. G. Gray, Finite element simulation of flow in deforming regions. J. Comput. Phys. 36, 135-153 (19S0). 6. B. Johns, S. K. Dube, U. C. Mohanty and P. C. Sinha, Numerical simulation of the surge generated by the 1977 Andhra cyclone. Q. J. R. Met. Soc. 107, 919-934 (1981). 7. S. K. Dube, P. C. Sinha and A. D. Rao, The response of different wind-stress forcings on the surges along the east coast of India. Mausam (Formerly Indian Journal of Meteorology, Hydrology and Geophysics) 32, 315-320 (1981).
0045.-79301821020105-12503.0010 Pergamon Press Ltd.
THE SIMULATION OF A CONTINUOUSLY DEFORMING LATERAL BOUNDARY IN PROBLEMS INVOLVING THE SHALLOW WATER EQUATIONS B. JOHNS Department of Meteorology, University of Reading, England and S. K. DUBE, P. C. SINHA, U. C. MOHANTYand A. D. RAO Centre for Atmospheric Sciences, Indian Institute of Technology, Delhi, India
(Received 5 August 1981) kbstract--A finite-differencemethod is given for the treatment of a continuously deforming lateral fluid boundary in computer simulations of systems involving the shallow water equations. The technique is described in the context of a numerical storm surge model for the east coast of India. Using a forcing wind-stressdistributionrepresentative of the 1977Andhra cyclone, a comparisonis made between simulations using models both with and without inland intrusion of water along the coast of Andhra Pradesh.
I. INTRODUCTION The treatment of deforming fluid boundaries is an essential part of computer simulations in many dynamical problems. A classical area of note is in the context of surface wave propagation where the position of the free surface is unknown a prior/. One of the difficulties here is that the free surface boundary conditions have to be applied at a level which may not coincide with any of the computational levels in a scheme of numerical solution. In the case of shallow water wave propagation, these problems may be overcome by firstmaking a coordinate transformation to ensure that the free surface always corresponds to a fixed computational level in the numerical scheme as described in[I,2]. A similar difficultyarises in problems involving the deformation of lateralfluid boundaries. An important example of this occurs in numerical storm surge prediction models. Frequently, the lateralboundaries in these are taken to be vertical side-walls through which there is no flux of water and many such studies have been undertaken. In actuality, however, the water will usually move continuously inland and the use of idealized vertical side-walls may lead to a misrepresentation of the surge development. In particular, a simplistic horizontal inland extrapolation of the predicted surge height at the fixed side-wall in order to estimate the distance of inland penetration does not in itselftake account of the dynamics of the intrusion process. The recognition of this has led several authors to develop numerical models that allow the flow of water across the initialposition of the model coastline. This has been accomplished in [3,4] by arranging for the advancing water front to move discontinuously from one grid-point to another according to a prescribed criterion. However, to obtain acceptable accuracy, this procedure requires an extremely fine grid-spacing which may not be otherwise necessary. More recently, a finiteelement simulation of shallow water flow in continuously deforming regions has been proposed by Lynch and Gray[5]. The attractions of their method are considerable because it accommodates a continuous motion of the shoreline. However, the computational overheads characteristic of finite element procedures may make the method expensive. In the present work, we describe a finite-difference method which also models a continuously moving lateral boundary. Although admitting of application in various contexts, we apply it to the numerical simulation of the surge generated by the 1977 Andhra cyclone which struck the east coast of India and led to disastrous inland flooding. An attempt to simulate this event numerically has previously been made by Johns et al. [6]. Their treatment incorporated fixed vertical side-wails through which there is no flux of water. However, they did employ a curvilinear representation of the lateral boundaries which allowed these to be modelled in a fairly realisticway. Using this framework, finite-differencecalculations were undertaken with wind-stress forcing data appropriate to the 1977 Andhra cyclone. 105
106
B. JOHNS
The method we employ here is a generalization of the curvilinear boundary technique described in[6] and is conceptually similar to the free surface treatment given in[I]. Using an upgraded representation of the data for the wind-stress forcing field, we have illustrated the method by comparing simulations of the surge with and without a deforming fluid boundary along the coast of Andhra Pradesh. Our conclusion is that the incorporation of a deforming coastline in storm surge models may lead to an important change in the prediction detail compared with that from a fixed coastline model. 2. FORMULATION All conditions are referred to a set of rectangular Cartesian coordinates in which the origin, O is within the equilibrium level of the free surface. Ox points towards the east, Oy towards the north and Oz is directed vertically upwards. The displaced position of the free surface is given by z = ~(x, y, t) and the position of the sea-floor by z = - h ( x , y). A western coastal boundary (the east coast of India) is situated at x = bl(y, t) and an eastern open-sea boundary is situated at x = b2(y). Southern and northern open-sea boundaries are at y = 0 and y = I, respectively. This configuration is shown in Fig. 1. The depth-averaged components of velocity, u and v, then satisfy
__aU+ u au + au 3~ + l_~ {,r~ _ kpu (u2 + v2),12} at -~x V ~ - y - [V = - g ax Hp
(2.1)
u2}.
(2.2)
and
_av+uav+ 3t
av
-~x
_
V ~y + [U =
a_~+ 1 {~.~_kpv(u2+v2)
gay
np
In (2.1) and (2.2), [ denotes the Coriolis parameter, the pressure is taken as hydrostatic and H is the total depth ~ + h. (~, t~) denotes the applied surface wind-stress and the bottom stress is parameterized in terms of a quadratic law. p denotes the water density and the friction coefficient, k, is taken as 2.6 x 10-3. The kinematical condition at the western coastal boundary requires either that the coastline consist of a vertical side-wall through which there is no flux of water or, if there is inland intrusion of water, that the shoreline move with the same velocity as that of its constituent water particles. In the first case, we require Obi u-v-~:--O at x = b t ( y , 0 )
(2.3)
_
oy
25*
y=t
BANGLADESH
-~ .-:...:f""~
INDIA/~//~/'-~" ' "v~!,. x : b,(y, t ) ~ /
". B U R M A
//
/i
f" iI
- - 20"
i"
!"~
c'"
-- 15"
X= bz(y)
I
,,
i'/,o %... '"
I
80°
....
°
d
"~ - - - ~ x
I
85" F i g . 1. T h e
I
[
90" analysis
95" area.
~uu" 5"
Problemsinvolvingthe shallowwaterequations
107
and, in the second Oy at x=b~(y,t).
(2.4)
Thus, in (2.3) the value of b~ is invariant in time and is always equal to its initial value. In (2.4), the position of the coastline changes and the kinematical condition must be applied at its a priori unknown location. At the eastern open-sea boundary, a radiation condition is applied that is identical to that used by Johns et al. [6] in their coastal zone model. This leads to
062_ ( ~ ) 1/2 u-v~y ~=0
at x=b2(y).
(2.5)
Conceptually similar radiation conditions are applied at the southern and northern open-sea boundaries to yield v+(~)
112 ,=0
at y=O
(2.6)
at y = l .
(2.7)
1/2 v-(~)
~=0
To accompany (2.1) and (2.2), the equation of continuity is applied in the vertically integrated form
o/-/+o__ Ot
Ox (Hu)+
(Hv) = O.
(2.8)
When bt is temporally invariant, the above formulation is the same as that of the coastal zone model described in [6] except that the northern boundary is now an open-sea boundary.
3. COORDINATE TRANSFORMATION A coordinate transformation is introduced which is a generalization of that used in[6] to facilitate the numerical treatment of an irregular boundary configuration. We now write _ x - b~(y, t)
-
b(y,t)
(3.1)
where b(y, t) = b2(y)- b~(y, t).
(3.2)
Thus, the western and eastern boundaries correspond, respectively, to ~ = 0 and ~ = 1. Equations (2.1) and (2.2) may firstly be written in flux form and then, with ~, y and t as new independent variables, may be transformed to Oti 0 . . ti)+ 0 ( v ~ ) _ f ~ = _ O~ b'r~_kfi . 2 ± . 2 , , / 2 ~-+-~(U 0), gH o!~+ P H ( u -,v j
(3.3)
and
p - -~ (u 2 + v2) 112,
(3.4)
106
B. JOHNS
where
(3.5) Hbu} = Hbv " =
(3.6)
Equation (2.8) is readily found to yield
o~(Hb)+
(Hb U)+
05)=0.
(3.7)
From (2.3), (2.4) and (3.5) we see that the kinematical condition at the curvilinear boundary is always satisfied provided that U=0
at
~=0.
(3.8)
The radiation condition (2.5) at the eastern boundary is satisfied if 1/2
As in [6], the numerical procedure will be to solve (3.3)-(3.7) in a rectangular computational area 0 ~<~ ~< 1, 0 ~< y ~< I. On satisfying (3.8) and (3.9), the conditions (2.3)-(2.5) on the western and eastern coastlines are automatically fulfilled. The conditions (2.6) and (2.7) are satisfied directly when implementing the numerical scheme. In the case of inland intrusion of water, it remains to establish a procedure for determining the changing position of the coastline. The kinematical condition (2.4) has already been satisfied by taking U = 0 at ~ = 0 and a further requirement is that the depth of the water be zero at the coastline. This leads to H=0
at
x=b~(y,t),
(3.10)
or, equivalently,
(,(~, = O, y, t) + h[x = b~(y, t), y] = O.
(3.11)
Depending on whether b,(y, t),~ ba(y,0), we may then either interpolate or use new inland orographical data to fix the value of h[x = b,(y, t), y]. For present purposes, however, the simplest procedure is to differentiate (3.11) with respect to t. This leads to O--~-(~= O, y, t) +
s=O
(3.12)
where S=
• =b,(y.,)"
(3.13)
If s is prescribed, (3.12) yields a prognostic equation for b,. The simplest case corresponds to a constant value of s when (3.12) immediately integrates to ba(y, t) = bl(y, 0) _ 1 ~(~ = 0, y, t).
(3.14)
Then, if ~(~ = 0, y, t ) < 0, the sea-surface at the shoreline is depressed and the shoreline has
Problemsinvolvingthe shallowwater equations
109
consequently receded from its initial position. If ~(~ = 0, y, t ) > 0, there is a positive surge; the elevation at the shoreline is raised above its equilibrium level and there is a corresponding inland penetration of water. The numerical procedure to be described in Section 4 leads to an updating scheme for ~(~ = 0, y, t). Our formulation may, therefore, be used to monitor the continuous movement of the shoreline. This procedure is based rigorously on the fulfilment of the physically correct coastline conditions (2.4) and (3.10) and allows full account to be taken of the dynamics of the intrusion (or recession) process. 4. NUMERICAL PROCEDURE We define discrete coordinate points in the ~-y plane by =
¢i
=
(i -
1 ) A s ¢,
i
=
1, 2 ....
m
;A¢ = ll(m ll(n
y = y~ = (j - l)Ay, j = I, 2 .... n ;Ay =
-
-
1)'~ I) J"
(4.1)
A sequence of time instants is defined by t = tp = p A t ,
p = O, 1 . . . . .
(4.2)
X~.
(4.3)
For any variable, X, we write
X(~i,
Y~,
tp)
=
In the ~-y plane we use a staggered grid where there are three distinct types of computational point. With i even and j odd, the point is a ~-point at which ~ is computed. If i is odd and j is odd the point is a u-point at which both u and U are computed. If i is even and j is even, the point is a v-point at which v is computed. Along s¢ = 0 we, therefore, have only u-points and at each of these we have U = 0 identically. We choose m to be even so that ~ = 1 consists of both ~-points and v-points. We also choose n to be odd thus ensuring that there are only ~-points and u-points along y = 0 and y = i. This arrangement of grid-points is shown in Fig. 2. With this scheme, it will be noted that computed elevations are not carried along the coastline ~ -- 0. In fact, an updating procedure is never invoked along ~ -- 0 where (in the case of inland intrusion) the depth is zero. This feature in our scheme avoids problems relating to the singular behavior of (2.1) and (2.2) when H = 0 to which reference has been made in[5]. In order to overcome this difficulty, it was suggested there that a physical re-formulation of the bottom friction may be required near the shoreline--a view to which we do not necessarily subscribe. In order to describe the finite-difference equations, we define difference operators by
A'X = (X~+' - X~)IAt
1 (4.4)
80~ = (X,P+~,j - X~-~,j)/(2AO t. 8yX = (X~+~ - X£~-d/(2,~y)J
.T:5
•
[
3=4 J=3
()
s-,
3=2
.T=I
Ay • I:1 Fig.
~p
I:2 1:3 1:4 1:5 I=6 2. Grid-pointarrangement.
110
B. JOHNS
Averaging operations are defined by
~ = ~(xr+,.j + xr-,.~) 1 ~?' =
'(xb+,
+
x~J-,) t
(4.5)
and a shift operator is defined by
EtX = vP'+l ,~,j .
(4.6)
The discretization of (3.7) is then based upon A,(Hb) + a~ (/4~b U) + ar 03) = 0.
(4.7)
Equation (4.7) is used in the updating procedure at points corresponding to i--2(2)m- 2 and j = 3(2)n - 2. The elevation at j = 1 is determined from (2.6) which, in practice, is applied at j = 2 and replaced by l(g~
'/2
~/
(W"+W')+vf2=o,
(4.8)
thus leading to an updating procedure for the elevation along the southern open-sea boundary. Likewise, the elevation at j = n is determined from (2.7) which is then replaced by 1( g ],,2 p+l X \z----.,,,.._,/ ( W ' + v~,,.-2,v~._, = O.
_
2
(4.9)
The elevation along ~ = 1 determined from (3.9) which is replaced by 1
"
g
xl12
-2 [-~\h,.-1,J
(~mjP+l+ ~m-2,j)p+* b Pj IJ,,-I,j p = O,
(4.10)
thus yielding updated elevations for i = m. In contrast with[6], the time-dependent nature of b should be noted in (4.10). Additionally, and most importantly, the application of (4.7)-(4.10) implies that the position of the grid points is changing with time in the physical x-y plane although they are fixed in the computational 6-y plane. This means that the grid-point values of h in the computational plane will be changing with time whereas they are independent of time at fixed points in the physical plane. Accordingly, after each updating of the prognostic variables it is necessary to re-compute new values of hit in the computational plane. This can only be done when the change has been determined in b~ corresponding to the current advance of the solution through a time At. This change may be computed from (3.12) by applying
(b,),~÷' = (b,),~ - ~ (¢V I - ¢~j).
(4.11)
The value of ~ 1 in (4.11) is not carried in our computational scheme. It is, therefore, determined by linear extrapolation from the adjacent ~-points by using ~¢j+1= ~(3~j+1_ ~j+l).
(4.12)
Values of ~,~+~ at the ~-points are determined from the prognostic variable Hb in (4.7) by application of
~+1 = _ hit + (Hb )~+llb~.
(4.13)
With the values of (bl)~ +1 and b~ +~ known, the x-coordinates of points in the physical plane at
Problems involving the shallow water equations
Ill
the advanced time level corresponding to the fixed computational points ~ = ~ are deduced from (3.1) in the form x
=(bDf+1+~bf +l.
(4.14)
Values of hij at the computational points in the ~-y plane corresponding to (4.14) may then be found by linear interpolation between values of h at the nearest adjacent points in the physical plane at which the bathymetry is known. After this procedure, the values of hij so determined are available for the next complete advance of the solution through a time At. However, before the next advance in ~, ti must be updated by applying a discretization of (3.3) based upon A,ti + Be(0e t1~) + 8y(~~t~r) - .f~'6Y= - gE,(H ~~ ) -~ br~x P _
k[u2 + (~y)2]1/2
Et(a)
E,(/-~ ~)
(4.15)
Equation (4.15) is used to update fi at the u-points. Updated values of u may then be deduced from fi and ~ by using u
fi b/~.
(4.16)
Equation (4.15) is applied for i - - 3 ( 2 ) m - 1 and j = 3(2)n-2. By virtue of (3.8), the boundary value of O at i = 1 (referenced by the averaging operator in (4.15)) is identically zero. When applied at j = 3 and j = n - 2, the averaging operator references values of u at j = 1 and j = n and to overcome this difficulty, an appropriate one-sided definition of 8y is then used. Similarly, when i = m - 1, the averaging operator in (4.15) references a value of O outside the analysis area and, at this position, we also use a one-sided definition of 8~. With this procedure, values of t~ (and hence u) may be updated at all the interior u-points. The discretization of (3.4) is based upon
A,~ + ,~(u' ~) + 8 ~ ( ~ ~) + E,q ~i~)
+ bTg _ k[(ueY)2 + v2] I/2 Et (1~)
p
(4.17)
E, (/-I y)
Equation (4.17) is used for updating ~ at the v-points after which values of v are deduced by applying
v
b/~y.
(4.18)
Equation (4.17) is applied for i = 2 ( 2 ) m - 2 and j = 2 ( 2 ) n - 1. One-sided definitions of ~y are used when j = 2 and j = n - 1. After the updating of u and v, the updated value of U may be obtained by applying (3.5) in the discretized form 1 O = ~ - [ u - ( 8 , bl+~8, b ) - ( S y b l + ~ S y b ) ~eY].
(4.19)
The ensuing values of U at the u-points are then available for the next complete advance of the solution through a time At. Details relating to the stability characteristics of this computational scheme are given in [6]. In fact, stability is only conditional upon the time step being limited by the space increment and CAFVol.10,No. 2-.-B
112
B.JOHNS
the gravity wave speed. However, with excessive shearing of the coastal boundary, that condition may become severely restrictive. It will, of course, be noted that the updated values of the elevation in the computational ~-y plane correspond to the elevation at moving points in the physical x-y plane. Accordingly, after a prescribed period of integration, the programmed procedure is designed to interpolate from the grid-point values of ~ in the computational plane to the initial positions of the grid-points in the physical plane. This then enables a comparison to be made of the time history of the sea-surface elevation between simulations with and without the incorporation of continuous inland intrusion. 5. NUMERICAL EXPERIMENTS Numerical experiments are performed using the analysis area depicted in Fig. 1. Thus, we consider that part of the Indian Ocean and Bay of Bengal lying between 6 and 22°N. On average, the model has its fixed eastern boundary at about 300 km from the east coast of India. The surgegenerating capacity of the cyclone is recorded for about 3 days before landfall at the coast of Andhra Pradesh. During this period, the cyclone is moving along an approximately north-westerly track (shown in Fig. 1) with a translation speed of between 13 and 14 km hr -I. In our simulations, the earliest recorded position of the centre of the cyclone (corresponding to t -- 0 when the system is at rest) is about 1000 km south-east of the point of landfall. The time of landfall then corresponds to t = 70 hr. An idealized cyclonic wind field is simulated by applying an empirically-based representation of the azimuthal wind speed in the form
i Vo(~-, ./" )312
for r ~
IV0 exp ( ~ - ~ )
for r > R
V =-
] •
(5.1)
In (5.1), V0 is the maximum wind speed, R the radius of maximum wind, r is the distance from the centre of the cyclone and ~ is a length scale fixing the areal extent of the far wind-field. Using reports from the India Meterological Department, the available wind field data may best be represented by choosing Vo = 70ms -l, R = 80km and a = 240km. The associated windstress is calculated by a quadratic law with a uniform friction coefficient of 2.8 Equation (5.1) should be contrasted with the corresponding representation used in [6]. It will be observed that (5.1) yields a reduced value for the wind speed at distances from the cyclone centre for which r ~>R. This is an upgraded representation compared with that in [6] and is in closer accord with estimated wind speeds at distances from the cylone centre in excess of 250 km. Dube et al. [7] have shown that by calculating the applied wind-stress forcing from (5.1), the bulk of the surge response is confined to the Andhra coast, its quantitative value there being little changed in comparison with the simulations given in[6]. The finite-difference grid in the computational plane is chosen with m = 10 and n = 49. With our choice of open-sea boundaries, the east-west space increment varies between 18 and 40 km and Ay--36.5 km. Tests reported in[6] for a non-deforming grid confirm the adequacy of this resolution. A generally representative bathymetry has been used with the depth increasing from between 0 and 10 m at the coastline to in excess of 500 m in the southerly regions near the open-sea boundary ~ = 1. In actuality, the maximum depth of water is greater than that used in the model. However, from the point of view of wind-generated Surges, 500 m is deep and an increase in the depth has little effect on the surge response, the only consequence being to make the computation more expensive by further limiting the maximum admissible value of At. We restrict the region of inland intrusion to the coastline of Andhra Pradesh. This implies that variations in the coastline position are confined to values of ] between 21 and 35 inclusive. If j is outside this range, inland intrusion does not occur and the coastline consists of a perfectly reflecting vertical side-wall at which we have to satisfy (2.3) or, equivalently, (2.4) with s = ~. Thus, with the initial configuration of the coastline we take hL~= 0 for 21 ~ j ~ 35. During the updating procedure this means that re-computed values of hij in the computational plane are only required for 21 ~<] ~ 35.
× 10-3.
Problemsinvolvingthe shallowwater equations
113
In the numerical experiments, we have used constant values of the inland gradient, s. Our experience has shown that these must be compatible with the sea-floor gradient in the off-shore region. We have successfully performed simulations with s = 2 x 10-4, 3 x lO-4 and 5 x 104, the equilibrium depth at the first off-shore grid-point being, respectively, 2.5, 5 and 10 m. Off-shore depths in excess of these values have been found to lead to numerical problems that are probably associated with an unacceptably large discontinuity in the bottom gradient between the inland and sea-floor regions. With s = 2 x 10-4, Fig. 3(a) gives the variation of the coastline displacement at two stations (Kavali and Divi) that are, respectively, south and north of the point of landfall in Andhra Pradesh. This is calculated from the deforming coastline model (MD). Negative values of the displacement correspond to inland intrusion and positive values correspond to a seaward recession of the coastline. We also give an estimate of the coastline displacement based upon a computation made with the fixed coastline model (MF). This estimate is derived from a horizontal inland (or seaward) extrapolation given by the value of g(~ = 0, y, t)/s. We note that MD predicts a maximum inland intrusion at Kavali of 36km at t = 60hr whereas MF yields a significantly greater maximum intrusion of 42 km at t = 55 hr. Consequently, calculations based on MF result in the predicted maximum inland inundation occurring some 5 hr earlier than in MD. The nature of the resurgence phase is also noteworthy. This commences for t > 70 hr when the centre of the forcing cyclone leaves the analysis area. At Kavali, we see that MD predicts that the coastline has returned to its initial position at t -- 74 hr after which it recedes seaward by about 17 km during the next 7 hr. Subsequently, the position of the coastline again moves landward and is within about 5 km of its initial position at t = 100 hr. This resurgence phase should be compared with the extrapolative calculation based on the use of MF which reveals that the maximum seaward recession is 37 km at t = 80 hr. During the 100 hr of model integration, there is only a slight recession of the coastline at Divi and comparable inundation is predicted by both MD and MF for t > 40 hr. The maximum intrusion in MD again occurs about 5 hr later than in MF. It is informative to compare the intrusion curves in Fig. 3(a) with the variation in time of the
_
I'~
8 --
,'X
/ I
I I I
6--
I
I
I
,,~ 2/,
IIi ~ I~. i
4 --
o
ii
"}
.....
_z 12
~o ~
I
-"-2:.'..:, /I _./if_..-'" -..'X
x"::, ' [ /'
~m 8
-
",
~\-~ If..."! l l: :
'\:"111 .F I
-24
--
i \
',J
iV I
-36
20
40
"./
60
~0
100
20
TIME IN HOURS
~ 60 TiME IN HOURS
80
100
Fig. 3. (a) Variationof coastlinedisplacementwith s = 2 × 10'4 calculatedfromMD and MF. Kavali(MD) ; Kavali(kiF). . . . . ; Divi (kiD)..... ; DiviOdF) .... Co)Variationof sea surfaceelevationat initial position of coastline with s = 2 × 10" calculated from MD and kiF. Kavali (MD) ; Kavali (MF) ..... ; Divi (kiD) ..... ; Divi (kiF) ....
114
B. JOHNS
sea-surface elevation calculated from MD and projected onto the initial position of the coastline. This information is given for s = 2 × 10-4 in Fig. 3(b) together with the variation of the coastal elevation computed from MF. At Kavali, we note that MD predicts a maximum sea-surface elevation of almost 7 m at t = 60 hr which coincides with the time of maximum inland intrusion. On the other hand, MF predicts a maximum elevation of approximately 8.5 m at t = 55 hr. Consequently, MF overpredicts in comparison with MD presumably because of the unrealistic piling up of water at the vertical side-wall boundary. In MD this piling up does not occur since the water may flow freely inland thus reducing the maximum elevation in MD at the initial position of the coastline. The resurgence phase at Kavali must be interpreted with due caution when using the derived value of the sea-surface elevation at the initial position of the coastline. This is because the initial position of the coastline will ultimately become dry during recession in consequence of which both MF and MD will lead to a fictitious negative depth of water. This feature is apparent for MF when t >/68 hr and for MD when t/> 75 hr. In the case of MD, we have already noted that the coastline has returned to its initial position at t-~ 74 hr which is consistent with our interpretation of the extrapolated, and therefore fictitious, negative depth of water at the initial position of the coastline. Of the two models, however, MD is clearly the more realistic during this phase since its formulation is not based upon the assumption of there always being water at the initial position of the coastline. A comparison of these predictions with observational estimates of the flooding in Andhra Pradesh documented in[6] is interesting. Undoubtedly, MD predicts an excessive distance of inland surge penetration (reckoned, in actuality, to be a maximum of 25 km) together with an overprediction of the maximum surge height at the initial position of the coastline. This is possibly a consequence of a misrepresentation of the off-shore bathymetry and a value of s that is too small. On the other hand, it is encouraging to note that MD predicts the maximum inland flooding to occur some 5 hr later than MF. This latter prediction is, indeed, in closer accord with estimates during the flooding in Andhra Pradesh than is the prediction based upon MF. In the next experiment we have taken s = 3 x 10-4. Figure 4(a) gives the intrusion and recession curves calculated from both MD and MF. The maximum inland penetration at Kavali computed from MD is about 23 km compared with 25 km in MF. These both occur at t -~ 55 hr. An examination of the entry into the resurgence phase reveals that MF yields a maximum seaward recession of 26 km compared with 16 km in MD. This overestimate in MF might be
8I
30
-!
2O
i
o z z
~o a. m
! i
N -tO
-20
:'. -'~'"
I
m
1 2O
"1
40 60 TIME IN HOURS
]
I
80
100
-8
I
20
I
40
60
~'-'1
TIME IN HOURS
Fig. 4. (a) s = 3 × 10-4: (b) s = 3 × l0 -4
80
J l
100
115
Problems involving the shallow water equations
r~
8 --
!1 z
"x
./".\
0
0 ¸
i
- - ~ \ .
,..J
/
-4 I
-
\~. :;: .i
.5
I l I l l i vl I
2-8 -12
I
20
I
40 60 TIME IN HOURS
I
80
I
100
I 20
l
40
IL}l
60
80
1 100
TIME IN HOURS
Fig. 5. (a) s = 5 × 10-4: (b) s = 5 x 10-4
attributed to the artificial reflection of the flood wave from a vertical side-wall compared with the more realistic simulation of the run-down of water in MD. It is noteworthy that the maximum inland penetration predicted by MD and MF at Divi is, respectively, 16 and 14 km. However, the time of maximum inundation is attained more than 5 hr later in MD compared with MF. The resurgence phase at Divi is not characterized by any significant seaward recession of the coastline from its initial position. A consequence of this is that the peak flood wave is longerqasting at Divi and inland inundation is predicted to persist for 55 hr compared with about 45 hr at Kavali. In Fig. 4(b), we give the corresponding variation of the sea-surface elevation projected onto the initial position of the coastline. It is seen that MD and MF yield respective peak elevations of about 6.7 and 7.4 m at Kavali. At Divi, MD produces a peak elevation of 4.8 m whilst MF leads to 4.3 m. The peak in MF occurs some 5 hr earlier than in MD. Finally, we consider the case of s = 5 × 10-4 with an off-shore depth of 10 m. The intrusion and recession curves are given in Fig. 5(a) from which we see that the differences between MD and MF are far less marked. With this parameter setting, it is clear that MF (relative to MD) will yield an acceptable estimate of the extent of inland flooding. At Kavali, for example, both models predict an inland penetration of 11 km at t = 55 hr. The corresponding variation of the sea-surface elevation at the initial position of the coastline is given in Fig. 5(b). As with the intrusion curves, the differences between MD and MF are not significant. As is to be expected the two models are tending to produce coincident results as the gradient, s, is increased. Acknowledgement--This work was supported by the U.K. Overseas Development Administration under the auspices of the Imperial College Delhi Committee.
REFERENCES 1. B. Johns, The modeling of tidal flow in a channel using a turbulence energy closure scheme. J. Phys. Oceanogr. 8, 1042-1049 (1978). 2. B. Johns and R. J. Jefferson, The numerical modeling of surface wave propagation in the surf zone. J. Phys. Oceanogr. 10, 1061-1069 (1980). 3. R. O. Reid and B. R. Bodine, Numerical model for storm surges in Galveston Bay. Proc. Am. S. Civil Engng, J. Waterways Harbors Div. 94, 33-5"}. (1968).
116
B. JOHNS
4. R. A. Flather and N. S. Heaps, Tidal computations in Morecambe Bay. Geophys. J.R. Astr. Soc. 42, 489-517 (1975). 5. D. R. Lynch and W. G. Gray, Finite element simulation of flow in deforming regions. J. Comput. Phys. 36, 135-153 (19S0). 6. B. Johns, S. K. Dube, U. C. Mohanty and P. C. Sinha, Numerical simulation of the surge generated by the 1977 Andhra cyclone. Q. J. R. Met. Soc. 107, 919-934 (1981). 7. S. K. Dube, P. C. Sinha and A. D. Rao, The response of different wind-stress forcings on the surges along the east coast of India. Mausam (Formerly Indian Journal of Meteorology, Hydrology and Geophysics) 32, 315-320 (1981).