New Features in GeoClaw

This notebook describes some new features developed in GeoClaw in the course of performing a tsunami hazard analysis for Whidbey Island and the Skagit River Delta for maps to be produced by the Washington State Department of Natural Resources. See the IslandWhidbeyTHA_2019 project webpage for that report and archive of the code used.

Some of these features were first introduced in earlier work, in particular for Whatcom County [Report, 2018] but have been improved further more recently.

In particular:

  1. See pcolorcells.html for a version of pcolormesh that works better for finite volume grid cells, found in visclaw.plottools.pcolorcells and used to improve geoclaw.Topography.plot. Functions added to geoclaw.kmltools allow making better png files and associated kml files for viewing on Google Earth.

  2. The introduction of Ruled Rectangles as a special class of polygons that can be used in various ways to specify a spatial region of interest. In particular, refinement regions in GeoClaw can now be specified as ruled rectangles rather than as simple rectangles. For an introduction to this type of polygon and the Python class region_tools.RuledRectangle, see RuledRectangles.html for a description of how Ruled Rectangles can be generated and used, using the region_tools.py module. Some examples are shown in the figure below from the IslandWhidbeyTHA_2019 project webpage.

  3. The manner in which AMR refinement regions are specified has been modified for regions that are Ruled Rectangles. A proposal to generalize the way all regions are specified, and to rename them as flagregions, can be found in the notebook FlagRegions.html.

  4. A new input format has been introduced for fgmax regions that allows providing a file in the format of a ASCII raster DEM (topotype == 3 in GeoClaw, see clawpack.org/topo.html), but in which the Z value takes only the values 1 or 0 to indicate whether each point should be selected as an fgmax point or not. Such a file is much smaller than an explicit list of the $(x,y)$ coordinates of each fgmax point (in cases where a large fraction of points over a large region are selected). It also has other advantages in terms of speeding up the code internally.

  5. The way fgmax grids are updated in the GeoClaw Fortran code has been sped up dramatically for problems with many fgmax points. This capability was originally designed for problems with a relatively small number of fgmax points, and was not written very efficiently. For problems with millions of fgmax points these improvements are critical.

  6. The user can now specify a Fortran function set_eta_init that is used to initialize the water surface eta when new grids at a fine resolution are first introduced during a simulation. In particular, this is used for coastal regions where subsidence or uplift takes place during the earthquake (e.g. in parts of Puget Sound for a Seattle Fault event). As the coast and sea floor subside, for example, the water offshore should also subside so that initially the shoreline is in the same location but the water elevation eta (measured relative to the vertical datum of the DEMs, typically MHW) has decreased. If it takes some time for the tsunami to arrive at this location then the finest 1/3 arcsecond grid patches may not be introduced until some time into the simulation. When they are first introduced, the coastline is better resolved so that a single coarse cell that is entirely wet (or dry) is split into many cells ($r^2$ cells if refining by $r$ in each direction), some of which are below MHW and others above. The default behavior of GeoClaw is to initialize cells so that those with cell-averaged topography value $B < 0$ are filled with water of depth $h = -B$, while those with $B\geq 0$ are initialized as dry. In other words, $h = \max(-B,0)$. The new capability specifies an eta_init value $\eta_i$ in each grid cell and the depth is initialized to $h = \max(\eta_i -B, 0)$. For the application of subsidence/uplift due to an earthquake described above, a version of the set_eta_init routine is provided that computes $\eta_i$ in each cell by interpolating from the dtopo file provided by the user (which gives the earth surface deformation at each point on a uniform grid). For other applications the user could provide their own custom set_eta_init file, e.g. to initialize an onshore lake to a higher surface elevation than the offshore ocean or Sound.

  7. A marching front algorithm has been implemented in the Python module marching_front.py, described in MarchingFront.html, that can be used to preprocess a topography DEM file and create a mask indicating which points in this file satisfy certain elevation requirements along with connectivity to the coast. This was originally developed to identify points in the DEM where the topography value $Z$ is below MHW but that should be initialized as dry land because they are in regions protected by dikes or levies. In this situation, the marching algorithm is used by initializing points well offshore (e.g. with $Z < -5$ meters) as wet and other points to unset. Then the front between wet/unset points is advanced by marking neighboring points as dry if $Z\geq0$ or as wet if $Z<0$. This is repeated iteratively for each new front until there are no more wet points with unset neighbors. At this point any points still unset are entirely buffered by dry points and must lie behind dikes, so these are also set to dry. There are several subtleties to how this is actually applied, see ForceDryMask.html for more details and examples.

  8. The marching front algorithm described above is generalized in the Python function select_by_flooding that can also be useful in at least two other ways. One is to select fgmax points that satisfy certain criteria, e.g. that have elevation $Z<Z_{\max}$ for some specified elevation (say 15 m) and that are also connected to the coast by points below this elevation. This is useful in choosing a minimal set of fgmax points that cover all coastal points that could possibly flood in a given tsunami simulation. See MarchingFront.html and FGmaxGrids.html for examples.

  9. Another application of the marching front algorithm is to select all nearshore points that lie in some range of elevations (say $-1000 < Z < 20$) that covers the continental shelf or some other topographic feature where a finer grid resolution must be used. The selected points can then be surrounded by a Ruled Rectangle that covers these points, which can in turn be used as a flag region for guiding the adaptive mesh refinement in GeoClaw, as described in point 1 above.

  10. A new gridtools.py module has some tools for reading in a solution as a set of grid patches computed using AMR and extracting data on a uniform 2d grid or along a 1d transect, by interpolating from the finest level patch available at each point. See gridtools.html.

In [1]:
%matplotlib inline
In [2]:
from pylab import *
from IPython.display import Image

Motivating examples

The figures below show a few examples to motivate and better explain some of the new features described above. For more details and example code, see some of the other notebooks cited above.

In [3]:
Image('figs/fgmax_polygons_bbox.png', width=350)
Out[3]:
In [4]:
Image('figs/fgmax_polygons_pcolor_wSE_labels.png', width=400)
Out[4]:

The figure above shows the study region for the recent project modeling Whidbey Island and the Skagit River Delta, as the union of the six polygons plotted and labelled. Each of the polygons is a Ruled Rectangle covering some portion of the coast and associated offshore water. Initially dry land is green (if above MHW) or pink (if below MHW). Note the following:

  • The Skagit region contains more than 150 km$^2$ of dry land below MHW that is protected by more than 200 km of dikes and levies. Initializing all of this pink area to be dry rather than wet is critical in order to properly model the extent of actual tsunami inundation. Several bays on Whidbey Island also exhibit regions of dry land below MHW protected by dikes.

  • The original study region was split into 6 subregions to limit the number of fgmax points (and the extent of the required 1/3" computational grids) in any single run. The Skagit region is the largest, where it was necessary to model this entire region in a single run to capture inundation properly. This region contains roughly 5.3 million fgmax points. Onshore points within each polygon are included in the set of fgmax points only if they are sufficiently low elevation and/or are close to the coast. See MarchingFront.html.

  • The polygons shown in the figure were designed so that they generally contain points on only one side of Whidbey Island (e.g. the polygon labelled E in the figure contains only points on the east side of Whidbey Island and the west side of Camano Island. The use of Ruled Rectangles to define these regions made this relatively easy. The marching front algorithm was then applied to choose points from within these polygons as the actual fgmax points. See ?? for an example.

Back to Index.html or go to the next notebook, pcolorcells.html.