Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

The numerical solution of a system of diffusion/convection equations describing coastal circulation Nicol, Thomas O. 1983-04-21

You don't seem to have a PDF reader installed, try download the pdf

Item Metadata


831-UBC_1983_A6_7 N52.pdf [ 3.65MB ]
JSON: 831-1.0095808.json
JSON-LD: 831-1.0095808-ld.json
RDF/XML (Pretty): 831-1.0095808-rdf.xml
RDF/JSON: 831-1.0095808-rdf.json
Turtle: 831-1.0095808-turtle.txt
N-Triples: 831-1.0095808-rdf-ntriples.txt
Original Record: 831-1.0095808-source.json
Full Text

Full Text

C . l THE NUMERICAL SOLUTION OF A SYSTEM OF EQUATIONS DESCRIBING COASTAL by THOMAS 0. NICOL B.Sc, The University Of British Columbia, 1976 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES Institute Of Applied Mathematics And Statistics We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA October 1983 DIFFUSION/CONVECTION CIRCULATION © Thomas 0. Nicol," 1983 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the Head of my Department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Department of Institute Of Applied Mathematics And Statistics The University of British Columbia 2075 Wesbrook Place Vancouver, Canada V6T 1W5 Date: October 14, 1983 i i Abstract A mathematical model describing nearshore ocean currents is examined. The motivation for the problem is discussed and a derivation of the model equations presented. The model equations consist of a pair of coupled, nonlinear, parabolic partial differential equations. A numerical method, the Crank-Nicolson finite difference scheme, is presented for solving the corresponding linear equations, and the convergence and stability of the method discussed. Methods for dealing with the nonlinear terms, and their effects on accuracy and stability, are examined. The initial and boundary conditions of the model equations present special problems, and we describe methods to solve them. A more efficient finite difference scheme for our equations, based on the Crank-Nicolson formula, is introduced, and its advantages discussed. A computer program incorporating these methods is developed to solve the model equations, and results for test data presented. We conclude with recommendations for additions to the program to model the currents more accurately. iii Table of Contents Abstract i i List of Figures v Acknowledgements v I. INTRODUCTION 1 11 . CHAPTER 1 4 - 2.1 Derivation Of The Model 4 2.2 The Domain Of The Model 8 2.3 Boundary Conditions2.4 Initial Conditions 9 2.5 How The Model Will Be Used 10 III. CHAPTER 2 13 3.1 Notation3.2 Finite Differences 14 3.3 Consistency, Stability And Convergence 19 3.4 Variable Coefficients 23.5 Stability Of Boundary Conditions 25 3.6 Nonlinear Terms 27 3.7 Solving The Equations 9 IV. CHAPTER 3 2 4.1 Introduction To ADI 33 4.2 Accuracy And Stability 5 4.3 Computational Details And Sample Problems 37 4.4 Applying ADI To The Model Equations 42 V. CHAPTER 4 4 5.1 Preliminary Tests 45.2 Initial Conditions At The Land Boundary 46 5.3 Finding H On The Boundaries 49 5.4 Boundary Condition For S 52 5.5 Tests Of The Model Equations 3 5.6 Interpretation Of Results 4 5.7 Recommendations For Proceeding 57 BIBLIOGRAPHY 78 iv List of Figures 1. Domain of the Coastal Circulation Model 11 2. Domain of the Numerical Model 12 3. Test 1 : S at time step 2 60 4. Test 1 : S at time step 55. Test 1 : S at time step 10 1 6. Test 1 : S at time step 15 2 7. Test 1 : S at time step 20 68. Test 1 : S at time step 30 3 9. Test 1 : S at time step 4010. Test 1: S at time step 50 4 11. Test 1 : H at time step 2 612. Test 1 : H at time step 5 5 13. Test 1 : H at time step 1014. Test 1 : H at time step 15 6 15. Test 1: H at time step 20 616. Test 1 : H at time step 30 7 17. Test 1 : H at time step 4018. Test 1 : H at time step 50 8 19. Test 1: H2S at time step 2 620. Test 1: H2S at time step 5 9 21. Test 1: H2S at time step 1022. Test 1: H2S at time step 15 70 23. Test 1: H2S at time step 2024. Test 1: H2S at time step 30 1 25. Test 1: H2S at time step 40 726. Test 1: H2S at time step 50 2 27. Test 2: S at time step 2028. Test 2: S at time step 30 3 29. Test 2: S at time step 40 730. Test 2: H at time step 20 4 31. Test 2: H at time step 30 5 32. Test 2: H at time step 4033. Test 2: H2S at time step 20 76 34. Test 2: H2S at time step 30 7 35. Test 2: H2S at time step 40V Acknowledgement The completion of this thesis has depended on the good advice, given over many hours, from Jim Varah, the multitude of suggestions from Paul Leblond, and the original idea of Bill Emery; and from all, a great deal of patience. 1 I. INTRODUCTION Studies of ocean currents in coastal regions have shown that fresh water run-off is an important factor in large-scale circulation (see [6] for several examples). While these studies consider the effect of large rivers on the coastal ocean, the influence of run-off over a large coastal region has not been closely examined. Dr. Paul LeBlond and Dr. William Emery of the Department of Oceanography at the University of British Columbia have undertaken a study of river run-off and salinity distribution off the coast of southern British Columbia and northern Washington State in an effort to understand their effects on alongshore currents. One aspect of the study has been the development of a mathematical model to describe those coastal currents created by river run-off. The model consists of a set of equations that describe nearshore hydrodynamics in terms of changes in salinity, and thus can be used to measure the impact of coastal run-off. The main purpose of this thesis is to determine a numerical method to solve the equations defining the model, and to develop a computer program implementing the method. As well as being of interest from a physical oceanographic point of view, it is hoped that the results from the model may be used to gain an understanding of offshore fisheries. In particular, it has been observed that in most years a large percentage of Fraser River salmon migrate via the Strait of Juan de Fuca, while the remainder return through Johnstone Strait at the northern end of Vancouver Island. In some years, however, a 2 majority return through Johnstone Strait. The reasons for this are not known, but it is likely that the large-scale coastal circulation plays a large part in which route is chosen. An understanding of these mechanisms may lead to an explanation for this behavior. Chapter One contains a discussion of the derivation of the problem, leading to the equations to be solved. Definitions of the variables and parameters follow, and the boundary and initial conditions are specified. A description of how the model is to be used concludes the chapter. In Chapter Two, the notation to be used to describe finite difference schemes is introduced. The concepts of consistency, accuracy and stability, as they relate to application of the Crank-Nicolson differencing of a linearized form of the equations, are discussed. The boundary conditions are then specified in terms of finite differences and their effects on stability examined. The nonlinear terms in the model equations are discussed and methods of dealing with them presented. In Chapter Three, we examine the systems of linear equations produced by the Crank-Nicolson scheme, and discuss methods to solve them. More efficient difference schemes for the equations are the alternating direction implicit, or ADI, methods. They are described and the concepts of Chapter Two are extended to one of these, the Peaceman-Rachford ADI scheme, for a linear equation. Methods of treating the additional problems presented by the model equations, including lower order space derivatives, variable coefficients and nonlinear terms, are 3 examined, and tested by solving sample problems with known solutions. Results from computer runs are shown and the accuracy and stability of the solutions discussed. Finally, we discretize the full equations in terms of the Peaceman-Rachford scheme. In Chapter Four, the ADI method is applied to the model equations with simplified boundaries to ensure the difference scheme satisfies basic physical properties. The discontinuous nature of the initial conditions is discussed, and we show how this problem is dealt with. The land boundary condition for the thickness layer is examined and the problems and methods we use to handle them are explained. Next the formula required by the ADI method for calculating the salinity defect boundary conditions is presented. Sample runs and their results for the full model equations are given, with brief physical interpretations. Recommendations for proceeding with the model conclude the chapter. 4 II. CHAPTER 1 2.1 Derivation Of The Model The primary concern of this thesis is the set of numerical problems presented by the mathematical model of coastal circulation. The derivation that follows is brief and includes little justification, from a physical oceanography viewpoint, for the assumptions and accomodations made in formulating the final equations. A complete derivation and discussion of the implications of the model is available in [6]. Our derivation is a condensation of that information. The model considers the nearshore ocean as having two layers; a lower layer of constant density, and hence salinity, and an upper layer of varying density and thickness. The motion of this upper layer is described by the two-dimensional momentum and conservation equations. If we consider a time scale greater than five days, local accelerations can be ignored, and the motion described is primarily geostrophic. In addition, we want only waves with a Rossby number much less than one; that is, U0f/Ax << 1, where U0 is the speed scale and f the Coriolis force. The space scale, Ax, must therefore be much greater than 1 km. We assume at this point that a value of Ax > 10 km. is suf f ic ient. Under these assumptions, the momentum equations for the upper layer are (1) -fv = (-1/p,)9p,/9x - Fx (2) fu = (-1/p,)9p,/9y ~ Fy 5 where: (u,v) = (u(x,y,t),v(x,y,t)) = components of velocity; p-, = pressure in the upper layer; p, = density of the upper layer; (Fx,Fy) = components of friction. Since the lower layer is assumed to be at rest, and using the assumption that the surface height is much less than the thickness of the upper layer and can therefore be neglected, (1) and (2) can be integrated over the upper layer to yield (3) -fV = -g/2 (H2Ap/p,)x - Fx (4) fU = -g/2 (H2Ap/p,)y - Fy where, (U,V) = (U(x,y,t),V(x,y,t)) = U = components of transport, where transport is velocity integrated over the upper layer; H = H(x,y,t) = the thickness of the upper layer; p2 = density of the lower layer; Ap = p2 - p, = density difference between layers; g = gravity; In the equations above, we have adopted the usual notation for partial derivatives; they are indicated here and throughout this thesis by 3 /3x (for example) as in equation (1) or by subscripts of the independent variable. In terms of transports, the volume and salinity conservation equations are (5) + VH-U = We (6) (HSU)+ + V^US = WeSb 6 where: SM = Sfc(x,Y/t), salinity in the upper layer; S^ = a constant lower layer salinity; We = the entrainment velocity; VK = the divergence operator. Density can be considered a function of salinity, with an equation of state given approximately by (7) Ap = p,a(T) S* where a is a function of temperature, T, and is assumed constant for the range of temperatures encountered. Introducing S(x,y,t) = Sfe(x,y,t) - SIA(x,y,t) , the difference in salinity between lower and upper layers (so S(x,y,t) > 0), and using (7) to relate density to salinity, equations (3) through (6) can be rewritten as (8) -fV = -(ga/2p,) (H2S)X - kU (9) fU = -(ga/2p,) (H2S)y - kv (10) Ht + VK.U= W€ (11) (HS)^ + VK-US = 0 Here we have taken friction components (Fx,Fy) to be proportional to the components of velocity (u,v), with k the friction coefficient, and integrated over the upper layer to produce transports. Multiplying (8) by k and (9) by f and combining gives (12) (f2 + k2) U = -(ga/2Pl) ( f(H2S)y + k(H2S)x ) (13) (f2 + k2 ) V = -(ga/2Pl) ( -f(H2S);< + k(H2S)y ) Substituting (12) and (13) into (10) yields, with C = gak/(2p,(f2+k2)) 7 (14) Ht = Wt + CV2( H2S ) . Similarly, substituting (12) and (13) into (11) gives (HS\ = C (kSV2(H2S) + fJ(S,H2S) + kVS-V(H2S)) or (15) St = C/H (kSV2(H2S) + fJ(S,H2S) + kVS V(H2S)) - (S/H)Ht where V2 is the Laplacian operator, J is the Jacobian operator in terms of x and y and V is the gradient operator. The function H2S represents the distribution of the density mass field; we expect fresh water run-off to create gradients in H2S, causing currents in the model domain. Both (14) and (15) are nonlinear, parabolic partial differential equations and are the governing equations for our model of the coastal circulation. Solving (14) yields H(x,y,t), the thickness of the upper layer, while solving (15) gives S(x,y,t), the salinity defect, the difference between S^, and Sb. The parameters and constants in (14) and (15) are defined and assigned values as follows: g = 9.8 m/sec2 f = 1.1E-4 sec"1 k = 0.25E-4 sec"1 a =0.77 ppt"1 p, = 103 ppt We = 6.0E-5 F2v where F2 = Froude number = (U2 + V2)/(gaSH3)/p, v = surface velocity = (U2 + V2 )''1*/H 8 Further discussion of parameters k, a, W and justification for the values chosen can be found in [6]. The components of transport, (U,V), and hence the velocity v, can be found algebraically from (12) and (13). 2.2 The Domain Of The Model The domain of the equations is an area bounded by open ocean to the north, west and south, and on the east by the coast from the Olympic Peninsula to the Queen Charlotte Islands (see Figure 1). The coast is broken between the Olympic Peninsula and Vancouver Island by the Strait of Juan de Fuca and north of Vancouver Island by Queen Charlotte Sound. It is through these two passages that river run-off enters coastal waters. For the purposes of a numerical model, the domain has been idealized to a square, with x designating the east-west axis and y the north-south axis, as illustrated by Figure 2. The origin for this domain is the southwest corner of the square. The land boundary is therefore taken to be parallel to the y-axis. 2 . 3 Boundary Conditions The equation for the salinity defect, (12), requires that the gradient on all boundaries be zero; that is, (16) 3S/3n = 0 where 9 /3n is the normal derivative. This is the usual form of the no-flow condition across a boundary. Where the land boundary is broken, representing fresh water flow into the domain, S is prescribed from historical data. 9 The boundary conditions for H are a little more difficult. It is assumed that near the open ocean boundaries, H does not change rapidly and can therefore be approximated from those values within the domain that are near the boundary. The actual method for finding H on these boundaries will be presented when the numerical scheme has been explained. On the land boundary, the x-direction transport must be zero, as there can be no flow across land, and therefore, from equation (12), f(H2S)y + k(H2S)x = 0 Expanding the derivative terms yields 2fSHHy + fH2Sy + 2kSHH;t + kH2S^ = 0 or (17) Hy = -(fHSy + 2kSHx + kHSx)/2fS To find H on the land boundary, we note that HA can be approximated numerically and Sy is known from the solution of (15). Thus, (17) is a first order initial value ordinary differential equation of H in y, with the initial value for the equation being the value of H given at the river discharge regions of the boundary. 2 . 4 Initial Conditions The solution of (14) and (15) requires initial conditions H(x,y,0) and S(x,y,0). Since the model is intended for use with historical river run-off data, some historic average of thickness and salinity data will eventually be considered as initial values. For the purposes of our initial computer runs, constant field values will be used for both H(x,y,0) and 1 0 S(x,y,0). 2.5 How The Model Will Be Used Equations (14) and (15), with boundary conditions as specified, describe the effect of river run-off on coastal circulation. The purpose of the model is to try to closely approximate actual currents and salinities by using collected data in boundary conditions. To that end, river discharge data has been collected for the coasts of Washington State, B.C. and Alaska. From these values, surface velocities (and hence transports) at the inflow regions were calculated and • a weekly time series produced to be used as boundary values for the model equations. Similarly, daily surface salinity readings, recorded at lighthouses along the coast, have been smoothed to a weekly series for use as salinity boundary values. Upper layer thickness values have been approximated from transports and surface salinities. A period of approximately twenty years of actual or interpolated weekly transport and surface salinity records are available to be used as boundary conditions for the model. The aim of the study is to approximate the coastal circulation for that length of time, and to compare the results with the known salmon migration data. From this it is hoped that some inferences can be drawn about the relationship between the two. 11 Figure 1 - Domain of the Coastal Circulation Model Figure 2 - Domain of the Numerical Model 1 3 III. CHAPTER 2 3.1 Notat ion To solve the partial differential equations (1) Ht = CV2(H2S) + We (2) St = C/H (kSV2(H2S)+fJ(S,H2S)+kVS V(H2S)) - (S/H)H_,. numerically, the domain is discretized by overlaying a rectangular grid. Initially, this domain will be a simplified version of that shown in Figure 2 by considering it to have just a single discharge point. The model equations are valid over the region represented by this grid, which has grid spacings of Ax and Ay in the x and y directions. Coordinates on the grid are specified by (iAx, jAy), with i=1,...,L and j=1,...,M, and LAx = MAy = 600 km. The time direction is also discretized by dividing the time domain in N equal time steps of size At. A specific time is represented by nAt. The solutions of (1) and (2) will be approximated at any grid point by the discrete solutions H(iAx,jAy,nAt) and S(iAx,jAy,nAt), respectively. For ease of notation, the approximate solutions are simplified to H-- and To specify each of (1) and (2) in discrete terms, the partial derivatives must be approximated in some way. Finite element or finite difference methods are used most often for finding numerical solutions to partial differential equations. There are a number of finite element methods for approximating (1) and (2), but they are for the most part difficult to program and in general are computationally expensive to solve. Finite 1 4 differences, in addition to being easier to program and less costly to solve, approximate derivatives using divided differences which are more suited to our rectangular domain. 3 . 2 Finite Differences We will use centred differences for the first and second order space derivatives and the first order time derivative. They will be defined using linear operators D^ and Djx for the first and second order x-direction derivatives, and Dey and Djy for the y-direction derivatives. For some function u, u=u(x,y,t), let D„X uij = (u„t|. - u-„(. )/2Ax and D*x ujj = ( ultt)- - 2UJJ + u;_(i. )/(Ax)2 From Taylor series expansion of u((i+1)Ax,jAy,t) and u((i-1)Ax,jAy,t), it is easily seen that uK = Dox u + 0((Ax)2) and uxx; = D„2X u + 0( (Ax) 2) where 0((Ax)2) is the order of the remainder term. The D^ and D2X linear operators are said to be second order accurate approximations. Similarly, Doyf and D^ are linear operators for the y-direction derivatives and are accurate to 0((Ay)2). Using the linear operators D^ , D|x , D^ and Djy , finite difference methods can be constructed to approximate the derivatives in (1) and (2). The solutions H and S are found by marching ahead in time in some fashion; the solutions at t=(n+l)At depend on values at lower time levels. 1 5 If H and S are found individually, using values from previous time levels, the method is said to be explicit. While computationally very simple, explicit methods are restricted to some maximum size of time step used to advance in time. Using a time step greater than this maximum, which depends on Ax and Ay, results in solutions that grow exponentially. Avoiding such instabilities by using a time step less than the maximum severely restricts the choice of Ax and Ay. As was mentioned in Chapter 1, various assumptions used in deriving the model equations put constraints on the size of Ax and Ay. An explicit method is therefore not appropriate for finding solutions to our equat ions. A finite difference method is implicit if the values of the dependent variable at points at the next time level are interrelated. Implicit methods thus require solution of a system of linear equations at each time level. While naturally taking more computing time than explicit methods, implicit methods generally do not place restrictions on the time step, beyond that required for accuracy. A well-known implicit method for solving parabolic equations is the Crank-Nicolson method. In terms of the linear operators defined above, the simple two-dimensional diffusion equat ion ut = UKX + uyy is approximated by the Crank-Nicolson (or CN) method at (iAx,jAy,(n+l/2)At) by - u|j = At(D£ ( + u{] )/2 + D/y ( u?j' + ujj )/2) 1 6 The difference scheme is centred about t=n+l/2 to preserve the second order accuracy in time, and is a two-level difference equation, solving for t=n+1 from values at t=n. The discussion in the remainder of the chapter which concerns finite difference equations and their application to (1) and (2) will be limited to the CN method. Before applying the CN scheme to equations (1) and (2), we must decide how the spatial derivatives are to be expanded. The equations have derivatives of products of the dependent variables which must be expressed as difference equations. Consider as an example the second order x-direction derivative , (H2S)^ , in equation (1). Applying the product rule in the usual way yields (3) 2HSHK)C + 2H2SX+ 4HHXSA + H2S^ while considering (H2S) as a product of functions (HS) and H gives (4) (HS)HX;< + 2(HS)XH^ + H(HS)XX Upon discretization, these expressions are clearly unequal, except for fortuitous values of H and S. Obviously, the solutions of both (1) and (2) will depend on which way the (H2S) derivatives are expanded. We are solving equation (1) for H, and it seems reasonable that those functions multiplying H should be considered as single, separate functions. That the derivatives should be expanded as (4) above was confirmed by experimental computer runs (as explained in Chapter 4). Expanding equation (1) therefore yields 1 7 (5) Ht = C((HS)V2H +2(HS);<H;t + 2(HS)/Hy +HV2(HS)) + We . The HV2(HS) term may be treated as the other terms in (5); that is, V2(HS) is considered the coefficient of H. It is usual ([12], p. 195) to assume that inhomogeneous terms involving the dependent variable are known. (5) becomes (6) Rt = C((HS)V2H + 2((HS)XHK + (HS)yHy)) + f, where f, = HV2(HS) +We. The equation for S also contains derivatives of (H2S), but since we are solving for S, their expansion must be different than that outlined for H. Treating H2 as a single entity was found to work best, and thus, for (H2^)^, the appropriate expansion is H2SXX+ S(H2)X;< + 2SX(H2)(X Expanding equation (2) gives (7) = C/H (kS(H2V2S + 2VH2-VS + SV2H2) + fSJ(S,H2) + k (SVS-V(H2)+H2(VS)2) - (S/H)Ht Just as for H, terms not involving derivatives of S are treated as forcing functions, and (7) becomes (8) St = C/H (kS(H2V2S + 2VH2* VS) + fSJ(S,H2) + k(SVS«V(H2)+H2(VS)2) + f2 where f2 = (CkS/H) (SV2H2) -(S/H)H<;. While the method of expanding the derivatives has been determined, the coefficients of the derivatives that result are obviously functions of the dependent variables; both equations are nonlinear. Discretization of such equations produces nonlinear systems of equations which would require some kind of iterative scheme for solution, such as Newton's method. Since 18 iterative methods for nonlinear equations are in general computationally very expensive, we approximate further by linearizing the equations. Although we leave the computational details until later in Chapter 2, the result of linearizing (6) and (8) are equations having variable coefficients. The equations can be simplified to and s-t = bisxx + b2Syy + b3Sx + b^Sy + f2 where aK , b^ , i = 1,..,4, are functions of x, y and t. In terms of the linear operators defined earlier, the CN discretization i s (9) Hf* - H-5 = At(a^H^+ a2D2yH?;Y-+ a3D^H^ + aaDoyH^yH f,) and (10) S?f ~ S% = At(b,DiS^+ b2DJySi^+ b3ELxS(^ + b«DoyS^+ f2). It should be noted here that evaluating the coefficients of each derivative of the equations shows that, for representative values of H and S, and using the constants defined in Chapter 1, both equations are strongly diffusive. This is important in determining the difference scheme because it is well known ([12], p. 26) that central differences are poor approximations to convective terms if their coefficients are large relative to those of the diffusion terms. While the convective terms are not insignificant, the primary moving force is produced by the 19 diffusion terms. 3.3 Consistency, Stability And Convergence Finite difference equations are naturally only approximations to differential equations. It is necessary to show that the finite difference solution is a 'good' approximation to the exact solution. Briefly, this means that we require the difference equation to tend to the differential equation for Ax,At -> 0; that the solution of the finite difference scheme will give a bounded solution (assuming the same for the exact solution); and that the finite difference solution will approach the exact solution as Ax,At -> 0. These properties are called consistency, stability and convergence. The standard theory for stability and convergence of initial value problems generally assumes constant coefficients, and does not consider the effect of either boundary conditions or variable coefficients. We will assume for the moment that the model problems meet these criteria, and consider later in the chapter the effect on our analysis of variable coefficients and boundary conditions. To simplify the presentation of these concepts for the CN scheme, we consider the one dimensional version of the model equations. What follows is easily extended to two dimensions. For u = u(x,t), let (11) ut = aux^ + bu„. + f(x,t) with a and b constant, and with initial condition u(x,0) = g(x). 20 Much of the analysis that follows is based on the notation and definitions of Kriess and Oliger ([5], p. 29). Applying the CN scheme to (11) yields (1 - aAtD^ - bAtDw)v(x,(n+1)At) = (1 + aAtD^ + bAtDC!X)v(xfnAt) + Atf(x,nAt) with initial condition v(x,0) = h(x), where v(x,t) is the solution to the finite difference equation. For any x, the discretization of (11) produces linear combinations of v(x,(n+l)At) and v(x,At). We can simplify these using matrix notation; (11) can be expressed as (12) Q,v(x,(n+1)At) = Q0v(x,nAt) + Atf(x,nAt) Definition The difference scheme is accurate of order (q(,q2) for the particular solution u(x,t), if there is a constant c and a function C(t), bounded on every finite interval [0,NAt], such that for all sufficiently small At, Ax | |Q,u(x,t+At) - Q0u(x,t) - Atf | | < C(t)(AxV + Atl* ) and | |h(x) - g (x) | | < c(Axtl + At?*" ) where ||» || is the 12 norm. We have seen that the linear operators, D,* and Dj^ , approximate 9/9x and 92/9x2 with accuracy 0((Ax)2), and that the CN approximation of 9/9t is accurate to 0((At)2). Assuming that h(x)=g(x), that the initial condition is known exactly, the scheme is accurate of order (2,2). The error in the scheme, caused by truncating the Taylor series in our approximations to the derivatives, is called the truncation error. For the CN scheme, the truncation error is At 21 times the accuracy, or 0(At3 +AtAx2). Note that this is the local error, the error in the solution from one time step to the next. If the difference equation converges to the differential equation -being approximated, the method is said to be consistent. For the CN scheme, consistency follows directly from the definition of accuracy. Since C(t) is assumed to be bounded, and qi,q2 > 0, C(t) ( (Ax)**' + (At)*1" ) -> 0 as At,Ax -> 0 and the scheme is consistent. While the difference solution approximates the exact solution with the order of accuracy given above, the exact error depends on the size of the partial derivatives multiplying the (Ax)2 and (At)2 terms in the Taylor series expansion. In practice, these of course are unknown. To minimize the error, several steps are taken when actually finding a numerical solution. First, the size of the time step is chosen to be about the same size as the space step; thus (At)3 will be roughly equal to At(Ax)2. Second, it is conventional to nondimensionalize the equations. This requires the equations to be scaled so that 0 < x,y < 1. The time scale now depends on the spatial resolution used. Stability requires that the solution to the difference equation, v(x,t), at t=nAt, is bounded as At,Ax -> 0. Def inition The difference approximation (11) is stable for a sequence At, Ax -> 0 if there are constants as, Kj such that for t ^ 0, the estimate 22 ||v(x,t)|| = | |S(t,0)v(x,0) | | < Ksexp(as t) | | v(x,0) | | holds. Here S(t,0) is a solution operator with S(0,0)=I, S(t2,0)=S(t2,t1)S(tl ,0); t2>t,>0. For t=nAt, (12) is (13) Q,v(x,(n+1)At) = Q0v(x,nAt) + Af(x,nAt) and S(nAt,0) can easily be found. For Q = Q7'Q0, v(x,(n+1)At) = Qv(x,nAt) + Qf1 Atf(x,nAt) It is evident that the right side of the stability condition in the definition above should include the contribution from the forcing function, f. That is, the bound should include upper limits on At^S (nAt, iAt) f (x , i At) . However, this only increases the bound, and for simplicity, we let the definition stand as is. It is clear that we have to show that | |Q| I**-1 is bounded for stability to hold. Assume first that a Fourier transform of v(x,t) exists. For every frequency, w, we must show that the amplification factor is not greater than one. That is, stability requires that from one time step to the next the amplitude of every frequency does not increase, except to allow truly increasing solutions. This is the von Neumann necessary condition for = Q(Qv(x,(n-1)At)+Q f1 Atf(x,(n-1)At) + Q ~1 Atf (x , nAt) 23 stability ([12], p. 70). Applying Fourier transforms to (13) requires only the transform of the linear operators D0)l and De2. As shown in ([9], p. 37), they are DPX = i sin(/3/Ax) (14) D.2, = -4 sin2(/3/2)/(Ax)2, and G, the amplification factor, is (from [8], p. 196), 1 - 2asin20 + i (b/2) (aAt/a) sin(2/3) + Atcf G = 1 + 2asin20 + i(b/2)(aAt/a) sin(20) with a = aAt/(Ax)2 and j3 = kAx. Simplifying G leads to an equation with (1 + 2asin2/3)2 + (b2/4)aAt/a sin2(20) in the denominator; it is clearly greater than 1 for any /3. So G < (1 - 2asin2/3 + i (b/2)(aAt/a) sin (20) + Atcf) (1 + 2asin2/3 - ib/2 ( aAt/a ) sin2/3 ) In modulus, the amplification factor can be given by (15) |G| < (l+h1(/3)At+h2d3)(At)2+h3(/3)(At)3+h«(0)(At)M = 1 + O(At), where h,(|3), h2(/3), h3(|3) and ha(j3) are bounded functions. The O(At) term allows exponential growth that may exist in the true solution because of the forcing function. The inequality above satisfies the von Neumann condition, which is necessary for stability; it is known ([12], p. 72) that for all two-level difference schemes for scalar equations, the von Neumann condition is sufficient as well. This shows that the solution from one time level to the 24 next is bounded. ||S(t,0)v(x,0)|| must therefore also be bounded, and the CN scheme meets the requirements of the stability definition. It is clear that |G| is bounded for any values of of At and Ax, and the scheme is said to be unconditionally stable. If the difference between the solutions to the differential and difference equations tends to zero as At, Ax -> 0, a difference scheme is said to be convergent. From [5], p. 31, we have the following: Theorem Assume that for fixed At and Ax the definitions for accuracy and stability are valid. Then for all t = nAt, n > 1, ||u(x,t)-v(x,t)|| < exp(a5 t) (Ax*' +At*1 )c . Since it was shown, for the CN scheme, that (q,,q2)=(2,2) and that stability is unconditional, there is no restriction on the way At, Ax -> 0. Therefore the right side tends to zero as At, Ax -> 0, and convergence is proven. 3.4 Variable Coefficients As was mentioned, the analysis above assumes a differential equation with constant coefficients, a condition the model equations clearly do not meet. While the accuracy is unaffected by variable coefficients, assuming of course that there are no errors in their calculation, it is quite possible that the stability condition may change. One method of ensuring stability is to look at the problem in terms of 'local stability'. That is, the coefficients are evaluated at each grid point in the domain and tested separately 25 for stability. While this does not guarantee convergence, it is usual that instablities in the approximate solution start locally ([12], p. 91). However, to apply von Neumann's method to a large number of points presents obvious difficulties. One practical way to test for stability of variable coefficient problems is to perform computer runs of sample problems. This was carried out for several problems, and the results are given in Chapter 3. 3.5 Stability Of Boundary Conditions To determine the effect of the boundary conditions on stability, we must examine the system of linear equations produced by the CN scheme for the initial value problem, and check that the introduction of boundary conditions do not add to the matrix Q any eigenvalues z with |z| > 1. The strategy is to assume that such an eigenvalue z is introduced by the boundary condition, and show that this leads to a contradiction. In particular, we want to show that this occurs for the boundary conditions of the model equations. As specified in Chapter 1, the boundary conditions for H are prescribed on all boundaries. The boundary conditions for S are prescribed at the inflow regions, and the no-flow condition elsewhere. Assuming smoothly changing functions of H and S (and thus no boundary layer problems), prescribed boundary conditions do not change the eigenvalues of the corresponding initial value problem, and so only the no-flow condition will be analysed. We consider the sample equation (12), with no forcing 26 function, and 0 ^ x < 1. It is known that for parabolic problems, stability of the left boundary condition implies stability for both boundaries [18], and so only the effect of the value at u(0,t) will be considered. In terms of u(x,t), the second order accurate approximation to the no-flow condition at the left boundary is u(0,nAt) = (4/3)u(Ax,nAt) - (1/3)u(2Ax,nAt) We have seen that the amplification matrix, Q = QpQor must have eigenvalues not greater than one for v(x,(n+l)At) to be a bounded solution. Let g = (g*,...,g«) be an eigenvector of Q, so that Qg = zg or zQi9 = Qo9 with z an eigenvalue. Q, and Q0 can be expressed as sums of their coefficients, so in terms of g and z, z TL c£ g; = X.a; g; or ( 16) jL(zCj - at )gj_ = 0 . From [18], (16) is a linear difference equation with constant coefficients, and the solution takes the form, 9 = £aj(T.(z)) where the Tj (z) are the roots of (16). .Substituting this into (16) gives (17) T2 - 2/3T + T/TJ = 0 where B= 1/(1+Ax)+((z-1)/(z+1))/(X(1+Ax)) r = X/2 - XAx/2 T? = X/2 + XAx/2 27 and X = At/(Ax)2 For |z| > 1, (17) can be shown to have roots T, and T2 with |T,(z)| < 1, |T2(z)| > 1. So (18) qi = a,T, (z) and, to be an eigenvector, gj must also satisfy the no-flow boundary condition. Thus, we require that ge = (4/3)g, - d/3)g2 . Substituting this expression into (18) gives (19) a,( 1 - (4/3)T,(z) + (l'/3)T^z) ) = 0 The boundary condition will produce an unstable solution if, for z, with |z|>1, an eigenvalue, the corresponding root of (17) is also a root of (19). The roots of (19), however, are 1 and 3, contradicting our finding that one of the roots of (17) is less than one. Thus z, with |z|>1, is not an eigenvalue; the no-flow boundary, approximated to second order accuracy, does not affect stability. 3.6 Nonlinear Terms To avoid having to solve our nonlinear equations directly, it is necessary to linearize the coefficients of the derivative terms. Linear equations are, of course, much easier to solve, but the linearization procedure yields only an approximation to the true coefficients. To examine this effect on the accuracy of the system, we note first the equations are quasilinear, since the coefficients of the second derivative terms involve only first derivatives or the function itself. Although the CN difference scheme can be 28 modified to produce predictor-corrector formulae [8], the truncation error estimate is 0(AtAx2+A\?fx-) rather than 0(AtAx2+At3) which the CN scheme (for linear equations) yields. Lees [8] has shown that by estimating the quasilinear coefficients using extrapolated values from the two previous time steps, the trucation error remains 0(AtAx2+At3). Since the CN scheme is centred about t=n+l/2, the coefficients of each term must be calculated at that time. In particular, if the HS terms in the coefficients of the second order partial derivatives of (9) are approximated by US"*1'7- = 1.5 HSw - 0.5 HS*",' an error of 0((At)3) is added to to the truncation error, which thus remains 0(AtAx2+At3). We see that the coefficients of 3H/3x and 3H/3y, which involve 3(HS)/3x and 3(HS)/3y terms, can be approximated in the same way. For (HS) at point (i,j), the derivative at time n+1/2 can be approximated by (HS)X = ( (HS)j^ - (HS)^- ) / (2Ax) and (HS)y = ( (HS)c;-„ - (HS)li>( ) / (2Ay) where all the (HS) values are calculated at t=n+l/2 as specified above. This method can also be applied to the S% and Sy2 terms that appear in the equation for S. One of each of the derivatives can be treated as the coefficient of the other, and approximated by extrapolating to time n+1/2 ([2], p. 141). Equation (10) has quasilinear coefficients which are functions of SH2, H2, Hx, Hy, Sx and Sy. All of these are estimated in the same fashion as outlined above, without altering the order 29 of the truncation error. The effect on stability of applying this extrapolation procedure to all the terms in each equation is not known. Varah [19], however, has shown that for an equation involving nonlinear terms in the first spatial derivative of the function, and the function itself, the stability of the CN scheme with extrapolation depends on the relative size of the coefficients of the diffusive and convective terms. For equations which are not highly convective, the CN scheme was found to be conditionally stable. Since both model equations are highly diffusive, it seems likely that stability will be maintained. While we can provide no stability analysis, the results of preliminary computer runs (described fully in Chapter 4) show that, for chosen space and time steps, this linearization produces results which appear to be conditionally stable. 3 . 7 Solving The Equations To solve the model equations, we note first that they are coupled, in that the two independent variables are present in the two equations. While the Crank-Nicolson difference scheme could be applied to solve for H and S simultaneously, the only new terms to be solved directly rather than approximately are the forcing functions. Since they are nonlinear as well, one of either H or S would have to be approximated in any case. Little would be gained from this approach. Instead, assume, as was mentioned in Chapter 1, that S varies slowly from one time level to the next. If equation (9) 30 is solved first, HS has to be estimated (by extrapolation). If S is in fact slowly varying, the error of approximation will be small (at least smaller than that caused by estimating H2, as solving equation (10) would require). As we have seen, the addition to the error is 0((At)3), which, will only cause a small increase in the overall error of approximating H. Once H has been calculated, the average of H^'and H* , H,wVa-, can used in evaluating the coefficients of equation (10), rather than an extrapolated value. Solving for SA+' , it is possible to repeat the process, using the new values of H*+v and then S**' to calculate new approximations for HM"''t and Sn+Vl-. Each step may be considered an iteration of a predictor-corrector method. We can lessen the effect of the interaction of H and S by alternating the order of solution from one iteration to the next. In practice, it is possible to continue iterating until the solutions from one iteration to the next change by only some arbitrary (small) amount. For smoothly changing values of H and S, the solutions will converge, it is hoped, in only a couple of iterations. Applying the CN difference scheme to the model equations produces two systems of linear equations that must be solved at each time step. If some method of iteration is used, which will change the coefficients of the matrices, each iteration will require the systems to be solved again. It seems likely that the majority of computation time will be spent in solving linear equations, particularly as we increase the resolution to obtain better accuracy. The consequences are discussed in the next 31 chapter. 32 IV. CHAPTER 3 We now consider the structure of the matrices we must solve, the possible methods of their solution and the computational effort, in terms of operation counts, each would need. The structure of the matrices for equations (9) and (10) of Chapter 2 depend on the linear operators and the way the H^ and S??' are ordered in the vector of unknowns. The linear operators involve only grid points one grid spacing away in each spatial direction; that is, to solve for (i, j) involves (i + 1,j), (i-1,j), (i,j + l), (i,j-l) as well as ( i , j ) . The conventional approach is to order the grid points by row with the resulting matrix being banded with half bandwidth L+1, where L is the number of grid points in the x-direction. For our equations, both matrices are nonsymmetric. A banded, nonsymmetric matrix can be solved using either a variant of Gaussian elimination or any number of iterative methods. The former has the advantage of reliability, but the latter can be considerably faster. The rate at which iterative methods converge, and thus the number of iterations they require, depends on the condition number of the matrix. For an ill-conditioned matrix iterative methods will often converge slowly, if at all. Noting that the coefficients of our matrix are the result of approximations to nonlinear terms, and vary with x, y and t, there is little we can say about how well-conditioned the matrices will be. (We assume of course-that they remain nonsingular.) This being the case, Gaussian elimination would appear to be the best method. 33 To calculate the amount of computational effort required, we need only look at the order and the half-bandwidth of the matrices. If L is the number of grid points in each direction (so Ax=Ay), the order is L2 and the half-bandwidth L+1. Although the matrices initially have zeros between the elements on the edge of the band and the diagonals, Gaussian elimination "fills in" these zeros. When calculating the amount of work, these elements must be considered as nonzeros. For such a matrix, having order L2, the operation count (the number of multiplications) necessary to reduce it to triangular form is 0(L"). However, difference schemes exist for equations with the same form of the model equations which produce matrices which can be solved in only 0(L2) operations. They are called alternating direction implicit methods, and are the subject of the next section. 4. 1 Introduction To ADI Alternating direction implicit, or ADI, methods are perturbations of difference methods formulated directly from the difference approximations to derivatives. ADI methods most often arise by adding difference operators to a multi-dimensional difference scheme in a way that does not affect the accuracy of the scheme but allows an easier solution. To illustrate the derivation of a particular ADI method, the CN scheme for ut = u^ + uyy , is, using the linear operators introduced in Chapter 2, (1) (1 - XD02x - XD^)u*" = (1 + XD2X + XD2y )u* 34 where X=At/2. If the operator X2D^ De2y is added to the difference operator on both sides, (1) becomes (2) (1 - XD^ - XDiy + X2D<2C D2„ )u^' = (1 + XD|X + XDe2y + X2D«2X D2y )u* . Factoring (2) in terms of D£x , D|y gives (3) (1 - XD2X )(1 - XD2y )u**> = (1 + XD,2, )(1 + XD2, )u* . If we introduce an intermediate time level between n and n+1, say n+*, (3) can be 'split' into two equations, namely (1 - XD2 )u*** = (1 + XD2 )un (4) V ( 1 - XDC2X ) u**» = (1 + XD2y ) u*** . The variable u*** is a dummy variable in the sense that the value of U"1** is not an approximation to the solution at any particular time. The method of splitting used in (4) is called the Peaceman-Rachford formula ([9], p. 60). Although there are many methods for splitting a two-dimensional parabolic equation like (4), this formula will hereafter be referred to as the ADI scheme. The advantage of (4) is that the solution at the next time level (either u**' or un**') , is found along lines parallel to either the x- or y-axis. The matrix to be solved thus only involves points immediately adjacent to the point being found; that is, the matrix is tridiagonal. The computational effort required to solve a tridiagonal system of equations is small compared to that for the CN scheme. For an L2 by L2 system, the tridiagonal matrix requires only 0(L2) multiplications, compared to 0(1/) required by the CN scheme. This considerable savings in computation time will not come 35 at the expense of accuracy or stability of the difference method. As we will show, the Peaceman-Rachford scheme, like the CN scheme, has a truncation error of 0((At)3+At(Ax)2) and is unconditionally stable, for a sample problem similar to that analyzed in Chapter 2. Note that in calculating the truncation error, we have assumed that Ax=Ay. For the rest of the thesis, this will be the case. 4.2 Accuracy And Stability To calculate the accuracy of the ADI scheme, we first consider an equation of the same form as our model equations, namely, for u=u(x,y,t), (5) ut = auXx + bUyy + cux + duv , where a,b,c and d are constant coefficients. Factoring (5) using the ADI scheme in the same manner as above, gives (6) (1 - XaD2x - XcD„x )(1 - XbDjy - XdDey )u^1 = (1 + XaD2x + XcD#Jf )(1 + XbD2y + XdD0j ) u* , where X is as above. Splitting (6) we get (1 - XbD2 - XdD*v )uM+* = (1 + XaD2 + XcD^ ) u* (7) * * (1 - XaD^ - XcDe;t )u*«" = (1 + XbDe2y + XdDcy )u*+* To calculate the error, we see first that the operators added to the CN scheme for (5) are (\2cd)Dox Dey , (X2bc)DoxD<2y , (X2ad)Doy D2^ and (X2ab)D^Dj^ . These terms do not result from Taylor series expansion of the derivatives, and the errors they produce can be added directly to the truncation error found for the CN scheme. Since each operator is applied to both uH and u"**' , the errors can easily be found by expanding u1**1 in a 36 Taylor series. Examining just one of these terms, (X2cd)D D , we have (X2cd)D D (u + Atu +0((At)2) = (X2cd)D D u Cancelling u on both sides, we see that the addition to the truncation error of this term is 0((At)3). The contribution to the error from the other operators is similarly found to be 0((At)3), so the truncation error for ADI is of the same order as the CN scheme discussed earlier; the added terms leave the overall accuracy unchanged. Consistency follows in the same way as shown for CN in the last chapter. Stability for the ADI scheme can be seen by again comparing with the CN scheme. Looking at the amplification factor, G, for CN applied to equation (5), we have 1-2(a,S,+a2S2)+i((b/2)(a,At/a)S3+(d/2)(a2At/d)S„) (8) G = 1+2(a1S1+a2S2)+i((b/2)(a,At/a)S3+(d/2)U2At/d)Sa) where S,=sin2/3, S2 = sin27j, S3 = sin2/3, Sa=sin2ry and a, = aAt/(Ax)2, a2 = bAt/(Ay)2, 0 = k,Ax and n = k2Ay ([12], p. 196) . It is easily shown that |G| < 1, for At,Ax->0. In calculating the amplification factor for ADI, we note that the additional terms of the ADI scheme are Fourier transformed to real valued terms involving products of sin2(0), sin2(7?), sin(2/3) and sin(2r?). But these terms are added to both the denominator and the numerator of (7); the inequality is unchanged. Thus, the amplification factor for ADI is also bounded by one and ADI is unconditionally stable for equation (5). Using the results from Chapter 2, convergence for ADI is guaranteed by the accuracy and unconditional stability of the 37 scheme. 4.3 Computational Details And Sample Problems We have assumed that the coefficients of (5) are constant, as required for the analysis to hold. To examine accuracy and stability for variable coefficients, and to test the effect of various extrapolation schemes for nonlinear terms, a series of equations with known solutions and of similar form to the model equations were solved using the ADI method. The Fortran code developed for this testing was also found to be useful when coding the program for the model equations. Before discussing the sample problems, we must determine how to evaluate variable coefficients and forcing functions so the ADI equations remain consistent. In addition, the requirement of the ADI scheme for boundary values at t=n+* must be considered, since our equations do not have boundary conditions defined at intermediate time levels. To see how to incorporate variable coefficients into the ADI scheme, consider equation (5), with a, b, c and d functions of x, y and t. The Peaceman-Rachford formula is the same as (6) and it is conventional([9], p. 68) to evaluate the coefficients at t=n+l/2. This is consistent with the CN difference scheme, which is centred about t=n+l/2. The accuracy of the ADI scheme is unchanged from the constant coefficient case, assuming the coeffients are known exactly. The addition of a forcing function to an equation poses a different problem. The forcing function itself must be split in 38 some way to maintain the consistency of the ADI formula. Considering as an example, ut= u^ + uyy + F(x,y,t) , we note first that F(x,y,t) must be evaluated at t=n+l/2. Using ( 1-XD2)uA+* = (l+XD2^ + F(x,y, (n+l/2)At)/2 (9) r (1-XD^)UwH = ( 1+XD.2 )un+*+ F(x,y, (n+l/2)At)/2 and solving for u*** , we find that (10) u*** = (1-XD^) "MO+XD^u* + F,WVV2) where FA+,/2- = F(x ,y, (n +1 /2) At) . Substituting (10) into the second split equation gives (l-XD.2y)(l-XD^)unH = O+XD^Ml+XD^Ju" + ((l-XDE2Y) + (l+XD<2Y))F^yV2 = ( 1+XD^) ( 1+XD^)u* + FN+^ which shows that (9) is consistent. When solving either Peaceman-Rachford equation, time-dependent boundary values at the intermediate time level, t=n+*, are needed. However, this time level does not correspond to any particular time. The values could be approximated by averaging values at t=n+1 and t=n, but the overall accuracy of the scheme would not be maintained ([15]). It is possible to find un+* explicitly in terms of the values at u" and u**' ; the formula follows from equation (7). Adding the split equations gives ( ( 1 - XDE2Y- XD^) + (1 + XD^+ \L\y) ) u«+* = ( 1 - XD£- XD^) u^1 + (l+XD^+XD^u'1 which is just u*** = ((1-XD^-XDeJ()un+, + (1+XD^+XELx)Um)/2 . 39 We require unt"* on the x-direction boundaries, and it is easily seen that the difference operators involve values only on those boundaries. This formula ensures that calculation of the boundary values preserves the order of accuracy of the method. We can now apply the ADI scheme to a series of sample problems. The problems we have chosen are intended to show that ADI has the accuracy predicted for the constant coefficient problem and that, although not as accurate for variable coefficient and nonlinear problems, the method appears to be at least conditionally stable. The sample problems all have the general form (11) u± = auw + buyy + cux + duy where the exact solution u and the coefficients are chosen to produce the various problems of interest. We use the exact solution on the boundaries except at the intermediate time levels, where the boundary conditions are calculated as explained above. The results, consisting of the maximum absolute error at the interior grid points at time = 1.0 for various values of At and Ax, are given in Table 1. For the constant coefficient problem, we set a = b = 0.5 and c = d = 1.0, which has exact solution u = exp(-x-y-t). The results are given in Table 1 under (I). The. errors vary according to the truncation error of the method, which is 0((At)3+At(Ax)2). We see that the errors are smaller than the truncation error would indicate, which is likely due to cancellation of error terms, whose coefficients are various higher derivatives of exp(-x-y-t). 40 To test a variable coefficient problem, (11) is solved with a=0.5/(y+.1)2, b=0.5/(x+.1)2, c=1/(y+.1) and d=l/(x+.1), which has true solution u=exp(-(x+.1)(y+.1)-t). The results are given under (II) and are much the same as the constant coefficient case. For both these cases the error becomes constant after the first ten to twenty time steps (for some cases, errors were checked for 200 time steps), indicating stability of the method for the time and space steps given in the table. One of the nonlinear terms both model equations contain is the square of the convective terms. As described in Chapter 2, these are approximated by discretizing one derivative and treating the other as its coefficient. The sample equation is solved with a and b as for the variable coefficient problem, c=-u</(uA(y+.1)) and d=-Uy/(uy(x+.1)), where u*, Uy are approximated by extrapolation and ux, Uy are known exactly. The exact solution is the same as for the variable coefficient problem. Case I'll in Table 1 contains the results found by approximating ux by D0<u" and Uy by D0yUw and case IV the results with approximations D0x (1 • 5un-0. 5uA'x ) and D0y(1.5u*-0.5urt"1). The two-point extrapolation is slightly more accurate, although not as much better as expected. For several runs with this extrapolation, however, the solution is unstable at later times, oscillating and finally blowing up. This is apparently caused by a parasitic solution introduced numerically by linking three time levels together. (That the one-point extrapolation did not blow up would seem to confirm this.) Since 41 the solution with At=.05, AX=.10 did not blow up, it appears this problem can be avoided by a careful choice of time and space steps. Since the two-point extrapolation requires solutions from two previous time levels, a solution for time=At is necessary. We compared the effect of approximating this solution by a Taylor series expansion with using the exact solution and found no difference in the errors for the test values of Ax and At. To test the effect of nonlinear coefficients of the diffusion terms, we solve (12) setting a=.5u/((y+.1)2u), b=.5u/((x+.1)2u) and c and d as for the variable coefficient problem. Here u = 1 . Su^-0. 5un"' and the exact solution is again u=exp(-(x+.1)(y+.1)-t). The solution at time=At was found approximately by Taylor series expansion. The results are given under (V). All the examples except At=.05, AX=.10 blow up at later times; the scheme is obviously no longer unconditionally stable. As for the previous example, this may be due to the presence of parasitic solutions. To find the convergence rate for this problem, several other runs were made with At and Ax chosen to produce stable results. These solutions are not only less accurate than the corresponding variable coefficient problem, but do not converge at the rate we expect; the errors changes only as At. The reasons for this are not clear but seem to be due to the form of the nonlinear problem. We conclude from these sample problems that the ADI scheme and the methods chosen for handling nonlinear terms produce 42 acceptable solutions. The larger errors and conditional stability for the nonlinear equations are due in part to the specific solution and the fact that the equations have coefficients which contain exact solutions or exact derivatives. The error that results from the approximations to these terms will always be in the same direction, which is certainly not necessarily the case for truly nonlinear equations. It should be noted as well that the error produced by extrapolating the coefficients of the convective terms is of little concern since the model equations are highly diffusive; the convective terms will have only a small contribution to the overall solution. At Ax (I) (II) (III) (IV) (V) .05 .05 1.5E-5 4.1E-6 5.5E-4 3.3E-4 1 .6E-4 .05 . 1 0 5.8E-5 1.5E-5 5.7E-4 3.4E-4 1.3E-4 . 1 0 .05 1.9E-5 5.9E-6 1.1E-3 1.3E-3 2.1E-3 .10 .10 6.2E-5 1.6E-5 1.1E-3 6.4E-4 2.9E-4 .20 .20 2.5E-4 6.7E-5 2.3E-3 4.8E-3 3.9E-4 Table 1. Results of ADI solutions of sample problems. 4.4 Applying ADI To The Model Equations The final form of the ADI scheme for the model equations is found by applying the procedures discussed in this chapter to equations (9) and (10) of Chapter 2. The split equations for (9) are 43 (l-Xa2D2 -Xa4D~ JH"** = ( 1 +Xa , D2 +Xa ^ )Hn + Atf,/2 (13) y ^ ^ ( l-Xa^ -Xa3D^ )U**[ = (l+Xa2D02y + Xa4D<,y )Hn1+*Atf, /2 and for (10) are O-XbzD,2 -Xb„D0y )S*+* = (l+Xb^2 + Xb3D01. ) S* + Atf2/2 (14) ^ y * (1-Xb^ -Xb3De>t )S*M = ( 1+Xb2D^ + Xb1)D^ JS*^ Atf 2/2 where the coefficients, f, and f2 are defined as in Chapter 2 and evaluated at time=n+1/2. 44 V. CHAPTER 4 In the first part of this chapter, several preliminary numerical tests of the model equations are discussed. They involve simplifying the boundary conditions to produce a closed system, allowing us to verify that the approximate solutions are conserving. We also discuss the computational details of both sets of initial and boundary conditions, neither of which, it was found, could be applied in a straightforward manner. The next section is devoted to testing the model equations with these initial and boundary conditions using various time series of input values. A .physical interpretation of the plotted results follows. We conclude with recommendations for continued work on the model, both with steps which might be taken to more accurately represent the physical system being modelled and with improvements to the numerical methods used. 5.1 Preliminary Tests The purpose of these tests is to ensure that applying the ADI scheme and the nonlinear approximations to the model equations conserve mass, as measured by H, and total salinity, as measured by the product of H and S. We consider the equations over the specified domain, but impose closed boundaries. By maintaining constant values of H and S on all boundaries, and setting the velocities u and v to zero on the y-and x-direction boundaries respectively, it is clear that the mass and total salinity in the system should remain constant. For initial conditions we use constant values of H and S 45 everywhere except at a single grid point. At this point, which is in the center of the domain, H or S or both is set to a value larger than the constant field. Three runs were made with constant values of 20 meters for H and 2 p.p.t. for S, and step sizes of Ax=Ay=At=0.1. The results are given in Table 2 below. Mass and HS are calculated simply by summing over the grid points. H(6,6) S(6,6) Mass HS Mass HS t = 0 t = 0 t = 0 t = 0 t=1 00 t=1 00 22.0 2.0 2422.0 4844.0 2420.2 4840.4 22.0 2.5 2422.0 4855.0 2422 . 9 4840.7 28.0 2.5 2428.0 4870.0 2421 . 4 4840.9 Table 2. Results of Conservation Tests The method appears to conserve both quantities very well, with an error over 100 time steps of less than .1 percent for the first case, and slightly larger for the other two. For the crude 11 by 11 grid used, this is very good. The results of the second and third tests also indicate, not surprisingly, that the larger the point sources, the less accurate the spatial derivatives are approximated, and the larger the overall error in terms of conservation. 46 5.2 Initial Conditions At The Land Boundary For purposes of testing, constant field values were assigned as initial conditions to both the height of the upper layer, H, and the salinity defect, S. The values chosen were 20 meters for H and 2 ppt for S, except at the river discharge point. Here the values we assigned to H and S depended on the time series we were testing. To test the problem created by a value of S at the discharge point larger than those nearby, we used a weekly series with S equal to five. The nearly-discontinuous nature of the function at the initial time level resulted in extremely unlikely values of the solutions at time At for both H and S at grid points near the discharge. These results persisted despite several iterations, and there was no indication that the solutions for H and S would converge. The problem stems from having to find the solution of equation (14), from Chapter 1, which is (1) Hy = -(fHSy + 2kSHx + kHS^)/2fS . Due to the large gradients of S created by a large input value, (1) produces large values of H. To avoid the effect of this near-discontinuity, we assume [7] that S at the discharge point has an immediate effect on points in the region; that is, the impact of large changes in salinity defect is smoothed over the region. We arbitrarily chose five grid points in either direction along the land boundary, two grid points into the domain, and varied S linearly over these points. We also want to find a value of H at the discharge point 47 such that H2S is initially in equilibrium with nearby values of H2S. To find such an H, we use the following procedure [7], It is assumed that the discharge region consists of only one grid point, but includes the region up to but not including the grid points immediately north and-south. The land boundary condition for U is U/C = -(f (H2S)y + k(H2S)x ) , with U given at the discharge point and U = 0 elsewhere. The gradient of H2S across the region is negative, due to the Coriolis force, and we approximate it using (H2SL. = (H2SL - 6(H2S) (2) (H2S)KrJ = (H2S)< + 5(H2S) where the discharge point (LAx,KAx) is denoted by the subscript K, the first grid point to the north by K+1, and the first to the south by K-1. 5(H2S) is the increase in H2S from point K to K-1. We see from (2) that a second order approximation to (H2S)y is (H2S)y = 5(H2S)/Ax . To find 8(H2S), we note that U = 0 at all grid points except point K, and so -f(H2S)y = k(H2S)x at these points. Denoting the initial constant field value by (H2S)©« , we can approximate -f(H2S) on the boundary north of the discharge point by -f((H2S)<*> - (H2S)^, )/5Ax . If we assume further that (H2S)X has the same value at the discharge point and at the five grid points north and south, we 48 have k(H2S)x = -f((H2S)eo " (H2S)(C+| )/5Ax and so at the discharge point, -U/C = f6(H2S)/Ax - f((H2S)co - (H2S)|C+) )/5Ax or (3) -UAx/(Cf) = 1.28(H2S) - 0.2((H2S)oo " (H2S)K) . We also require that the grid point inside the domain adjacent to the discharge point be equal to (H2S)o© . Since this point is ((L-1)Ax,KAx), H2S can be expanded in a Taylor series in the x-direction about LAx, which is, to O(Ax), H2S((L-1)Ax,KAx) = H2S(LAx,KAx) - Ax(H2S(LAx,KAx)) or (4) (H2S)oo = (H2SV + (f/5k) ( (H2S)o* " (H2S)IC - 6(H2S)) . We know S everywhere along the boundary, U at the discharge point, and (H2S)oe , our initial condition, so the only unknowns in equations (3) and (4) are 6(H2S) and (H2S)K. Two equations with two unknowns are easily solved, and once 6(H2S) and (H2S)^ are known, H(LAx,(K+1)Ax) and H(LAx,(K-1)Ax) can be determined from (1). H(LAx,(K±i)Ax), i=1,...,5 , are found by assuming H2S is linear from (H2S)(LAx,(K±1)Ax) to (H2S) . Finally, using a first order difference approximation for (H2S) , the values for H((L-1)Ax,(K±i)Ax), i=1,...,5, are easily determined from (2). The result of all this is simply to smooth large initial values of S, and to provide a method of finding H at the discharge region so that initially an equilibrium in (H2S) is maintained. When incorporated into the model, this procedure smoothed the results sufficiently to produce reasonable results. 49 5.3 Finding H On The Boundaries In Chapter 1, we mentioned only that the boundary conditions for H were prescribed on both the ocean and land boundaries. The computational details are given in this section. On the open ocean boundaries H is unknown. It is reasonable to assume, however, that H does not change rapidly, and that any changes that do occur are due to gradients of H (and hence S) within the domain. Except at the first time level, where H is set to the values used initially, we calculate H on the ocean boundaries by extrapolating from values at the two previous time levels using the extrapolation formula discussed in Chapter 2. After the model equation for H has been solved at the current time level, H on the boundary is further modified by extrapolating from the two nearest interior points on a grid line. That is, we use at point (iAx,0), for example. The values on the ocean boundary thus reflect more closely what is occurring within the domain. For H on the land boundary, we must solve equation (1) from the previous section with as an initial value the input value of H at the river discharge point. Since the discharge point is in the middle of the land boundary, (l) must be solved in each direction. When solving it for the southern portion of the boundary, the signs of the y-derivatives must be reversed. In the process of finding a solution to equation (1), a number of difficulties were encountered. One was the effect of 50 the salinity defect values chosen as input at the river discharge point. For test purposes, we assigned input values to S that were large compared to adjacent boundary and interior points. Even with the smoothing procedure described in the last section, this led to a large gradient, Sy, in each direction and produced a stiff equation that was difficult to solve near the discharge point. Integration methods using Ax as a step size only produced physically impossible results. It is necessary therefore to use an integration step size much smaller than Ax, which requires values for S and U at points between our grid points. For S, this is accomplished by a fitting a cubic spline through S at the grid points and evaluating the curve at the required points. For U, which is zero at the grid points corresponding to the land boundary, but nonzero between the discharge point (JAx,KAx) and grid points (JAx,(K+1)Ax) to the north and (JAx,(K-1)Ax) to the south, we need some function which describes its behavior. We chose an exponential curve with a large negative exponent, such that U(JAx,KAx) is the input discharge value and U(JAx,(K+1)Ax), U(JAx,(K-1)Ax) are very nearly zero. This, it is hoped, is a reasonable approximation to U at the "river mouth". To include these nonzero values of U in equation (1), we must reintroduce the term -u/(2fSCH) on the right side of (1), and the ODE becomes (5) Hy = -(U/(CH) + fHSy +2kSHx + kHSx)/(2fS) A further problem arises when evaluating the x-derivative terms. Approximating H and S with the second order difference 51 equation detailed earlier involves values at grid points inside the domain. Unfortunately, these interior points must be at the same time level (say time=n+l) we are solving equation (5) for; values either at the interior points or on the land boundary must be estimated in some way. Rather than extrapolating from previous time levels to provide approximate values at the interior points, we first solve the model equations for H and S. The boundary values for H on all boundaries are extrapolated from the two previous levels in the usual way. Once solved, we have interior values for H and S at time n+1, H and S can now be approximated, and equation (5) solved. This results in new boundary values for H which replace the extrapolated values. A solution to equation (5) is found using the UBC GEARB subroutine, which is an iterative method using a backward differentiation formula as a predictor and the chord method as a corrector. The reader should consult [10] for a more detailed explanat ion. We note here that since H must be extrapolated, and since H can change rapidly making the extrapolation less accurate, we wish to avoid having to find land boundary values at the intermediate time level, which does not correspond to an actual time. We do this by splitting the equation such that we solve for H first in the y-direction, then in the x-direction. Land boundary values are thus only required for times n and n+1. 52 5.4 Boundary Condition For S The boundary condition for the salinity defect is the no-flow condition, which requires setting the first derivative to zero. Using a first order accurate approximation to the derivative has one major advantage; when it is discretized, the matrix equations remain tridiagonal. As we have noted, the salinity defect changes rapidly on the land boundary, particularly near the river discharge region (depending on the input values used). For such cases, it is very important to approximate the boundary condition here with at least the same accuracy as the ADI scheme. A second order accurate approximation to the no-flow condition is usually centred about the grid point on the boundary, requiring a value for S outside the domain. On the open ocean boundaries, this is a reasonable approach, but on the land boundary, such a value has no physical meaning. Instead we use a difference formula involving only interior points; for the western boundary, for example, this is (6) S.j - 4St>j /3 + S^ /3 = 0 , j = 1,...,L. Similar equations apply on the other boundaries. Incorporating (6) into the ADI matrix equations for S, we see that each row corresponding to the boundary contains two elements to the right of the diagonal; the tridiagonality of the system is lost. Rather than try to solve a system with a half-bandwidth of three instead of two, a row reduction is performed on each of these rows before solving the whole system. For example, the first two rows of the matrix equation for S at 53 t=n+* are 1 -4/3 1/3 0 a 21 a23 0 This can be reduced to (1-a2,/(3a23)) (-4/3-a22/(3a23)) 0 A21 3 2 2 A23 Reducing each such row for both left and right boundaries returns the system to tridiagonality, and increases the operation count by only 0(L). 5.5 Tests Of The Model Equations Incorporating these methods for initial and boundary conditions into the FORTRAN program tested in the first section of this chapter, we can proceed with tests of the full model equations. Input to the program consists of test data at the inflow, which was considered to be a single grid point located at the midpoint of the x=LAx boundary. For each time step, the salinity defect, S, and x-direction velocity, u, corresponding to this point are input. The thickness of the upper layer, H, is calculated as explained in section 4.3 so as to be in equilibrium with these values of S and u at time t=0. Input for the first experimental run consists of a pulse inflow, with S=5.0 ppt, u=0.3 m/sec, and the corresponding equilibrium value for H, which is H=17.66. Initial conditions are S=2.0 ppt and H=20.0 m throughout the domain. The inflow values are introduced into the model at time t=0, and remain the same until the final time step, t=50. 54 Based on the results of the first run, the second experimental run is identical with the first until t=l5. At this point the boundary conditions at the inflow are changed to u=0 and dS/dx=0, corresponding to shutting off the supply of fresh water into the model. H is found by extrapolating from the two previous time steps in the usual way. This run is continued until time step 40. 5.6 Interpretation Of Results The results are presented in Figures 3 to 35 and consist of contour plots of H and S and a contour plot of H2S (representing the pressure field) overlayed with a vector plot of velocity. The vector plot is normalized to the largest velocity at the time level, whose value in m/sec is printed on the figure. The point source of fresh water is at the midpoint of the right y-axis, at x=600 km, y=300 km. Due to the algorithm used for contouring, a number of plots have regions of unusual looking contours. The areas where this occurs have values that are all approximately the same, and the contouring program cannot produce accurate contours. On all plots the results for the western half of the domain should be ignored; the accurate contours are always smooth curves. The Figures for t=2 reflect numerical difficulties created by the initial large gradients of S. Despite the smoothing procedure used, these gradients induce unlikely values of H around the inflow point. By using a large number of iterations over the next few time steps, values for H are smoothed, as can 55 be seen by t=5. The region of low salinity moves slowly north, due to Coriolis force, growing larger as more fresh water flows in the system. This is accompanied by an increase in the thickness of the upper layer in this region; to the south, H decreases. The figures for the pressure term H2S and velocity show the flow to be primarily northwards. That the velocities do not lie exactly along the H2S contours reflect the effects of friction; the flow is pushed to the left of these contours. By time step 15 this pattern is well established and Figures for later time steps show only an increase in the upper layer thickness and a spreading by advection of the fresh water along the surface. As noted in [6], increasing amounts of fresh water move west and south as well due to diffusive effects. By time step 30, the limits of the domain and the boundary conditions in force there begin to have an affect on both upper layer thickness and flow. It appears that the flow is being "reflected" away from the northeast corner, causing flow back into the domain. A possible cause is the extrapolation scheme for H, which may give too large a value for H on the boundary in this area. While the results are not shown, the same boundary conditions were used until time step 100. The reflecting flow increases along the northern boundary, but otherwise the results change very little from those at time step 40. Being identical to the first test until t=l5, the results for the second experimental run begin at time step 20 (Figures 56 27, 30 and 33). At t=16, the fresh water flow into the model is "turned off"; we would expect the values to slowly return to their initial conditions. As can be seen in Figures 27-35 this is in fact what happens, as the gradients of S and H decrease due to advection and entrainment of saltier lower layer water into the upper layer. These values decrease from one time step to the next until at t=40 the maximum value of S is about 3.6 ppt while that of H has decreased from a high of about 28 m. to about 22 m. The region of highest salinity defect moves south pushed by the flow, whose speeds (the maximum of which is given on the H2S figures) diminish rapidly. The flow continues to the left of the H2S contours. There appears to be an eddy (Figure 35) just to the south of the inflow point, perhaps due to the "reflection" in the northeast corner of the domain. Both these experimental runs show that the numerical model produces reasonable results with no sign of any instabilies. The second run also shows that a given physical situation, the system returning to rest, is adequatedly modelled. Two other tests were carried out which serve to illustrate further difficulties with solving the model equations. Both produced unstable solutions, so no figures are available. The first used the boundary conditions of the first test discussed above until time step 40. Thereafter, the velocity u remained the same, but S was slowly decreased over time to S=2. The values of H were those found to be in equilibrium. This models the effect on the system of a flow of varying salinity. 57 The salinity defect was decreased at the rate of as little as 0.1 ppt per time step; even at this rate the model quickly developed instabilities in the values of S (soon followed by H). Iterating a large number of times showed the solution for S to be actually diverging. The other unsuccessful test involved the usual boundary conditions, but with a spatial mesh size one quarter of that indicated above. Instabilities resulted within two or three time steps. Increasing the number of iterations again led to a diverging solution. 5.7 Recommendations For Proceeding There are several shortcomings in the model equations and the numerical methods we have chosen to solve them that should initially be considered when improving the model. The numerical difficulties must be dealt with before any changes can be made in the equations, so they will be discussed first. The most obvious problem is illustrated by the two tests discussed at the end of the last section. In both cases, the instabilities seem to be caused by large gradients near the coastal boundary. In the first , it resulted from reducing the input salinity defect from 5.0 ppt after reaching a near-steady state solution. The input salinity defect dropped below S at surrounding grid points producing gradients that apparently created instabilities in equation (5), the land boundary ODE; the resulting values of S on the coast increase to values larger than any entering the model at the inflow. 58 Results from the second indicate that the same problem results from a decrease in mesh size. Since the right-hand side of (5) depends on x-direction derivatives of H and S which are approximated by difference equations, these values become larger (assuming adjacent values of H and S are approximately the same) and (5) gets even stiffer as the resolution increases. Again the solution blows up along the coastal boundary, even when increasing the accuracy of the derivative approximations to fourth order. A method of solving this problem is not readily apparent, although a variation of the averaging scheme used for the initial conditions might be considered. Another numerical problem as yet unresolved is created when the true domain is considered; the land boundary now contains two segments of the coast where run-off enters the model. To the north of the northern segment and to the south of the southern segment, equation (5) can be solved as usual. Between these segments, however, the initial value problem (IVP) becomes a boundary value problem (BVP), as the values for H at both ends are fixed. Considering the difficulties encountered in solving (5) as an IVP, it seems likely that it will present further problems as a BVP, probably having boundary layers at both ends. To more accurately reflect the physical processes of the model, a number of additions should be considered. First, a realistic set of initial conditions for H and S are needed, involving perhaps some sort of historical average. We should also overlay a geostrophic current field, both initially (the model is presently considered at rest), and at each time level. 59 As noted in the previous section, some method of more accurately reflecting behavior on the open-ocean boundaries is required, particularly at the intersection of the coastal boundary. Last, the model equations should include advection current terms and wind forcing, and wind-induced mixing should be added to W , the entrainment term. The effects of these changes on the numerical methods for the simplified model would have to be determined. 60 Figure 4 - Test 1 : S at time step 5 61 0.0 100.C 200.0 300.0 400.0 soo.o soo.o KM . Figure 5 - Test 1: S at time step 10 KM Figure 6 - Test 1: S at time step 15 62 KM -200.0 too. a soo.o Figure 8 - Test 1: S at time step 30 63 Figure 10 - Test 1: S at time step 50 64 ^P-3 200.3 300.3 400.3 SOO.O SOof'3 i I 1 1 - I i ; 0 0.3 100.3 200.3 300.3 400.0 SOO.O SOO.O KM Figure 11 - Test 1: H at time step 2 t I ! : i 1 1 ~ 0 0.0 100..3 230.0 300.0 400.0 SOO.O S00.3 KM Figure 12 - Test 1: H at time step 5 65 Figure 14 - Test 1: H at time step 15 66 Figure 15 - Test 1: H at time step 20 Figure 16 - Test H at time step 30 67 Figure 18 - Test 1: H at time step 50 68 Figure 20 - Test 1: H2S at time step 5 69 Figure 22 - Test 1: H2S at time step 15 70 71 Figure 26 - Test 1: H2S at time step 50 72 Figure' 28 - Test 2: S'at time step 30 73 I I i I l iV o 3.3 lttD.3 2M.3 330.3 40D.3 Sflfl.3 503.3 KM Figure 29 - Test 2: S at time step 40 74 n i ! 1 1 1 r o 0.0 100.0 200.0 300.0 400.0 SOO.O 500.0 KM Figure 30 - Test 2: H at time step 20 KM Figure 3 1 - Test 2: H at time step 30 Figure 32 - Test 2: H at time step 40 76 Figure 34 - Test 2: H2S at time step 30 77 sac.s saol2 9.3 •200.3 300.3 KM 300.3 / 600.3 Figure 35 - Test 2: H2S at time step 40 78 BIBLIOGRAPHY 1. Bo Pedersen, F., A Monograph on Turbulent Entrainment and Friction in Two Layer Stratified Flow, Inst. Hydrodynamics and Hydraulic Engeneering, Tech. Univ. Denmark, Lyngby. Series Paper No. 25. (1980) 2. Forsythe, G.E. and Wasow, W.R., Finite Difference Methods  for Partial Differential Equations, John Wiley & Sons, New York (1960). 3. Gear, C.W., Numerical Initial Value Problems in Ordinary  Differential Equations, Prentice-Hall, Englewood Cliffs, W\~J~. ( 1 971 ) . 4. Gottlieb, D. and Gustafsson, B., Generalized Du Fort-Frankel Methods for Parabolic Initial-Boundary Value Problems, SIAM J. Numer. Anal., 13 (1976) 129-144. 5. Kriess, H. and Oliger, J., Methods of the Approximate  Solution of Time Dependent Problems, GARP Publications Series No. 10 (1973) . 6. LeBlond, P.H., Emery, W.J., Nicol, T.O., A Model of Runoff Driven Coastal Circulation, Submitted for publication to Coastal, Estuarine and Shelf Science, (1983) 7. LeBlond, P.H., Department of Oceanography, University of British Columbia, Vancouver, B.C., personal communication. 8. Lees, M., An Extrapolated Crank-Nicolson Difference Scheme for Quasilinear Parabolic Equations, Numerical Solution of  Partial Differential Equations, Ames, W.F. (Ed.), Barnes & Noble, New York (1969), 193-201. 9. Mitchell, A.R. and Griffiths, D.F., The Finite Difference  Method in Partial Differential Equations, John Wiley & Sons, Chichester (1980) 10. Moore, C.F., Integration of Stiff Systems of Ordinary  Differential Equations Having a Banded Jacobian, Computing Centre, University of British Columbia (1979)., 11. Morris, J.L., On the Numerical Solution of a Heat Equation Associated with a Thermal Print-Head, J. Comp. Phys., 5 (1970), 208-228. 12. Richtmyer, R.D. and Morton, K.W., Difference Methods for  Initial-Value Problems, Interscience Publishers, New York (1967) . 79 13. Royer, T.C, On the Effect of Precipitation and Runoff on Coastal Circulation in the Gulf of Alaska, J. Phys. Oceanogr., 9 (1979), 555-563. 14. Royer, T.C., Hansen, D.V. and Paslinski, D.J., Coastal Flow in the Northern Gulf of Alaska as Observed by Dynamic Topography and Satellite-Tracked Drogued Drift Buoys, J. Phys. Oceanogr., 9 (1979), 785-801. 15. Royer, T.C, Coastal Fresh Water Discharge in the Northeast Pacific, J. Geophys. Res., 83, C3 (1982) 2017-2021. 16. Sommeijer, B.P., Van Der Houwen, P.J. and Verwer, J.G., On the Treatment of Time-Dependent Boundary Conditions in Splitting Methods for Parabolic Differential Equations, Int. J. Numer. Meth. Eng. 17 (1981> 335-346. 17. Varah, J.M., Stability of Difference Approximations to the Mixed Initial Boundary Value Problems for Parabolic Systems, SIAM J. Numer. Anal., 8 (1971) 598-615. 18. Varah, J.M., On the Stability of Boundary Conditions for Separable Difference Approximations to Parabolic Equations, SIAM J. Numer. Anal., 14 (1977) 1114-1125. 19. Varah, J.M., Stability Restrictions on Second Order, Three Level Finite Difference Schemes for Parabolic Systems, SIAM J. Numer. Anal., 17 (1980) 300-309. 


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items