POP simulation in an eddying regime

1
POP simulation in an eddying regime
Matthew Hecht, Mathew Maltrud Los Alamos National
Laboratory and collaborators
2
How did Ocean Modelingcome to the high desert?
3

• 1979 Bob Malone builds support for a new
  endeavor at the Lab
• Build a climate modeling capability
• Bob goes to Boulder
• 5 consecutive summers at NCAR
• Contributed to first version of NCAR Community
  Climate Model
• earliest version of the atmospheric component of
  the Community Climate System Model
• Led by Eric Pitcher, V. Ramanathan

4
(No Transcript)
5
(No Transcript)
6

• Climate research at LANL goes dormant for a few
  years,
• but a base of experience has been gained

7

• Support from the research arm of the Dept of
  Energy comes in 1990
• DOE CHAMMP program
• Bob Malone coordinates model development
• First LANL code projects
• parallelization of 2 models, atm and ocn
• Phil Jones works on GFDLs Skyhigh,
• Rick Smith takes on the Bryan-Cox code
• John Dukowicz is interested in efficient solution
  of the 2-D barotropic mode.

8
Los Alamos ocean model POP

• Its a Bryan Cox code
• Z-coordinate
• Branched off from the Semtner-Chervin Code
• P for Parallel (Ocean Program)
• Back then, the big models used maybe 8 pes
• POP ran on 100s

9
Physical science, along with model development

• Model development is more successful with
  in-house expertise in the physical science
• High resolution simulation allows
• Oceanography that even an oceanographer might
  find interesting,
• along with computational challenge

10
First ocean modeling study
(0.28º global)
11
Barotropic stream function
12
From analysis of Fu and Smith (1996)

• From analysis of Fu and Smith (1996)
• Overall energy level of model lower than obs
• Spectral slope in frequency steeper than observed
• At wavelengths less than 200km,
• Obsd spectrum is white
• Model spectrum is red

13
Aside Another opportunity in model development

• Elizabeth Hunke and John Dukowicz establish new
  rheology
• 1997 first public release
• 1998 Model adopted as sea ice component by
  Community Climate System Model
• Our ocean model is also in the CCSM

14
Second ocean modeling study
15
1996-97, on a Connection Machine
16

• State of the field at this time (1996-97)
• Models still eddy-admitting
• TOPEX/Poseidon altimetry available
• Several excellent reviews in Warm-watersphere of
  the North Atlantic Ocean, 1996
• Stammer and Boning Boning and Bryan Dengg et al.

17

• Dengg et al. 1996 reviewed the Gulf Stream
  problem
• No single mechanism entirely dominant
• Downplayed role of DWBC
• Note New reviews now from Chassignet and
  Marshall, Marshall and Kiss.
• Boning and Bryan 1996 raise likelihood that
  threshold in resolution yet to be crossed

18
Eddy-admitting era

• Community Modeling Effort (CME, 1989)
• 1º, 1/3º and 1/6º North Atlantic
• Better with resolution, in many respects
• but Gulf Stream System and some other features
  still badly biased
• Chao et al. (1996) had a GS that was nearly
  separating, using modified CME configuration, at
  1/6, in POP
• Fine Resolution Antarctic Model (FRAM, 1991)

19

• POP North Atlantic modeling
• First presentations in 1997
• 1998 WOCE Newsletter
• 2000 J. Physical Oceanography
• Smith, Maltrud, Bryan and Hecht
• Note also Paiva et al. 1999 simulation with
  MICOM at 1/12

20
Forcing of Atlantic simulations(unless otherwise
specified)

• Barnier monthly heat flux climatology
• with penetrative solar
• E/P from 30 day (over 10 m) restoring to Levitus
  monthly climatology
• Restore SST to -2º under diagnosed ice
• Daily ECMWF winds, 1986--2000
• Restoring at lateral boundaries
• North of NA Sill, S Atlantic, Sicily Channel

21

• Domain
• Mercator grids (?x ?y, so ?y also scales w/
  cos(y))
• Mostly 0.1º res, but also 0.2 º, 0.4 º
• 40 (42) Levels
• 10 m in upper ocn, 250 m in deep, 5500 (6000) m
  max
• 20º S to 73º N
• Including Gulf of Mexico, Western Mediterranean
• Mixing (unless otherwise specified)
• Vertical
• Pacanowski and Philander Ri-number-based mixing
• Lateral viscosity and diffusion
• Biharmonic form

22
This had been the Gulf Stream problem
Sea Surface Height Variability
0.1º
0.2º
0.4º
Obs
23
Modern day obs (Iselin, 1936)
24
Modern day obs (Iselin, 1936)
Last Glacial Max, from Robinson et al. 1995
25
Last Glacial Max, from Robinson et al. 1995
Climate model simulation (CCSM.3)
26
Some improvement as resolution had increased
Sea Surface Height Variability
0.1º
0.2º
0.4º
Obs
27
0.1º
0.2º
Note also that the North Atlantic Current formed
at 0.1º
0.4º
Obs
28
Increase in heat transport with resolution

• Max heat transport between 25 and 30ºN, all cases
• 50 rise from lowest to highest
• (all without GM)

Dotted for 0.4º dashed for 0.2º solid for 0.1º
(BryanSmith 98)
29
Model Gulf Stream core at 10 day intervals, along
with Synop array observations
Stream separates at Cape Hatteras, if too
zonally, and with excess variability just off of
the Cape
30
Gulf Stream in stream coordinates, at 68º W
Observed data from Johns et al (1995)
Model Gulf Stream
31
Azores current in POP

• Upper panel is time-mean SSH, lower is
  instantaneous
• The variability not small in comparison with the
  mean!

32
Front had been weak at 0.28º
0.1º North Atlantic 0.28º Global
33

• Deep western boundary current off of Abaco was
  stronger, narrower and deeper than in previous
  simulation
• Transport compared well with obs of Lee et al.
  (1996)

34
First fully global model at 0.1
Note Masumoto et al. 2004 in J. of Earth
Simulator reported on nearly global 0.1
simulation
35
But 0.1º resolution not necessarily
enough. GS/NAC was not so good in first
fully-global 0.1º simulation. (Note, also,
Kuroshio meander).
From Maltrud and McClean, Ocean Modelling, 2005.
36
Whats preventing Gulf Stream System from forming?

• In terms of GS separation
• In real world, separation point is fixed at
  Hatteras.
• Almost certainly due to topography (what else is
  so unchanging?)
• In models
• Controlling influence of topography would appear
  to be easily overwhelmed by bias(es).
• Certainly so if Rossby radius not adequately
  resolved.
• Evidently so with slightly non-optimal
  configuration even at Order(0.1º)

37

• Fig 1 from Smith, Maltrud, Bryan and Hecht, JPO
  2000
• 1rst Rossby Radius ranges from gt100km to 5km, in
  zonal avg
• Stammer and Boning 1996 had reported size of
  eddies as 1.7R1 86km, from altimetry

38

• Rapid rise of Eddy Kinetic Energy with resolution
• Cases from Smith et al. 2000, Bryan and Smith
  1998
• This plot from Hecht and Smith 2008

Resolution
39
Perhaps some insight on poor Gulf Stream System
here?
Here, sensitivity is primarily with respect to
lateral dissipation.
40
Sensitivity to lateral dissipation?scaling of
biharmonic coeficients

• As standard practice we scale horizontal mixing
  coefs with (area)3/2 so as to maintain constant
  grid-Reynolds (-Peclet) number
• Grid-Re U (area/area0)3/2/?, for biharmonic
  form
• Scaled within simulation, as grid cell area
  varies, and between simulation, by default
• In these experiments we also include an
  adjustable factor C so that
• ? C ?0 (area/area0)3/2
• ?0-2.7x1010m4/s,
• area0(11.2 km)2 (equatorial grid cell area)
• C1 for standard set of simulations

41
8 runs at 3 resolutions, 5 different levels of
horizontal dissipation (5 values of C)

• C ?0 (area/area0)3/2
• where C1/8, 1/4, 1/2, 1 or 8.

42
Basin Mean Kinetic Energy
50
C1/4
40
C1/2
SMBH
C1
g(cm/s)2
0.2º, C1/8
same dissipative parameters
C8
same dissipative parameters
0.2º, C1
10
0.4º, C1
0
2000
1996
1992
1988

• 15 year integration gives reasonable
  equilibration of KE
• Time series of three 0.1º cases shown
• Averages over last 3 years shown for other cases

43
Gulf Stream paths from intersection of 12C
isotherm and 400m, 1998-2000
C1, 0.1º
Decrease dissipation

• By this point in time (years 13-15) the more
  viscous (C1) case has a more realistic
  separation (obs in green, model mean, 1 ? and
  extreme envelope in blue)

C1/2, 0.1º
C1/4, 0.1º
44
Decrease dissipation
45
C1
C1/2
Decrease dissipation
C1/4
Sea Surface Height Variability (1998-2000)
46
C1
C1
C1/2
Decrease dissipation
C1/4
C1/4
PT at 729 m
Sea Surface Height Variability (1998-2000)
47
C1/4 (lowest dissipation)
-20.0
32.5

• Peak velocities in NAC match pretty well
• but see how 10 cm/s isotach is at 1600 m in
  model, 3500 m in obs

48
Eddy Kinetic Energies, 55º W

• Deeper penetration of NA Current in less viscous
  case
• Ozgokmen (97) found jet needed to be highly
  inertial with low eddy activity to separate and
  cross f/H, in process study.
• We dont see relation between low eddy activity
  and separation at Hatteras
• but maybe high eddy activity for reattachment at
  Grand Banks
• and Ozgokmen may be correct, even at Hatteras
  Deep Western Boundary Current provides greater
  isolation of GS from topography in less viscous
  cases

C8
Decrease dissipation
C1
C1/4
49
Findings

• With basic biharmonic form of horizontal
  dissipation in POP at 0.1º, best compromise
  between GS separation and NW Corner, Azores
  Current is in moderate viscosity (moderate C1/2)
  case.

50
Findings

• With basic biharmonic form of horizontal
  dissipation in POP at 0.1º, best compromise
  between GS separation and NW Corner, Azores
  Current is in moderate viscosity (moderate C1/2)
  case.
• Some evidence for convergence at undesirably high
  level of dissipation.
• As expected, not converged at more desirable
  level of dissipation

51
Findings

• With basic biharmonic form of horizontal
  dissipation in POP at 0.1º, best compromise
  between GS separation and NW Corner, Azores
  Current is in moderate viscosity (moderate C1/2)
  case.
• Some evidence for convergence at undesirably high
  level of dissipation.
• As expected, not converged at more desirable
  level of dissipation
• Deep penetration of Stream, KE, appears to be
  required as Gulf Stream Extensions encounters
  topography
• At the Southeast Newfoundland Rise (where NAC and
  Azores branch off)

52
Findings

• With basic biharmonic form of horizontal
  dissipation in POP at 0.1º, best compromise
  between GS separation and NW Corner, Azores
  Current is in moderate viscosity (moderate C1/2)
  case.

This was with biharmonic forms of viscosity and
diffusion. Possible also to find this good
balance between GS sep and strong penetration of
GSExt with joint use of Laplacian and Biharmonic
Chassignet and Garraffo (2001) Hecht et al.
(2008).
53
NA sector of global

• Just lowering dissipation in 0.1 did not allow
  for good GS separation.
• Sensitivity to domain
• Lateral boundary conditions pose a constraint in
  regional models
• Are these lateral boundary conditions enough to
  explain the poorer behavior in global domain?

54
But 0.1º resolution not necessarily
enough. GS/NAC was not so good in first
fully-global 0.1º simulation. (Note, also,
Kuroshio meander).
From Maltrud and McClean, Ocean Modelling, 2005.
55
Cutting North Atlantic sector out of global
grid, in order to diagnose role of lateral
boundary conditions. Note lat/lon in S.
hemisphere, but distorted in North.
56
Sector

• Lateral boundary conditions sufficient to
  explain inability to produce Gulf Stream
  Separation at Hatteras, lack of penetration of N.
  Atlantic Current into NW Corner.
• Nice that we understand this, so whats one to do?

Global
SSH
57
Other sensitivities

• Anisotropic Gent-McWilliams Parameterization for
  Ocean Models
• As suggested by Roberts and Marshall (JPO, 1998),
  adiabatic closure beneficial, even at 0.1º
  resolution
• Anisotropic GM shown to allow energy levels to
  remain high, along with use of anisotropic
  horizontal viscosity, in Smith and Gent, JPO 2004
• Treatment of topography
• partial bottom cells
• Moderate smoothing of topography
• First in North Atlantic,
• then more limited testing in global

58
Partial Bottom Cellsvertical cross section
through the grid
See also Adcroft et al. 1997, MWR.
59
(No Transcript)
60
Full cells Partial bottom cells
61
Full cells Partial bottom cells Add anisotropic
dissipation
62
Full cells Partial bottom cells Add anisotropic
dissipation Add smoothing of topography
63
Sea Surface Height
Full cells
Partial bottom cells
Add smoothing of topography
Add anisotropic dissipation
64
SSH Variability
Full cells
Partial bottom cells
AVISO
Add smoothing of topography
Add anisotropic viscosity
65
2-3º Temp Transport
full
pbc
aniso
smooth
66
Improved global model

• Tripole grid
• Same number of points (3600x2400), but more
  uniform gridding of the ocean
• Smoother representation of bottom topography
• Partial bottom cells

67
Dipole Grid (puts pole in Greenland, resolution
focused around Greenland)
Tripole Grid (relatively uniform resolution)
68
Improved circulation in newer version of 0.1º POP
(SSH field)
older dipole grid configuration
newer tripole grid configuration
TOPEX/ERS1 altimetry, years 1994-2001 (courtesy
of AVISO)
69
North Atlantic Current turns north around Grand
Banks
older dipole
newer tripole
TOPEX/ERS1 (AVISO)
70
Brazil-Malvinas Confluence and Zappiola
Anticyclone improved
older dipole
newer tripole
TOPEX/ERS1 (AVISO)
71
snapshot of relative vorticity at 15m depth from
0.1 tripole
72
surface chlorophyl concentration after 2 years of
simulation using the 24-component ecosystem model
of Moore etal.
73
cfc11 concentration (after several decades) on a
plane sloping from near surface in the north to
the abyss in the south.
74
In strongly-eddying models

• Some real successes
• GS/NAC speed, depth-profile, path are far better
• Similar improvements in Kuroshio, Agulhas, ACC,
  etc.
• Room for improvement
• Kuroshio may set up standing meander
• Agulhas eddies may follow too regular a path
• GS/NAC still sensitive to model setup

75
Does it matter?

• Will these less adequately modeled features
  matter more to 21rst Century response than to
  20th Century control?

76
Hallberg and Gnanadesikan showed that
parameterized eddies dont provide correct
response
77
Does the path taken by warm, salty North Atlantic
waters matter
to stability of poleward circulation?
IPCC Assessment Report III.
78
Department of Energy supporting climate modeling
with a strongly eddying ocean

• Using the tripolar grid version of POP, as shown
  above, to be shown in detail later today
• Again, in collaboration with scientists outside
  of DOE, particularly with NCAR

79
And ocean response to fresh water fluxes from
Greenland

• Weijer, Maltrud and Hecht (Los Alamos), Henk
  Dijkstra and Michael Kliphuis (The Netherlands)

80
Resolution also matters to vertical mixing
response to wind forcing
Wind energy input to oceanic near-inertial
motions comparable to the work done by the
steady largescale winds on the general
circulation of the ocean (this comment based,
in turn, on research of Wunsch Ferrari, and
Alford)
81
Surface vorticity field (top) and wind energy
input (bottom)
82

• Vertical velocity as a function of depth, in 0.1º
  POP NA model, with 12-hourly wind forcing and
  daily averaged forcing

83
Same plot from POP 0.1º on left, from higher
resolution idealized Earth Simulator run on
right. Note strong secondary maximum at a few
hundred meters.
84
POP-alpha for climate
85
The test problem Idealization of Antarctic
Circumpolar Current
up
solid boundary
N
12ºC
surface thermal forcing
periodic bndry
zonal wind
E
periodic bndry
solid boundary
2ºC
deep-sea ridge
test problem invokes baroclinic instability.
86
low resolution 0.8º
0.4º
0.2º
high resolution 0.1º
note SST and thus heat flux is main influence
on atmosphere for climate
87
Test Problem Results Baroclinic Instability
Vertical temperature profile
6C isotherm
0.1(high res)
0.2
0.4
0.8(low res)
0.2
0.1
0.4
0.8(low res)
88
Lagrangian-Averaged Navier-Stokes Equation
(LANS-?)
Two ways to take averages Lagrangian and Eulerian
Lagrangian averaged trajectory
Lagrangian averaged velocity
actual velocity
particle trajectory
Eulerian averaged velocity, u, is average at xo
Lagrangian averaged velocity
u Eulerian averaged velocity
rough
smooth
Helmholtz operator
advection
Coriolis
pressure gradient
extra nonlinear term
diffusion
89
Test Problem Results Baroclinic Instability
Vertical temperature profile
6C isotherm
0.1(high res)
0.2
0.4
0.8(low res)
0.2
0.1
0.4
0.8(low res)
90
Test Problem Results Baroclinic Instability
Vertical temperature profile
6C isotherm
0.1(high res)
0.2 POP-?
0.2
0.4 POP-?
0.4
0.8 POP-?
0.8(low res)
0.2
0.1
0.4
0.8(low res)
0.2 POP-?
0.4 POP-?
0.8 POP-?
91
More model development
92
POP scaling to 10,000 processors, on Cray
Pat Worley, Oak Ridge National Lab
93
POP scaling beyond 10,000 processors, on blue
gene architecture
John Dennis, NCAR
94
Implicit POP

• Explicit methods (POPs leap frog scheme)
• Time step limited by stability concerns (CFL,
  etc.)
• Momentum eqns are split
• Implicit methods (Crank-Nicholson, Backward
  Euler, etc.)
• Unconditionally stable, time step chosen to
  resolve relevant processes
• No barotropic/baroclinic mode splitting
• Demonstrated capability to determine THC
  stability as a function of forcing, at lower
  resolution
• Parameter continuation
• Can it be made efficient enough to also be
  advantageous when the flow is variable?

95
imPOP

• How?
• Jacobian-Free Newton-Krylov (JFNK) methods
• Tailored preconditioning
• Couple POP to scalable and portable solver
  packages Trilinos
• Who?
• Wilbert Weijer (LANL)
• with contributions from Kate Evans (Oak Ridge),
  Erik Bernsen and Henk Dijkstra (Utrech), Jonas
  Thies and Fred Wubs (Groningen), and Trilinos
  group (Sandia)

96
test of imPOP code infrastructure and
preconditioning

• Wind-forced periodic channel with bump
• Leap frog (10 min)
• Crank Nicolson (10 hr)

97
Ice sheets in climate modeling -- a community
effort (including UK)
98
Conformal, Variable-Resolution Meshes
A conformal mesh is a mesh with no hanging nodes.

• Developing algorithms that target conformal
  meshes, for ocean, atmosphere and glacial ice
• Arbitrary Voronoi tessellations and Delaunay
  triangulations
• as well as Lat/Lon and conformally-mapped cubed
  sphere.

99
Began with Todd Ringler, LANL Max Gunzburger,
Florida State U. Lili Ju, U.of South
Carolina Now with Joe Klemp, NCAR Phil Jones,
LANL Jamal Mohd-Yosuf, LANL Bill Skamarock,
NCAR John Thuburn, Exeter U.
Research Topics Mesh Generation Numerical
Methods Analysis of Large-Scale Dynamics Analysis
of Cloud-Resolving Scales of Motion Implementation
on hybrid CPU/GPU architectures
100
Modus Operandi of Climate Modeling
There is absolutely nothing wrong with this
approach. Clearly it works. The fact that it
works should not preclude other complementary
approaches.
101
havent abandoned modis operandi -- RoadRunner

• IBMs heterogeneous architecture
• AMD Opterons cell processors
• Game rendering technology
• First to reach a petaflop
• Successful port of DNS code
• Attempt now to apply that experience to POP

102
AGU volume -- First collection of papers on ocean
modeling in an eddying regime
103
Parting thoughts

• Resolution can solve a lot of problems.

104
Parting thoughts

• Resolution can solve a lot of problems.
• You can never have enough resolution

105
Parting thoughts

• Resolution can solve a lot of problems.
• You can never have enough resolution
• So lets do keep thinking about smart
  parameterizations, techniques

106
Parting thoughts

• Resolution can solve a lot of problems.
• You can never have enough resolution
• So lets do keep thinking about smart
  parameterizations, techniques
• And well also see if we can get our codes to run
  on the next biggest computer.