Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

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

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Notice for Google Chrome users:
If you are having trouble viewing or searching the PDF with Google Chrome, please download it here instead.

Item Metadata

Download

Media
831-UBC_1983_A6_7 N52.pdf [ 3.65MB ]
Metadata
JSON: 831-1.0095808.json
JSON-LD: 831-1.0095808-ld.json
RDF/XML (Pretty): 831-1.0095808-rdf.xml
RDF/JSON: 831-1.0095808-rdf.json
Turtle: 831-1.0095808-turtle.txt
N-Triples: 831-1.0095808-rdf-ntriples.txt
Original Record: 831-1.0095808-source.json
Full Text
831-1.0095808-fulltext.txt
Citation
831-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 S^ , 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:

        

Citations by CSL (citeproc-js)

Usage Statistics

Share

Embed

Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                        
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            src="{[{embed.src}]}"
                            data-item="{[{embed.item}]}"
                            data-collection="{[{embed.collection}]}"
                            data-metadata="{[{embed.showMetadata}]}"
                            data-width="{[{embed.width}]}"
                            data-media="{[{embed.selectedMedia}]}"
                            async >
                            </script>
                            </div>
                        
                    
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:
https://iiif.library.ubc.ca/presentation/dsp.831.1-0095808/manifest

Comment

Related Items