VIC Frozen Soils Simulation in Permafrost Regions: Modifications, Improvements, and Testing

1. VIC Frozen Soils Simulation in Permafrost Regions: Modifications, Improvements, and Testing Jennifer Adam Amanda Tan May 2, 2007 Hydro Group Seminar

2. Talk Overview Motivation Model Description Model Improvements Study Domain Testing of Model Improvements Conclusions

3. Motivation

4. 1940 1960 1980 2000 Trend, mm/year-2 1. Precipitation/Streamflow Trend Inconsistency in Permafrost Basins (e.g. Aldan River) Streamflow Precipitation Adam and Lettenmaier (2007)

5. cm km cm m Permafrost Definition (loose) Seasonally Frozen Ground Permafrost Active Layer Permafrost Layer Winter Summer Winter Summer

6. 2. Active Layer Depth Simulated Poorly • Snow cover extent, lake freeze and break-up dates, streamflow climatology – satisfactory • Permafrost active layer depth – unsatisfactory Su et al. (2005)

7. Yenisei Ob’ Q, 103 m3s-1 Q, 103 m3s-1 Lena 1940 1960 1980 1940 1960 1980 2000 Observed Simulated 3. Historical Streamflow Trends Not Captured in Permafrost Basins Permafrost Basins Non-Permafrost Basin

8. Model Description

10. Frozen Soils Simulation: Heat Equation Term 2: Vertical Heat Conduction (Space) Term 1: Heat Storage (Time) Term 3: Latent Heat (Time) (Cherkauer et al. 1999) With parameterizations for Cs, κ, and θi

11. VIC Frozen Soils Algorithm • Cherkauer and Lettenmaier (1999) finite difference algorithm • solving of thermal fluxes through soil column • infiltration/runoff response adjusted to account for effects of soil ice content • parameterization for spatial distribution of frost • tracks multiple freeze/thaw layers • can use either “zero flux” or “constant temperature” bottom boundary (Cherkauer et al. 1999)

12. Current Implementation (Cherkauer et al. 1999)

13. Spatial Frost Algorithm • To produce a spatial distribution of ice content and subsequently soil moisture drainage. (Cherkauer et al. 2003)

14. Frozen Soils Options… WOW! • NODES 18 # number of soil thermal nodes • FULL_ENERGY TRUE # run full energy balance mode • GRND_FLUX TRUE # solve surface energy balance • FROZEN_SOIL TRUE # run frozen soils • QUICK_FLUX FALSE # Liang et al. 1999, otherwise Cherkauer et al. 1999 • QUICK_SOLVE FALSE # Cherkauer et al. 1999 for final step only, Liang et al. 1999 for rest • NOFLUX TRUE # zero flux bottom boundary • IMPLICIT TRUE # uses implicit solver (NEW) • EXP_TRANS TRUE # exponential grid transform (NEW) • QUICK_FS FALSE # linear equations for max unfrozen water content (not tested?) • QUICK_FS_TEMPS 7 • SPATIAL_FROST TRUE # sub-grid frost distribution • FROST_SUBAREAS 10 • Dp 15 # damping depth (m) • Tave -4.2 # temp. at damping depth (°C)

15. Changes and Improvements Zero Flux Bottom Boundary Parameterization Exponential Grid Transformation Implicit Solver Patch for the “Cold Nose” Problem

16. Constant T BB Dp = 4 m Tb,init = -12 °C Zero Flux BB Dp = 4 m Tb,init = -12 °C Zero Flux BB Dp = 15 m Tb,init = -3 °C 10 Soil Temperature, °C 0 -10 -20 1930 1931 1932 1930 1931 1932 1930 1931 1932 Soil Surface (Top Boundary) Soil Bottom Boundary 1. Soil Temperature Sensitivity to Bottom Boundary Specification

17. Annual Soil Temperature @ 3.2 m, C -6 -4 -2 0 2 4 6 BB Temperature Initialization A) Frauenfield et al. 2004 station data B) Interpolated station data

18. i=0 1 m 3 m 5 m 7 m 9 m 11 m 13 m 15 m i=Nnodes-1 2. Exponential Grid Transformation • Parameters: Dp=15m, Nnodes=18 Linear Distribution Exponential Distribution z1=soil depth1 z2=2·soil depth1 zi=linearly distributed to Dp Problem: discontinuity in Δz between nodes 2 and 3 zi=exp(b·i)+c Problem: discontinuity in Δz betweenall nodes

19. Exponential Grid Transformation, continued • Transform spatial derivatives only (temporal derivatives are unaffected) • Expand heat conduction term (chain rule), because κ varies with z

20. Exponential Grid Transformation, continued • Introduce new space variable, η (T will vary exponentially with z, but linearly with η) • Develop transform function, η=f(z) (chain rule)

21. Exponential Grid Transformation, continued • Determine partial derivatives (in η) • Substitute partial derivatives with η into heat conduction terms new elements

22. Exponential Grid Transformation, continued • Determine constants, b and c • Boundary condition 1: z(η=Nnodes-1) = Dp • Boundary condition 2: z(η=0) = 0 • Therefore: c = -1 • Solve in η-space (because the η-nodes are linearly distributed, the finite difference assumptions are not compromised) • Map temperatures in η-space back to z-space • Recalculate Cs, κ, and θi as a function of T for each z-node

23. 0 0 1 1 2 2 Space (Thermal Nodes) Space (Thermal Nodes) 3 3 4 4 5 5 6 6 t-1 t t-1 t Time Time unknown known 3. Implicit Solvers Explicit (forward in time) Implicit (backward in time) 7 explicit equations solved independently 7 implicit equations solved simultaneously

24. Stability Issues • Stability of convergence • Implicit: unconditionally stable • Explicit: satisfy the Courant-Friedrichs-Lewy Condition λ≤1/4 for no oscillating errors λ=1/6 to minimize truncation error (solution: make Δt≤1hr or Δz≥0.2m ) • Ability to find a solution • Explicit: not an issue with physically reasonable values (root_brent is very robust!) • Implicit: often is unable to find a solution at the initial formation of ice in the soil column =1 to 1.5 for Δt=3hr, Δz=0.1m

25. Implicit VIC Frozen Soils Algorithm • Newton-Raphson method to solve non-linear system of simultaneous equations (Ming Pan) • New Functions (mainly from Numerical Recipes) • solve_T_profile_implicit (replaces solve_T_profile) • fda_heat_eqn (replaces soil_thermal_eqn) • newt_raph (replaces root_brent) • fdjac3 (approximation for Jacobian) • tridiag (solves tridiagonal linear system) • Modifications from original code • Merged to VIC 4.1.0 from VIC 4.0.3 • Added NOFLUX and EXP_TRANS options • Nodal updating of Cs, κ, and θi as a function of T during iteration • Allowed for time-varying Cs: • When unable to find a solution, defers to explicit solver for that time-step

26. Temperature, °C Time Step Soil Top Boundary Soil Bottom Boundary 4. “Cold Nose” Problem • Run-away temperatures in near-surface thermal nodes • The coldest node becomes colder and breaks the 2nd law of thermodynamics, e.g. “Heat cannot of itself pass from a colder to a hotter body” • VIC crashes – error statement is “increase SOIL_DT” or “increase SURF_DT” • Occurs for all versions of VIC when using the finite-difference scheme, with all modes (implicit/explicit, noflux, exp/linear, etc…)

27. Heat Equation: Conduction: Term1 Term2 Explanation chain rule • As T decreases, κ increases (especially if θi increases) • Therefore, at a “cold nose”: but • If |Term1|>|Term2|, heat flows from the cold node • T at that node escape towards -∞

28. TU 0 i=0 T 10 i=1 Depth, cm Finite Difference Approx. TL i=2 20 Actual i=3 30 -10 -5 0 -15 Temperature, °C The “Cold Nose” Patch • Is calculation being made for the two near-surface nodes? • Is the node fall on a “cold nose”? • Is Term1 greater than Term2 in absolute value? • THEN Term1 = 0 • To allow for some lenience, can also check that • |TL-TU| > 5

29. Summary of Changes • Zero Flux Bottom Boundary Parameterization • Change in implementation, not in code • Involves also increasing Dp and Nnodes • Necessary for climate change studies in permafrost regions • Exponential Grid Transformation • Allows for closer node spacing near surface • Solves problem of discontinuity in Δz • Implicit Solver • Should give more accurate solution • No convergence instability (no wildly wrong results) • Nodal updating, change of Cs with time • Defers to explicit when no solution is found • Should give lower simulation time, but doesn’t… • Patch for the “Cold Nose” Problem • Inelegant, but it works (exists in implicit and explicit modes) • Does this problem go away if Δz becomes infinitesimally small?

30. Arctic Ocean Study Area: Aldan River Basin

31. Forest Shrubland Savanna Grassland Wetland Cropland Urban Barren Tundra Lena River Basin Aldan Tributary 700,000 km2 Revenga et al. 1998

32. Continuous , 90-100% Isolated, <10% Discontinuous, 50-90% Seasonally Frozen Ground Sporadic, 10-50% Permafrost Distribution Aldan: 89% continuous 10% discontinuous 1% sporadic Brown et al. 1998

33. Model Testing

34. Soil Moisture Soil Temperature Snow Depth Precipitation Air Temperature In-Situ Observations

35. Simulations conducted

36. Comparison of simulation time

37. Soil Moisture Comparison

38. Soil Moisture Comparison, Optimized Run (Run 1) --- Observed (Liquid) --- Total Soil Moisture --- Liquid Water --- Ice Content

39. Soil Temperature Comparison

40. Soil Temperature Comparison, Optimized Run (Run 1)

41. Streamflow comparison

42. Streamflow Comparison

43. Conclusions

44. Conclusions • Optimum set of parameters: • Damping depth = 15m • Number of nodes = 18 • Implicit solution • Exponential distribution of thermal nodes • No flux bottom boundary • Traditional setup • Constant flux, Dp=4m, 5 nodes • Gives lowest simulation time • Not suitable for climate change/permafrost studies • When using the optimized parameters, calibration for streamflow is required