Alaska Tsunami Sources for WA State modeling projects¶
The Statement of Work for the 2019 UW modeling study (Island County and Skagit River Delta) refers to using a model of the 1964 Mw 9.2 earthquake, and refers to several past studies. However, the past studies mentioned have also/instead used a model sometimes called AKmax, which not displaced somewhat from the 1964 event, and that has magnitude 9.15 or 9.2 depending on the version used.
This notebook summarizes our exploration of this source.
%pylab inline
import os,sys
from scipy.interpolate import griddata
from IPython.display import Image
from clawpack.geoclaw import topotools,dtopotools
Load topo and define background map¶
etopo1 = topotools.read_netcdf('etopo1', extent=[-160,-140,53,63], coarsen=3)
def fig_background():
figure(figsize=(12,5))
subplot(121)
contourf(etopo1.X, etopo1.Y, etopo1.Z, [0,10000], colors=[[.8,1,.8]])
contour(etopo1.X, etopo1.Y, etopo1.Z, [0], colors='g', lw=0.5)
plot([-149.9],[61.2],'ro')
text(-150.2,61.4,'Anchorage',color='r',fontsize=14)
plot([-147.339],[60.908],'k*',markersize=15)
axis([-158, -143, 54, 62])
gca().set_aspect(1/cos(57*pi/180))
grid(True)
In the figures below, the black star is the epicenter of the 1964 event, at 60.908°N 147.339°W according to https://earthquake.usgs.gov/earthquakes/eventpage/official19640328033616_30/origin/detail
fig_background()
Scope of Work statement¶
Here is the statement from the SOW:
Image('sow.png', width=650)
Reviewing the Witter et al. (2011) study of Bandon, OR, it seems that they only used the 15 T-shirt size CSZ realizations in that study and no Alaska events. They do however refer to the paper of Johnson et al. (1996), where a slip model for the Alaska 1964 event is presented. This paper is also referenced in the Allen et al. (2018) DOGAMI report, and is presumably the AK64 slip model discussed further at the end of this this document.
DOGAMI Report on tsunami modeling in the Columbia River¶
In Allan et al., 2018 two Alaska sources were used, called AK64 and AKMax, as described in the paragraphs below.
Image('DOGAMI_ColumbiaRiver.png', width=650)
Seaside OR paper¶
The Gonzalez et al. 2009 paper on PTHA for Seaside, OR has Table 1 below. Note that the slip values for Source 3 agree with the values listed in the DOGAMI report above, 15, 20, 25, and 30m.
Image('GonzalezGeist_Table1.png', width=800)
AASZ03 from original Crescent City PTHA study¶
For the Crescent City study of 2014 we used a source AASZ03 that was generated using the AASZ03.csv file printed below. Note that the penultimate column contains the slip values, which are 19.5, 24.5, 29.5, and 34.5. Each is 4.5m higher than the corresponding value from the references above. I am not sure where we got these files.
Longitude,Latitude,Depth(km),Length,Width,Strike,Dip,Rake,Slip3,UnitSrc
204.89500,55.97000,17.94,100.0,50.0,236.000,15.000,90.0,19.5000,acsza31
205.34000,55.59800,5.00,100.0,50.0,236.000,15.000,90.0,19.5000,acszb31
206.20800,56.47300,17.94,100.0,50.0,236.000,15.000,90.0,24.5000,acsza32
206.65800,56.10000,5.00,100.0,50.0,236.000,15.000,90.0,24.5000,acszb32
207.53700,56.97500,17.94,100.0,50.0,236.000,15.000,90.0,29.5000,acsza33
207.99300,56.60300,5.00,100.0,50.0,236.000,15.000,90.0,29.5000,acszb33
208.93710,57.51240,17.94,100.0,50.0,236.000,15.000,90.0,34.5000,acsza34
209.40000,57.14000,5.00,100.0,50.0,236.000,15.000,90.0,34.5000,acszb34
210.25970,58.04410,17.94,100.0,50.0,230.000,15.000,90.0,34.5000,acsza35
210.80000,57.70000,5.00,100.0,50.0,230.000,15.000,90.0,34.5000,acszb35
211.32490,58.65650,17.94,100.0,50.0,218.000,15.000,90.0,34.5000,acsza36
212.00000,58.38000,5.00,100.0,50.0,218.000,15.000,90.0,34.5000,acszb36It is possible the slips were adjusted to give a Mw 9.2 event, which this does:
sift_slip_large = {'acsza31':19.5,'acszb31':19.5, 'acsza32':24.5,'acszb32':24.5, \
'acsza33':29.5,'acszb33':29.5, 'acsza34':34.5,'acszb34':34.5, \
'acsza35':34.5,'acszb35':34.5, 'acsza36':34.5,'acszb36':34.5}
AKmax_sift_large = dtopotools.SiftFault(sift_slip_large)
print('Magnitude = %.2f' % AKmax_sift_large.Mw())
Confirm that the parameters agree with the csv file values:
print('Longitude Latitude Depth Strike Dip Rake Slip')
for s in AKmax_sift_large.subfaults:
print ('%.4f %.4f %4.1f %.1f %.1f %.1f %.1f' \
% (s.longitude, s.latitude, s.depth/1e3, s.strike, s.dip, s.rake, s.slip))
We shift the longitude by -360 for consistency with recent work:
AKmax_sift_large = dtopotools.SiftFault(sift_slip_large, longitude_shift=-360.)
Apply the Okada model to compute the seafloor deformation¶
We use the GeoClaw utility for generating seafloor deformation from NOAA SIFT unit sources.
x,y = AKmax_sift_large.create_dtopo_xy(dx=1/60., buffer_size=0.5)
print('Will create dtopo with shape %i by %i' % (len(x),len(y)))
dtopo_AKmax_sift_large = AKmax_sift_large.create_dtopography(x,y, times=[1.],verbose=True)
fig_background()
clines_neg = arange(-10,0.,1.)
clines_pos = arange(1,16,1.)
dtopo = dtopo_AKmax_sift_large
dZ = dtopo.dZ[-1,:,:]
contour(dtopo.X, dtopo.Y, dZ, clines_neg, colors='b', linestyles='-')
contour(dtopo.X, dtopo.Y, dZ, clines_pos, colors='r', linestyles='-')
dZmin = dZ.min()
dZmax = dZ.max()
title('AKmax SIFT large\n min(dZ) = %.1f, max(dZ) = %.1f' % (dZmin,dZmax))
print('Blue contours at ', [c for c in clines_neg if c>=dZmin])
print('Red contours at ', [c for c in clines_pos if c<=dZmax])
subplot(122)
j = find(dtopo.x < -150).max()
plot(dtopo.y, dZ[:,j],'b')
axis([55,62,-10,20])
xlabel('Latitude')
ylabel('meters')
title('Transect at Longitude -150')
grid(True);
Note that although the magnitude agrees with the 1964 event, the deformation is not centered around the epicenter of that event, shown by the black star in the figure above, at 60.908°N 147.339°W according to https://earthquake.usgs.gov/earthquakes/eventpage/official19640328033616_30/origin/detail
SIFT source with slip values from Table 1.¶
Recreating the source with the smaller slip values from Table 1 of the Seaside paper gives a similar seafloor deformation but with smaller maxima, and also a smaller magnitude Mw 9.15:
sift_slip = {'acsza31':15.,'acszb31':15., 'acsza32':20.,'acszb32':20., \
'acsza33':25.,'acszb33':25., 'acsza34':30.,'acszb34':30., \
'acsza35':30.,'acszb35':30., 'acsza36':30.,'acszb36':30.}
AKmax_sift = dtopotools.SiftFault(sift_slip, longitude_shift=-360.)
print('Magnitude = %.2f' % AKmax_sift.Mw())
#x,y = AKmax_sift.create_dtopo_xy(dx=1/60., buffer_size=0.5) # same as above
print('Will create dtopo with shape %i by %i' % (len(x),len(y)))
dtopo_AKmax_sift = AKmax_sift.create_dtopography(x,y, times=[1.],verbose=True)
fig_background()
clines_neg = arange(-10,0.,1.)
clines_pos = arange(1,16,1.)
dtopo = dtopo_AKmax_sift
dZ = dtopo.dZ[-1,:,:]
contour(dtopo.X, dtopo.Y, dZ, clines_neg, colors='b', linestyles='-')
contour(dtopo.X, dtopo.Y, dZ, clines_pos, colors='r', linestyles='-')
dZmin = dZ.min()
dZmax = dZ.max()
title('AKmax SIFT\n min(dZ) = %.1f, max(dZ) = %.1f' % (dZmin,dZmax))
print('Blue contours at ', [c for c in clines_neg if c>=dZmin])
print('Red contours at ', [c for c in clines_pos if c<=dZmax])
subplot(122)
j = find(dtopo.x < -150).max()
plot(dtopo.y, dZ[:,j],'b')
axis([55,62,-10,20])
xlabel('Latitude')
ylabel('meters')
title('Transect at Longitude -150')
grid(True);
AKmax.txt source¶
During work on the Crescent City project (ca. 2012) Frank purchased a CD that contained files AK64.txt and AKmax.txt in a directory CSZ_Bandon/def_grids. These are not in GeoClaw format and do not appear to have been used in our Crescent City modeling work. They are on a coarse grid with 15-20 arcsecond resolution. Plotting AKmax.txt shows it is similar to the above, with min and max of dZ (-6.2m and 13.4m respectively) that agree well with the dZ generated using the SIFT unit sources with the slip values from Table 1. (And hence correspond to the smaller Mw 9.15 event.)
fname = '/Users/rjl/git/CCptha/Sources/CSZ_Bandon/def_grids/AKmax.txt'
d = loadtxt(fname,skiprows=1)
data = d[:,1:] # ignore first column
data[:,0] = data[:,0] #+ 360
points = data[:,:2]
values = data[:,2]
xpoints = points[:,0]
ypoints = points[:,1]
# Interpolate onto finer grid for plotting:
mx = 600
my = 400
xlower = xpoints.min()
xupper = xpoints.max()
ylower = ypoints.min()
yupper = ypoints.max()
x = linspace(xlower,xupper,mx)
y = linspace(ylower,yupper,my)
X_AKmax,Y_AKmax = meshgrid(x,y)
dZ_AKmax = griddata(points, values, (X_AKmax,Y_AKmax), method='nearest', fill_value=0.)
fig_background()
clines_neg = arange(-10,0.,1.)
clines_pos = arange(1,15,1.)
contour(X_AKmax,Y_AKmax,dZ_AKmax, clines_neg, colors='b', linestyles='-')
contour(X_AKmax,Y_AKmax,dZ_AKmax, clines_pos, colors='r', linestyles='-')
dZmin = dZ_AKmax.min()
dZmax = dZ_AKmax.max()
title('AKmax\n min(dZ) = %.1f, max(dZ) = %.1f' % (dZmin,dZmax))
print('Blue contours at ', [c for c in clines_neg if c>=dZmin])
print('Red contours at ', [c for c in clines_pos if c<=dZmax])
subplot(122)
j = find(x < -150).max()
plot(y, dZ_AKmax[:,j],'b')
axis([55,62,-10,20])
xlabel('Latitude')
ylabel('meters')
title('Transect at Longitude -150')
grid(True);
AK64.txt source¶
The AK64.txt file (of unknown origin) contains a similarly coarse resolution deformation that is centered more around the 1964 epicenter.
fname = '/Users/rjl/git/CCptha/Sources/CSZ_Bandon/def_grids/AK64.txt'
d = loadtxt(fname,skiprows=1)
data = d[:,1:] # ignore first column
data[:,0] = data[:,0] #+ 360
points = data[:,:2]
values = data[:,2]
xpoints = points[:,0]
ypoints = points[:,1]
# Interpolate onto finer grid for plotting:
mx = 600
my = 400
xlower = xpoints.min()
xupper = xpoints.max()
ylower = ypoints.min()
yupper = ypoints.max()
x = linspace(xlower,xupper,mx)
y = linspace(ylower,yupper,my)
X_AK64,Y_AK64 = meshgrid(x,y)
dZ_AK64 = griddata(points, values, (X_AK64,Y_AK64), method='nearest', fill_value=0.)
fig_background()
clines_neg = arange(-10,0.,1.)
clines_pos = arange(1,15,1.)
contour(X_AK64,Y_AK64,dZ_AK64, clines_neg, colors='b', linestyles='-')
contour(X_AK64,Y_AK64,dZ_AK64, clines_pos, colors='r', linestyles='-')
#axis([-155, -140, 55, 61])
dZmin = dZ_AK64.min()
dZmax = dZ_AK64.max()
title('AK64\n min(dZ) = %.1f, max(dZ) = %.1f' % (dZmin,dZmax))
print('Blue contours at ', [c for c in clines_neg if c>=dZmin])
print('Red contours at ', [c for c in clines_pos if c<=dZmax])
subplot(122)
j = find(x < -150).max()
plot(y, dZ_AK64[:,j],'b')
axis([55,62,-10,20])
xlabel('Latitude')
ylabel('meters')
title('Transect at Longitude -150')
grid(True);
Johnson, Satake, et al. 1996¶
The paper on the Bandon, OR study by Witter et al. refers to this paper for the Alaska 1964 source,
The 1964 Prince William Sound earthquake: Joint inversion of tsunami and geodetic data, by Jean M. Johnson Kenji Satake Sanford R. Holdahl Jeanne Sauber, JGR 1996. DOI 10.1029/95JB02806
The figure below from this paper shows subfaults that may have been used to derive the source above.
This paper has a table of subfault parameters so this could be used to generate a better version of the sea floor deformation.
Image('JohnsonSatakeFig9.png', width=400)
Source from PMEL tech report¶
Further discussion led us to consider the source constructed in
MODELING TSUNAMI INUNDATION IMPACT ON THE WASHINGTON COAST FROM DISTANT SEISMIC SOURCES
Christopher D. Chamberlin, Vasily V. Titov, and Diego Arcas
PMEL Tech report 2009. (Available online?)In this report a source consisting of 20 unit sources with uniform slip is constructed that gives the maximum amplitude on the Washington coast. It is stated that a uniform slip of 20m on each of 20 subfaults gives Mw 9.2. However, this is not consistent with this statement from the report:
"In the database, a single unit source, a 50 x 100 km region, with a scaling factor (slip) of 1.0 represents a Mw 7.5 earthquake."
We believe this is correct, but that one must place a uniform slip of roughly 17.7m on each subfault to obtain Mw 9.2. Here are our calculations, using
$$ M_w = \frac 2 3 (\log_{10}(M_o) - 9.05), $$
where the seismic moment $M_o$ is in N-m, using a rigidity $\mu = 4\times 10^{10}$ Pa or 40 GPa.
mu = 40e9 # rigidity in Pascals
L = 100e3 # length in meters
W = 50e3 # width in meters
s = 1. # unit slip 1 meter
Mo = L*W*s*mu # seismic moment in N-m
Mw = 2/3. * (log10(Mo) - 9.05)
print('Using mu = %i GPa, a unit source with slip 1m has Mo = %g and Mw = %g' % (mu/1e9,Mo,Mw))
This agrees with the paper, but then if we use $N=20$ subfaults with slip $s=20$m we get
N = 20
s = 20
print('This uniform slip %g on %i subfaults gives Mw = %g' \
%(s, N, 2/3. * (log10(N*s*Mo) - 9.05)))
We can solve for how much slip we need with $N$ unit sources:
Mw = 9.2
print('Uniform slip required for Mw = %g:' % Mw)
for N in [16,18,20,22,24]:
s = 10**(1.5*Mw + 9.05) / (Mo*N)
print(' for %s unit sources, s = %g' % (N,s))
These values agree with the values given in Table 1 of the Seaside paper, where sources AASZ01 and 2 have 20 sources and slip 17.7, while sources 4 and 5 have 24 sources and slip 14.8. (But their source 3 has the wrong seismic moment unless the slips are scaled upward as discussed above. The total of s over all unit sources should be 354 to give Mw 9.2 but it is only 300.)
Rigidity value?¶
Above we used a rigidity (shear modulus) of $\mu = $40 GPa. We believe PMEL uses the value 40 GPa in constructing deformation from unit sources.
However, in some work a value of 30 GPa has been used instead, e.g. for the Seattle Fault model we have been using from
S. Koshimura, H. O. Mofjeld, F. I. Gonzalez, and A. L. Moore, Modeling the 1100 BP paleotsunami in Puget Sound, Washington, Geophysical Research Letters, 29 (2002), p. 1948, DOI 10.1029/2002GL015170.
Note that reducing $\mu$ by 3/4 would require increasing the slip by 4/3 to get the same value of Mw. Re-doing the above calculations with 30 GPa gives:
mu = 30e9 # rigidity in Pascals
L = 100e3 # length in meters
W = 50e3 # width in meters
s = 1. # unit slip 1 meter
Mo = L*W*s*mu # seismic moment in N-m
Mw = 2/3. * (log10(Mo) - 9.05)
print('With mu = %i GPa, a unit source with slip 1m has Mo = %g and Mw = %g' % ((mu/1e9),Mo,Mw))
N = 20
s = 20
print('With mu = %i GPa, this uniform slip %g on %i subfaults gives Mw = %g' \
%((mu/1e9),s, N, 2/3. * (log10(N*s*Mo) - 9.05)))
Mw = 9.2
print('With mu = %i GPa, uniform slip required for Mw = %g:' % ((mu/1e9),Mw))
for N in [16,18,20,22,24]:
s = 10**(1.5*Mw + 9.05) / (Mo*N)
print(' for %s unit sources, s = %g' % (N,s))
This does not give values that agree any better with what is in some of the references, so we still believe 40 GPa is the correct value to use for CSZ.
Also note that the value of $\mu$ used comes into the calculation of Mw, but does not come into the Okada model except indirectly through the definition of the Poisson ratio. We set the Poisson ratio to $0.25$ in our Okada model and believe this is standard, so changing the rigidity $\mu$ for a given set of slip values affects the magnitude Mw calculated but not the sea floor deformation.
Create a fault for PMEL maximum scenario 1¶
Use slip 17.7m rather than 20m, and rigidity 40 GPa.
sift_slip_pmelmax = {}
slip_uniform = 17.7
for usrc in range(29,39):
k = 'acsza%s' % usrc
sift_slip_pmelmax[k] = slip_uniform
k = 'acszb%s' % usrc
sift_slip_pmelmax[k] = slip_uniform
AK_pmelmax = dtopotools.SiftFault(sift_slip_pmelmax, longitude_shift=-360.)
print('Magnitude = %.2f' % AK_pmelmax.Mw())
Check that the rigidity is set to 40 GPa:
AK_pmelmax.subfaults[0].mu/1e9
Create dtopo:¶
x,y = AK_pmelmax.create_dtopo_xy(dx=1/60., buffer_size=0.5)
print('Will create dtopo with shape %i by %i' % (len(x),len(y)))
dtopo_AK_pmelmax = AK_pmelmax.create_dtopography(x,y, times=[0.],verbose=True)
fig_background()
clines_neg = arange(-10,0.,1.)
clines_pos = arange(1,16,1.)
dtopo = dtopo_AK_pmelmax
dZ = dtopo.dZ[-1,:,:]
contour(dtopo.X, dtopo.Y, dZ, clines_neg, colors='b', linestyles='-')
contour(dtopo.X, dtopo.Y, dZ, clines_pos, colors='r', linestyles='-')
dZmin = dZ.min()
dZmax = dZ.max()
axis([-160,-143,52,62])
title('AK pmelmax\n min(dZ) = %.1f, max(dZ) = %.1f' % (dZmin,dZmax))
print('Blue contours at ', [c for c in clines_neg if c>=dZmin])
print('Red contours at ', [c for c in clines_pos if c<=dZmax])
subplot(122)
j = find(dtopo.x < -150).max()
plot(dtopo.y, dZ[:,j],'b')
axis([55,62,-10,20])
xlabel('Latitude')
ylabel('meters')
title('Transect at Longitude -150')
grid(True);
The plot on the left is similar to what is in the top plot of Figure 3 of the PMEL report.
We will call this AK_PMEL_2009_Scenario1_slip_177.
dtopo.write('AK_PMEL_2009_Scenario1_slip_177.tt3', dtopo_type=3)
Comments on other sources in the PMEL report¶
Scenario 3 in the PMEL report was created from source #2 of the Seaside paper. However this 2009 report did referred to the original 2006 Seaside technical report "Seaside, Oregon, Tsunami Pilot Study : Modernization of FEMA Flood Hazard Maps", available at https://pdxscholar.library.pdx.edu/geology_fac/16/, rather than to the 2009 Seaside paper. Table 1 of the technical report contained typos that were fixed in the 2009 paper. In the paper Source #2 is a Mw 9.2 created using 17.7 m of slip on each of 20 unit sources, which we believe is correct. The incorrect table the 2006 report stated that a slip of 18.1 m was used on each of 22 unit sources, which would give a Mw 9.234 event. This is of no concern to the present work since Scenario 3, even with this incorrect larger value, was found to have much smaller impact on the WA coast than the other scenarios considered due to the orientation of the fault.
Scenario 2 of the PMEL report was a model of 1964 generated by Tang et al (2006) by starting with the Johnson et al. model and with additional inversion of tide gauge records. It is stated that it has a maximum slip value of 38 m but the unit source solution is not provided in the 2009 report (. Diego later provided it as
4.000*ac33a+5.000*ac33z+28.000*ac34a+5.000*ac34x+37.000*ac34y+38.000*ac34z
+2.000*ac35a+3.000*ac35x+14.000*ac35y+29.000*ac35z+10.000*ac36x`.This is a different unit source labeling system than we are currently using.
This gives a Mw 9.0 event:
total_slip = 4+5+28+5+37+38+2+3+14+29+10
mu = 4e10 # rigidity in Pascals
L = 100e3 # length in meters
W = 50e3 # width in meters
Mo = L*W*total_slip*mu # seismic moment in N-m
Mw = 2/3. * (log10(Mo) - 9.05)
print('A total slip of %g has Mo = %g and Mw = %g' % (total_slip,Mo,Mw))
Preliminary Conclusion¶
What should be used for our Alaska Source? There are at least 5 possibilities:
- The Mw 9.2 AKmax deformation as defined above using the SIFT subfaults,
- The Mw 9.15 AKmax deformation as defined above using the SIFT subfaults,
- A fine resolution version of the 1964 source model from Johnson, Satake, et al. 1996,
- A better source model for the 1964 event developed more recently, e.g from PMEL,
- The PMEL Scenario 1 event with uniform slip 20 m.
- The PMEL Scenario 1 event with uniform slip 17.7 m.
We now think we should choose option #2 (which we call AASZ03) or #6, which we will refer to as AK_PMEL_2009_Scenario1_slip_177.
Comparison of GeoClaw results with different sources¶
It probably doesn't matter very much which we use. One fault is longer than the other and has less uniform slip as a result, but they are both in the same region. Once the slips are adjusted to have roughly the same magnitude, Mw 9.2, they give very similar tsunamis.
fig_background()
clines_neg = arange(-10,0.,1.)
clines_pos = arange(1,16,1.)
dtopo = dtopo_AKmax_sift_large
dZ = dtopo.dZ[-1,:,:]
contour(dtopo.X, dtopo.Y, dZ, clines_neg, colors='b', linestyles='-')
contour(dtopo.X, dtopo.Y, dZ, clines_pos, colors='r', linestyles='-')
dZmin = dZ.min()
dZmax = dZ.max()
axis([-160,-143,52,62])
title('AKmax SIFT large\n min(dZ) = %.1f, max(dZ) = %.1f' % (dZmin,dZmax))
print('Blue contours at ', [c for c in clines_neg if c>=dZmin])
print('Red contours at ', [c for c in clines_pos if c<=dZmax])
subplot(122)
j = find(dtopo.x < -150).max()
plot(dtopo.y, dZ[:,j],'b')
axis([55,62,-10,20])
xlabel('Latitude')
ylabel('meters')
title('Transect at Longitude -150')
grid(True);
fig_background()
clines_neg = arange(-10,0.,1.)
clines_pos = arange(1,16,1.)
dtopo = dtopo_AK_pmelmax
dZ = dtopo.dZ[-1,:,:]
contour(dtopo.X, dtopo.Y, dZ, clines_neg, colors='b', linestyles='-')
contour(dtopo.X, dtopo.Y, dZ, clines_pos, colors='r', linestyles='-')
dZmin = dZ.min()
dZmax = dZ.max()
axis([-160,-143,52,62])
title('AK pmelmax\n min(dZ) = %.1f, max(dZ) = %.1f' % (dZmin,dZmax))
print('Blue contours at ', [c for c in clines_neg if c>=dZmin])
print('Red contours at ', [c for c in clines_pos if c<=dZmax])
subplot(122)
j = find(dtopo.x < -150).max()
plot(dtopo.y, dZ[:,j],'b')
axis([55,62,-10,20])
xlabel('Latitude')
ylabel('meters')
title('Transect at Longitude -150')
grid(True);
Tsunami simulation results¶
Below are a few initial tests on a relatively coarse grid (6 arcminutes finest resolution), looking at several gauges as indicated in this figure:
from IPython.display import Image
Image(filename='AKplots/gauges.png', width=800)
