Atmospheric Transport Models - Numerical Integration

Atmospheric Physics & Chemistry Solving the Basic Equations for the Atmosphere in 3-D Mass Continuity Equations of Motion (momentum continuity) ∂ρ d (uρ ) d (vρ ) d (wρ ) =− − − ∂t dx dy dz d (uu) d (vu) d (wu) ∂u =− − − + ... ∂t dx dy dz ∂v d (uv) d (vv) d (wv) =− − − + ... ∂t dx dy dz ∂w d (uw) d (vw) d (ww) + ... =− − − ∂t dx dy dz Thermodynamic ∂T d (uT ) d (vT ) d (wT ) 1 ⎛ =− − − + ⎜J Equation ∂t dx dy dz cv ⎝ (energy continuity) J : radiation, conduction, latent heat release, etc + pressure gradient + Coriolis force + gravity + friction −p D(1/ ρ ) ⎞ ⎟ Dt ⎠ D(1/ρ) / Dt : conversion between thermal and mechanical energy in fluid system Chemical Continuity Equation ∂χ d (uχ ) d (vχ ) d (wχ ) + Chemical Production =− − − – Chemical Loss ∂t dx dy dz 3. Spatial representations a. Finite difference schemes (truncated Taylor expansion at J grid-points) b. Spectral techniques (express variables using truncated series of N orthogonal harmonic functions and solve for N coefficients of expansion;) see c. Interpolation schemes (interpolates between grid points e.g. using a polynomial) d. Finite element schemes (minimizes error between actual and approximate solutions using a “basis function”, good for irregular geometries, c.f. (b) above which is good for regular geometries) 4. Explicit and Implicit time stepping t+Δt t Explicit: ( )x = f ⎡⎣...., ( )x* ,...⎤⎦ t+Δt t t+Δt Implicit: ( )x = f ⎡⎣..., ( )x* , ( )x* ,....⎤⎦ (Implicit methods more stable (but often less accurate) than explicit methods for longer time steps) Time stepping and stability In the numerical model, time is treated in discrete units and the time intervals chosen depend on the size of the model grid boxes. Intuitively, don’t want to transport across more than a grid cell over a time step. General Rule for stability: the CFL condition u∆t ≤1 ∆x ∆x u∆t e.g. Typically in the atmosphere, max u = 100m/s & grid spacing = 200 km, so constraint is Δt < 2000 seconds (33min) 5. Example Finite Difference Schemes ∆t n-1 n t n+1 ∆x x j-1 j j+1 (A) Advection (a) Forward / Upstream (explicit, conditionally* stable) (time) (space) ( ∂( ) → ∂t )j n +1 −( Δt ⎧ ⎧( ⎪u nj ⎪⎨ ⎪ ∂ ( ) ⎪ ⎩⎪ u →⎨ ∂x ⎪ n ⎪⎧ ( ⎪u j ⎨ ⎪⎩ ⎩⎪ )j n (first order accurate) ) j − ( ) j−1 ⎫⎪⎫⎪ n n Δx ) j+1 − ( n Δx ⎬ ⎪⎭⎪⎪ ⎬ n ) j ⎫⎪⎪ ⎬⎪ ⎭⎪⎭⎪ *For stability need Courant No. (u j > 0) (first order) (u j ≤ 0) uΔt ≤1 Δx (b) Centered / Centered or Leap-frog (explicit, neutrally** stable) ( ∂( ) → ∂t )j n +1 −( )j n−1 (second order) 2Δt ⎧ ( )n − ( )n ⎫ ∂( ) j+1 j−1 ⎪ n ⎪ (second order) → uj ⎨ u ⎬ ∂x 2Δx ⎪⎩ ⎪⎭ **No amplitude dissipation but still need Courant No. ≤ 1 (Note that forward/centered is linearly unstable to small perturbations) (B) Diffusion (i) Forward / Centered (explicit, conditionally*** stable) ∂( ) as in (a) above ∂t ⎧( ∂ ⎛ ∂( ) ⎞ ⎪ → K K ⎜ ⎟ 1 ⎨ j+ ∂x ⎝ ∂x ⎠ 2 ⎪ ⎩ ⎧( ) − ( ) ⎫ ) j+1 − ( ) j ⎫⎪ ⎪ j j−1 ⎪ − K ⎬ ⎬ 1 ⎨ 2 2 j− Δx ( Δx ) ⎪⎭ 2 ⎪ ⎪⎭ ⎩ ( ) n n ***For stability need Fourier No. n KΔt ( Δx ) 2 ≤ n 1 4 (ii) Backward / Centered (implicit, unconditionally stable) ∂( ) ∂ ⎛ ∂( ) ⎞ as in (a) above and ⎟ as in (i) above but replace n by n + 1 ⎜K ∂t ∂x ⎝ ∂x ⎠ (iii) Backward-forward / Centered (semi-implicit, stable for all Δt , Crank-Nicholson Scheme) ∂( ) as in (a) above ∂t ∂ ⎛ ∂( ) ⎞ ⎟ is AVERAGE of (i) and (ii) above ⎜K ∂x ⎝ ∂x ⎠ (Note that this method second-order accurate in space and time since both centered at 1 time n + ) 2 (C) Chemistry (involves a non-linear vectorial equation) (a) Explicit forward G G X nj +1 − X nj Δt G G = R X nj , t n requires Δt ≤ ( 2 λ max ) where λ max G ∂R is the largest eigenvalue of G (often ∂X impractical) G (b) Implicit backward (replace n in R above by n +1 ) and semi-implicit (average of explicit forward / implicit backward) allow larger Δt and are preferred. (c) For higher order accuracy than above schemes (which are first order accurate) use predictor-corrector methods (e.g. Gear) or generalized Runga-Kutta (e.g. KapsRentrop). These methods are however generally too computationally demanding for 3D chemical transport models. (d) Hybrid schemes. For 3D models can use: diagnostic equations (steady-state) for [i ] < Δt ; (a) or (b) above for τ > 100Δt and the analytical species i with τi = i Li 10 solution assuming constant Pi and τi : ⎡ Δt ⎤ n n Δt < τi < 100Δt . X ijn +1 = ( Pijτij ) + ⎡⎢ X ijn − ( Pijτij ) ⎤⎥ exp ⎢ − n ⎥ for ⎦ ⎣ 10 ⎣⎢ τij ⎦⎥ Note that the hybrid method is inherently non-conservative so corrections required. (e) Negative mixing ratios/concentrations. Whenever Xi < 0 (nonphysical!) then need to replace by Xi = 0 and lower either adjacent (grid-point) or global (spectral) Xi to compensate (sometimes called “borrow and fill”) 6. Surface fluxes ( φij ) Consider chemical species i and surface grid points j (a) specified natural or anthropogenic emissions 1 1 n+ n+ ⎛ [i ]n +1 − [i ]n ⎞ φij 2 ΔxΔy φij 2 j j ⎜ ⎟ = = Δz ⎜ Δt ⎟ ΔxΔyΔz ⎝ ⎠emissions (b) interactive fluxes (ocean) (two-way) w φij = p [i ]eq, j − [i ] j HL ( ) ( [i ]eq, j is the concentration in surface air if gas is in equilibrium with surface ocean) where w p = “piston velocity” (monotonically increases with surface wind speed, determined empirically) H L = dimensionless Henry’s Law coefficient (measured in laboratory) (c) deposition fluxes (one way) φij = − w s [i ] j ( w s = dry deposition velocity; empirical, depends on gas and surface type) Wet deposition important for soluble species: φij = −P ⎡⎣i( aq ) ⎤⎦ (P = precipitation rate; cm/sec) − P [i ] j j (equilibrium assumed between gas and raindrop) HL 7. Upper boundary conditions • Specified [i ] j (including [i ] = 0 ). Specification can be direct or specify τi  0 in uppermost layer. • Specified φij (including φij = 0 ). Can usually obtain φij = 0 by equating Xij values in top 2-3 layers (depending on finite difference scheme in vertical), or assuming K j or w j = 0. Recall: φij  − K j [ M ] j φij = w j [i ] j ∂X ij ∂z (diffusion), or (advection) 8. Subgridscale parameterizations a. Eddy diffusion coefficients ∂ ln θ • K zz (RR i ) ; i = g ∂z ⎛ ⎛ ∂u ⎞ 2 ⎛ ∂v ⎞ 2 ⎞ + ⎜⎜ ⎝⎜ ∂z ⎟⎠ ⎜⎝ ∂z ⎠⎟ ⎟⎟ ⎝ ⎠ (gradient Richardson no.) 1 R i > 0 → stable (if R i > get laminar flow) 4 R i < 0 → unstable (if R i ≤ 1 then forced convection and if R i > 1 then free convection) R i  0 → neutral ( K zz  kzu* , u * = friction velocity) ⎛ ∂T ∂T ∂u ∂u ⎞ K xx , K yy ⎜ , , , , etc. ⎟ (e.g. due to baroclinic eddies) ⎠ ⎝ ∂x ∂y ∂x ∂y b. moist convection ∂θE • ≤ 0 where θE = equivalent potential temperature implies convective ∂z instability • treatments range from simple “convective adjustment” (transport heat/mass ∂θ necessary to restore E = 0 ), to more complex “process-resolving” models ∂z (Kuo, Arakawa-Schubert, Hack (NCAR,CCM2), Tiedtke (ECMWF and ECHAM3), and Emanuel). For chemical models we want the mass fluxes not to the energy fluxes from these various treatments. 9. Chemical rate constants Consider the simplified ozone layer chemical reactions: 1 O 2 + hν ⎯J⎯ →O + O • l O + O 2 + M ⎯⎯ → O3 + M 1 O3 + hν ⎯J⎯ → O2 + O O + O3 ⎯k⎯ → O2 + O2 ( catalysed!) The relevant chemical reaction rates are expressed using first ( Ji ), second (k) and third (l) order rate constants: ⎧ -1 −3 ⎪−J i [i ] ( sec ⋅ molecule ⋅ cm ) 2 d [i ] ⎛ molecule ⎞ ⎪⎪ sec-1 ⋅ cm3 ⋅ molecule−1 ⋅ ( molecule ⋅ cm −3 ) ⎜ ⎟ = ⎨−k ij [i ] [ j] 3 dt ⎝ cm sec ⎠ ⎪ 3 ⎪−l [i ][ j] [ M ] sec-1 ⋅ cm 6 ⋅ molecule −2 ⋅ ( molecule ⋅ cm −3 ) ijM ⎪⎩ The chemical rate constants (k,l) are measured in the laboratory. ( ( ) ) Some typical expressions for their dependence on temperature (T) and density ([ M ]) are: ⎛ B⎞ k = A exp ⎜ − ⎟ ⎝ T⎠ ⎛ T ⎞ l = l ( Tref , [ M ]) ⎜ ⎟ ⎝ Tref ⎠ −α ( measure A and B ) ( measure l ( T [ M]) and α ) ref, The rate constant for photodissociation is calculated in a non-scattering atmosphere using: λ2 ⎡ N M j (z) ⎤ J i = ∫ σi ( λ ) φi ( λ ) I ( ∞ ) exp ⎢ −∑ σ j ( λ ) ⎥ dλ cos θ ⎦ λ1 ⎣ j=1 where σi ( λ ) = absorption cross-section at wavelength λ ( cm 2 ⋅ molecule −1 ) φi ( λ ) = photodissociation yeild ( dimensionless ) λ 2 − λ1 = width of electronic absorption band I ( ∞ ) = solar photon flux at altitude z = ∞ ( photon ⋅ cm -2 ⋅sec −1 ) N = number of gases ( j) absorbing at wavelength λ M j ( z ) = molecules of j per unit area above z ( molecule ⋅ cm -2 ) θ = solar zenith angle Sun θ z Surface 10. Some essential chemistry and radiation components a. UV fluxes for photodissociation rates b. For all species involving OH in their chemistry need to include: 1. O3 , O 2 , O ( 1 D ) 2. H, OH, HO2 , H 2 O2 , with latter 3 in gas and aqueous phase 3. NO, NO2 , NO3 , N 2 O5 , HNO3 with latter 2 in gas and aqueous phase 4. CH 4 , CH3 , CH3O2 , CH3O, CH3O2 H, CH 2 O, CHO, CO (also selected heavier hydrocarbons such as isoprene and terpenes in forested areas and anthropogenic hydrocarbons in urban areas) 11. Spectral (spherical harmonic) models Ynm (λ, μ) ≡ Pnm (μ)eimλ natural since eigensolutions of barotropic wave equation grid value: ξ(λ i , μ j ) = N( m ) M ∑ ∑ξ m = − M n =|m| spectral coefficient: ξ nm = P (μ j )eimλi m m n n ⎧ 1 2M ⎫ ξ(λ l , μ j )eimλl ⎬Pnm (μ j )w j ⎨ ∑ ∑ j =1 ⎩ 2M l =1 ⎭ J (wj - Gaussian weights) where λ and μ respectively denote the zonal and meridional independent variables, Pnm are the associated Legendre functions for which m denotes zonal wavenumber, n-m denotes a form of meridional wavenumber, M and N are the spectral truncation limits, and J is the order of the north-south transform grid, a function of the truncation parameters. 30 n = m Rhomboidal (R 15) Triangular (T 21) 21 N n Insert the image on page 11 of lec11.pdf. 15 M 0 0 15 21 30 m (zonal wave number) Figure by MIT OCW. Note: Physics and chemistry are computed on the equivalent grid points and then transformed to spherical harmonic representation. 12. Examples (a) MIT (1989) Chemistry-Dynamics Model CHClF2 simulations Model details: Horizontal – R6 (spectral), 15 lat x 16 long (grid) Vertical – ln P (surface to ≈ 72 km), 26 levels Explicit chemistry and radiation “balance-type” dynamics, parameterized trop. hating time – 4 x 1 hour, Lorenz n-cycle.