Thunderstorm computer simulations

[MClarkson]

I have been working on a project to create a fairly high resolution 2 dimensional model covering a domain that is roughly the same size as a common set of thunderstorm cells. The 2 dimensional nature should allow me to capture most thunderstorm structure, but I will not be able to resolve 3d structure such as mesocyclones. Each pixel in the following image represents a 50m by 50m grid cell, of which there are 400 in the vertical and 800 in the horizontal (20km high by 40km wide at the moment.) The model is fully compressible, non-hydrostatic, and accounts for warm phase precipitation only (ice will be added soon). The equations of motions are solved using forward upstream differencing on a staggered c-grid. The time step is 0.05 seconds in order to remain numerically stable.

My first full resolution run that stayed stable just completed! This run took 18 hours on a single computer core, completing about 50 trillion calculations. The model itself runs in C++ and I created the graphics in java. Anyway, here is an animation of the cloud field. The simulated time between frames is 75 seconds. The atmosphere was initialized based on a very moist, unstable, and low vertical shear tropical sounding. I then triggered convection with a small warm bubble. I'll post some other fields in a bit. If you guys see anything that looks wrong, especially if you have an idea why, definitely point it out! This model can input any radiosonde observation, so if you guys see anything that you want a closer look at, give me a shout. Now that the model is mostly complete, I am looking for good test cases for the verification process.

 

[MClarkson]

Heres the simulated radar reflectivity and rain fields from the same run. Units for rain is in g/kg

 

[mkavulich]

I know you're probably doing this for fun (and I admit, I'd probably enjoy playing with a toy model like this for a while), but you're not really going to get anything close to reality with a 2-dimensional model. Any plans to extend it to 3-d?

As far as things that look wrong, it looks fairly good for a simple 2-d convection model, but your boundary conditions might need some tweaking: that convection on the right boundary seems to have been caused by the outflow from the first storm. Also, the reflectivity (rain) that isn't obviously being held up by an updraft takes way too long to fall out. Might want to adjust that parameterization.

 

[MClarkson]

3d would require either a massive reduction in grid spacing or me learning parallel processing, so short answer... no. Not going to be 3d. I think it should do a decent job for simulating linear features like squall lines. Trying now to get some initialization soundings to cooperate with my boundary conditions to see exactly how well it does for squall lines.

I agree with the right boundary convection. The forced vertical velocity from the outflow off the first cell seems to be magnifying at the right boundary when that boundary should be letting everything through. Boundary conditions are really annoying me lately. My rain has a terminal velocity of 7 m/s at the surface, increased for decreased density aloft by 1/sqrt(density/sfc density). I am pretty sure I tested that. I think thats right. Let me go check to see if that is turned on. From initialization to first precip at the ground takes 37 minutes here, and for reference one of the pioneering thunderstorm models (Klemp 78) took 32 minutes.

 

[MClarkson]

Well now that I think about maybe not so massive of a reduction. If I went from 50m to 200m thats 4^3 times an extra 2 for my CFL number being better... thats 128 times faster, and if I add 100 grid squares in a 3rd dimension, the overall runtime should be similar. Hmmm... maybe I will get around to that.

 

[David Wolfson]

I'm not sure why you need a .05 second dt for stability? So long as your time-step is short enough that your fastest-moving parcel would traverse significantly less than half the dimension of a grid-cell, I think you're basically ok. With a 50m grid, 1/2 second dt should be plenty. Perhaps that's what you meant?

That also points to another possible way of saving CPU cycles -- nesting your grid temporally. You might have a "t-scale" field in the grid-cell record that is set to the modulus with which you process the cell. For the many cells where nothing is happening very fast in the area you set it to, say 1 <= t-scale <= 16. Your code, then, would increment a cycling ring flag in the range of 1 to 16 at the start of a time-step and only process a cell if the flag MOD t-scale is 0. Otherwise you would skip all the calculations for the cell.

 

[Jeff Duda]

This looks like pretty good work and all, so if you were just building this model to say you could do it, congratulations. That's quite a feat. However, the WRF package has several cases designed explicitly for idealized 2D and 3D convective modeling. If you really want to investigate convective processes at such high resolutions, I would recommend you just download and compile em_squall2d_x, em_squall2d_y, or em_quarter_ss.

 

[MClarkson]

"I'm not sure why you need a .05 second dt for stability?"

Excellent question. I wish I knew. That would save a lot of compute time if my CFL criteria was 1 instead of 0.1 But I can't keep it stable at these small grid spacings. If I recall Robert Edmonds, another poster here, had similarly low CFL on a similar scale model he made.

"nesting your grid temporally."

ya, I might try something like that. Sort of an adaptive time stepping applied to only part of the domain? I gave up on the idea of domain-wide adaptive dt because of the stability problems...

Ive run em_quarter_ss. I suppose I should look at the 2d squall lines. This was a semester project, built from scratch in part because I wanted to see if I could do it, and in part because I am tired of using code that I didnt develop that can be hard to keep track of. At least now, when my code breaks, I usually know what its doing. Usually...

 

[David Wolfson]

I think what you may be referring to is this thread: http://www.stormtrack.org/forum/sho...rrent-Modeling&p=251135&viewfull=1#post251135

Robert is far more expert than I, and I've been out of the modeling biz for some years now. But... considering the eight grid cells neighboring a 2-d target cell, how do you diffuse parameter values from the neighbors into the target during a time-step? I think if you're just taking some sort of average, that may be trouble. Most of the values associated with a parcel diffuse according to the vector velocity of the parameter. So downwind cells should contribute nothing, whereas upwind cells should contribute proportionate to how much of the cell parcel "penetrates" the target during the time step. Unless I'm way off base this filtered approach both reflects physical reality and should damp many noise issues that are forcing a very short dt.

 

[MClarkson]

Heres a kindof cool simulation... Thunderstorms forming in response to surface heating over an island. The three images, 2d cross sections each, are simulated radar reflectivity, clouds, and temperature anomaly. The sounding is a bit drier than typical tropics, so the cloud bases are a bit higher.

Radar:

Clouds:

Delta-T

 

[Jeff Duda]

Is the land-surface temperature the same as the sea surface temperature? I didn't see a sea-breeze front in those plots.

 

[MClarkson]

Heating was applied only to the land and not to the sea. There is a sea breeze front that forms in the early-mid portion of the simulation. It then breaks down under the outflow of the first thunderstorm, and perhaps reforms later and is blown off to the east side of the "island" This is looking mainly at the temp anomaly field.

Also here is a u wind image that shows the sea breeze fronts about to converge in the center of the island. Winds are a little bit higher than is realistic, but keep in mind the actual surface values are not shown. The first layer displayed is 50m winds. 10m winds would be about 70% of this.

U wind (m/s)