Governing Equations and Model Reduction in Continuous Stratified Models

Detailed analysis of the governing equations in a continuous stratified model with vertical mode decomposition and Sturm-Liouville problem. The density perturbation and pressure perturbation are examined to derive a reduced 3D system. The Boussinesq approximation and hydrostatic balance are applied, popular in ocean general circulation models.

• Stratified Model
• Vertical Mode Decomposition
• Sturm-Liouville Problem
• Boussinesq Approximation
• Ocean Circulation

cha_k Follow

Uploaded on Sep 11, 2024 | 0 Views

Download Presentation

Please find below an Image/Link to download the presentation.

The content on the website is provided AS IS for your information and personal use only. It may not be sold, licensed, or shared on other websites without obtaining consent from the author. Download presentation by click this link. If you encounter any issues during the download, it is possible that the publisher has removed the file from their server.

Presentation Transcript

1. Eigenvalue problem ex) continuous stratified model Keywords vertical mode decomposition Sturm-Liouville problem

2. 2 Governing equations Primitive equations Du f u Dt r r r r r g r r 1 D Dt + = , p g = ( , ) T S + = 0, u Here, we do not assume the constant density However, the strategy is again reduction of the dimension (4D => 3D) First, the density perturbation is assumed as follows ( , , , ) x y z t = ? ( , , , ) x y z t + % % Boussinesq approximation 0 0

3. 3 Governing equations This density is substituted into the density equation D Dt r % r + + % g = ( ) 0 u 0 r r g = 0 u The largest term is 0 The remaining terms are r % % r D Dt D Dt + % g = = 0 u Hence, we obtain r r g = % D Dt = 0 u 0 and

4. 4 Governing equations The vertical momentum equation is w w u t x 1 + w y w z p z + + + = v w g % 0 We introduce p0, which is related to 0 p = ( , , , ) p x y z t + % 0( ) p z dp dz = 0 g 0 Thus, % % 1 w t w x w y w z p z + + + = u v w g 0 0

5. 5 Governing equations The zonal pressure term is + p % 1 1 + p p x = 0 % % x 0 1 p x 0 Similarly, the meridional pressure term is % 1 1 p y p y = 0

6. 6 Governing equations + + v v u v t x % 1 u t u x u y v y u z v z p x p y + = u v w fv 1 0 % + + + + = w fu 0 % % 1 w t w x w y w z p z + + + = u v w g 0 0 = % r r g D Dt 0 = 0 u Boussinesq equations We have 5 equations and 5 unknowns

7. 7 Governing equations We also assumed the hydrostatic balance % 1 u t u x u y u z p x + + + = u v w fv 0 % 1 v t v x v y v z p y + + + + = u v w fu 0 r r g = % D Dt % p z = % = 0 g 0 u These equations are popular for ocean general circulation models (OGCMs)

8. 8 Governing equations 1. 2. We consider the linear system The density perturbation is assumed as follows ( , , , ) x y z t = % ? ( , , , ) x y z t + ( ) z B where B The corresponding pressure perturbation is also ( , , , ) p x y z t % ( , , , ) p x y z t = + ( ) p z B Therefore, p = g B B z

9. 9 Governing equations A caution is needed for the density equation D Dt t = = % % % % % = + + + u v w x d y z + w B t dz N g 2 0 0 w t where d g buoyancy frequency (Brunt-V is l frequency) 2 N B dz 0

10. 10 Governing equations Thus, the equations are 1 u t v t p x p y = fv 1 0 + = fu 0 r r g p z = = g 0 u 2 0N g = w t We still have 5 equations and 5 unknowns

11. 11 Vertical mode decomposition To reduce the dimension, it is assumed that u, v, and p have the same vertical structures ( , , , ) ( , , , ) ( , , , )/ p x y z t = = = p x y t F z ( , , ) ( ) ( , , ) ( ) ( , , ) ( ) u x y z t v x y z t u x y t F z v x y t F z 0 Similarly, w and have the same vertical structures ( , , , ) ( , , , ) x y z t = = ( , , ) ( ) ( , , ) ( ) x y t G z w x y z t w x y t G z

12. 12 Vertical mode decomposition The equations become as follows u t v t p x p y = fv (1) + = fu (2) dF p dz + g (3) = G 0 u x v y t dG dz + = (4) 0 F w 2 N g = 0 w (5)

13. 13 Vertical mode decomposition By using the equations (3) and (5), we can eliminate the density perturbation p dF t dz (6) = 2 N Gw We also use the equation (4) to eliminate w and F t p + u v G N G separation of variables Note that it is assumed that G 0 x y = = zz 2 where is a constant value We have obtained two differential equations

14. 14 Vertical mode decomposition One equation is u x v y p t + = We also have two momentum equations fv t v t p y p x u + = fu = This set of the equations is equivalent to the linear shallow water equation one-layer model: if = gh and = (gH)-1 p

15. 15 Vertical mode decomposition The other differential equation is 2 d G dz N G = 2 2 If N is constant, this equation is easily solved The boundary conditions are that G(0) = G( H) = 0 (i) If > 0 = N z + N z 2 2 sin( ) cos( ) G A B (ii) If < 0 = N z + N z 2 2 exp( ) exp( ) G A B (iii) If = 0 = + G Az B

16. 16 Vertical mode decomposition From the boundary conditions, the cases of 0 are inadequate In the case of > 0, to meet the boundary conditions, n H B = 0 and = 2 N where n is integer Therefore, 2 n NH = is not an arbitrary value

17. 17 Vertical mode decomposition vertical structures of G n=1 n=2 n=3 0 0 0 -200 -200 -200 depth [m] -400 -400 -400 -600 -600 -600 -800 -800 -800 -1000 -1000 -1000 -1 0 1 -1 0 1 -1 0 1

18. 18 Vertical mode decomposition dG dz F can be obtained from the equation (4) as, F vertical structures of F n=1 0 n=2 n=3 0 0 -100 -100 -100 -200 -200 -200 -300 -300 -300 depth [m] -400 -400 -400 -500 -500 -500 -600 -600 -600 -700 -700 -700 -800 -800 -800 -900 -900 -900 -1000 -1000 -1000 -1 0 1 -1 0 1 -1 0 1

19. 19 Vertical mode decomposition However, N generally depends on the depth at 30 N, 180 E depth [m] N2 [s 2] potential density [kg/m3] In this case, does this equation have solutions?

20. 20 Sturm-Liouville problem If the differential equation has the following form d df p x g x f dx dx + ( ) r x f = ( ) ( ) 0 (0 1) x In addition, if the boundary conditions are ( ) A f x df dx df dx + = = 0 (at 0) B x 1 1 + = = ( ) 0 (at 1) A f x B x 2 2 This problem is called the Sturm-Liouville problem

21. 21 Sturm-Liouville problem 1. If a eigenvalue is a certain real value, the problem has a non-trivial answer (f 0) 2. Each n(n is integer) has a eigenfunction fn, which meets the original differential equation 3. The eigenvalues can ordered as follows: L L 1 2 n The eigenfunction fnhas n 1 zeros (node) in 0 < x < 1 If m n, the eigenfunctions are orthogonal 1 ( ) ( ) n r x f x f 4. 5. ( ) x dx = 0 m 0

22. 22 Sturm-Liouville problem In our case, p(x) = 1, g(x) = 0, and r(x) = N2 In the example that N is constant, 1. The eigenvalues are 2. The eigenfunction fnis a sine function 3. The eigenfunction fnhas n 1 zeros in H < x < 0 4. If m n, the eigenfunctions are orthogonal ( ) 2 = / n NH In general, we have to numerically solve

23. 23 How to solve If is known, the equation is an elliptic type By using the finite difference method, the equation can be represented by this form A = F We know how to solve this set of the linear equations However, is unknown

24. 24 How to solve In this case, we need to solve the eigenvalue problem The form of the original equation is changed 1 d G N dz 2 = G 2 2 By using the finite difference method, the LHS is 2 1 1 d G N dz N + 2 dz G G G = + 1 1 i i i 2 2 2 2 G G G G + 1 i = 2 2 2 2 2 2 1/ 2/ 1/ N dz N dz N dz i 1 i A (known) (unknown)

25. 25 How to solve Then, the problem becomes = AG G eigenvalue problem There are many program packages to solve eigenvalue problems e.g) Fortran: LAPACK, Python: np.linalg.eig MATLAB has the function eig [lambda, G] = eig(A) gives eigenvalue diagonal matrix lambda and matrix G whose columns are eigenvectors

26. 26 Example F G depth [m] 1= 0.2051, 2= 0.9427

27. 27 Physical interpretation The first mode resembles to the baroclinic mode that is obtained from the two layer model The first mode is called the first baroclinic mode Similarly, the second mode is called the second baroclinic mode n is the inverse number of the gravity speed of the corresponding mode In addition, 1/ n f gives the Rossby deformation radius of the corresponding mode

28. 28 Where is a barotropic mode? We assumed that G is not zero If G = 0, the equation (6) becomes p dF t dz = = 2 0 N Gw Hence, dF dz = = 0 and F const This indicates that u, v, and p are vertically uniform barotropic mode

29. 29 Orthogonality = ( , , , ) u x y z t ( , , ) u x y t F z ( ) n n n Multiplying Fm(z) and r(z) = N2and integrating from the bottom to the surface, we obtain 0 2 ( ) z F z u x y z t dz ( ) ( , , , ) N m H 0 = 2 ( ) z F z ( ) ( , , ) u x y t F z dz ( ) N m n n n H = ( , , ) m u x y t

30. 30 Orthogonality 1 u t p x = fv 0 By using the same procedure for the governing equations, we can obtain u p = m t m fv m x Therefore, the each mode independently obeys the shallow water equations By solving the equations of each mode and by summing the solutions, we can obtain the motion of the original governing equations

31. 31 Summary Boussinesq equations vertical mode decomposition Shallow water equations the sum of all modes Boussinesq equations

32. 32 Miscellaneous The barotropic and first baroclinic modes generally play major roles in the formation of ocean circulation The higher baroclinic modes ( 2nd) sometimes play an important role To apply the vertical mode decomposition has several limitations, such as linearity and flat ocean bottom In real ocean, there is interaction between different modes

33. 33 Sample program of Matlab 1/2 % read data load matlab.mat % caluculate N dz=depth(2:end)-depth(1:end-1); drhodz=(pden(:,2:end)-pden(:,1:end-1))./dz; N2=-9.8.*drhodz./1000; depth2=0.5*(depth(1:end-1)+depth(2:end)); % calculate A L=length(depth2); dz1=depth2(1,2:end-1)-depth2(1,1:end-2); dz2=depth2(1,3:end)-depth2(1,2:end-1); dZ=dz1+dz2; A=zeros(L,L); for ii=2:L-1 A(ii,ii-1:ii+1)=[2/dz1(ii-1)/dZ(ii-1) -2/dz1(ii-1)/dz2(ii-1)... 2/dz2(ii-1)/dZ(ii-1)]./(N2(ii)); end

34. 34 Sample program of Matlab 2/2 %%% solve A [V,q]=eig(A); q=diag(q); [q,idx]=sort(q,'descend'); for ii=1:L G(:,ii)=V(:,idx(ii)); end % calculate F dz=depth2(3:end)-depth2(1:end-2); dz1=depth2(2)-depth2(1); dz2=depth2(end)-depth2(end-1); F=NaN*ones(L); F(1,:)=(G(2,:)-G(1,:))./dz1; F(2:end-1,:)=(G(3:end,:)-G(1:end-2,:))./(dz'*ones(1,L)); F(end,:)=(G(end,:)-G(end-1,:))./dz2;

35. 35 Sample program of Python 1/2 import numpy as np from scipy import io import matplotlib as mpl import matplotlib.pyplot as plt matdata = io.loadmat('matlab.mat') pden= matdata['pden'] depth= matdata['depth'] zsz = depth.size # caluculate N dz = depth[0,1:]-depth[0,0:zsz-1] drhodz = (pden[0,1:]-pden[0,0:zsz-1])/dz N2 = -9.8*drhodz/1000 depth2 = 0.5*(depth[0,1:]+depth[0,0:zsz-1]) # calculate A L=depth2.size dz1 = depth2[1:L-1]-depth2[0:L-2] dz2 = depth2[2:]-depth2[1:L-1] dZ = dz1+dz2

36. 36 Sample program of Python 2/2 A = np.zeros((L, L)) for ii in range(1,L-1): A[ii,ii-1] = 2/dz1[ii-1]/dZ[ii-1]/N2[ii] A[ii,ii] = -2/dz1[ii-1]/dz2[ii-1]/N2[ii] A[ii,ii+1] = 2/dz2[ii-1]/dZ[ii-1]/N2[ii] # solve A q, G=np.linalg.eig(A) idx = np.argsort(q)[::-1] q = q[idx] G = G[:,idx] # calculate F dz = depth2[2:]-depth2[1:L-1] dz=np.reshape(dz,(L-2,1)) dz1 = depth2[1]-depth2[0] dz2 = depth2[-1]-depth2[-2] F = np.zeros((L, L)) F[0,:]=(G[1,:]-G[0,:])/dz1 F[1:L-1,:]=(G[2:L,:]-G[1:L-1,:])/dz F[L-1,:]=(G[L-1,:]-G[L-2,:])/dz2