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 1983

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

Item Metadata

Download

Media
UBC_1983_A6_7 N52.pdf
UBC_1983_A6_7 N52.pdf [ 3.65MB ]
Metadata
JSON: 1.0095808.json
JSON-LD: 1.0095808+ld.json
RDF/XML (Pretty): 1.0095808.xml
RDF/JSON: 1.0095808+rdf.json
Turtle: 1.0095808+rdf-turtle.txt
N-Triples: 1.0095808+rdf-ntriples.txt
Citation
1.0095808.ris

Full Text

C . l THE NUMERICAL SOLUTION OF A SYSTEM OF EQUATIONS DESCRIBING COASTAL by THOMAS 0. NICOL B . S c , The U n i v e r s i t y Of B r i t i s h C o l u m b i a , 1976 A THESIS SUBMITTED IN PARTIAL FULFILMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE i n THE FACULTY OF GRADUATE STUDIES I n s t i t u t e Of A p p l i e d M a t h e m a t i c s And S t a t i s t i c s We a c c e p t t h i s t h e s i s a s c o n f o r m i n g t o t h e r e q u i r e d s t a n d a r d THE UNIVERSITY OF BRITISH COLUMBIA O c t o b e r 1983 DIFFUSION/CONVECTION CIRCULATION © Thomas 0. N i c o l , " 1983 In p r e s e n t i n g t h i s t h e s i s in p a r t i a l f u l f i l m e n t of the requirements f o r an advanced degree at the U n i v e r s i t y of B r i t i s h C o l u m b i a , I agree tha t the L i b r a r y s h a l l make i t f r e e l y a v a i l a b l e f o r r e f e r e n c e and s t u d y . I f u r t h e r agree tha t p e r m i s s i o n f o r e x t e n s i v e c o p y i n g of t h i s t h e s i s f o r s c h o l a r l y purposes may be granted by the Head of my Department or by h i s or her r e p r e s e n t a t i v e s . I t i s unders tood tha t c o p y i n g or p u b l i c a t i o n of t h i s t h e s i s for f i n a n c i a l g a i n s h a l l not be a l l o w e d wi thout my w r i t t e n p e r m i s s i o n . Department of I n s t i t u t e Of A p p l i e d Mathematics And S t a t i s t i c s The U n i v e r s i t y of B r i t i s h Columbia 2075 Wesbrook P l a c e Vancouver , Canada V6T 1W5 Date : October 14, 1983 i i A b s t r a c t A m a t h e m a t i c a l model d e s c r i b i n g n e a r s h o r e ocean c u r r e n t s i s examined. The m o t i v a t i o n f o r the problem i s d i s c u s s e d and a d e r i v a t i o n of the model e q u a t i o n s p r e s e n t e d . The model e q u a t i o n s c o n s i s t of a p a i r of c o u p l e d , n o n l i n e a r , p a r a b o l i c p a r t i a l d i f f e r e n t i a l e q u a t i o n s . A n u m e r i c a l method, the C r a n k - N i c o l s o n f i n i t e d i f f e r e n c e scheme, i s p r e s e n t e d f o r s o l v i n g the c o r r e s p o n d i n g l i n e a r e q u a t i o n s , and the convergence and s t a b i l i t y of the method d i s c u s s e d . Methods f o r d e a l i n g w i t h t h e n o n l i n e a r terms, and t h e i r e f f e c t s on a c c u r a c y and s t a b i l i t y , a r e examined. The i n i t i a l and boundary c o n d i t i o n s of the model e q u a t i o n s p r e s e n t s p e c i a l problems, and we d e s c r i b e methods t o s o l v e them. A more e f f i c i e n t f i n i t e d i f f e r e n c e scheme f o r our e q u a t i o n s , based on the C r a n k - N i c o l s o n f o r m u l a , i s i n t r o d u c e d , and i t s advantages d i s c u s s e d . A computer program i n c o r p o r a t i n g t h e s e methods i s developed t o s o l v e the model e q u a t i o n s , and r e s u l t s f o r t e s t d a t a p r e s e n t e d . We co n c l u d e w i t h recommendations f o r a d d i t i o n s t o the program t o model the c u r r e n t s more a c c u r a t e l y . i i i T a b l e of C o n t e n t s A b s t r a c t i i L i s t of F i g u r e s i v Acknowledgements v I . INTRODUCTION 1 11 . CHAPTER 1 4 - 2.1 D e r i v a t i o n Of The Model 4 2.2 The Domain Of The Model 8 2.3 Boundary C o n d i t i o n s 8 2.4 I n i t i a l C o n d i t i o n s 9 2.5 How The Model W i l l Be Used 10 I I I . CHAPTER 2 13 3.1 N o t a t i o n 13 3.2 F i n i t e D i f f e r e n c e s 14 3.3 C o n s i s t e n c y , S t a b i l i t y And Convergence 19 3.4 V a r i a b l e C o e f f i c i e n t s 24 3.5 S t a b i l i t y Of Boundary C o n d i t i o n s 25 3.6 N o n l i n e a r Terms 27 3.7 S o l v i n g The E q u a t i o n s 29 IV. CHAPTER 3 32 4.1 I n t r o d u c t i o n To ADI 33 4.2 A c c u r a c y And S t a b i l i t y 35 4.3 C o m p u t a t i o n a l D e t a i l s And Sample Problems 37 4.4 A p p l y i n g ADI To The Model E q u a t i o n s 42 V. CHAPTER 4 44 5.1 P r e l i m i n a r y T e s t s 44 5.2 I n i t i a l C o n d i t i o n s At The Land Boundary 46 5.3 F i n d i n g H On The Bou n d a r i e s 49 5.4 Boundary C o n d i t i o n For S 52 5.5 T e s t s Of The Model E q u a t i o n s 53 5.6 I n t e r p r e t a t i o n Of R e s u l t s 54 5.7 Recommendations For P r o c e e d i n g 57 BIBLIOGRAPHY 78 i v L i s t of F i g u r e s 1. Domain of the C o a s t a l C i r c u l a t i o n Model 11 2. Domain of the N u m e r i c a l Model 12 3. Test 1 : S a t time s t e p 2 60 4. Test 1 : S a t time s t e p 5 61 5. T e s t 1 : S a t time s t e p 10 61 6. Test 1 : S a t time s t e p 15 62 7. Test 1 : S a t time s t e p 20 62 8. Test 1 : S a t time s t e p 30 63 9. Test 1 : S a t time s t e p 40 63 10. Test 1: S a t time s t e p 50 64 11. Test 1 : H a t time s t e p 2 64 12. Test 1 : H at time s t e p 5 65 13. Test 1 : H a t time s t e p 10 65 14. Test 1 : H a t time s t e p 15 66 15. Test 1: H a t time s t e p 20 66 16. Test 1 : H a t time s t e p 30 67 17. Test 1 : H a t time s t e p 40 67 18. T e s t 1 : H a t time s t e p 50 68 19. Test 1: H 2S a t time s t e p 2 68 20. T e s t 1: H 2S a t time s t e p 5 69 21. Test 1: H 2S a t time s t e p 10 69 22. Test 1: H 2S a t time s t e p 15 70 23. Test 1: H 2S a t time s t e p 20 70 24. Test 1: H 2S a t time s t e p 30 71 25. Test 1: H 2S a t time s t e p 40 71 26. Test 1: H 2S a t time s t e p 50 72 27. Test 2: S a t time s t e p 20 72 28. Test 2: S a t time s t e p 30 73 29. Test 2: S a t time s t e p 40 73 30. Test 2: H a t time s t e p 20 74 31. T e s t 2: H a t time s t e p 30 75 32. Test 2: H at time s t e p 40 75 33. Test 2: H 2S a t time s t e p 20 76 34. Test 2: H 2S a t time s t e p 30 77 35. Test 2: H 2S a t time s t e p 40 77 V Acknowledgement The c o m p l e t i o n of t h i s t h e s i s has depended on the good a d v i c e , g i v e n over many h o u r s , from J im V a r a h , the m u l t i t u d e of sugges t ions from Paul L e b l o n d , and the o r i g i n a l idea of B i l l Emery; and from a l l , a grea t d e a l of p a t i e n c e . 1 I . INTRODUCTION S t u d i e s of ocean c u r r e n t s i n c o a s t a l r e g i o n s have shown t h a t f r e s h water r u n - o f f i s an i m p o r t a n t f a c t o r i n l a r g e - s c a l e c i r c u l a t i o n (see [6] f o r s e v e r a l e x amples). W h i l e these s t u d i e s c o n s i d e r the e f f e c t of l a r g e r i v e r s on the c o a s t a l ocean, the i n f l u e n c e of r u n - o f f over a l a r g e c o a s t a l r e g i o n has not been c l o s e l y examined. Dr. P a u l LeBlond and Dr. W i l l i a m Emery of the Department of Oceanography a t the U n i v e r s i t y of B r i t i s h Columbia have undertaken a study of r i v e r r u n - o f f and s a l i n i t y d i s t r i b u t i o n o f f the c o a s t of s o u t h e r n B r i t i s h Columbia and n o r t h e r n Washington S t a t e i n an e f f o r t t o u n d e r s t a n d t h e i r e f f e c t s on a l o n g s h o r e c u r r e n t s . One a s p e c t of the study has been the development of a ma t h e m a t i c a l model t o d e s c r i b e those c o a s t a l c u r r e n t s c r e a t e d by r i v e r r u n - o f f . The model c o n s i s t s of a s e t of e q u a t i o n s t h a t d e s c r i b e nearshore hydrodynamics i n terms of changes i n s a l i n i t y , and thus can be used t o measure the impact of c o a s t a l r u n - o f f . The main purpose of t h i s t h e s i s i s t o determine a n u m e r i c a l method to s o l v e the e q u a t i o n s d e f i n i n g the model, and to d e v e l o p a computer program implementing the method. As w e l l as b e i n g of i n t e r e s t from a p h y s i c a l o c e a n o g r a p h i c p o i n t of view, i t i s hoped t h a t the r e s u l t s from the model may be used t o g a i n an u n d e r s t a n d i n g of o f f s h o r e f i s h e r i e s . In p a r t i c u l a r , i t has been observed t h a t i n most y e a r s a l a r g e p e r c e n t a g e of F r a s e r R i v e r salmon m i g r a t e v i a t h e S t r a i t of Juan de Fuca, w h i l e the remainder r e t u r n through Johnstone S t r a i t a t the n o r t h e r n end of Vancouver I s l a n d . In some y e a r s , however, a 2 m a j o r i t y r e t u r n t h rough Johnstone S t r a i t . The reasons f o r t h i s a r e not known, but i t i s l i k e l y t h a t the l a r g e - s c a l e c o a s t a l c i r c u l a t i o n p l a y s a l a r g e p a r t i n which r o u t e i s chosen. An u n d e r s t a n d i n g of these mechanisms may l e a d t o an e x p l a n a t i o n f o r t h i s b e h a v i o r . Chapter One c o n t a i n s a d i s c u s s i o n of the d e r i v a t i o n of the problem, l e a d i n g t o the e q u a t i o n s t o be s o l v e d . D e f i n i t i o n s of the v a r i a b l e s and parameters f o l l o w , and the boundary and i n i t i a l c o n d i t i o n s a re s p e c i f i e d . A d e s c r i p t i o n of how the model i s t o be used c o n c l u d e s the c h a p t e r . In Chapter Two, the n o t a t i o n t o be used t o d e s c r i b e f i n i t e d i f f e r e n c e schemes i s i n t r o d u c e d . The c o n c e p t s of c o n s i s t e n c y , a c c u r a c y and s t a b i l i t y , as they r e l a t e t o a p p l i c a t i o n of the C r a n k - N i c o l s o n d i f f e r e n c i n g of a l i n e a r i z e d form of the e q u a t i o n s , a re d i s c u s s e d . The boundary c o n d i t i o n s a r e then s p e c i f i e d i n terms of f i n i t e d i f f e r e n c e s and t h e i r e f f e c t s on s t a b i l i t y examined. The n o n l i n e a r terms i n the model e q u a t i o n s a r e d i s c u s s e d and methods of d e a l i n g w i t h them p r e s e n t e d . In Chapter Three, we examine the systems of l i n e a r e q u a t i o n s produced by the C r a n k - N i c o l s o n scheme, and d i s c u s s methods t o s o l v e them. More e f f i c i e n t d i f f e r e n c e schemes f o r the e q u a t i o n s a r e the a l t e r n a t i n g d i r e c t i o n i m p l i c i t , or ADI, methods. They a r e d e s c r i b e d and the c o n c e p t s of Chapter Two are extended t o one of t h e s e , the Peaceman-Rachford ADI scheme, f o r a l i n e a r e q u a t i o n . Methods of t r e a t i n g the a d d i t i o n a l problems p r e s e n t e d by the model e q u a t i o n s , i n c l u d i n g lower o r d e r space d e r i v a t i v e s , v a r i a b l e c o e f f i c i e n t s and n o n l i n e a r terms, are 3 examined, and t e s t e d by s o l v i n g sample problems w i t h known s o l u t i o n s . R e s u l t s from computer runs a r e shown and the a c c u r a c y and s t a b i l i t y of the s o l u t i o n s d i s c u s s e d . F i n a l l y , we d i s c r e t i z e the f u l l e q u a t i o n s i n terms of the Peaceman-Rachford scheme. In Chapter F o u r , the ADI method i s a p p l i e d t o the model e q u a t i o n s w i t h s i m p l i f i e d b o u n d a r i e s t o ensure the d i f f e r e n c e scheme s a t i s f i e s b a s i c p h y s i c a l p r o p e r t i e s . The d i s c o n t i n u o u s n a t u r e of the i n i t i a l c o n d i t i o n s i s d i s c u s s e d , and we show how t h i s problem i s d e a l t w i t h . The l a n d boundary c o n d i t i o n f o r the t h i c k n e s s l a y e r i s examined and the problems and methods we use t o h a n d l e them a r e e x p l a i n e d . Next the f o r m u l a r e q u i r e d by the ADI method f o r c a l c u l a t i n g the s a l i n i t y d e f e c t boundary c o n d i t i o n s i s p r e s e n t e d . Sample runs and t h e i r r e s u l t s f o r the f u l l model e q u a t i o n s are g i v e n , w i t h b r i e f p h y s i c a l i n t e r p r e t a t i o n s . Recommendations f o r p r o c e e d i n g w i t h the model c o n c l u d e the c h a p t e r . 4 I I . CHAPTER 1 2.1 D e r i v a t i o n Of The Model The p r i m a r y c o n c e r n of t h i s t h e s i s i s the s e t of n u m e r i c a l problems p r e s e n t e d by the m a t h e m a t i c a l model of c o a s t a l c i r c u l a t i o n . The d e r i v a t i o n t h a t f o l l o w s i s b r i e f and i n c l u d e s l i t t l e j u s t i f i c a t i o n , from a p h y s i c a l oceanography v i e w p o i n t , f o r the assumptions and accomodations made i n f o r m u l a t i n g the f i n a l e q u a t i o n s . A complete d e r i v a t i o n and d i s c u s s i o n of the i m p l i c a t i o n s of the model i s a v a i l a b l e i n [ 6 ] . Our d e r i v a t i o n i s a c o n d e n s a t i o n of t h a t i n f o r m a t i o n . The model c o n s i d e r s the n e a r s h o r e ocean as h a v i n g two l a y e r s ; a lower l a y e r of c o n s t a n t d e n s i t y , and hence s a l i n i t y , and an upper l a y e r of v a r y i n g d e n s i t y and t h i c k n e s s . The motion of t h i s upper l a y e r i s d e s c r i b e d by the t w o - d i m e n s i o n a l momentum and c o n s e r v a t i o n e q u a t i o n s . I f we c o n s i d e r a time s c a l e g r e a t e r than f i v e days, l o c a l a c c e l e r a t i o n s can be i g n o r e d , and the motion d e s c r i b e d i s p r i m a r i l y g e o s t r o p h i c . In a d d i t i o n , we want o n l y waves w i t h a Rossby number much l e s s than one; t h a t i s , U 0 f / A x << 1, where U 0 i s the speed s c a l e and f the C o r i o l i s f o r c e . The space s c a l e , Ax, must t h e r e f o r e be much g r e a t e r than 1 km. We assume a t t h i s p o i n t t h a t a v a l u e of Ax > 10 km. i s suf f i c i e n t . Under these a s s u m p t i o n s , the momentum e q u a t i o n s f o r the upper l a y e r are (1) - f v = (-1/p , )9p,/9x - Fx (2) f u = (-1/p,)9p,/9y ~ Fy 5 where: (u,v) = ( u ( x , y , t ) , v ( x , y , t ) ) = components of v e l o c i t y ; p-, = p r e s s u r e i n the upper l a y e r ; p, = d e n s i t y of the upper l a y e r ; (Fx,Fy) = components of f r i c t i o n . S i n c e the lower l a y e r i s assumed t o be a t r e s t , and u s i n g the assumption t h a t the s u r f a c e h e i g h t i s much l e s s than the t h i c k n e s s of the upper l a y e r and can t h e r e f o r e be n e g l e c t e d , (1) and (2) can be i n t e g r a t e d over the upper l a y e r t o y i e l d (3) - f V = -g/2 ( H 2 A p / p , ) x - Fx (4) fU = -g/2 ( H 2 A p / p , ) y - Fy where, (U,V) = ( U ( x , y , t ) , V ( x , y , t ) ) = U = components of t r a n s p o r t , where t r a n s p o r t i s v e l o c i t y i n t e g r a t e d over the upper l a y e r ; H = H ( x , y , t ) = the t h i c k n e s s of the upper l a y e r ; p 2 = d e n s i t y of the lower l a y e r ; Ap = p 2 - p, = d e n s i t y d i f f e r e n c e between l a y e r s ; g = g r a v i t y ; In the e q u a t i o n s above, we have adopted the u s u a l n o t a t i o n f o r p a r t i a l d e r i v a t i v e s ; they a r e i n d i c a t e d here and throughout t h i s t h e s i s by 3 /3x ( f o r example) as i n e q u a t i o n (1) or by s u b s c r i p t s of the independent v a r i a b l e . In terms of t r a n s p o r t s , the volume and s a l i n i t y c o n s e r v a t i o n e q u a t i o n s a r e (5) + VH-U = We (6) ( H S U ) + + V^US = W eS b 6 where: SM = S f c ( x , Y / t ) , s a l i n i t y in the upper l a y e r ; S^ = a cons tant lower l a y e r s a l i n i t y ; W e = the entra inment v e l o c i t y ; V K = the d i v e r g e n c e o p e r a t o r . D e n s i t y can be c o n s i d e r e d a f u n c t i o n of s a l i n i t y , wi th an e q u a t i o n of s t a t e g iven a p p r o x i m a t e l y by (7) Ap = p , a ( T ) S* where a i s a f u n c t i o n of t e m p e r a t u r e , T , and i s assumed c o n s t a n t for the range of temperatures e n c o u n t e r e d . I n t r o d u c i n g S ( x , y , t ) = S f e ( x , y , t ) - S I A ( x , y , t ) , the d i f f e r e n c e i n s a l i n i t y between lower and upper l a y e r s (so S ( x , y , t ) > 0 ) , and u s i n g (7) to r e l a t e d e n s i t y to s a l i n i t y , e q u a t i o n s (3) through (6) can be r e w r i t t e n as (8) - f V = - ( g a / 2 p , ) ( H 2 S ) X - kU (9) fU = - ( g a / 2 p , ) ( H 2 S ) y - kv (10) H t + VK.U= W € (11) (HS)^ + VK-US = 0 Here we have taken f r i c t i o n components ( F x , F y ) to be p r o p o r t i o n a l to the components of v e l o c i t y ( u , v ) , w i th k the f r i c t i o n c o e f f i c i e n t , and i n t e g r a t e d over the upper l a y e r to produce t r a n s p o r t s . M u l t i p l y i n g (8) by k and (9) by f and combining g i v e s (12) ( f 2 + k 2 ) U = - ( g a / 2 P l ) ( f ( H 2 S ) y + k ( H 2 S ) x ) (13) ( f 2 + k 2 ) V = - ( g a / 2 P l ) ( - f ( H 2 S ) ; < + k ( H 2 S ) y ) S u b s t i t u t i n g (12) and (13) i n t o (10) y i e l d s , w i th C = g a k / ( 2 p , ( f 2 + k 2 ) ) 7 (14) H t = Wt + CV 2( H 2S ) . S i m i l a r l y , s u b s t i t u t i n g (12) and (13) i n t o (11) g i v e s (HS\ = C ( k S V 2 ( H 2 S ) + f J ( S , H 2 S ) + kVS-V(H 2S)) or (15) S t = C/H ( k S V 2 ( H 2 S ) + f J ( S , H 2 S ) + kVS V ( H 2 S ) ) - ( S / H ) H t where V 2 i s the L a p l a c i a n o p e r a t o r , J i s the J a c o b i a n o p e r a t o r i n terms of x and y and V i s the g r a d i e n t o p e r a t o r . The f u n c t i o n H 2S r e p r e s e n t s the d i s t r i b u t i o n of the d e n s i t y mass f i e l d ; we expect f r e s h water r u n - o f f t o c r e a t e g r a d i e n t s i n H 2S, c a u s i n g c u r r e n t s i n the model domain. Both (14) and (15) are n o n l i n e a r , p a r a b o l i c p a r t i a l d i f f e r e n t i a l e q u a t i o n s and are the g o v e r n i n g e q u a t i o n s f o r our model of the c o a s t a l c i r c u l a t i o n . S o l v i n g (14) y i e l d s H ( x , y , t ) , the t h i c k n e s s of the upper l a y e r , w h i l e s o l v i n g (15) g i v e s S ( x , y , t ) , the s a l i n i t y d e f e c t , the d i f f e r e n c e between Ŝ , and S b. The parameters and c o n s t a n t s i n (14) and (15) are d e f i n e d and a s s i g n e d v a l u e s as f o l l o w s : g = 9.8 m/sec 2 f = 1.1E-4 s e c " 1 k = 0.25E-4 s e c " 1 a =0.77 p p t " 1 p, = 10 3 ppt We = 6.0E-5 F 2 v where F 2 = Froude number = ( U 2 + V 2 ) / ( g a S H 3 ) / p , v = s u r f a c e v e l o c i t y = ( U 2 + V 2 )''1*/H 8 F u r t h e r d i s c u s s i o n of parameters k, a, W and j u s t i f i c a t i o n f o r the v a l u e s chosen can be found i n [ 6 ] . The components of t r a n s p o r t , (U,V), and hence the v e l o c i t y v, can be found a l g e b r a i c a l l y from (12) and ( 1 3 ) . 2.2 The Domain Of The Model The domain of the e q u a t i o n s i s an a r e a bounded by open ocean t o the n o r t h , west and s o u t h , and on the e a s t by the c o a s t from the Olympic P e n i n s u l a t o the Queen C h a r l o t t e I s l a n d s (see F i g u r e 1). The c o a s t i s broken between the Olympic P e n i n s u l a and Vancouver I s l a n d by the S t r a i t of Juan de Fuca and n o r t h of Vancouver I s l a n d by Queen C h a r l o t t e Sound. I t i s t h rough th e s e two passages t h a t r i v e r r u n - o f f e n t e r s c o a s t a l w a t e r s . For t h e purposes of a n u m e r i c a l model, the domain has been i d e a l i z e d t o a square, w i t h x d e s i g n a t i n g the e a s t - w e s t a x i s and y the n o r t h - s o u t h a x i s , as i l l u s t r a t e d by F i g u r e 2. The o r i g i n f o r t h i s domain i s the southwest c o r n e r of the s q u a r e . The l a n d boundary i s t h e r e f o r e taken t o be p a r a l l e l t o the y - a x i s . 2 . 3 Boundary C o n d i t i o n s The e q u a t i o n f o r the s a l i n i t y d e f e c t , ( 1 2 ) , r e q u i r e s t h a t the g r a d i e n t on a l l b o u n d a r i e s be z e r o ; t h a t i s , (16) 3S/3n = 0 where 9 /3n i s the normal d e r i v a t i v e . T h i s i s the u s u a l form of the n o - f l o w c o n d i t i o n a c r o s s a boundary. Where the l a n d boundary i s broken, r e p r e s e n t i n g f r e s h water f l o w i n t o the domain, S i s p r e s c r i b e d from h i s t o r i c a l d a t a . 9 The boundary c o n d i t i o n s f o r H are a l i t t l e more d i f f i c u l t . I t i s assumed t h a t near the open ocean b o u n d a r i e s , H does not change r a p i d l y and can t h e r e f o r e be approximated from those v a l u e s w i t h i n t h e domain t h a t a r e near the boundary. The a c t u a l method f o r f i n d i n g H on these b o u n d a r i e s w i l l be p r e s e n t e d when the n u m e r i c a l scheme has been e x p l a i n e d . On the l a n d boundary, the x - d i r e c t i o n t r a n s p o r t must be z e r o , as t h e r e can be no f l o w a c r o s s l a n d , and t h e r e f o r e , from e q u a t i o n ( 1 2 ) , f ( H 2 S ) y + k ( H 2 S ) x = 0 Expanding the d e r i v a t i v e terms y i e l d s 2fSHH y + fH 2Sy + 2kSHH;t + kH2S^ = 0 or (17) Hy = -(fHSy + 2kSH x + k H S x ) / 2 f S To f i n d H on t h e l a n d boundary, we note t h a t H A can be a p p r o x i m a t e d n u m e r i c a l l y and Sy i s known from the s o l u t i o n of ( 1 5 ) . Thus, (17) i s a f i r s t o r d e r i n i t i a l v a l u e o r d i n a r y d i f f e r e n t i a l e q u a t i o n of H i n y, w i t h the i n i t i a l v a l u e f o r the e q u a t i o n b e i n g the v a l u e of H g i v e n a t the r i v e r d i s c h a r g e r e g i o n s of the boundary. 2 . 4 I n i t i a l C o n d i t i o n s The s o l u t i o n of (14) and (15) r e q u i r e s i n i t i a l c o n d i t i o n s H(x,y,0) and S ( x , y , 0 ) . S i n c e the model i s i n t e n d e d f o r use w i t h h i s t o r i c a l r i v e r r u n - o f f d a t a , some h i s t o r i c average of t h i c k n e s s and s a l i n i t y data w i l l e v e n t u a l l y be c o n s i d e r e d as i n i t i a l v a l u e s . For the purposes of our i n i t i a l computer r u n s , c o n s t a n t f i e l d v a l u e s w i l l be used f o r both H(x,y,0) and 1 0 S ( x , y , 0 ) . 2.5 How The Model W i l l Be Used E q u a t i o n s (14) and ( 1 5 ) , w i t h boundary c o n d i t i o n s as s p e c i f i e d , d e s c r i b e the e f f e c t of r i v e r r u n - o f f on c o a s t a l c i r c u l a t i o n . The purpose of the model i s t o t r y t o c l o s e l y a pproximate a c t u a l c u r r e n t s and s a l i n i t i e s by u s i n g c o l l e c t e d d a t a i n boundary c o n d i t i o n s . To t h a t end, r i v e r d i s c h a r g e data has been c o l l e c t e d f o r the c o a s t s of Washington S t a t e , B.C. and A l a s k a . From the s e v a l u e s , s u r f a c e v e l o c i t i e s (and hence t r a n s p o r t s ) a t the i n f l o w r e g i o n s were c a l c u l a t e d and • a weekly time s e r i e s produced t o be used as boundary v a l u e s f o r the model e q u a t i o n s . S i m i l a r l y , d a i l y s u r f a c e s a l i n i t y r e a d i n g s , r e c o r d e d at l i g h t h o u s e s a l o n g the c o a s t , have been smoothed t o a weekly s e r i e s f o r use as s a l i n i t y boundary v a l u e s . Upper l a y e r t h i c k n e s s v a l u e s have been appr o x i m a t e d from t r a n s p o r t s and s u r f a c e s a l i n i t i e s . A p e r i o d of a p p r o x i m a t e l y twenty y e a r s of a c t u a l or i n t e r p o l a t e d weekly t r a n s p o r t and s u r f a c e s a l i n i t y r e c o r d s are a v a i l a b l e t o be used as boundary c o n d i t i o n s f o r the model. The aim of the study i s t o approximate the c o a s t a l c i r c u l a t i o n f o r t h a t l e n g t h of t i m e , and t o compare the r e s u l t s w i t h the known salmon m i g r a t i o n d a t a . From t h i s i t i s hoped t h a t some i n f e r e n c e s can be drawn about the r e l a t i o n s h i p between the two. 11 F i g u r e 1 - Domain of the C o a s t a l C i r c u l a t i o n Model F i g u r e 2 - Domain of the Numerica l Model 1 3 I I I . CHAPTER 2 3.1 Notat i o n To s o l v e the p a r t i a l d i f f e r e n t i a l e q u a t i o n s (1) H t = C V 2 ( H 2 S ) + We (2) St = C/H ( k S V 2 ( H 2 S ) + f J ( S , H 2 S ) + k V S V ( H 2 S ) ) - (S/H)H_,. n u m e r i c a l l y , the domain i s d i s c r e t i z e d by o v e r l a y i n g a r e c t a n g u l a r g r i d . I n i t i a l l y , t h i s domain w i l l be a s i m p l i f i e d v e r s i o n of t h a t shown i n F i g u r e 2 by c o n s i d e r i n g i t t o have j u s t a s i n g l e d i s c h a r g e p o i n t . The model e q u a t i o n s a r e v a l i d over the r e g i o n r e p r e s e n t e d by t h i s g r i d , which has g r i d s p a c i n g s of Ax and Ay i n the x and y d i r e c t i o n s . C o o r d i n a t e s on the g r i d a r e s p e c i f i e d by ( i A x , j A y ) , w i t h i=1,...,L and j=1,...,M, and LAx = MAy = 600 km. The time d i r e c t i o n i s a l s o d i s c r e t i z e d by d i v i d i n g the time domain i n N e q u a l time s t e p s of s i z e A t . A s p e c i f i c time i s r e p r e s e n t e d by nAt. The s o l u t i o n s of (1) and (2) w i l l be a p p r o x i m a t e d a t any g r i d p o i n t by the d i s c r e t e s o l u t i o n s H ( i A x , j A y , n A t ) and S ( i A x , j A y , n A t ) , r e s p e c t i v e l y . For ease of n o t a t i o n , the approximate s o l u t i o n s are s i m p l i f i e d t o H-- and To s p e c i f y each of (1) and (2) i n d i s c r e t e terms, the p a r t i a l d e r i v a t i v e s must be approximated i n some way. F i n i t e element or f i n i t e d i f f e r e n c e methods a r e used most o f t e n f o r f i n d i n g n u m e r i c a l s o l u t i o n s t o p a r t i a l d i f f e r e n t i a l e q u a t i o n s . There are a number of f i n i t e element methods f o r a p p r o x i m a t i n g (1) and ( 2 ) , but they a r e f o r the most p a r t d i f f i c u l t t o program and i n g e n e r a l a r e c o m p u t a t i o n a l l y e x p e n s i v e t o s o l v e . F i n i t e 1 4 d i f f e r e n c e s , i n a d d i t i o n t o b e i n g e a s i e r t o program and l e s s c o s t l y t o s o l v e , approximate d e r i v a t i v e s u s i n g d i v i d e d d i f f e r e n c e s which a r e more s u i t e d t o our r e c t a n g u l a r domain. 3 . 2 F i n i t e D i f f e r e n c e s We w i l l use c e n t r e d d i f f e r e n c e s f o r the f i r s t and second o r d e r space d e r i v a t i v e s and the f i r s t o r d e r time d e r i v a t i v e . They w i l l be d e f i n e d u s i n g l i n e a r o p e r a t o r s D^ and D j x f o r the f i r s t and second o r d e r x - d i r e c t i o n d e r i v a t i v e s , and D e y and D j y f o r the y - d i r e c t i o n d e r i v a t i v e s . For some f u n c t i o n u, u = u ( x , y , t ) , l e t D„X uij = (u„t|. - u-„(. )/2Ax and D*x ujj = ( u l t t )- - 2UJJ + u;_(i. ) / ( A x ) 2 From T a y l o r s e r i e s e x p a n s i o n of u ( ( i + 1 ) A x , j A y , t ) and u ( ( i - 1 ) A x , j A y , t ) , i t i s e a s i l y seen t h a t uK = D o x u + 0 ( ( A x ) 2 ) and u x x ; = D„2X u + 0( (Ax) 2 ) where 0 ( ( A x ) 2 ) i s the o r d e r of the remainder term. The D^ and D 2 X l i n e a r o p e r a t o r s are s a i d t o be second o r d e r a c c u r a t e a p p r o x i m a t i o n s . S i m i l a r l y , Doyf and D^ a r e l i n e a r o p e r a t o r s f o r the y - d i r e c t i o n d e r i v a t i v e s and are a c c u r a t e t o 0 ( ( A y ) 2 ) . U s i n g t h e l i n e a r o p e r a t o r s D^ , D|x , D^ and Djy , f i n i t e d i f f e r e n c e methods can be c o n s t r u c t e d t o approximate the d e r i v a t i v e s i n ( 1 ) and ( 2 ) . The s o l u t i o n s H and S a r e found by marching ahead i n time i n some f a s h i o n ; the s o l u t i o n s a t t = ( n + l ) A t depend on v a l u e s a t lower time l e v e l s . 1 5 I f H and S are found i n d i v i d u a l l y , u s i n g v a l u e s from p r e v i o u s time l e v e l s , the method i s s a i d t o be e x p l i c i t . While c o m p u t a t i o n a l l y v e r y s i m p l e , e x p l i c i t methods a r e r e s t r i c t e d t o some maximum s i z e of time s t e p used t o advance i n t i m e . U s i n g a time s t e p g r e a t e r than t h i s maximum, which depends on Ax and Ay, r e s u l t s i n s o l u t i o n s t h a t grow e x p o n e n t i a l l y . A v o i d i n g such i n s t a b i l i t i e s by u s i n g a time s t e p l e s s than the maximum s e v e r e l y r e s t r i c t s the c h o i c e of Ax and Ay. As was mentioned i n Chapter 1, v a r i o u s assumptions used i n d e r i v i n g the model e q u a t i o n s put c o n s t r a i n t s on the s i z e of Ax and Ay. An e x p l i c i t method i s t h e r e f o r e not a p p r o p r i a t e f o r f i n d i n g s o l u t i o n s t o our equat i o n s . A f i n i t e d i f f e r e n c e method i s i m p l i c i t i f the v a l u e s of the dependent v a r i a b l e at p o i n t s a t the next time l e v e l a r e i n t e r r e l a t e d . I m p l i c i t methods thus r e q u i r e s o l u t i o n of a system of l i n e a r e q u a t i o n s a t each time l e v e l . W h i l e n a t u r a l l y t a k i n g more computing time than e x p l i c i t methods, i m p l i c i t methods g e n e r a l l y do not p l a c e r e s t r i c t i o n s on the time s t e p , beyond t h a t r e q u i r e d f o r a c c u r a c y . A well-known i m p l i c i t method f o r s o l v i n g p a r a b o l i c e q u a t i o n s i s the C r a n k - N i c o l s o n method. In terms of the l i n e a r o p e r a t o r s d e f i n e d above, the s i m p l e t w o - d i m e n s i o n a l d i f f u s i o n equat i o n u t = UKX + u y y i s a p p r o x i m a t e d by the C r a n k - N i c o l s o n (or CN) method a t ( i A x , j A y , ( n + l / 2 ) A t ) by - u|j = At(D£ ( + u{] )/2 + D/y ( u?j' + ujj ) / 2 ) 1 6 The d i f f e r e n c e scheme i s c e n t r e d about t=n+l/2 t o p r e s e r v e the second o r d e r a c c u r a c y i n t i m e , and i s a t w o - l e v e l d i f f e r e n c e e q u a t i o n , s o l v i n g f o r t=n+1 from v a l u e s at t=n. The d i s c u s s i o n i n the remainder of the c h a p t e r which c o n c e r n s f i n i t e d i f f e r e n c e e q u a t i o n s and t h e i r a p p l i c a t i o n t o (1) and (2) w i l l be l i m i t e d t o the CN method. B e f o r e a p p l y i n g the CN scheme t o e q u a t i o n s (1) and ( 2 ) , we must d e c i d e how the s p a t i a l d e r i v a t i v e s are t o be expanded. The e q u a t i o n s have d e r i v a t i v e s of p r o d u c t s of the dependent v a r i a b l e s which must be e x p r e s s e d as d i f f e r e n c e e q u a t i o n s . C o n s i d e r as an example the second o r d e r x - d i r e c t i o n d e r i v a t i v e , ( H 2 S ) ^ , i n e q u a t i o n ( 1 ) . A p p l y i n g the p r o d u c t r u l e i n the u s u a l way y i e l d s (3) 2HSH K ) C + 2H 2S X+ 4HH XS A + H 2 S ^ w h i l e c o n s i d e r i n g (H 2S) as a p r o d u c t of f u n c t i o n s (HS) and H g i v e s (4) (HS)H X ; < + 2(HS) XH^ + H ( H S ) X X Upon d i s c r e t i z a t i o n , t h e s e e x p r e s s i o n s a r e c l e a r l y u n e q u a l , except f o r f o r t u i t o u s v a l u e s of H and S. O b v i o u s l y , the s o l u t i o n s of both (1) and (2) w i l l depend on which way the (H 2S) d e r i v a t i v e s a r e expanded. We a r e s o l v i n g e q u a t i o n (1) f o r H, and i t seems r e a s o n a b l e t h a t those f u n c t i o n s m u l t i p l y i n g H s h o u l d be c o n s i d e r e d as s i n g l e , s e p a r a t e f u n c t i o n s . That the d e r i v a t i v e s s h o u l d be expanded as (4) above was c o n f i r m e d by e x p e r i m e n t a l computer runs (as e x p l a i n e d i n Chapter 4 ) . Expanding e q u a t i o n (1) t h e r e f o r e y i e l d s 1 7 (5) H t = C ( ( H S ) V 2 H +2(HS) ; <H ; t + 2 ( H S ) / H y +HV 2(HS)) + We . The HV 2(HS) term may be t r e a t e d as the o t h e r terms i n ( 5 ) ; t h a t i s , V 2(HS) i s c o n s i d e r e d the c o e f f i c i e n t of H. I t i s u s u a l ( [ 1 2 ] , p. 195) t o assume t h a t inhomogeneous terms i n v o l v i n g the dependent v a r i a b l e a r e known. (5) becomes (6) Rt = C ( ( H S ) V 2 H + 2 ( ( H S ) X H K + ( H S ) y H y ) ) + f, where f, = HV 2(HS) +We. The e q u a t i o n f o r S a l s o c o n t a i n s d e r i v a t i v e s of ( H 2 S ) , but s i n c e we are s o l v i n g f o r S, t h e i r e x p a n s i o n must be d i f f e r e n t than t h a t o u t l i n e d f o r H. T r e a t i n g H 2 as a s i n g l e e n t i t y was found t o work b e s t , and t h u s , f o r ( H 2 ^ ) ^ , the a p p r o p r i a t e e x p a n s i o n i s H 2 S X X + S ( H 2 ) X ; < + 2 S X ( H 2 ) ( X Expanding e q u a t i o n (2) g i v e s (7) = C/H ( k S ( H 2 V 2 S + 2VH 2-VS + SV 2H 2) + f S J ( S , H 2 ) + k ( S V S - V ( H 2 ) + H 2 ( V S ) 2 ) - ( S /H)H t J u s t as f o r H, terms not i n v o l v i n g d e r i v a t i v e s of S a r e t r e a t e d as f o r c i n g f u n c t i o n s , and (7) becomes (8) S t = C/H ( k S ( H 2 V 2 S + 2VH 2* VS) + f S J ( S , H 2 ) + k(SVS«V(H 2)+H 2(VS) 2) + f 2 where f 2 = (CkS/H) ( S V 2 H 2 ) -(S/H)H < ;. Whil e the method of expanding the d e r i v a t i v e s has been d e t e r m i n e d , the c o e f f i c i e n t s of the d e r i v a t i v e s t h a t r e s u l t a r e o b v i o u s l y f u n c t i o n s of the dependent v a r i a b l e s ; both e q u a t i o n s a r e n o n l i n e a r . D i s c r e t i z a t i o n of such e q u a t i o n s produces n o n l i n e a r systems of e q u a t i o n s which would r e q u i r e some k i n d of i t e r a t i v e scheme f o r s o l u t i o n , such as Newton's method. S i n c e 18 i t e r a t i v e methods for n o n l i n e a r e q u a t i o n s are in g e n e r a l c o m p u t a t i o n a l l y very e x p e n s i v e , we approximate f u r t h e r by l i n e a r i z i n g the e q u a t i o n s . A l though we l eave the c o m p u t a t i o n a l d e t a i l s u n t i l l a t e r in Chapter 2, the r e s u l t of l i n e a r i z i n g (6) and (8) are equat ions hav ing v a r i a b l e c o e f f i c i e n t s . The e q u a t i o n s can be s i m p l i f i e d to and s-t = b i s x x + b2Syy + b 3 S x + b^Sy + f 2 where aK , b̂  , i = 1 , . . , 4 , are f u n c t i o n s of x, y and t . In terms of the l i n e a r o p e r a t o r s d e f i n e d e a r l i e r , the CN d i s c r e t i z a t i o n i s (9) Hf* - H-5 = At(a^H^+ a 2 D 2 y H ? ; Y - + a 3 D ^ H ^ + a a D o y H ^ y H f , ) and (10) S?f ~ S% = A t ( b , D i S ^ + b 2 D J y S i ^ + b 3 E L x S ( ^ + b « D o y S ^ + f 2 ) . I t shou ld be noted here that e v a l u a t i n g the c o e f f i c i e n t s of each d e r i v a t i v e of the e q u a t i o n s shows t h a t , f or r e p r e s e n t a t i v e v a l u e s of H and S, and u s i n g the c o n s t a n t s d e f i n e d in Chapter 1, both equat ions are s t r o n g l y d i f f u s i v e . T h i s i s important in d e t e r m i n i n g the d i f f e r e n c e scheme because i t i s w e l l known ( [ 1 2 ] , p . 26) tha t c e n t r a l d i f f e r e n c e s are poor a p p r o x i m a t i o n s to c o n v e c t i v e terms i f t h e i r c o e f f i c i e n t s are l a r g e r e l a t i v e to those of the d i f f u s i o n terms . Whi le the c o n v e c t i v e terms are not i n s i g n i f i c a n t , the p r i m a r y moving f o r c e i s produced by the 19 d i f f u s i o n terms. 3 . 3 C o n s i s t e n c y , S t a b i l i t y And Convergence F i n i t e d i f f e r e n c e e q u a t i o n s are n a t u r a l l y o n l y a p p r o x i m a t i o n s t o d i f f e r e n t i a l e q u a t i o n s . I t i s n e c e s s a r y t o show t h a t the f i n i t e d i f f e r e n c e s o l u t i o n i s a 'good' a p p r o x i m a t i o n t o the e x a c t s o l u t i o n . B r i e f l y , t h i s means t h a t we r e q u i r e the d i f f e r e n c e e q u a t i o n t o tend t o the d i f f e r e n t i a l e q u a t i o n f o r Ax,At -> 0; t h a t the s o l u t i o n of the f i n i t e d i f f e r e n c e scheme w i l l g i v e a bounded s o l u t i o n (assuming the same f o r the exact s o l u t i o n ) ; and t h a t the f i n i t e d i f f e r e n c e s o l u t i o n w i l l approach the e x a c t s o l u t i o n as Ax,At -> 0. These p r o p e r t i e s a r e c a l l e d c o n s i s t e n c y , s t a b i l i t y and convergence. The s t a n d a r d t h e o r y f o r s t a b i l i t y and convergence of i n i t i a l v a l u e problems g e n e r a l l y assumes c o n s t a n t c o e f f i c i e n t s , and does not c o n s i d e r the e f f e c t of e i t h e r boundary c o n d i t i o n s or v a r i a b l e c o e f f i c i e n t s . We w i l l assume f o r the moment t h a t the model problems meet the s e c r i t e r i a , and c o n s i d e r l a t e r i n the c h a p t e r the e f f e c t on our a n a l y s i s of v a r i a b l e c o e f f i c i e n t s and boundary c o n d i t i o n s . To s i m p l i f y the p r e s e n t a t i o n of these c o n c e p t s f o r the CN scheme, we c o n s i d e r the one d i m e n s i o n a l v e r s i o n of the model e q u a t i o n s . What f o l l o w s i s e a s i l y extended to two d i m e n s i o n s . For u = u ( x , t ) , l e t ( 1 1 ) u t = a u x ^ + bu„. + f ( x , t ) w i t h a and b c o n s t a n t , and w i t h i n i t i a l c o n d i t i o n u(x,0) = g ( x ) . 20 Much of the a n a l y s i s t h a t f o l l o w s i s based on the n o t a t i o n and d e f i n i t i o n s of K r i e s s and O l i g e r ( [ 5 ] , p. 2 9 ) . A p p l y i n g the CN scheme t o (11) y i e l d s (1 - a A tD^ - b A t D w ) v ( x , ( n + 1 ) A t ) = (1 + a A t D ^ + b A t D C ! X ) v ( x f n A t ) + A t f ( x , n A t ) w i t h i n i t i a l c o n d i t i o n v ( x , 0 ) = h ( x ) , where v ( x , t ) i s the s o l u t i o n t o the f i n i t e d i f f e r e n c e e q u a t i o n . For any x, the d i s c r e t i z a t i o n of (11) produces l i n e a r c o m b i n a t i o n s of v ( x , ( n + l ) A t ) and v ( x , A t ) . We can s i m p l i f y t h e s e u s i n g m a t r i x n o t a t i o n ; (11) can be e x p r e s s e d as (12) Q,v(x,(n+1)At) = Q 0 v ( x , n A t ) + A t f ( x , n A t ) D e f i n i t i o n The d i f f e r e n c e scheme i s a c c u r a t e of o r d e r ( q ( , q 2 ) f o r the p a r t i c u l a r s o l u t i o n u ( x , t ) , i f t h e r e i s a c o n s t a n t c and a f u n c t i o n C ( t ) , bounded on e v e r y f i n i t e i n t e r v a l [ 0 , N A t ] , such t h a t f o r a l l s u f f i c i e n t l y s m a l l A t , Ax | |Q,u(x,t+At) - Q 0 u ( x , t ) - A t f | | < C ( t ) ( A x V + A t l * ) and | |h(x) - g (x) | | < c ( A x t l + At?*" ) where ||» || i s the 1 2 norm. We have seen t h a t the l i n e a r o p e r a t o r s , D,* and Dj^ , a p proximate 9/9x and 9 2 / 9 x 2 w i t h a c c u r a c y 0 ( ( A x ) 2 ) , and t h a t the CN a p p r o x i m a t i o n of 9/9t i s a c c u r a t e t o 0 ( ( A t ) 2 ) . Assuming t h a t h ( x ) = g ( x ) , t h a t the i n i t i a l c o n d i t i o n i s known e x a c t l y , the scheme i s a c c u r a t e of o r d e r ( 2 , 2 ) . The e r r o r i n the scheme, caused by t r u n c a t i n g the T a y l o r s e r i e s i n our a p p r o x i m a t i o n s t o the d e r i v a t i v e s , i s c a l l e d the t r u n c a t i o n e r r o r . For the CN scheme, the t r u n c a t i o n e r r o r i s At 2 1 t i m e s the a c c u r a c y , or 0 ( A t 3 + A t A x 2 ) . Note t h a t t h i s i s the l o c a l e r r o r , the e r r o r i n the s o l u t i o n from one time s t e p t o the n e x t . I f the d i f f e r e n c e e q u a t i o n converges t o the d i f f e r e n t i a l e q u a t i o n -being a p p r o x i m a t e d , the method i s s a i d t o be c o n s i s t e n t . For the CN scheme, c o n s i s t e n c y f o l l o w s d i r e c t l y from the d e f i n i t i o n of a c c u r a c y . S i n c e C ( t ) i s assumed t o be bounded, and q i , q 2 > 0, C ( t ) ( (Ax)**' + (At)* 1" ) -> 0 as At,Ax -> 0 and the scheme i s c o n s i s t e n t . W h i l e the d i f f e r e n c e s o l u t i o n a p proximates the e x a c t s o l u t i o n w i t h the o r d e r of a c c u r a c y g i v e n above, the ex a c t e r r o r depends on the s i z e of the p a r t i a l d e r i v a t i v e s m u l t i p l y i n g the ( A x ) 2 and ( A t ) 2 terms i n the T a y l o r s e r i e s e x p a n s i o n . In p r a c t i c e , t h e s e of c o u r s e a re unknown. To m i n i m i z e the e r r o r , s e v e r a l s t e p s a r e t a k e n when a c t u a l l y f i n d i n g a n u m e r i c a l s o l u t i o n . F i r s t , the s i z e of the time s t e p i s chosen t o be about the same s i z e as the space s t e p ; thus ( A t ) 3 w i l l be r o u g h l y e q u a l t o A t ( A x ) 2 . Second, i t i s c o n v e n t i o n a l t o n o n d i m e n s i o n a l i z e the e q u a t i o n s . T h i s r e q u i r e s the e q u a t i o n s t o be s c a l e d so t h a t 0 < x,y < 1 . The time s c a l e now depends on the s p a t i a l r e s o l u t i o n used. S t a b i l i t y r e q u i r e s t h a t the s o l u t i o n t o the d i f f e r e n c e e q u a t i o n , v ( x , t ) , at t=nAt, i s bounded as At,Ax -> 0. Def i n i t i o n The d i f f e r e n c e a p p r o x i m a t i o n ( 1 1 ) i s s t a b l e f o r a sequence A t , Ax -> 0 i f t h e r e a re c o n s t a n t s a s , Kj such t h a t f o r t ^ 0, the e s t i m a t e 22 | | v ( x , t ) | | = | | S ( t , 0 ) v ( x , 0 ) | | < K s e x p ( a s t ) | | v( x , 0 ) | | h o l d s . Here S ( t , 0 ) i s a s o l u t i o n o p e r a t o r w i t h S ( 0,0)=I, S ( t 2 , 0 ) = S ( t 2 , t 1 ) S ( t l , 0 ) ; t 2>t,>0. For t=nAt, (12) i s (13) Q,v(x,(n+1)At) = Q 0 v ( x , n A t ) + A f ( x , n A t ) and S(nAt,0) can e a s i l y be found. For Q = Q7'Q 0, v( x , ( n + 1 ) A t ) = Qv(x,nAt) + Q f 1 A t f ( x , n A t ) I t i s e v i d e n t t h a t the r i g h t s i d e of the s t a b i l i t y c o n d i t i o n i n the d e f i n i t i o n above s h o u l d i n c l u d e the c o n t r i b u t i o n from the f o r c i n g f u n c t i o n , f . That i s , the bound s h o u l d i n c l u d e upper l i m i t s on A t ^ S (nAt, i A t ) f (x , i At) . However, t h i s o n l y i n c r e a s e s the bound, and f o r s i m p l i c i t y , we l e t the d e f i n i t i o n s t a n d as i s . I t i s c l e a r t h a t we have t o show t h a t | |Q| I**-1 i s bounded f o r s t a b i l i t y t o h o l d . Assume f i r s t t h a t a F o u r i e r t r a n s f o r m of v ( x , t ) e x i s t s . For e v e r y f r e q u e n c y , w, we must show t h a t the a m p l i f i c a t i o n f a c t o r i s not g r e a t e r than one. That i s , s t a b i l i t y r e q u i r e s t h a t from one time s t e p t o the next the a m p l i t u d e of every f r e q u e n c y does not i n c r e a s e , e x c e p t t o a l l o w t r u l y i n c r e a s i n g s o l u t i o n s . T h i s i s the von Neumann n e c e s s a r y c o n d i t i o n f o r = Q(Qv(x,(n-1)At)+Q f 1 A t f ( x , ( n - 1 ) A t ) + Q ~1 A t f (x , nAt) 23 s t a b i l i t y ( [ 1 2 ] , p . 70) . A p p l y i n g F o u r i e r t r a n s f o r m s to (13) r e q u i r e s on ly the t r a n s f o r m of the l i n e a r o p e r a t o r s D0)l and D e 2 . As shown in ( [ 9 ] , p . 37) , they are D P X = i s in( /3 /Ax) (14) D.2, = -4 s i n 2 ( / 3 / 2 ) / ( A x ) 2 , and G , the a m p l i f i c a t i o n f a c t o r , i s (from [ 8 ] , p . 196), 1 - 2 a s i n 2 0 + i (b/2) ( a A t / a ) sin(2/3) + A t c f G = 1 + 2 a s i n 2 0 + i ( b / 2 ) ( a A t / a ) sin(20) w i t h a = a A t / ( A x ) 2 and j3 = kAx. S i m p l i f y i n g G l e a d s to an e q u a t i o n wi th (1 + 2 a s i n 2 / 3 ) 2 + ( b 2 / 4 ) a A t / a s i n 2 ( 2 0 ) in the denominator; i t i s c l e a r l y g r e a t e r than 1 for any /3. So G < (1 - 2asin 2 /3 + i ( b / 2 ) ( a A t / a ) s i n (20) + A t c f ) (1 + 2asin 2 /3 - i b / 2 ( a A t / a ) sin2/3 ) In modulus, the a m p l i f i c a t i o n f a c t o r can be g i v e n by (15) | G | < ( l + h 1 ( / 3 ) A t + h 2 d 3 ) ( A t ) 2 + h 3 ( / 3 ) ( A t ) 3 + h « ( 0 ) ( A t ) M = 1 + O ( A t ) , where h , ( |3 ) , h 2 ( /3 ) , h 3(|3) and h a(j3) are bounded f u n c t i o n s . The O(At) term a l l o w s e x p o n e n t i a l growth that may e x i s t in the t r u e s o l u t i o n because of the f o r c i n g f u n c t i o n . The i n e q u a l i t y above s a t i s f i e s the von Neumann c o n d i t i o n , which i s necessary for s t a b i l i t y ; i t i s known ( [ 1 2 ] , p . 72) tha t for a l l t w o - l e v e l d i f f e r e n c e schemes for s c a l a r e q u a t i o n s , the von Neumann c o n d i t i o n i s s u f f i c i e n t as w e l l . T h i s shows that the s o l u t i o n from one time l e v e l to the 24 next i s bounded. | | S ( t , 0 ) v ( x , 0 ) | | must t h e r e f o r e a l s o be bounded, and the CN scheme meets the r e q u i r e m e n t s of the s t a b i l i t y d e f i n i t i o n . I t i s c l e a r t h a t |G| i s bounded f o r any v a l u e s of of At and Ax, and the scheme i s s a i d t o be u n c o n d i t i o n a l l y s t a b l e . I f the d i f f e r e n c e between the s o l u t i o n s t o the d i f f e r e n t i a l and d i f f e r e n c e e q u a t i o n s tends t o z e r o as A t , Ax -> 0, a d i f f e r e n c e scheme i s s a i d t o be c o n v e r g e n t . From [ 5 ] , p. 31, we have the f o l l o w i n g : Theorem Assume t h a t f o r f i x e d At and Ax the d e f i n i t i o n s f o r a c c u r a c y and s t a b i l i t y a r e v a l i d . Then f o r a l l t = nAt, n > 1, | | u ( x , t ) - v ( x , t ) | | < e x p ( a 5 t ) (Ax*' +At*1 )c . S i n c e i t was shown, f o r the CN scheme, t h a t (q,,q 2)=(2,2) and t h a t s t a b i l i t y i s u n c o n d i t i o n a l , t h e r e i s no r e s t r i c t i o n on the way A t , Ax -> 0. T h e r e f o r e the r i g h t s i d e tends t o z e r o as A t , Ax -> 0, and convergence i s p r o v e n . 3.4 V a r i a b l e C o e f f i c i e n t s As was mentioned, the a n a l y s i s above assumes a d i f f e r e n t i a l e q u a t i o n w i t h c o n s t a n t c o e f f i c i e n t s , a c o n d i t i o n the model e q u a t i o n s c l e a r l y do not meet. While the a c c u r a c y i s u n a f f e c t e d by v a r i a b l e c o e f f i c i e n t s , assuming of c o u r s e t h a t t h e r e are no e r r o r s i n t h e i r c a l c u l a t i o n , i t i s q u i t e p o s s i b l e t h a t the s t a b i l i t y c o n d i t i o n may change. One method of e n s u r i n g s t a b i l i t y i s t o l o o k a t the problem i n terms of ' l o c a l s t a b i l i t y ' . That i s , the c o e f f i c i e n t s are e v a l u a t e d a t each g r i d p o i n t i n the domain and t e s t e d s e p a r a t e l y 25 f o r s t a b i l i t y . W h i l e t h i s does not guarantee convergence, i t i s u s u a l t h a t i n s t a b l i t i e s i n the approximate s o l u t i o n s t a r t l o c a l l y ( [ 1 2 ] , p. 9 1 ) . However, t o a p p l y von Neumann's method t o a l a r g e number of p o i n t s p r e s e n t s o b v i o u s d i f f i c u l t i e s . One p r a c t i c a l way t o t e s t f o r s t a b i l i t y of v a r i a b l e c o e f f i c i e n t problems i s t o p e r f o r m computer runs of sample problems. T h i s was c a r r i e d out f o r s e v e r a l problems, and the r e s u l t s are g i v e n i n C hapter 3. 3.5 S t a b i l i t y Of Boundary C o n d i t i o n s To determine the e f f e c t of the boundary c o n d i t i o n s on s t a b i l i t y , we must examine the system of l i n e a r e q u a t i o n s produced by the CN scheme f o r the i n i t i a l v a l u e problem, and check t h a t the i n t r o d u c t i o n of boundary c o n d i t i o n s do not add t o the m a t r i x Q any e i g e n v a l u e s z w i t h |z| > 1. The s t r a t e g y i s t o assume t h a t such an e i g e n v a l u e z i s i n t r o d u c e d by the boundary c o n d i t i o n , and show t h a t t h i s l e a d s t o a c o n t r a d i c t i o n . In p a r t i c u l a r , we want t o show t h a t t h i s o c c u r s f o r the boundary c o n d i t i o n s of the model e q u a t i o n s . As s p e c i f i e d i n Chapter 1, the boundary c o n d i t i o n s f o r H are p r e s c r i b e d on a l l b o u n d a r i e s . The boundary c o n d i t i o n s f o r S are p r e s c r i b e d a t the i n f l o w r e g i o n s , and the n o - f l o w c o n d i t i o n e l s e w h e r e . Assuming smoothly changing f u n c t i o n s of H and S (and thus no boundary l a y e r p r o b l e m s ) , p r e s c r i b e d boundary c o n d i t i o n s do not change the e i g e n v a l u e s of the c o r r e s p o n d i n g i n i t i a l v a l u e problem, and so o n l y the n o-flow c o n d i t i o n w i l l be a n a l y s e d . We c o n s i d e r the sample e q u a t i o n ( 1 2 ) , w i t h no f o r c i n g 26 f u n c t i o n , and 0 ^ x < 1. I t i s known t h a t f o r p a r a b o l i c problems, s t a b i l i t y of the l e f t boundary c o n d i t i o n i m p l i e s s t a b i l i t y f o r both b o u n d a r i e s [ 1 8 ] , and so o n l y the e f f e c t of the v a l u e a t u ( 0 , t ) w i l l be c o n s i d e r e d . In terms of u ( x , t ) , the second ord e r a c c u r a t e a p p r o x i m a t i o n t o the no - f l o w c o n d i t i o n a t the l e f t boundary i s u(0,nAt) = (4/3)u(Ax,nAt) - (1/3)u(2Ax,nAt) We have seen t h a t the a m p l i f i c a t i o n m a t r i x , Q = Q p Q o r must have e i g e n v a l u e s not g r e a t e r than one f o r v ( x , ( n + l ) A t ) t o be a bounded s o l u t i o n . L e t g = (g*,...,g«) be an e i g e n v e c t o r of Q, so t h a t Qg = zg or zQi9 = Qo9 w i t h z an e i g e n v a l u e . Q, and Q 0 can be e x p r e s s e d as sums of t h e i r c o e f f i c i e n t s , so i n terms of g and z, z TL c £ g; = X.a; g; or ( 16) jL(zCj - a t )gj_ = 0 . From [ 1 8 ] , (16) i s a l i n e a r d i f f e r e n c e e q u a t i o n w i t h c o n s t a n t c o e f f i c i e n t s , and the s o l u t i o n t a k e s the form, 9 = £aj(T.(z)) where the Tj (z) a r e the r o o t s of ( 1 6 ) . . S u b s t i t u t i n g t h i s i n t o (16) g i v e s (17) T 2 - 2/3T + T / T J = 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 = A t / ( A x ) 2 For |z| > 1, (17) can be shown t o have r o o t s T, and T 2 w i t h | T , ( z ) | < 1, | T 2 ( z ) | > 1. So (18) qi = a,T, (z) and, t o be an e i g e n v e c t o r , gj must a l s o s a t i s f y the no-flow boundary c o n d i t i o n . Thus, we r e q u i r e t h a t g e = ( 4 / 3 ) g , - d / 3 ) g 2 . S u b s t i t u t i n g t h i s e x p r e s s i o n i n t o (18) g i v e s (19) a,( 1 - ( 4 / 3 ) T , ( z ) + ( l ' / 3 ) T ^ z ) ) = 0 The boundary c o n d i t i o n w i l l produce an u n s t a b l e s o l u t i o n i f , f o r z, w i t h |z|>1, an e i g e n v a l u e , the c o r r e s p o n d i n g r o o t of (17) i s a l s o a r o o t of ( 1 9 ) . The r o o t s of ( 1 9 ) , however, are 1 and 3, c o n t r a d i c t i n g our f i n d i n g t h a t one of the r o o t s of (17) i s l e s s than one. Thus z, w i t h |z|>1, i s not an e i g e n v a l u e ; the no-flow boundary, approximated t o second o r d e r a c c u r a c y , does not a f f e c t s t a b i l i t y . 3.6 N o n l i n e a r Terms To a v o i d h a v i n g t o s o l v e our n o n l i n e a r e q u a t i o n s d i r e c t l y , i t i s n e c e s s a r y to l i n e a r i z e the c o e f f i c i e n t s of the d e r i v a t i v e terms. L i n e a r e q u a t i o n s a r e , of c o u r s e , much e a s i e r t o s o l v e , but the l i n e a r i z a t i o n p r o c e d u r e y i e l d s o n l y an a p p r o x i m a t i o n t o the t r u e c o e f f i c i e n t s . To examine t h i s e f f e c t on the a c c u r a c y of the system, we note f i r s t the e q u a t i o n s a r e q u a s i l i n e a r , s i n c e the c o e f f i c i e n t s of the second d e r i v a t i v e terms i n v o l v e o n l y f i r s t d e r i v a t i v e s or the f u n c t i o n i t s e l f . A l t h o u g h the CN d i f f e r e n c e scheme can be 28 m o d i f i e d t o produce p r e d i c t o r - c o r r e c t o r formulae [ 8 ] , the t r u n c a t i o n e r r o r e s t i m a t e i s 0(AtAx 2+A\? f x-) r a t h e r than 0 ( A t A x 2 + A t 3 ) which the CN scheme ( f o r l i n e a r e q u a t i o n s ) y i e l d s . Lees [8] has shown t h a t by e s t i m a t i n g the q u a s i l i n e a r c o e f f i c i e n t s u s i n g e x t r a p o l a t e d v a l u e s from the two p r e v i o u s time s t e p s , the t r u c a t i o n e r r o r remains 0 ( A t A x 2 + A t 3 ) . S i n c e the CN scheme i s c e n t r e d about t=n+l/2, the c o e f f i c i e n t s of each term must be c a l c u l a t e d a t t h a t t i m e . In p a r t i c u l a r , i f the HS terms i n the c o e f f i c i e n t s of the second o r d e r p a r t i a l d e r i v a t i v e s of (9) a r e approximated by US"*1'7- = 1.5 HS w - 0.5 HS*",' an e r r o r of 0 ( ( A t ) 3 ) i s added to t o the t r u n c a t i o n e r r o r , which thus remains 0 ( A t A x 2 + A t 3 ) . We see t h a t the c o e f f i c i e n t s of 3H/3x and 3H/3y, which i n v o l v e 3(HS)/3x and 3(HS)/3y terms, can be appro x i m a t e d i n the same way. For (HS) a t p o i n t ( i , j ) , t he d e r i v a t i v e a t time n+1/2 can be appro x i m a t e d by (HS) X = ( ( H S ) j ^ - (HS)^- ) / (2Ax) and (HS) y = ( (HS)c;-„ - (HS) l i > ( ) / (2Ay) where a l l t h e (HS) v a l u e s a re c a l c u l a t e d a t t=n+l/2 as s p e c i f i e d above. T h i s method can a l s o be a p p l i e d t o the S% and Sy2 terms t h a t appear i n the e q u a t i o n f o r S. One of each of the d e r i v a t i v e s can be t r e a t e d as the c o e f f i c i e n t of the o t h e r , and approximated by e x t r a p o l a t i n g t o time n+1/2 ( [ 2 ] , p. 141). E q u a t i o n (10) has q u a s i l i n e a r c o e f f i c i e n t s which a r e f u n c t i o n s of SH 2, H 2, H x, Hy, S x and Sy. A l l of the s e a re e s t i m a t e d i n the same f a s h i o n as o u t l i n e d above, w i t h o u t a l t e r i n g the o r d e r 29 of the t r u n c a t i o n e r r o r . The e f f e c t on s t a b i l i t y of a p p l y i n g t h i s e x t r a p o l a t i o n procedure t o a l l the terms i n each e q u a t i o n i s not known. Varah [ 1 9 ] , however, has shown t h a t f o r an e q u a t i o n i n v o l v i n g n o n l i n e a r terms i n the f i r s t s p a t i a l d e r i v a t i v e of the f u n c t i o n , and the f u n c t i o n i t s e l f , the s t a b i l i t y of the CN scheme w i t h e x t r a p o l a t i o n depends on the r e l a t i v e s i z e of the c o e f f i c i e n t s of the d i f f u s i v e and c o n v e c t i v e terms. For e q u a t i o n s which are not h i g h l y c o n v e c t i v e , the CN scheme was found t o be c o n d i t i o n a l l y s t a b l e . S i n c e both model e q u a t i o n s a r e h i g h l y d i f f u s i v e , i t seems l i k e l y t h a t s t a b i l i t y w i l l be m a i n t a i n e d . W h i l e we can p r o v i d e no s t a b i l i t y a n a l y s i s , the r e s u l t s of p r e l i m i n a r y computer runs ( d e s c r i b e d f u l l y i n Chapter 4) show t h a t , f o r chosen space and time s t e p s , t h i s l i n e a r i z a t i o n produces r e s u l t s which appear to be c o n d i t i o n a l l y s t a b l e . 3 . 7 S o l v i n g The E q u a t i o n s To s o l v e the model e q u a t i o n s , we note f i r s t t h a t they a r e c o u p l e d , i n t h a t the two independent v a r i a b l e s a r e p r e s e n t i n the two e q u a t i o n s . W h i l e the C r a n k - N i c o l s o n d i f f e r e n c e scheme c o u l d be a p p l i e d t o s o l v e f o r H and S s i m u l t a n e o u s l y , the o n l y new terms t o be s o l v e d d i r e c t l y r a t h e r than a p p r o x i m a t e l y are the f o r c i n g f u n c t i o n s . S i n c e they a r e n o n l i n e a r as w e l l , one of e i t h e r H or S would have t o be approximated i n any c a s e . L i t t l e would be g a i n e d from t h i s approach. I n s t e a d , assume, as was mentioned i n Chapter 1, t h a t S v a r i e s s l o w l y from one time l e v e l t o the n e x t . I f e q u a t i o n (9) 30 i s s o l v e d f i r s t , HS has t o be e s t i m a t e d (by e x t r a p o l a t i o n ) . I f S i s i n f a c t s l o w l y v a r y i n g , the e r r o r of a p p r o x i m a t i o n w i l l be s m a l l ( a t l e a s t s m a l l e r than t h a t caused by e s t i m a t i n g H 2, as s o l v i n g e q u a t i o n (10) would r e q u i r e ) . As we have seen, the a d d i t i o n t o the e r r o r i s 0 ( ( A t ) 3 ) , which, w i l l o n l y cause a s m a l l i n c r e a s e i n the o v e r a l l e r r o r of a p p r o x i m a t i n g H. Once H has been c a l c u l a t e d , the average of H^'and H* , H,wVa-, can used i n e v a l u a t i n g the c o e f f i c i e n t s of e q u a t i o n ( 1 0 ) , r a t h e r than an e x t r a p o l a t e d v a l u e . S o l v i n g f o r S A +' , i t i s p o s s i b l e t o repeat the p r o c e s s , u s i n g the new v a l u e s of H*+v and then S**' t o c a l c u l a t e new a p p r o x i m a t i o n s f o r HM"''t and S n + V l-. Each s t e p may be c o n s i d e r e d an i t e r a t i o n of a p r e d i c t o r - c o r r e c t o r method. We can l e s s e n the e f f e c t of the i n t e r a c t i o n of H and S by a l t e r n a t i n g the o r d e r of s o l u t i o n from one i t e r a t i o n t o the n e x t . In p r a c t i c e , i t i s p o s s i b l e t o c o n t i n u e i t e r a t i n g u n t i l the s o l u t i o n s from one i t e r a t i o n t o the next change by o n l y some a r b i t r a r y ( s m a l l ) amount. For smoothly changing v a l u e s of H and S, the s o l u t i o n s w i l l c onverge, i t i s hoped, i n o n l y a c o u p l e of i t e r a t i o n s . A p p l y i n g the CN d i f f e r e n c e scheme t o the model e q u a t i o n s produces two systems of l i n e a r e q u a t i o n s t h a t must be s o l v e d a t each time s t e p . I f some method of i t e r a t i o n i s used, which w i l l change the c o e f f i c i e n t s of the m a t r i c e s , each i t e r a t i o n w i l l r e q u i r e the systems t o be s o l v e d a g a i n . I t seems l i k e l y t h a t the m a j o r i t y of c o m p u t a t i o n time w i l l be spent i n s o l v i n g l i n e a r e q u a t i o n s , p a r t i c u l a r l y as we i n c r e a s e the r e s o l u t i o n t o o b t a i n b e t t e r a c c u r a c y . The consequences a r e d i s c u s s e d i n the next 3 1 c h a p t e r . 32 IV. CHAPTER 3 We now c o n s i d e r the s t r u c t u r e of the m a t r i c e s we must s o l v e , the p o s s i b l e methods of t h e i r s o l u t i o n and the c o m p u t a t i o n a l e f f o r t , i n terms of o p e r a t i o n c o u n t s , each would need. The s t r u c t u r e of the m a t r i c e s f o r e q u a t i o n s (9) and (10) of Chapter 2 depend on the l i n e a r o p e r a t o r s and the way the H ^ and S??' are o r d e r e d i n the v e c t o r of unknowns. The l i n e a r o p e r a t o r s i n v o l v e o n l y g r i d p o i n t s one g r i d s p a c i n g away i n each s p a t i a l d i r e c t i o n ; t h a t i s , t o s o l v e f o r ( i , j ) i n v o l v e s ( i + 1 , j ) , ( i - 1 , j ) , ( i , j + l ) , ( i , j - l ) as w e l l as ( i , j ) . The c o n v e n t i o n a l approach i s t o o r d e r the g r i d p o i n t s by row w i t h the r e s u l t i n g m a t r i x b e i n g banded w i t h h a l f bandwidth L+1, where L i s the number of g r i d p o i n t s i n the x - d i r e c t i o n . For our e q u a t i o n s , b o t h m a t r i c e s a r e nonsymmetric. A banded, nonsymmetric m a t r i x can be s o l v e d u s i n g e i t h e r a v a r i a n t of G a u s s i a n e l i m i n a t i o n or any number of i t e r a t i v e methods. The former has the advantage of r e l i a b i l i t y , but the l a t t e r can be c o n s i d e r a b l y f a s t e r . The r a t e a t which i t e r a t i v e methods converge, and thus the number of i t e r a t i o n s they r e q u i r e , depends on the c o n d i t i o n number of the m a t r i x . For an i l l - c o n d i t i o n e d m a t r i x i t e r a t i v e methods w i l l o f t e n converge s l o w l y , i f a t a l l . N o t i n g t h a t the c o e f f i c i e n t s of our m a t r i x a r e the r e s u l t of a p p r o x i m a t i o n s t o n o n l i n e a r terms, and vary w i t h x, y and t , t h e r e i s l i t t l e we can say about how w e l l - c o n d i t i o n e d the m a t r i c e s w i l l be. (We assume of course- t h a t they remain n o n s i n g u l a r . ) T h i s b e i n g the c a s e , G a u s s i a n e l i m i n a t i o n would appear t o be the best method. 33 To c a l c u l a t e the amount of c o m p u t a t i o n a l e f f o r t r e q u i r e d , we need o n l y l o o k a t the o r d e r and the h a l f - b a n d w i d t h of the m a t r i c e s . I f L i s the number of g r i d p o i n t s i n each d i r e c t i o n (so Ax=Ay), the o r d e r i s L 2 and the h a l f - b a n d w i d t h L+1. A l t h o u g h the m a t r i c e s i n i t i a l l y have z e r o s between the elements on the edge of the band and the d i a g o n a l s , G a u s s i a n e l i m i n a t i o n " f i l l s i n " these z e r o s . When c a l c u l a t i n g the amount of work, t h e s e elements must be c o n s i d e r e d as nonze r o s . F o r such a m a t r i x , h a v i n g o r d e r L 2 , the o p e r a t i o n count ( t h e number of m u l t i p l i c a t i o n s ) n e c e s s a r y t o reduce i t t o t r i a n g u l a r form i s 0 ( L " ) . However, d i f f e r e n c e schemes e x i s t f o r e q u a t i o n s w i t h the same form of the model e q u a t i o n s which produce m a t r i c e s which can be s o l v e d i n o n l y 0 ( L 2 ) o p e r a t i o n s . They a r e c a l l e d a l t e r n a t i n g d i r e c t i o n i m p l i c i t methods, and are the s u b j e c t of the next s e c t i o n . 4. 1 I n t r o d u c t i o n To ADI A l t e r n a t i n g d i r e c t i o n i m p l i c i t , or ADI, methods are p e r t u r b a t i o n s of d i f f e r e n c e methods f o r m u l a t e d d i r e c t l y from the d i f f e r e n c e a p p r o x i m a t i o n s t o d e r i v a t i v e s . ADI methods most o f t e n a r i s e by add i n g d i f f e r e n c e o p e r a t o r s t o a m u l t i - d i m e n s i o n a l d i f f e r e n c e scheme i n a way t h a t does not a f f e c t the a c c u r a c y of the scheme but a l l o w s an e a s i e r s o l u t i o n . To i l l u s t r a t e the d e r i v a t i o n of a p a r t i c u l a r ADI method, the CN scheme f o r u t = u ^ + u y y , i s , u s i n g the l i n e a r o p e r a t o r s i n t r o d u c e d i n Chapter 2, (1) (1 - XD02x - XD^)u*" = (1 + XD 2 X + XD 2 y ) u * 34 where X=At/2. I f the o p e r a t o r X 2D^ De2y i s added t o the d i f f e r e n c e o p e r a t o r on both s i d e s , (1) becomes (2) (1 - XD^ - XDiy + X 2D < 2 C D2„ ) u ^ ' = (1 + XD|X + XDe2y + X2D«2X D 2 y ) u * . F a c t o r i n g (2) i n terms of D£ x , D|y g i v e s (3) (1 - XD 2 X )(1 - XD 2 y )u**> = (1 + XD,2, )(1 + XD2, )u* . I f we i n t r o d u c e an i n t e r m e d i a t e time l e v e l between n and n+1, say n+*, (3) can be ' s p l i t ' i n t o two e q u a t i o n s , namely (1 - XD 2 )u*** = (1 + XD 2 ) u n (4) V ( 1 - XDC2X ) u**» = (1 + XD 2 y ) u*** . The v a r i a b l e u*** i s a dummy v a r i a b l e i n the sense t h a t the v a l u e of U"1** i s not an a p p r o x i m a t i o n t o the s o l u t i o n a t any p a r t i c u l a r t i m e . The method of s p l i t t i n g used i n (4) i s c a l l e d the Peaceman-Rachford f o r m u l a ( [ 9 ] , p. 6 0 ) . A l t h o u g h t h e r e a r e many methods f o r s p l i t t i n g a t w o - d i m e n s i o n a l p a r a b o l i c e q u a t i o n l i k e ( 4 ) , t h i s f o r m u l a w i l l h e r e a f t e r be r e f e r r e d t o as the ADI scheme. The advantage of (4) i s t h a t the s o l u t i o n a t the next time l e v e l ( e i t h e r u**' or un**') , i s found a l o n g l i n e s p a r a l l e l t o e i t h e r the x- or y - a x i s . The m a t r i x t o be s o l v e d thus o n l y i n v o l v e s p o i n t s i m m e d i a t e l y a d j a c e n t t o the p o i n t b e i n g found; t h a t i s , the m a t r i x i s t r i d i a g o n a l . The c o m p u t a t i o n a l e f f o r t r e q u i r e d t o s o l v e a t r i d i a g o n a l system of e q u a t i o n s i s s m a l l compared t o t h a t f o r the CN scheme. For an L 2 by L 2 system, the t r i d i a g o n a l m a t r i x r e q u i r e s o n l y 0 ( L 2 ) m u l t i p l i c a t i o n s , compared t o 0(1/) r e q u i r e d by the CN scheme. T h i s c o n s i d e r a b l e s a v i n g s i n com p u t a t i o n time w i l l not come 35 a t the expense of a c c u r a c y or s t a b i l i t y of the d i f f e r e n c e method. As we w i l l show, the Peaceman-Rachford scheme, l i k e the CN scheme, has a t r u n c a t i o n e r r o r of 0 ( ( A t ) 3 + A t ( A x ) 2 ) and i s u n c o n d i t i o n a l l y s t a b l e , f o r a sample problem s i m i l a r t o t h a t a n a l y z e d i n Chapter 2. Note t h a t i n c a l c u l a t i n g the t r u n c a t i o n e r r o r , we have assumed t h a t Ax=Ay. For the r e s t of the t h e s i s , t h i s w i l l be the c a s e . 4.2 Accur a c y And S t a b i l i t y To c a l c u l a t e the a c c u r a c y of the ADI scheme, we f i r s t c o n s i d e r an e q u a t i o n of the same form as our model e q u a t i o n s , namely, f o r u = u ( x , y , t ) , (5) u t = a u X x + b U y y + c u x + d u v , where a,b,c and d a r e c o n s t a n t c o e f f i c i e n t s . F a c t o r i n g (5) u s i n g the ADI scheme i n the same manner as above, g i v e s (6) (1 - XaD 2 x - XcD„x )(1 - XbDj y - XdDey ) u ^ 1 = (1 + XaD 2 x + XcD # J f )(1 + XbD 2 y + XdD 0j ) u* , where X i s as above. S p l i t t i n g (6) we get (1 - XbD 2 - XdD*v ) u M + * = (1 + XaD 2 + XcD^ ) u* (7) * * (1 - XaD^ - XcDe;t )u*«" = (1 + XbDe2y + XdD c y ) u * + * To c a l c u l a t e the e r r o r , we see f i r s t t h a t the o p e r a t o r s added t o the CN scheme f o r (5) are (\ 2cd)D o x D e y , ( X 2 b c ) D o x D < 2 y , ( X 2 a d ) D o y D2^ and ( X 2 a b ) D ^ D j ^ . These terms do not r e s u l t from T a y l o r s e r i e s e x p a n s i o n of the d e r i v a t i v e s , and the e r r o r s they produce can be added d i r e c t l y t o the t r u n c a t i o n e r r o r found f o r the CN scheme. S i n c e each o p e r a t o r i s a p p l i e d t o both u H and u"**' , the e r r o r s can e a s i l y be found by expanding u1**1 i n a 36 T a y l o r s e r i e s . Examining j u s t one of these terms, ( X 2 c d ) D D , we have ( X 2 c d ) D D (u + Atu + 0 ( ( A t ) 2 ) = ( X 2 c d ) D D u C a n c e l l i n g u on both s i d e s , we see t h a t the a d d i t i o n t o the t r u n c a t i o n e r r o r of t h i s term i s 0 ( ( A t ) 3 ) . The c o n t r i b u t i o n t o the e r r o r from the o t h e r o p e r a t o r s i s s i m i l a r l y found t o be 0 ( ( A t ) 3 ) , so the t r u n c a t i o n e r r o r f o r ADI i s of the same o r d e r as the CN scheme d i s c u s s e d e a r l i e r ; the added terms l e a v e the o v e r a l l a c c u r a c y unchanged. C o n s i s t e n c y f o l l o w s i n the same way as shown f o r CN i n the l a s t c h a p t e r . S t a b i l i t y f o r the ADI scheme can be seen by a g a i n comparing w i t h the CN scheme. L o o k i n g at the a m p l i f i c a t i o n f a c t o r , G, f o r CN a p p l i e d t o e q u a t i o n ( 5 ) , we have 1-2(a,S,+a 2S 2)+i((b/2)(a,At/a)S 3+(d/2)(a 2At/d)S„) (8) G = 1 + 2 ( a 1 S 1 + a 2 S 2 ) + i ( ( b / 2 ) ( a , A t / a ) S 3 + ( d / 2 ) U 2 A t / d ) S a ) where S , = s i n 2 / 3 , S 2 = s i n 2 7 j , S 3 = s i n 2 / 3 , S a=sin2ry and a, = a A t / ( A x ) 2 , a 2 = b A t / ( A y ) 2 , 0 = k,Ax and n = k 2Ay ( [ 1 2 ] , p. 196) . I t i s e a s i l y shown t h a t |G| < 1, f o r At,Ax->0. In c a l c u l a t i n g the a m p l i f i c a t i o n f a c t o r f o r ADI, we note t h a t the a d d i t i o n a l terms of the ADI scheme are F o u r i e r t r a n s f o r m e d to r e a l v a l u e d terms i n v o l v i n g p r o d u c t s of s i n 2 ( 0 ) , s i n 2 ( 7 ? ) , s i n ( 2 / 3 ) and s i n ( 2 r ? ) . But t h e s e terms a r e added t o both the denominator and the numerator of ( 7 ) ; the i n e q u a l i t y i s unchanged. Thus, the a m p l i f i c a t i o n f a c t o r f o r ADI i s a l s o bounded by one and ADI i s u n c o n d i t i o n a l l y s t a b l e f o r e q u a t i o n ( 5 ) . U s i n g the r e s u l t s from Chapter 2, convergence f o r ADI i s guaranteed by the a c c u r a c y and u n c o n d i t i o n a l s t a b i l i t y of the 37 scheme. 4 . 3 C o m p u t a t i o n a l D e t a i l s And Sample Problems We have assumed t h a t the c o e f f i c i e n t s of ( 5 ) are c o n s t a n t , as r e q u i r e d f o r the a n a l y s i s t o h o l d . To examine a c c u r a c y and s t a b i l i t y f o r v a r i a b l e c o e f f i c i e n t s , and t o t e s t the e f f e c t of v a r i o u s e x t r a p o l a t i o n schemes f o r n o n l i n e a r terms, a s e r i e s of e q u a t i o n s w i t h known s o l u t i o n s and of s i m i l a r form t o the model e q u a t i o n s were s o l v e d u s i n g the ADI method. The F o r t r a n code d e v e l o p e d f o r t h i s t e s t i n g was a l s o found t o be u s e f u l when c o d i n g the program f o r the model e q u a t i o n s . B e f o r e d i s c u s s i n g the sample problems, we must det e r m i n e how t o e v a l u a t e v a r i a b l e c o e f f i c i e n t s and f o r c i n g f u n c t i o n s so the ADI e q u a t i o n s remain c o n s i s t e n t . In a d d i t i o n , the r e quirement of the ADI scheme f o r boundary v a l u e s at t=n+* must be c o n s i d e r e d , s i n c e our e q u a t i o n s do not have boundary c o n d i t i o n s d e f i n e d a t i n t e r m e d i a t e time l e v e l s . To see how t o i n c o r p o r a t e v a r i a b l e c o e f f i c i e n t s i n t o the ADI scheme, c o n s i d e r e q u a t i o n ( 5 ) , w i t h a, b, c and d f u n c t i o n s of x, y and t . The Peaceman-Rachford f o r m u l a i s the same as ( 6 ) and i t i s c o n v e n t i o n a l ( [ 9 ] , p. 6 8 ) t o e v a l u a t e the c o e f f i c i e n t s a t t = n + l / 2 . T h i s i s c o n s i s t e n t w i t h the CN d i f f e r e n c e scheme, which i s c e n t r e d about t = n + l / 2 . The a c c u r a c y of the ADI scheme i s unchanged from the c o n s t a n t c o e f f i c i e n t c a s e , assuming the c o e f f i e n t s are known e x a c t l y . The a d d i t i o n of a f o r c i n g f u n c t i o n t o an e q u a t i o n poses a d i f f e r e n t problem. The f o r c i n g f u n c t i o n i t s e l f must be s p l i t i n 38 some way t o m a i n t a i n the c o n s i s t e n c y of the ADI f o r m u l a . C o n s i d e r i n g as an example, u t= u ^ + u y y + F ( x , y , t ) , we note f i r s t t h a t F ( x , y , t ) must be e v a l u a t e d a t t=n+l/2. U s i n g ( 1 - X D 2 ) u A + * = ( l + X D 2 ^ + F ( x , y , (n+l/2)At)/2 (9) r ( 1 - X D ^ ) U w H = ( 1+XD.2 )u n +*+ F ( x , y , ( n + l / 2 ) A t ) / 2 and s o l v i n g f o r u*** , we f i n d t h a t (10) u*** = (1-XD^) "MO+XD^u* + F , W V V 2 ) where F A + , / 2 - = F ( x ,y, (n +1 /2) At) . S u b s t i t u t i n g (10) i n t o the second s p l i t e q u a t i o n g i v e s (l-XD. 2 y)(l-XD^)u n H = O+XD^Ml+XD^Ju" + ((l-XDE2Y) + (l+XD < 2 Y))F^ yV2 = ( 1+XD^) ( 1+XD^)u* + F N + ^ which shows t h a t (9) i s c o n s i s t e n t . When s o l v i n g e i t h e r Peaceman-Rachford e q u a t i o n , time-dependent boundary v a l u e s a t the i n t e r m e d i a t e time l e v e l , t=n+*, are needed. However, t h i s time l e v e l does not c o r r e s p o n d t o any p a r t i c u l a r t i m e . The v a l u e s c o u l d be approximated by a v e r a g i n g v a l u e s a t t=n+1 and t=n, but the o v e r a l l a c c u r a c y of the scheme would not be m a i n t a i n e d ( [ 1 5 ] ) . I t i s p o s s i b l e t o f i n d u n + * e x p l i c i t l y i n terms of the v a l u e s a t u" and u**' ; the form u l a f o l l o w s from e q u a t i o n ( 7 ) . Adding the s p l i t e q u a t i o n s g i v e s ( ( 1 - XDE2Y- XD^) + (1 + XD^+ \L\y) ) u« +* = ( 1 - XD£- XD^) u ^ 1 + (l+XD^+XD^u' 1 which i s j u s t u*** = (( 1 - X D ^ - X D e J ( ) u n + , + (1+XD^+XEL x)U m)/2 . 39 We r e q u i r e u n t"* on the x - d i r e c t i o n b o u n d a r i e s , and i t i s e a s i l y seen t h a t the d i f f e r e n c e o p e r a t o r s i n v o l v e v a l u e s o n l y on those b o u n d a r i e s . T h i s f o r m u l a e n s u r e s t h a t c a l c u l a t i o n of the boundary v a l u e s p r e s e r v e s the o r d e r of a c c u r a c y of the method. We can now a p p l y the ADI scheme t o a s e r i e s of sample problems. The problems we have chosen are i n t e n d e d t o show t h a t ADI has the a c c u r a c y p r e d i c t e d f o r the c o n s t a n t c o e f f i c i e n t problem and t h a t , a l t h o u g h not as a c c u r a t e f o r v a r i a b l e c o e f f i c i e n t and n o n l i n e a r problems, the method appears t o be a t l e a s t c o n d i t i o n a l l y s t a b l e . The sample problems a l l have the g e n e r a l form (11) u± = a u w + b u y y + c u x + duy where the e x a c t s o l u t i o n u and the c o e f f i c i e n t s a r e chosen t o produce the v a r i o u s problems of i n t e r e s t . We use the e x a c t s o l u t i o n on the b o u n d a r i e s except a t the i n t e r m e d i a t e time l e v e l s , where the boundary c o n d i t i o n s a r e c a l c u l a t e d as e x p l a i n e d above. The r e s u l t s , c o n s i s t i n g of the maximum a b s o l u t e e r r o r a t the i n t e r i o r g r i d p o i n t s a t time = 1.0 f o r v a r i o u s v a l u e s of At and Ax, a r e g i v e n i n Table 1. For the c o n s t a n t c o e f f i c i e n t problem, we s e t a = b = 0.5 and c = d = 1.0, which has e x a c t s o l u t i o n u = e x p ( - x - y - t ) . The r e s u l t s are g i v e n i n T a b l e 1 under ( I ) . The. e r r o r s v a r y a c c o r d i n g t o the t r u n c a t i o n e r r o r of the method, which i s 0 ( ( A t ) 3 + A t ( A x ) 2 ) . We see t h a t the e r r o r s are s m a l l e r than the t r u n c a t i o n e r r o r would i n d i c a t e , which i s l i k e l y due t o c a n c e l l a t i o n of e r r o r terms, whose c o e f f i c i e n t s a r e v a r i o u s h i g h e r d e r i v a t i v e s of e x p ( - x - y - t ) . 40 To t e s t a v a r i a b l e c o e f f i c i e n t problem, (11) i s s o l v e d w i t h a=0.5/(y+.1) 2, b=0.5/(x+.1) 2, c=1/(y+.1) and d=l/(x+.1), which has t r u e s o l u t i o n u = e x p ( - ( x + . 1 ) ( y + . 1 ) - t ) . The r e s u l t s a r e g i v e n under ( I I ) and a r e much the same as the c o n s t a n t c o e f f i c i e n t c a s e . For both t h e s e cases the e r r o r becomes c o n s t a n t a f t e r the f i r s t t en t o twenty time s t e p s ( f o r some c a s e s , e r r o r s were checked f o r 200 time s t e p s ) , i n d i c a t i n g s t a b i l i t y of the method f o r the time and space s t e p s g i v e n i n the t a b l e . One of the n o n l i n e a r terms both model e q u a t i o n s c o n t a i n i s the square of the c o n v e c t i v e terms. As d e s c r i b e d i n Chapter 2, these are a p p r o x i m a t e d by d i s c r e t i z i n g one d e r i v a t i v e and t r e a t i n g the o t h e r as i t s c o e f f i c i e n t . The sample e q u a t i o n i s s o l v e d w i t h a and b as f o r the v a r i a b l e c o e f f i c i e n t problem, c=-u</(u A(y+.1)) and d = - U y / ( u y ( x + . 1 ) ) , where u*, Uy a r e a p p r o x i m a t e d by e x t r a p o l a t i o n and u x , Uy are known e x a c t l y . The e x a c t s o l u t i o n i s the same as f o r the v a r i a b l e c o e f f i c i e n t problem. Case I ' l l i n Table 1 c o n t a i n s the r e s u l t s found by a p p r o x i m a t i n g u x by D 0 < u " and Uy by D 0 y U w and case IV the r e s u l t s w i t h a p p r o x i m a t i o n s D 0x (1 • 5u n-0. 5uA'x ) and D 0y(1.5u*-0.5u r t" 1). The t w o - p o i n t e x t r a p o l a t i o n i s s l i g h t l y more a c c u r a t e , a l t h o u g h not as much b e t t e r as e x p e c t e d . For s e v e r a l runs w i t h t h i s e x t r a p o l a t i o n , however, the s o l u t i o n i s u n s t a b l e at l a t e r t i m e s , o s c i l l a t i n g and f i n a l l y b l o w i n g up. T h i s i s a p p a r e n t l y caused by a p a r a s i t i c s o l u t i o n i n t r o d u c e d n u m e r i c a l l y by l i n k i n g t h r e e time l e v e l s t o g e t h e r . (That the o n e - p o i n t e x t r a p o l a t i o n d i d not blow up would seem t o c o n f i r m t h i s . ) S i n c e 41 the s o l u t i o n w i t h At=.05, AX=.10 d i d not blow up, i t appears t h i s problem can be a v o i d e d by a c a r e f u l c h o i c e of time and space s t e p s . S i n c e the t w o - p o i n t e x t r a p o l a t i o n r e q u i r e s s o l u t i o n s from two p r e v i o u s time l e v e l s , a s o l u t i o n f o r time=At i s n e c e s s a r y . We compared the e f f e c t of a p p r o x i m a t i n g t h i s s o l u t i o n by a T a y l o r s e r i e s e x p a n s i o n w i t h u s i n g the e x a c t s o l u t i o n and found no d i f f e r e n c e i n the e r r o r s f o r the t e s t v a l u e s of Ax and A t . To t e s t the e f f e c t of n o n l i n e a r c o e f f i c i e n t s of the d i f f u s i o n terms, we s o l v e (12) s e t t i n g a = . 5 u / ( ( y + . 1 ) 2 u ) , b=.5u/((x+.1) 2u) and c and d as f o r the v a r i a b l e c o e f f i c i e n t problem. Here u = 1 . Su^-0. 5u n"' and the e x a c t s o l u t i o n i s a g a i n u = e x p ( - ( x + . 1 ) ( y + . 1 ) - t ) . The s o l u t i o n a t time=At was found a p p r o x i m a t e l y by T a y l o r s e r i e s e x p a n s i o n . The r e s u l t s are g i v e n under ( V ) . A l l the examples except At=.05, AX=.10 blow up a t l a t e r t i m e s ; the scheme i s o b v i o u s l y no l o n g e r u n c o n d i t i o n a l l y s t a b l e . As f o r the p r e v i o u s example, t h i s may be due t o the presence of p a r a s i t i c s o l u t i o n s . To f i n d the convergence r a t e f o r t h i s problem, s e v e r a l o t h e r runs were made w i t h At and Ax chosen t o produce s t a b l e r e s u l t s . These s o l u t i o n s a r e not o n l y l e s s a c c u r a t e than the c o r r e s p o n d i n g v a r i a b l e c o e f f i c i e n t problem, but do not converge a t the r a t e we e x p e c t ; the e r r o r s changes o n l y as A t . The reasons f o r t h i s a r e not c l e a r but seem t o be due t o the form of the n o n l i n e a r problem. We c o n c l u d e from t h e s e sample problems t h a t the ADI scheme and the methods chosen f o r h a n d l i n g n o n l i n e a r terms produce 42 a c c e p t a b l e s o l u t i o n s . The l a r g e r e r r o r s and c o n d i t i o n a l s t a b i l i t y f o r the n o n l i n e a r e q u a t i o n s a re due i n p a r t t o the s p e c i f i c s o l u t i o n and the f a c t t h a t the e q u a t i o n s have c o e f f i c i e n t s which c o n t a i n e x a c t s o l u t i o n s or e x a c t d e r i v a t i v e s . The e r r o r t h a t r e s u l t s from the a p p r o x i m a t i o n s to the s e terms w i l l always be i n the same d i r e c t i o n , which i s c e r t a i n l y not n e c e s s a r i l y the case f o r t r u l y n o n l i n e a r e q u a t i o n s . I t s h o u l d be noted as w e l l t h a t the e r r o r produced by e x t r a p o l a t i n g the c o e f f i c i e n t s of the c o n v e c t i v e terms i s of l i t t l e c o n c e r n s i n c e the model e q u a t i o n s a re h i g h l y d i f f u s i v e ; the c o n v e c t i v e terms w i l l have o n l y a s m a l l c o n t r i b u t i o n t o the o v e r a l l s o l u t i o n . At Ax ( I ) ( I I ) ( I I I ) (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. R e s u l t s of ADI s o l u t i o n s of sample problems. 4.4 A p p l y i n g ADI To The Model E q u a t i o n s The f i n a l form of the ADI scheme f o r the model e q u a t i o n s i s found by a p p l y i n g the p r o c e d u r e s d i s c u s s e d i n t h i s c h a p t e r t o e q u a t i o n s (9) and (10) of Chapter 2. The s p l i t e q u a t i o n s f o r (9) a re 43 ( l - X a 2 D 2 -Xa 4D~ JH"** = ( 1 +Xa , D 2 +Xa ̂  ) H n + A t f , / 2 (13) y ^ ^ ( l - X a ^ -Xa 3D^ )U**[ = (l+Xa 2D 0 2 y + Xa4D<,y )H n 1+*Atf, /2 and f o r (10) a r e O-XbzD, 2 -Xb„D 0 y )S*+* = ( l + X b ^ 2 + Xb3D01. ) S* + A t f 2 / 2 (14) ^ y * ( 1 - X b ^ -Xb 3D e > t ) S * M = ( 1+Xb 2D^ + Xb 1 )D^ J S * ^ A t f 2/2 where the c o e f f i c i e n t s , f , and f 2 a r e d e f i n e d as i n Chapter 2 and e v a l u a t e d at time=n+1/2. 44 V. CHAPTER 4 In the f i r s t p a r t of t h i s c h a p t e r , s e v e r a l p r e l i m i n a r y n u m e r i c a l t e s t s of the model e q u a t i o n s a r e d i s c u s s e d . They i n v o l v e s i m p l i f y i n g the boundary c o n d i t i o n s t o produce a c l o s e d system, a l l o w i n g us t o v e r i f y t h a t the approximate s o l u t i o n s are c o n s e r v i n g . We a l s o d i s c u s s the c o m p u t a t i o n a l d e t a i l s of both s e t s of i n i t i a l and boundary c o n d i t i o n s , n e i t h e r of which, i t was found, c o u l d be a p p l i e d i n a s t r a i g h t f o r w a r d manner. The next s e c t i o n i s d e v o t e d t o t e s t i n g the model e q u a t i o n s w i t h these i n i t i a l and boundary c o n d i t i o n s u s i n g v a r i o u s time s e r i e s of i n p u t v a l u e s . A . p h y s i c a l i n t e r p r e t a t i o n of the p l o t t e d r e s u l t s f o l l o w s . We c o n c l u d e w i t h recommendations f o r c o n t i n u e d work on the model, both w i t h s t e p s which might be taken t o more a c c u r a t e l y r e p r e s e n t the p h y s i c a l system b e i n g m o d e l l e d and w i t h improvements t o the n u m e r i c a l methods used. 5.1 P r e l i m i n a r y T e s t s The purpose of t h e s e t e s t s i s t o ensure t h a t a p p l y i n g the ADI scheme and the n o n l i n e a r a p p r o x i m a t i o n s t o the model e q u a t i o n s c o n s e r v e mass, as measured by H, and t o t a l s a l i n i t y , as measured by the p r o d u c t of H and S. We c o n s i d e r the e q u a t i o n s over the s p e c i f i e d domain, but impose c l o s e d b o u n d a r i e s . By m a i n t a i n i n g c o n s t a n t v a l u e s of H and S on a l l b o u n d a r i e s , and s e t t i n g the v e l o c i t i e s u and v t o z e r o on the y- and x - d i r e c t i o n b o u n d a r i e s r e s p e c t i v e l y , i t i s c l e a r t h a t the mass and t o t a l s a l i n i t y i n the system s h o u l d remain c o n s t a n t . For i n i t i a l c o n d i t i o n s we use c o n s t a n t v a l u e s of H and S 45 everywhere except a t a s i n g l e g r i d p o i n t . At t h i s p o i n t , which i s i n the c e n t e r of the domain, H or S or both i s s e t t o a v a l u e l a r g e r than the c o n s t a n t f i e l d . Three runs were made w i t h c o n s t a n t v a l u e s of 20 meters f o r H and 2 p . p . t . f o r S, and s t e p s i z e s of Ax=Ay=At=0.1. The r e s u l t s a re g i v e n i n T a b l e 2 below. Mass and HS a r e c a l c u l a t e d s i m p l y by summing over the g r i d p o i n t s . 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. R e s u l t s of C o n s e r v a t i o n T e s t s The method appears t o c o n s e r v e both q u a n t i t i e s v e r y w e l l , w i t h an e r r o r over 100 time s t e p s of l e s s than .1 p e r c e n t f o r the f i r s t c a s e , and s l i g h t l y l a r g e r f o r the o t h e r two. For the crude 11 by 11 g r i d used, t h i s i s v e r y good. The r e s u l t s of the second and t h i r d t e s t s a l s o i n d i c a t e , not s u r p r i s i n g l y , t h a t the l a r g e r the p o i n t s o u r c e s , the l e s s a c c u r a t e the s p a t i a l d e r i v a t i v e s a r e approx i m a t e d , and the l a r g e r the o v e r a l l e r r o r i n terms of c o n s e r v a t i o n . 46 5.2 I n i t i a l C o n d i t i o n s At The Land Boundary For purposes of t e s t i n g , c o n s t a n t f i e l d v a l u e s were a s s i g n e d as i n i t i a l c o n d i t i o n s t o both the h e i g h t of the upper l a y e r , H, and the s a l i n i t y d e f e c t , S. The v a l u e s chosen were 20 meters f o r H and 2 ppt f o r S, except a t the r i v e r d i s c h a r g e p o i n t . Here the v a l u e s we a s s i g n e d t o H and S depended on the time s e r i e s we were t e s t i n g . To t e s t the problem c r e a t e d by a v a l u e of S a t the d i s c h a r g e p o i n t l a r g e r than those nearby, we used a weekly s e r i e s w i t h S e q u a l t o f i v e . The n e a r l y - d i s c o n t i n u o u s n a t u r e of the f u n c t i o n a t the i n i t i a l time l e v e l r e s u l t e d i n e x t r e m e l y u n l i k e l y v a l u e s of the s o l u t i o n s a t time At f o r both H and S a t g r i d p o i n t s near the d i s c h a r g e . These r e s u l t s p e r s i s t e d d e s p i t e s e v e r a l i t e r a t i o n s , and t h e r e was no i n d i c a t i o n t h a t the s o l u t i o n s f o r H and S would converge. The problem stems from h a v i n g t o f i n d the s o l u t i o n of e q u a t i o n ( 1 4 ) , from Chapter 1, which i s (1) H y = - ( f H S y + 2kSH x + kHS^)/2fS . Due t o the l a r g e g r a d i e n t s of S c r e a t e d by a l a r g e i n p u t v a l u e , (1) produces l a r g e v a l u e s of H. To a v o i d the e f f e c t of t h i s n e a r - d i s c o n t i n u i t y , we assume [7] t h a t S a t the d i s c h a r g e p o i n t has an immediate e f f e c t on p o i n t s i n the r e g i o n ; t h a t i s , the impact of l a r g e changes i n s a l i n i t y d e f e c t i s smoothed over the r e g i o n . We a r b i t r a r i l y chose f i v e g r i d p o i n t s i n e i t h e r d i r e c t i o n a l o n g the l a n d boundary, two g r i d p o i n t s i n t o the domain, and v a r i e d S l i n e a r l y over these p o i n t s . We a l s o want t o f i n d a v a l u e of H at the d i s c h a r g e p o i n t 47 such tha t H 2 S i s i n i t i a l l y in e q u i l i b r i u m wi th nearby va lues of H 2 S . To f i n d such an H , we use the f o l l o w i n g procedure [ 7 ] , I t i s assumed that the d i s c h a r g e r e g i o n c o n s i s t s of o n l y one g r i d p o i n t , but i n c l u d e s the r e g i o n up to but not i n c l u d i n g the g r i d p o i n t s immediate ly n o r t h a n d - s o u t h . The l and boundary c o n d i t i o n for U i s U / C = - ( f ( H 2 S ) y + k ( H 2 S ) x ) , w i th U g i v e n at the d i s c h a r g e p o i n t and U = 0 e l s ewhere . The g r a d i e n t of H 2 S a c r o s s the r e g i o n i s n e g a t i v e , due to the C o r i o l i s f o r c e , and we approximate i t u s i n g ( H 2 S L . = ( H 2 S L - 6 (H 2 S) (2) ( H 2 S ) K r J = (H 2 S)< + 5 (H 2 S) where the d i s c h a r g e p o i n t (LAx,KAx) i s denoted by the s u b s c r i p t K, the f i r s t g r i d p o i n t to the n o r t h by K+1, and the f i r s t to the south by K - 1 . 5 (H 2 S) i s the i n c r e a s e in H 2 S from p o i n t K to K - 1 . We see from (2) that a second order a p p r o x i m a t i o n to ( H 2 S ) y i s ( H 2 S ) y = 5 ( H 2 S ) / A x . To f i n d 8 ( H 2 S ) , we note tha t U = 0 at a l l g r i d p o i n t s except p o i n t K, and so - f ( H 2 S ) y = k ( H 2 S ) x at these p o i n t s . Denot ing the i n i t i a l cons tant f i e l d va lue by ( H 2 S ) © « , we can approximate - f ( H 2 S ) on the boundary n o r t h of the d i s c h a r g e p o i n t by -f((H2S)<*> - ( H 2 S ) ^ , ) /5Ax . I f we assume f u r t h e r that ( H 2 S ) X has the same v a l u e at the d i s c h a r g e p o i n t and at the f i v e g r i d p o i n t s n o r t h and south , we 48 have k ( H 2 S ) x = - f ( ( H 2 S ) e o " ( H 2 S ) ( C + | )/5Ax and so a t the d i s c h a r g e p o i n t , -U/C = f 6 ( H 2 S ) / A x - f ( ( H 2 S ) c o - ( H 2 S ) | C + ) )/5Ax or (3) -UAx/(Cf) = 1.28 (H 2S) - 0.2 ( (H 2S ) oo " ( H 2 S ) K ) . We a l s o r e q u i r e t h a t t he g r i d p o i n t i n s i d e the domain a d j a c e n t t o the d i s c h a r g e p o i n t be e q u a l t o (H 2S )o© . S i n c e t h i s p o i n t i s ( ( L - 1 ) A x , K A x ) , H 2 S can be expanded i n a T a y l o r s e r i e s i n t he x - d i r e c t i o n about LAx, which i s , t o O(Ax), H 2 S ((L-1)Ax,KAx) = H 2 S(LAx,KAx) - Ax(H 2 S ( L A x , K A x ) ) or (4) (H 2 S )oo = ( H 2 S V + ( f / 5 k ) ( (H 2 S )o* " ( H 2 S ) I C - 6 ( H 2 S ) ) . We know S everywhere a l o n g the boundary, U a t the d i s c h a r g e p o i n t , and ( H 2 S ) o e , our i n i t i a l c o n d i t i o n , so the o n l y unknowns i n e q u a t i o n s (3) and (4) a r e 6 (H 2 S ) and ( H 2 S ) K . Two e q u a t i o n s w i t h two unknowns are e a s i l y s o l v e d , and once 6 (H 2 S ) and (H 2 S )^ are known, H(LAx,(K+1)Ax) and H (LAx,(K-1)Ax) can be deter m i n e d from ( 1 ) . H(LAx,(K±i)Ax), i=1,...,5 , a r e found by assuming H 2 S i s l i n e a r from (H 2 S)(LAx,(K±1)Ax) t o ( H 2 S ) . F i n a l l y , u s i n g a f i r s t o r d e r d i f f e r e n c e a p p r o x i m a t i o n f o r (H 2 S ) , the v a l u e s f o r H((L-1)Ax,(K±i)Ax), i=1,...,5, a r e e a s i l y d e t e r m i n e d from ( 2 ) . The r e s u l t of a l l t h i s i s s i m p l y t o smooth l a r g e i n i t i a l v a l u e s of S, and t o p r o v i d e a method of f i n d i n g H a t the d i s c h a r g e r e g i o n so t h a t i n i t i a l l y an e q u i l i b r i u m i n (H 2 S ) i s m a i n t a i n e d . When i n c o r p o r a t e d i n t o the model, t h i s p r o c e d u r e smoothed the r e s u l t s s u f f i c i e n t l y t o produce r e a s o n a b l e r e s u l t s . 49 5.3 F i n d i n g H On The B o u n d a r i e s In Chapter 1, we mentioned o n l y t h a t the boundary c o n d i t i o n s f o r H were p r e s c r i b e d on b oth th e ocean and l a n d b o u n d a r i e s . The c o m p u t a t i o n a l d e t a i l s a r e g i v e n i n t h i s s e c t i o n . On the open ocean b o u n d a r i e s H i s unknown. I t i s r e a s o n a b l e t o assume, however, t h a t H does not change r a p i d l y , and t h a t any changes t h a t do occur a r e due t o g r a d i e n t s of H (and hence S) w i t h i n the domain. Except at the f i r s t time l e v e l , where H i s s e t t o the v a l u e s used i n i t i a l l y , we c a l c u l a t e H on the ocean b o u n d a r i e s by e x t r a p o l a t i n g from v a l u e s a t the two p r e v i o u s time l e v e l s u s i n g the e x t r a p o l a t i o n f o r m u l a d i s c u s s e d i n Chapter 2. A f t e r the model e q u a t i o n f o r H has been s o l v e d a t the c u r r e n t time l e v e l , H on the boundary i s f u r t h e r m o d i f i e d by e x t r a p o l a t i n g from the two n e a r e s t i n t e r i o r p o i n t s on a g r i d l i n e . That i s , we use a t p o i n t ( i A x , 0 ) , f o r example. The v a l u e s on the ocean boundary thus r e f l e c t more c l o s e l y what i s o c c u r r i n g w i t h i n the domain. For H on the l a n d boundary, we must s o l v e e q u a t i o n ( 1 ) from the p r e v i o u s s e c t i o n w i t h as an i n i t i a l v a l u e the i n p u t v a l u e of H a t the r i v e r d i s c h a r g e p o i n t . S i n c e the d i s c h a r g e p o i n t i s i n the m i d d l e of the l a n d boundary, ( l ) must be s o l v e d i n each d i r e c t i o n . When s o l v i n g i t f o r the s o u t h e r n p o r t i o n of the boundary, the s i g n s of the y - d e r i v a t i v e s must be r e v e r s e d . In the p r o c e s s of f i n d i n g a s o l u t i o n t o e q u a t i o n ( 1 ) , a number of d i f f i c u l t i e s were e n c o u n t e r e d . One was the e f f e c t of 50 the s a l i n i t y d e f e c t v a l u e s chosen as i n p u t a t the r i v e r d i s c h a r g e p o i n t . For t e s t p u r p o s e s , we a s s i g n e d i n p u t v a l u e s t o S t h a t were l a r g e compared t o a d j a c e n t boundary and i n t e r i o r p o i n t s . Even w i t h the smoothing pr o c e d u r e d e s c r i b e d i n the l a s t s e c t i o n , t h i s l e d t o a l a r g e g r a d i e n t , Sy, i n each d i r e c t i o n and produced a s t i f f e q u a t i o n t h a t was d i f f i c u l t t o s o l v e near the d i s c h a r g e p o i n t . I n t e g r a t i o n methods u s i n g Ax as a s t e p s i z e o n l y produced p h y s i c a l l y i m p o s s i b l e r e s u l t s . I t i s n e c e s s a r y t h e r e f o r e t o use an i n t e g r a t i o n s t e p s i z e much s m a l l e r than Ax, which r e q u i r e s v a l u e s f o r S and U a t p o i n t s between our g r i d p o i n t s . For S, t h i s i s a c c o m p l i s h e d by a f i t t i n g a c u b i c s p l i n e t h rough S a t the g r i d p o i n t s and e v a l u a t i n g the cu r v e a t the r e q u i r e d p o i n t s . For U, which i s z e r o a t the g r i d p o i n t s c o r r e s p o n d i n g t o the l a n d boundary, but nonzero between the d i s c h a r g e p o i n t (JAx,KAx) and g r i d p o i n t s (JAx,(K+1)Ax) t o the n o r t h and (JAx,(K-1)Ax) t o the s o u t h , we need some f u n c t i o n which d e s c r i b e s i t s b e h a v i o r . We chose an e x p o n e n t i a l c u r v e w i t h a l a r g e n e g a t i v e exponent, such t h a t U(JAx,KAx) i s the i n p u t d i s c h a r g e v a l u e and U(JAx,(K+1)Ax), U(JAx,(K-1)Ax) a r e v e r y n e a r l y z e r o . T h i s , i t i s hoped, i s a r e a s o n a b l e a p p r o x i m a t i o n t o U a t the " r i v e r mouth". To i n c l u d e these nonzero v a l u e s of U i n e q u a t i o n ( 1 ) , we must r e i n t r o d u c e the term - u/(2fSCH) on the r i g h t s i d e of ( 1 ) , and the ODE becomes (5) Hy = -(U/(CH) + fH S y +2kSH x + k H S x ) / ( 2 f S ) A f u r t h e r problem a r i s e s when e v a l u a t i n g the x - d e r i v a t i v e terms. A p p r o x i m a t i n g H and S w i t h the second o r d e r d i f f e r e n c e 51 e q u a t i o n d e t a i l e d e a r l i e r i n v o l v e s v a l u e s a t g r i d p o i n t s i n s i d e the domain. U n f o r t u n a t e l y , these i n t e r i o r p o i n t s must be a t the same time l e v e l (say time=n+l) we a r e s o l v i n g e q u a t i o n (5) f o r ; v a l u e s e i t h e r at the i n t e r i o r p o i n t s or on the l a n d boundary must be e s t i m a t e d i n some way. Rather than e x t r a p o l a t i n g from p r e v i o u s time l e v e l s t o p r o v i d e approximate v a l u e s a t the i n t e r i o r p o i n t s , we f i r s t s o l v e the model e q u a t i o n s f o r H and S. The boundary v a l u e s f o r H on a l l b o u n d a r i e s are e x t r a p o l a t e d from the two p r e v i o u s l e v e l s i n the u s u a l way. Once s o l v e d , we have i n t e r i o r v a l u e s f o r H and S a t time n+1, H and S can now be a p p r o x i m a t e d , and e q u a t i o n (5) s o l v e d . T h i s r e s u l t s i n new boundary v a l u e s f o r H which r e p l a c e the e x t r a p o l a t e d v a l u e s . A s o l u t i o n t o e q u a t i o n (5) i s found u s i n g the UBC GEARB s u b r o u t i n e , which i s an i t e r a t i v e method u s i n g a backward d i f f e r e n t i a t i o n f o r m u l a as a p r e d i c t o r and the c h o r d method as a c o r r e c t o r . The reader s h o u l d c o n s u l t [10] f o r a more d e t a i l e d e x p l a n a t i o n . We note here t h a t s i n c e H must be e x t r a p o l a t e d , and s i n c e H can change r a p i d l y making the e x t r a p o l a t i o n l e s s a c c u r a t e , we w ish t o a v o i d h a v i n g t o f i n d l a n d boundary v a l u e s a t the i n t e r m e d i a t e time l e v e l , which does not c o r r e s p o n d to an a c t u a l t i m e . We do t h i s by s p l i t t i n g the e q u a t i o n such t h a t we s o l v e f o r H f i r s t i n the y - d i r e c t i o n , then i n the x - d i r e c t i o n . Land boundary v a l u e s a r e thus o n l y r e q u i r e d f o r times n and n+1. 52 5.4 Boundary C o n d i t i o n For S The boundary c o n d i t i o n f o r the s a l i n i t y d e f e c t i s the n o - f l o w c o n d i t i o n , which r e q u i r e s s e t t i n g the f i r s t d e r i v a t i v e to z e r o . U s i n g a f i r s t o r d e r a c c u r a t e a p p r o x i m a t i o n t o the d e r i v a t i v e has one major advantage; when i t i s d i s c r e t i z e d , the m a t r i x e q u a t i o n s remain t r i d i a g o n a l . As we have n o t e d , the s a l i n i t y d e f e c t changes r a p i d l y on the l a n d boundary, p a r t i c u l a r l y near the r i v e r d i s c h a r g e r e g i o n (depending on the i n p u t v a l u e s u s e d ) . For such c a s e s , i t i s v e r y i m p o r t a n t t o approximate the boundary c o n d i t i o n here w i t h a t l e a s t the same a c c u r a c y as the ADI scheme. A second o r d e r a c c u r a t e a p p r o x i m a t i o n t o the no-flow c o n d i t i o n i s u s u a l l y c e n t r e d about the g r i d p o i n t on the boundary, r e q u i r i n g a v a l u e f o r S o u t s i d e the domain. On the open ocean b o u n d a r i e s , t h i s i s a r e a s o n a b l e approach, but on the l a n d boundary, such a v a l u e has no p h y s i c a l meaning. I n s t e a d we use a d i f f e r e n c e f o r m u l a i n v o l v i n g o n l y i n t e r i o r p o i n t s ; f o r the western boundary, f o r example, t h i s i s ( 6 ) S.j - 4St>j /3 + S^ /3 = 0 , j = 1,...,L. S i m i l a r e q u a t i o n s a p p l y on the o t h e r b o u n d a r i e s . I n c o r p o r a t i n g (6) i n t o the ADI m a t r i x e q u a t i o n s f o r S, we see t h a t each row c o r r e s p o n d i n g t o the boundary c o n t a i n s two elements t o the r i g h t of the d i a g o n a l ; the t r i d i a g o n a l i t y of the system i s l o s t . Rather than t r y t o s o l v e a system w i t h a h a l f - b a n d w i d t h of t h r e e i n s t e a d of two, a row r e d u c t i o n i s performed on each of t h e s e rows b e f o r e s o l v i n g the whole system. For example, the f i r s t two rows of the m a t r i x e q u a t i o n f o r S a t 53 t=n+* a re 1 -4/3 1/3 0 a 21 a 2 3 0 T h i s can be reduced t o (1-a 2 , / (3a 2 3 ) ) (-4/3-a 2 2/(3a 2 3 )) 0 A 2 1 3 2 2 A 2 3 Reducing each such row f o r both l e f t and r i g h t b o u n d a r i e s r e t u r n s the system t o t r i d i a g o n a l i t y , and i n c r e a s e s the o p e r a t i o n count by o n l y 0 ( L ) . 5.5 T e s t s Of The Model E q u a t i o n s I n c o r p o r a t i n g these methods f o r i n i t i a l and boundary c o n d i t i o n s i n t o the FORTRAN program t e s t e d i n the f i r s t s e c t i o n of t h i s c h a p t e r , we can pro c e e d w i t h t e s t s of the f u l l model e q u a t i o n s . Input t o the program c o n s i s t s of t e s t d a t a a t the i n f l o w , which was c o n s i d e r e d t o be a s i n g l e g r i d p o i n t l o c a t e d a t the m i d p o i n t of the x=LAx boundary. For each time s t e p , the s a l i n i t y d e f e c t , S, and x - d i r e c t i o n v e l o c i t y , u, c o r r e s p o n d i n g t o t h i s p o i n t a r e i n p u t . The t h i c k n e s s of the upper l a y e r , H, i s c a l c u l a t e d as e x p l a i n e d i n s e c t i o n 4.3 so as t o be i n e q u i l i b r i u m w i t h these v a l u e s of S and u a t time t=0. Input f o r the f i r s t e x p e r i m e n t a l run c o n s i s t s of a p u l s e i n f l o w , w i t h S=5.0 p p t , u=0.3 m/sec, and the c o r r e s p o n d i n g e q u i l i b r i u m v a l u e f o r H, which i s H=17.66. I n i t i a l c o n d i t i o n s are S=2.0 ppt and H=20.0 m throughout the domain. The i n f l o w v a l u e s a re i n t r o d u c e d i n t o the model a t time t=0, and remain the same u n t i l the f i n a l time s t e p , t=50. 54 Based on the r e s u l t s of the f i r s t r u n , the second e x p e r i m e n t a l run i s i d e n t i c a l w i t h the f i r s t u n t i l t = l 5 . At t h i s p o i n t the boundary c o n d i t i o n s a t the i n f l o w a re changed t o u=0 and dS/dx=0, c o r r e s p o n d i n g t o s h u t t i n g o f f the s u p p l y of f r e s h water i n t o the model. H i s found by e x t r a p o l a t i n g from the two p r e v i o u s time s t e p s i n the u s u a l way. T h i s run i s c o n t i n u e d u n t i l time s t e p 40. 5.6 I n t e r p r e t a t i o n Of R e s u l t s The r e s u l t s a r e p r e s e n t e d i n F i g u r e s 3 t o 35 and c o n s i s t of c o n t o u r p l o t s of H and S and a c o n t o u r p l o t of H 2S ( r e p r e s e n t i n g the p r e s s u r e f i e l d ) o v e r l a y e d w i t h a v e c t o r p l o t of v e l o c i t y . The v e c t o r p l o t i s n o r m a l i z e d t o the l a r g e s t v e l o c i t y a t the time l e v e l , whose v a l u e i n m/sec i s p r i n t e d on the f i g u r e . The p o i n t source of f r e s h water i s a t the m i d p o i n t of the r i g h t y - a x i s , a t x=600 km, y=300 km. Due t o the a l g o r i t h m used f o r c o n t o u r i n g , a number of p l o t s have r e g i o n s of unusual l o o k i n g c o n t o u r s . The a r e a s where t h i s o c c u r s have v a l u e s t h a t a r e a l l a p p r o x i m a t e l y the same, and the c o n t o u r i n g program cannot produce a c c u r a t e c o n t o u r s . On a l l p l o t s the r e s u l t s f o r the western h a l f of the domain s h o u l d be i g n o r e d ; the a c c u r a t e c o n t o u r s a r e always smooth c u r v e s . The F i g u r e s f o r t=2 r e f l e c t n u m e r i c a l d i f f i c u l t i e s c r e a t e d by the i n i t i a l l a r g e g r a d i e n t s of S. D e s p i t e the smoothing procedure used, these g r a d i e n t s induce u n l i k e l y v a l u e s of H around the i n f l o w p o i n t . By u s i n g a l a r g e number of i t e r a t i o n s over the next few time s t e p s , v a l u e s f o r H are smoothed, as can 55 be seen by t=5. The r e g i o n of low s a l i n i t y moves s l o w l y n o r t h , due t o C o r i o l i s f o r c e , growing l a r g e r as more f r e s h water f l o w s i n the system. T h i s i s accompanied by an i n c r e a s e i n the t h i c k n e s s of the upper l a y e r i n t h i s r e g i o n ; to the s o u t h , H d e c r e a s e s . The f i g u r e s f o r the p r e s s u r e term H 2S and v e l o c i t y show the flow t o be p r i m a r i l y northwards. That the v e l o c i t i e s do not l i e e x a c t l y a l o n g the H 2S c o n t o u r s r e f l e c t the e f f e c t s of f r i c t i o n ; the f l o w i s pushed t o the l e f t of t h e s e c o n t o u r s . By time s t e p 15 t h i s p a t t e r n i s w e l l e s t a b l i s h e d and F i g u r e s f o r l a t e r time s t e p s show o n l y an i n c r e a s e i n the upper l a y e r t h i c k n e s s and a s p r e a d i n g by a d v e c t i o n of the f r e s h water a l o n g the s u r f a c e . As noted i n [ 6 ] , i n c r e a s i n g amounts of f r e s h water move west and south as w e l l due t o d i f f u s i v e e f f e c t s . By time s t e p 30, the l i m i t s of the domain and the boundary c o n d i t i o n s i n f o r c e t h e r e b e g i n to have an a f f e c t on both upper l a y e r t h i c k n e s s and f l o w . I t appears t h a t the f l o w i s b e i n g " r e f l e c t e d " away from the n o r t h e a s t c o r n e r , c a u s i n g f l o w back i n t o the domain. A p o s s i b l e cause i s the e x t r a p o l a t i o n scheme f o r H, which may g i v e too l a r g e a v a l u e f o r H on the boundary i n t h i s a r e a . W h i l e the r e s u l t s a r e not shown, the same boundary c o n d i t i o n s were used u n t i l time s t e p 100. The r e f l e c t i n g f l o w i n c r e a s e s a l o n g the n o r t h e r n boundary, but o t h e r w i s e the r e s u l t s change very l i t t l e from those a t time s t e p 40. B e i n g i d e n t i c a l t o the f i r s t t e s t u n t i l t = l 5 , the r e s u l t s f o r the second e x p e r i m e n t a l run begin a t time s t e p 20 ( F i g u r e s 56 27, 30 and 3 3 ) . At t=16, the f r e s h water f l o w i n t o the model i s " t u r n e d o f f " ; we would expect the v a l u e s t o s l o w l y r e t u r n t o t h e i r i n i t i a l c o n d i t i o n s . As can be seen i n F i g u r e s 27-35 t h i s i s i n f a c t what happens, as the g r a d i e n t s of S and H decrease due t o a d v e c t i o n and e n t r a i n m e n t of s a l t i e r lower l a y e r water i n t o the upper l a y e r . These v a l u e s d e c r e a s e from one time s t e p to the next u n t i l a t t=40 the maximum v a l u e of S i s about 3.6 ppt w h i l e t h a t of H has d e c r e a s e d from a h i g h of about 28 m. t o about 22 m. The r e g i o n of h i g h e s t s a l i n i t y d e f e c t moves s o u t h pushed by the f l o w , whose speeds (the maximum of which i s g i v e n on the H 2S f i g u r e s ) d i m i n i s h r a p i d l y . The f l o w c o n t i n u e s t o the l e f t of the H 2S c o n t o u r s . There appears t o be an eddy ( F i g u r e 35) j u s t t o the s o u t h of the i n f l o w p o i n t , perhaps due t o the " r e f l e c t i o n " i n the n o r t h e a s t c o r n e r of the domain. Both t h e s e e x p e r i m e n t a l runs show t h a t the n u m e r i c a l model produces r e a s o n a b l e r e s u l t s w i t h no s i g n of any i n s t a b i l i e s . The second run a l s o shows t h a t a g i v e n p h y s i c a l s i t u a t i o n , the system r e t u r n i n g t o r e s t , i s a d e q u a t e d l y m o d e l l e d . Two o t h e r t e s t s were c a r r i e d out which s e r v e t o i l l u s t r a t e f u r t h e r d i f f i c u l t i e s w i t h s o l v i n g the model e q u a t i o n s . Both produced u n s t a b l e s o l u t i o n s , so no f i g u r e s a r e a v a i l a b l e . The f i r s t used the boundary c o n d i t i o n s of the f i r s t t e s t d i s c u s s e d above u n t i l time s t e p 40. T h e r e a f t e r , the v e l o c i t y u remained the same, but S was s l o w l y d e c r e a s e d over time t o S=2. The v a l u e s of H were those found t o be i n e q u i l i b r i u m . T h i s models the e f f e c t on the system of a f l o w of v a r y i n g s a l i n i t y . 57 The s a l i n i t y d e f e c t was d e c r e a s e d a t the r a t e of as l i t t l e as 0.1 ppt per time s t e p ; even a t t h i s r a t e the model q u i c k l y d e v e l o p e d i n s t a b i l i t i e s i n the v a l u e s of S (soon f o l l o w e d by H). I t e r a t i n g a l a r g e number of t i m e s showed the s o l u t i o n f o r S t o be a c t u a l l y d i v e r g i n g . The o t h e r u n s u c c e s s f u l t e s t i n v o l v e d the u s u a l boundary c o n d i t i o n s , but w i t h a s p a t i a l mesh s i z e one q u a r t e r of t h a t i n d i c a t e d above. I n s t a b i l i t i e s r e s u l t e d w i t h i n two or t h r e e time s t e p s . I n c r e a s i n g the number of i t e r a t i o n s a g a i n l e d t o a d i v e r g i n g s o l u t i o n . 5.7 Recommendations For P r o c e e d i n g There a r e s e v e r a l s h o r t c o m i n g s i n the model e q u a t i o n s and the n u m e r i c a l methods we have chosen t o s o l v e them t h a t s h o u l d i n i t i a l l y be c o n s i d e r e d when i m p r o v i n g the model. The n u m e r i c a l d i f f i c u l t i e s must be d e a l t w i t h b e f o r e any changes can be made i n the e q u a t i o n s , so they w i l l be d i s c u s s e d f i r s t . The most o b v i o u s problem i s i l l u s t r a t e d by the two t e s t s d i s c u s s e d a t the end of the l a s t s e c t i o n . In b o t h c a s e s , the i n s t a b i l i t i e s seem t o be caused by l a r g e g r a d i e n t s near the c o a s t a l boundary. In the f i r s t , i t r e s u l t e d from r e d u c i n g the i n p u t s a l i n i t y d e f e c t from 5.0 ppt a f t e r r e a c h i n g a n e a r - s t e a d y s t a t e s o l u t i o n . The i n p u t s a l i n i t y d e f e c t dropped below S a t s u r r o u n d i n g g r i d p o i n t s p r o d u c i n g g r a d i e n t s t h a t a p p a r e n t l y c r e a t e d i n s t a b i l i t i e s i n e q u a t i o n ( 5 ) , the l a n d boundary ODE; the r e s u l t i n g v a l u e s of S on the c o a s t i n c r e a s e t o v a l u e s l a r g e r than any e n t e r i n g the model at the i n f l o w . 58 R e s u l t s from the second i n d i c a t e t h a t the same problem r e s u l t s from a decrease i n mesh s i z e . S i n c e the r i g h t - h a n d s i d e of (5) depends on x - d i r e c t i o n d e r i v a t i v e s of H and S which a r e approximated by d i f f e r e n c e e q u a t i o n s , t h e s e v a l u e s become l a r g e r (assuming a d j a c e n t v a l u e s of H and S a r e a p p r o x i m a t e l y the same) and (5) g e t s even s t i f f e r as the r e s o l u t i o n i n c r e a s e s . Again the s o l u t i o n blows up a l o n g the c o a s t a l boundary, even when i n c r e a s i n g the a c c u r a c y of the d e r i v a t i v e a p p r o x i m a t i o n s t o f o u r t h o r d e r . A method of s o l v i n g t h i s problem i s not r e a d i l y a p p a r e n t , a l t h o u g h a v a r i a t i o n of the a v e r a g i n g scheme used f o r the i n i t i a l c o n d i t i o n s might be c o n s i d e r e d . Another n u m e r i c a l problem as y e t u n r e s o l v e d i s c r e a t e d when the t r u e domain i s c o n s i d e r e d ; the l a n d boundary now c o n t a i n s two segments of the c o a s t where r u n - o f f e n t e r s the model. To the n o r t h of the n o r t h e r n segment and t o the so u t h of the so u t h e r n segment, e q u a t i o n (5) can be s o l v e d as u s u a l . Between thes e segments, however, the i n i t i a l v a l u e problem (IVP) becomes a boundary v a l u e problem (BVP), as t h e v a l u e s f o r H a t both ends ar e f i x e d . C o n s i d e r i n g the d i f f i c u l t i e s e n c o untered i n s o l v i n g (5) as an IVP, i t seems l i k e l y t h a t i t w i l l p r e s e n t f u r t h e r problems as a BVP, p r o b a b l y h a v i n g boundary l a y e r s a t both ends. To more a c c u r a t e l y r e f l e c t the p h y s i c a l p r o c e s s e s of the model, a number of a d d i t i o n s s h o u l d be c o n s i d e r e d . F i r s t , a r e a l i s t i c s e t of i n i t i a l c o n d i t i o n s f o r H and S a r e needed, i n v o l v i n g perhaps some s o r t of h i s t o r i c a l average. We s h o u l d a l s o o v e r l a y a g e o s t r o p h i c c u r r e n t f i e l d , both i n i t i a l l y ( t he model i s p r e s e n t l y c o n s i d e r e d a t r e s t ) , and a t each time l e v e l . 59 As noted i n the p r e v i o u s s e c t i o n , some method of more a c c u r a t e l y r e f l e c t i n g b e h a v i o r on the open-ocean b o u n d a r i e s i s r e q u i r e d , p a r t i c u l a r l y a t the i n t e r s e c t i o n of the c o a s t a l boundary. L a s t , the model e q u a t i o n s s h o u l d i n c l u d e a d v e c t i o n c u r r e n t terms and wind f o r c i n g , and wind-induced m i x i n g s h o u l d be added t o W , the entr a i n m e n t term. The e f f e c t s of t h e s e changes on the n u m e r i c a l methods f o r the s i m p l i f i e d model would have t o be d e t e r m i n e d . 60 F i g u r e 4 - T e s t 1 : S at time s tep 5 61 0 . 0 1 0 0 . C 2 0 0 . 0 3 0 0 . 0 4 0 0 . 0 soo.o soo.o K M . F i g u r e 5 - T e s t 1 : S at t ime s t ep 10 KM F i g u r e 6 - T e s t 1 : S at time s tep 15 62 KM -200.0 too. a soo.o F i g u r e 8 - Tes t 1: S at time s tep 30 63 F i g u r e 10 - Tes t 1: S at t ime s tep 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 F i g u r e 11 - Tes t 1: H a t t ime s tep 2 t I ! : i 1 1 ~ 0 0.0 100..3 230.0 300.0 400.0 SOO.O S00.3 KM F i g u r e 12 - T e s t 1: H a t time step 5 65 F i g u r e 14 - Tes t 1: H at time s tep 15 66 F i g u r e 15 - Tes t 1: H at time s t ep 20 F i g u r e 16 - Tes t H at time s tep 30 67 F i g u r e 18 - Tes t 1: H at t ime s tep 50 68 F i g u r e 20 - T e s t 1: H 2 S at time s tep 5 69 F i g u r e 22 - Tes t 1: H 2 S at time s tep 15 70 71 F i g u r e 26 - T e s t 1: H 2 S at time s tep 50 72 F i g u r e ' 28 - T e s t 2: S ' a t t i m e s t e p 30 73 I I i I l i V o 3.3 lttD.3 2 M . 3 330.3 40D.3 Sflfl.3 503.3 KM F i g u r e 29 - Tes t 2: S at t ime s tep 40 74 n i ! 1 1 1 r o 0 . 0 100.0 200 . 0 300 . 0 400 . 0 SOO.O 500 . 0 KM F i g u r e 30 - Tes t 2: H at t ime s tep 20 KM F i g u r e 3 1 - Tes t 2: H at time s tep 30 F i g u r e 32 - Tes t 2: H at t ime s tep 40 76 F i g u r e 34 - Tes t 2: H 2S at time s tep 30 77 sac.s saol2 9 . 3 • 2 0 0 . 3 3 0 0 . 3 KM 3 0 0 . 3 / 6 0 0 . 3 F i g u r e 35 - T e s t 2: H 2 S at time s tep 40 78 BIBLIOGRAPHY 1. Bo Pe d e r s e n , F., A Monograph on T u r b u l e n t E n t r a i n m e n t and F r i c t i o n i n Two Layer S t r a t i f i e d Flow, I n s t . Hydrodynamics and H y d r a u l i c E n g e n e e r i n g , Tech. U n i v . Denmark, Lyngby. S e r i e s Paper No. 25. (1980) 2. F o r s y t h e , G.E. and Wasow, W.R., F i n i t e D i f f e r e n c e Methods f o r P a r t i a l D i f f e r e n t i a l E q u a t i o n s , John W i l e y & Sons, New York (1960). 3. Gear, C.W., N u m e r i c a l I n i t i a l V a l u e Problems i n O r d i n a r y D i f f e r e n t i a l E q u a t i o n s , P r e n t i c e - H a l l , Englewood C l i f f s , W\~J~. ( 1 971 ) . 4. G o t t l i e b , D. and G u s t a f s s o n , B., G e n e r a l i z e d Du F o r t - F r a n k e l Methods f o r P a r a b o l i c I n i t i a l - B o u n d a r y V a l u e Problems, SIAM J . Numer. A n a l . , 13 (1976) 129-144. 5. K r i e s s , H. and O l i g e r , J . , Methods of the Approximate S o l u t i o n of Time Dependent Problems, GARP P u b l i c a t i o n s S e r i e s No. 10 (1973) . 6. L e B l o n d , P.H., Emery, W.J., N i c o l , T.O., A Model of Runoff D r i v e n C o a s t a l C i r c u l a t i o n , S u b m i t t e d f o r p u b l i c a t i o n t o C o a s t a l , E s t u a r i n e and S h e l f S c i e n c e , (1983) 7. L e B l o n d , P.H., Department of Oceanography, U n i v e r s i t y of B r i t i s h Columbia, Vancouver, B.C., p e r s o n a l communication. 8. Lees, M., An E x t r a p o l a t e d C r a n k - N i c o l s o n D i f f e r e n c e Scheme f o r Q u a s i l i n e a r P a r a b o l i c E q u a t i o n s , N u m e r i c a l S o l u t i o n of P a r t i a l D i f f e r e n t i a l E q u a t i o n s , Ames, W.F. ( E d . ) , Barnes & Noble, New York (1969), 193-201. 9. M i t c h e l l , A.R. and G r i f f i t h s , D.F., The F i n i t e D i f f e r e n c e Method i n P a r t i a l D i f f e r e n t i a l E q u a t i o n s , John W i l e y & Sons, C h i c h e s t e r (1980) 10. Moore, C.F., I n t e g r a t i o n of S t i f f Systems of O r d i n a r y D i f f e r e n t i a l E q u a t i o n s Having a Banded J a c o b i a n , Computing C e n t r e , U n i v e r s i t y of B r i t i s h Columbia (1979) . , 11. M o r r i s , J . L . , On the N u m e r i c a l S o l u t i o n of a Heat E q u a t i o n A s s o c i a t e d w i t h a Thermal P r i n t - H e a d , J . Comp. Phys., 5 (1970), 208-228. 12. R i c h t m y e r , R.D. and Morton, K.W., D i f f e r e n c e Methods f o r I n i t i a l - V a l u e Problems, I n t e r s c i e n c e P u b l i s h e r s , New York (1967) . 79 13. Royer, T . C , On the E f f e c t of P r e c i p i t a t i o n and Runoff on C o a s t a l C i r c u l a t i o n i n the G u l f of A l a s k a , J . Phys. Oceanogr., 9 (1979), 555-563. 14. Royer, T.C., Hansen, D.V. and P a s l i n s k i , D.J., C o a s t a l Flow i n the N o r t h e r n G u l f of A l a s k a as Observed by Dynamic Topography and S a t e l l i t e - T r a c k e d Drogued D r i f t Buoys, J . Phys. Oceanogr., 9 (1979), 785-801. 15. Royer, T . C , C o a s t a l F r e s h Water D i s c h a r g e i n the N o r t h e a s t P a c i f i c , 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 C o n d i t i o n s i n S p l i t t i n g Methods f o r P a r a b o l i c D i f f e r e n t i a l E q u a t i o n s , I n t . J . Numer. Meth. Eng. 17 (1981> 335-346. 17. Varah, J.M., S t a b i l i t y of D i f f e r e n c e A p p r o x i m a t i o n s t o the Mixed I n i t i a l Boundary V a l u e Problems f o r P a r a b o l i c Systems, SIAM J . Numer. A n a l . , 8 (1971) 598-615. 18. Varah, J.M., On the S t a b i l i t y of Boundary C o n d i t i o n s f o r S e p a r a b l e D i f f e r e n c e A p p r o x i m a t i o n s t o P a r a b o l i c E q u a t i o n s , SIAM J . Numer. A n a l . , 14 (1977) 1114-1125. 19. Varah, J.M., S t a b i l i t y R e s t r i c t i o n s on Second Order, Three L e v e l F i n i t e D i f f e r e n c e Schemes f o r P a r a b o l i c Systems, SIAM J . Numer. A n a l . , 17 (1980) 300-309.

Cite

Citation Scheme:

    

Usage Statistics

Country Views Downloads
China 13 0
Japan 8 0
United States 3 0
City Views Downloads
Beijing 13 0
Tokyo 8 0
Unknown 2 0
Ashburn 1 0

{[{ mDataHeader[type] }]} {[{ month[type] }]} {[{ tData[type] }]}

Share

Share to:

Comment

Related Items