Are oscillations in gauge data unavoidable?

[tovogt]

For testing purposes, I run GeoClaw with a synthetic tropical cyclone over open waters (constant bathymetry at -1000, no obstacles). It seems like even this simplistic example generates high oscillations in gauge output: Is this kind of picture unavoidable from a mathematical (numerical) point of view? I also experience these kinds of oscillations in GeoClaw runs with real-world data and it doesn't seem to be very realistic.

Here are the parameters (data files) I use to run this example: adjoint.data amr.data claw.data dtopo.data fgmax_grids.data fixed_grids.data flagregions.data friction.data gauges.data geoclaw.data multilayer.data qinit.data refinement.data regions.data surge.data topo.data

[mandli]

That's pretty surprising data! Can you provide the setrun.py and some of the 2D plots by change that you have made for this? The plots with the storm track and wind field relative to the entire domain would be particularly interesting to see.

[tovogt]

I tried to reduce the example as much as possible in order to track the reason down. Now I get the following gauge output: For that, I use the following set of files (for the tide gauge plot I included the script plot.py): example.zip Can you reproduce this output on your machine?

As for the storm track and wind field: The track runs straight through the rectangular domain from south to north with constant wind speed (18 m/s), pressure (985 hPa) and size. Here is a plot (orange: domain boundary, blue: eye_location, black: storm_radius):

[mjberger]

By any chance are you gauges right on a grid line, so that sometimes it's in one cell and sometimes another, or at the edge of a patch, where only piecewise constant interpolation is done instead of linear?

I know it's a long shot, but worth asking,

— Marsha

On Jul 15, 2021, at 4:02 AM, Thomas Vogt @.***> wrote:

I tried to reduce the example as much as possible in order to track the reason down. Now I get the following gauge output: https://user-images.githubusercontent.com/57705593/125749357-0b86f4ee-8bab-4bdd-817e-b2960443c3f4.png For that, I use the following set of files (for the tide gauge plot I included the script plot.py): example.zip https://github.com/clawpack/geoclaw/files/6821280/example.zip Can you reproduce this output on your machine?

As for the storm track and wind field: The track runs straight through the rectangular domain from south to north with constant wind speed (18 m/s), pressure (985 hPa) and size. Here is a plot (orange: domain boundary, blue: eye_location, black: storm_radius): https://user-images.githubusercontent.com/57705593/125751655-36e8c167-7225-4ec0-a54f-2b075da5e29b.png — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/clawpack/geoclaw/issues/521#issuecomment-880484925, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAGUGC5HOOXDACOYLXEJPI3TX2I2RANCNFSM5ALX5UTA.

[tovogt]

@mjberger Thanks for the suggestion! Unfortunately, that doesn't seem to be the case. I wiggled the gauge coordinates a bit using np.random.rand and the result doesn't change visually.

[mandli]

While have not seen the code produce oscillations like that the code definitely does not react well when a storm is not contained well away from domain boundaries. Can you try and expand the domain by at least a factor of 3 and see if the issue goes away or see if expanding it to some amount will remove the issue?

[tovogt]

I expanded the domain by a factor of 3 to make sure that the storm is contained well (blue dots are gauge locations): There are still oscillations, but notice how the absolute values on the y-axis shrinked by a factor of 3, maximum surge height is now down to 0.2 m: The result is similar when increasing storm intensity by factor 3, or when expanding the domain by another factor 2.

The code to reproduce these results is in this gist: https://gist.github.com/tovogt/b9bcd48c3fc46556b98cbc0ffb4f50e9 And I will update the gist when you have further suggestions about what to try.

[mandli]

Do you have 2D fields available for the surface and winds? Also are you ramping up the storm in any way? My guess is that you are getting waves bouncing around in the domain as you put a storm down "hard" onto an ocean at rest and the boundary conditions are not quite able to handle the resulting waves. I would note that after about two days those oscillations are being damped to what looks like a steady state negative value.

[tovogt]

It's true that I'm not ramping up the storm intensity in the above examples: The storm's pressure and wind speed are constant from the beginning until the end. However, the solver can solve this kind of problem within 1-2 minutes. Compared to that, when ramping up the intensity from 0, the solver will take hours. I wanted to create some toy examples for testing purposes (the bug #514 has been discovered using such a toy example) and it seems exaggerated to have such a long run time just for testing basic functionality, right?

Here are the gauge plots when ramping up the storm intensity from 0 (and ramping down to 0 at the end): Notice the surge heights of several hundred meters (!).

[tovogt]

Here is a plot from a single tide gauge from the ramping-up phase of the storm: And here is a plot from the intermediate stage, zoomed-in so that oscillations at the centimeter-scale become visible: I can very well imagine that this kind of behavior over open water is perfectly fine. I just want to be sure that it's not a result of any misconfiguration on my side.

[mandli]

That blow up does not seem good. I would again guess that was the storm getting too close to the boundaries but without seeing the 2D surfaces it is hard to tell. I am going to try and recreate your toy example myself and see what I can see. Going through it and putting in the parameters myself might make clear what the issue is or illuminate the issue you are seeing.

[tovogt]

I'm sorry for not posting the surface plots right away. So far, I'm not very experienced with clawpack's plotting functionality, but used the fgmax and gauge output only. Here are the first 10 frames of the surface plot: Note that the colorbar is ranging from -50 till +50 meters (!). This is with a ramping up of storm intensity from 0 at the beginning. Note that this is not a strong storm, but at 18 m/s barely crosses the threshold of being called a "tropical storm".

[mandli]

Wow, there's all sorts of problems it looks like with that. The waves for one should not be that big to begin with and those 4 vortices are also pretty interesting. My guess at this point is that the synthetic storm you generated is much too strong.

[tovogt]

According to the storm definition in the input file track.storm, the wind speed should never exceed 18 m/s. Translational speed is at 5 m/s so that shouldn't make too much of a difference. However, the wind field plots show that the wind speed starts at 0 in the first frame, reaches almost 50 m/s (!) in the third frame and comes back to 18 m/s after that: Maybe there is a bug in the Holland 1980 wind model for low mws? The pressure field is well-contained between 985 and 1013 mbar:

[mandli]

That pulse is pretty interesting. I will have to see what seems to be generating it.

[mandli]

So the pulse is coming entirely from the y-velocity of the wind field, which never actually becomes symmetric. My guess at this point is there is a bug in the code that added southern hemisphere storms.

[mandli]

[tovogt]

Well, that's easy to check. I have rerun the example on commit 01c9a8e07cbb65b6bb21d8e5f79a09abecb1c4ae (v5.8.0) and can report that I get the same high wind speeds, only the asymmetry is reversed, so the southern hemisphere correction is clearly not the cause of this behavior:

The surge height still amounts to more than 50 meters:

Note that the wind field remains asymmetric the whole time because the wind speeds are very low, allowing the storm's translational speed to play a dominant role.

[mandli]

It definitely has something to do with the adjustment to the wind due to translational speed of the storm. I am not sure why it's only causing this for this storm though. What is the effective forward velocity of this storm?

[tovogt]

I agree, it definitely has something to do with the translational speed adjustment. The translational velocity is 0.18 degrees (latitude) per hour, which is effectively 6 meters per second. The maximum wind speed of the storm ramps up from 0 to 18 meters per second within the first 8 hours. There seems to be some kind of singularity at the point where the translational speed equals the maximum wind speed in the ramp-up phase. Unfortunately, I won't have time to look into the details until next week or so.

[tovogt]

The problem is in these lines: https://github.com/clawpack/geoclaw/blob/5f94c15b9366dd6226e12a8f5dd6c497b2d347e9/src/2d/shallow/surge/model_storm_module.f90#L456-L458 Before dividing by mod_mws, we need to make sure that it's not only greater zero, but also that it's not close to zero. A quick fix is to replace mod_mws > 0 by mod_mws > 1. This means that the difference between translational speed and max. wind speed is at least 1 m/s.

However, I would even question whether this is really an appropriate way to adjust the influence of the translational velocity. Can we find out on which paper or study this method is based? I would rather use something like the "absorbing factor" in this paper: https://unepgrid.ch/en/resource/19B7D302 (Fig. 7)

Edit: The division by mod_mws has been introduced by @bolliger32 in b7111744cb14831698b7c7fab235c300a1059764 two years ago.

[tovogt]

I checked the historical storm tracks in IBTrACS and found that, since 1980, there are 342 storms for which the translational speed equals the maximum wind speed at some point in time. For all of those 342 storms, GeoClaw will produce the unrealistically high surge heights discussed in this thread, unless the storm tracks are truncated somehow before feeding into GeoClaw. This is really not an artificial test case, but it potentially affects the work with real-world input data.

Furthermore, this only explains the issue I ran into once I started to ramp up the storm intensity from 0 (see https://github.com/clawpack/geoclaw/issues/521#issuecomment-882794318). The original issue at the beginning of this thread is not related to this.

[bolliger32]

@tovogt Thanks for the mention and apologies for this bug. If I remember correctly, this adjustment approach is based on the ADCIRC implementation. I would imagine there are better ways of doing this if we wanted to explore those. That was one reason for separating this "post-processing" functionality into its own subroutine.

IIRC, the intent here was that wind was always going to be < mod_mws. If this is not happening and we're getting high translational speeds, I think there is probably something I didn't implement correctly upstream of the application of these lines. I'll take a look

[bolliger32]

I think the wind calculation is implemented correctly according to the holland equations (here's a ADCIRC powerpoint about their implementation about it). I think the issue might be that mod_mws is related to the holland B parameter assuming that at the RMW point, the coriolis force is negligable. With smallmod_mws, this is probably not true and thus this Holland equation is giving some estimated wind values that are higher thanmod_mws? One approach might be to find the max of this field, determine if it is greater than mod_mws and if so, scale the entire wind field by this ratio such that max(wind) = mod_mws?

[tovogt]

Thanks for your quick response, @bolliger32!

Here is the current wind model used by ADCIRC: https://wiki.adcirc.org/wiki/Generalized_Asymmetric_Holland_Model It's still loosely based on Holland 1980 - moreover, it includes comments about handling of translational velocity:

Occasionally, we have to deal with situations where Vmax < Vr, which violate (13) so Rmax couldn't be calculated. These situations mostly happen in the right hand quadrants (in the Northern Atmosphere) of a weak storm with a relatively high translational speed. For cases like this, we assign Vmax = Vr, which is equivalent to assigning Rmax = Rr.

That's basically what you suggested in your last post.

[mandli]

It's frustrating that things like this are not better documented. I guess they are there but just kind of hidden...

Anyway, we need to decide what the best way to handle this is. Suggestions?

[bolliger32]

Perhaps scale such that max(V_r) = mod_mws always? And/or then also implement this GAHM model as well, which seems to do better with low rossby number systems (i.e. those that seem to be causing these problems)?

[tovogt]

As already mentioned by @bolliger32, the value of wind should (from a mathematical point of view) never be larger than mod_mws, so this bug should never happen - unless there is a bug in the computation of wind further upstream. Apparently, we didn't look far enough upstream. There is an error in the computation of wind: It's in the computation of the coriolis force term in the file geoclaw_module.f90, a file that I never looked into at any occasion before. Here is the fix: #523

[bolliger32]

@tovogt good find! I don't think changing that function is the right way to go b/c the coriolis force SHOULD be signed depending on hemisphere, right? Maybe we just need to just take abs(f) in this line? And then make sure that the sign of the coriolis force is properly handled in other vortex functions?: https://github.com/clawpack/geoclaw/blob/5f94c15b9366dd6226e12a8f5dd6c497b2d347e9/src/2d/shallow/surge/model_storm_module.f90#L596

[tovogt]

@bolliger32 I agree. The parameter is usually defined to be signed. However, this is a good opportunity to clean up all the other wind field models - because, actually, only Holland 1980 and CLE do depend on the Coriolis parameter. All the others never made use of it.

Are you familiar enough with the CLE model to say for sure whether or not the signed Coriolis parameter is to be used there?

[bolliger32]

@tovogt not 100% positive but I think the same should apply here (i.e. it should be abs(coriolis)). Probably could use someone else to confirm however.

Maybe if we do want to change the function itself we just make a new one called unsigned_coriolis() which just equals abs(coriolis())? Just in case a future wind model comes along that wants to use the actual signed coriolis parameter. Just worried that we'll forget about this change :)

[tovogt]

@bolliger32 Good idea! I added a function unsigned_coriolis in my PR.

[tovogt]

This is what the result looks like after the fix: There are still oscillations, but much less and they look much more natural to me.

So, I think we can close this once we have merged #523

[mandli]

I looked at this as well and with these changes in the 2D fields looks good to me. I am wondering if these types of issues might be worthwhile documenting somewhere but maybe having them here on github is at least enough to have people aware of these issues. Great work guys!

[mandli]

Thanks to @tovogt for all your work on this and to @bolliger32 for also helping to close the loop on the causes!