Numerical Simulation of the Vortex on Jupiter
PURPOSE
The purpose of our study is to model the Great Red Spot (GRS) on Jupiter using the principles of fluid dynamics and a super-computer. Once the initial model was completed, we changed the initial size of the vortex by altering the Gaussian Hill (Valley) Function. This enabled us to compare the different outputs of the simulation for structure and stability.
SIGNIFICANCE
Fluid dynamics, the study of motion and equilibrium of fluids, is instrumental in the study of earth’s dynamical systems. Our problem uses the principles of fluid dynamics and is therefore a significant problem and of value to Earth Science. The GRS is a storm system that is very similar to storm systems on earth, such as hurricanes and tornadoes. By simulating the GRS on Jupiter, it is possible to simulate a similar storm on Earth using our code.
BACKGROUND INFORMATION
Jupiter, the fifth planet from the sun, is the site of the GRS. For more than three hundred years the storm has raged. Jupiter is the largest planet, with a diameter of 142,984 kilometers, and has the fastest rotation period of any planet in our solar system, making one rotation every nine hours and fifty-five minutes. The gravitational pull of Jupiter is 0.0249 kilometers/second2, or approximately 2.5 times greater than that of earth. Jupiter has a mean orbital radius of 778,400,000 kilometers from the sun, and completes a revolution in twelve earth years. The atmosphere of Jupiter is made up of hydrogen (82%), helium (15%), and the remaining three- percent methane, ammonia, carbon monoxide, ethane, acetylene, phosphate, and water vapor. The temperature at the top of Jupiter’s clouds is about –157o C, but increases greatly with depth.
The GRS was a mystery for many years, and was recently (within 100 years) discovered to be a storm system. It is found about eight kilometers above the clouds of Jupiter, at a latitude of 22.5o south. It is a large, shallow feature that measures about 40,200 kilometers long and 32,000 kilometers wide. The vortex varies in depth from 20-40 kilometers. The periphery of the GRS rotates counterclockwise, once every six days, which indicates that it is a high-pressure zone. The GRS is between opposing zonal winds, which reach speeds up to five hundred kilometers per hour. The red color of the spot is due mostly to phosphorus, which is produced when phosphine in the lower atmosphere is broken up by sunlight.
MATERIALS AND PROCEDURE
During the course of our project, we used a C++ compiler to create the program. We used both personal and super-computers to test and run the program. Due to unforeseen circumstances, Dr. Jones ran our output files in Ferret, a graphics program to plot vector simulations. Dr. Jones then placed .gif files on our Pi account. We also used a text editor to write the final report.
To begin our project, we met with Dr. Philip Jones of Los Alamos National Labs to learn the basics of fluid dynamics. By looking at scientific papers written about the GRS, we determined the phenomena that would be factors in the evolution of the GRS in our simulation. We calculated the acceleration due to gravity on Jupiter, as well as the Coriolis force due to Jupiter’s rotation. We referred to scientific papers to obtain the values of the densities of the gases on Jupiter, the Rossby deformation radius, and the zonal velocities of the underlying layers.
DERIVATION OF EQUATIONS
For the derivation of the equations, we had copious amounts of assistance from Dr. Philip Jones. He came to our high school to explain all the basics of fluid dynamics so that we would have the essentials for our project.
To derive the equations that were used to model the vortex, we studied several different aspects of fluid dynamics including: conservation of mass, conservation of momentum, density of fluids, conservation of vorticity, and the coriolis force. Each of these important principles was needed to define the aspects of our problem.
Mass Conservation: Consider a cube in which fluid can flow freely through each side. We will use x to describe the east-west direction and u as the velocity in that direction. Similarly y, v will be the north- south direction and z, w will be the up- down direction. The mass inside the box must be conserved. If there are no other forces acting on the box, the change in mass will only be due to the amount of fluid entering and leaving the box. The GRS is like this box, having particles entering and leaving, but the total mass remains constant.
Lagrangian Form: Using the example of the box again, we can say that the mass of the box will still remain the same, even if the density changes due to change in volume.
Conservation of Momentum: The conservation of momentum and the conservation of mass follow a similar pattern. Referring to the box again, the momentum can change due to the velocity of particles flowing from other sides of the box. Other forces acting on the box, gravity and pressure forces can also affect the momentum. In our simulation, we assume that no forces act on the box to change the momentum, thereby simplifying the problem.
Coriolis Force: This is, the rotation of a fluid due to the rotation of a planet. Particles rotate faster near the equator than at the poles. The Coriolis parameter is defined as f=2Wsinq where W is the rotation of the planet (in radians/ second) and q is the latitude.
Geostrophic/ Quasigeostrophic Balance: Geostrophic balance means that fluid flows along lines of constant pressure. We use a geostrophic balance in our initial conditions to start the vortex. Since geostrophic conditions do not usually exist in nature, the simulation quickly becomes quasigeostrophic. The quasigeostrophic approximation assumes that the major forces affecting momentum are the Coriolis force and the pressures due to density differences between layers.
Vorticity: Vorticity is the measure of the rotation of a fluid. Vorticity only changes due to a change in latitude or depth of the layer. Vorticity changes with latitude because the Coriolis force changes with latitude, and consequently the speed at which the vortex will rotate. A change in depth will stretch the vortex, which will change the angular velocity of the particles in the vortex.
In our simulation of the GRS, we are assuming a closed system, which allows the principle of conservation of mass to be used in simplifying our problem. The opposing winds, flowing east –west, south of the GRS, and west –east, north of the GRS, cause the vortex to rotate counterclockwise. Within these winds, particles are entering and exiting the system. In our simulation, those particles that leave on the west, reappear in the east. If an imaginary band was placed around Jupiter at the latitude of the GRS, this same phenomenon would occur. We use boundary conditions in the north and south directions, setting the velocity to zero, to assure that no matter leaves the system in these directions. To conserve mass, we have allowed for our simulation of the GRS to have a fluctuating thickness. The winds that cause the GRS to rotate also add to a fundamental principle of physics, conservation of momentum. The GRS will continue to rotate due to angular momentum, making it more stable than other systems. The conservation of vorticity, or angular velocity of a fluid at a point, is related to the conservation of momentum.
The GRS is a very shallow layer of gases sitting on top of a more dense, deeper layer of gas, therefore, the density of fluids is important to our simulation. The thickness and pressure gradient of the GRS and the layer below are important to our simulation. Using shallow water equations, which say that the GRS is a relatively shallow feature compared to its other dimensions, we were able to put a term into our equation that involved the different densities of the separate layers of gases. This term included the Rossby deformation radius, a measure of the vertical stratification of the gases. The density of the lower layer of gas, below the GRS, is not significantly larger than that of the GRS. Therefore the thickness of the GRS does not fluctuate much. Hydrostatic balance dictates that the less dense gases will be in a higher layer (or strata) than the more dense gases.
Because a plane rotates like a rigid body, the velocity of a particle due to the planet’s rotation varies with its position on the planet. At the equator the velocity is faster because the surface of the planet is farther from the axis of rotation. The coriolis force deflects the motion of a particle to the right in the Northern Hemisphere and to the left in the Southern Hemisphere.
Using these principles, we were able to use common fluid dynamics equations to model the GRS. The fluid dynamics equations are:
Boundary Conditions are:
East-West: periodic (cyclic)
North-South:
Initial Conditions are:
The equations for fluid flow can be difficult to use and are frequently chaotic, requiring a numerical solution. With the help of Dr. Philip Jones, we were able to transform the differential equations into numerical equations. The equations that we derived are on the next page. The boundary conditions for our numerical equations are:
The variables in each equation represent the different principles and factors, which we have mentioned previously. F0 is the coriolis parameter, at the reference latitude of 22.5° south. b results from the beta plane approximation, which approximates the Coriolis parameter by a linear function. Its units are radians/ second/ kilometer, and it represents the change in rotation per change in latitude. The east/ west velocities are represented by u, and the north/south velocities are represented by v. H is the thickness of the GRS, found to be about 30 kilometers, which fluctuates slightly (Marcus, p.524). The gravitational pull of Jupiter (.0249 km/s2) is represented by g. r is the density in the top layer of the GRS, and r2 is the density of the layer below the GRS. r2 is only slightly larger than r, having a ratio close to one. D is a diffusion coefficient, which was added to keep the values generated from the grid from spreading.
Equation one simulates the east/west velocity and equation two the north/south velocity. The third equation represents the thickness of the GRS. In the first two equations, terms one and two are advection terms. This means that particles are carried with the fluid flow, from one place to another, advecting momentum. Advecting momentum dictates that the particles will keep moving in a circular pattern. The third term, which contains the f0 and by, is the coriolis term, which results from the rotation rate of the GRS. The fourth term takes into account the density and thickness of the GRS. The last term in the brackets includes D, the diffusion term. We added this term to keep the solution stable and prevent checkerboard instability. Checkerboard instability results from a combination of leapfrog and centered differences, which allows the solution to separate on even and odd points.
INITIAL CONDITIONS
Our initial thickness is based on an average of 30 kilometers, although there is some variation. This variation is caused by a Gaussian hill, which changes the depths to make the initial shape of the vortex. The diameter of this hill was varied in our program to change the initial diameter of the vortex.
Using our boundary condition that sets our simulation in geostrophic balance from terms in the partial differential equation for the v velocity, we solved the equation for u, the east-west velocity. Solving this for u gives us our initial geostrophic conditions for the vortex.
We also needed the initial conditions for v to be geostrophic, so from terms in the partial differential equation for the u velocity, we solved for v.
BOUNDARY CONDITIONS
We had to include boundary conditions inside our code to conserve mass within the vortex. We made all the north/south velocity at j =1 and j=y_grid (southern most point in the model) equal to 0. Also we made all the north/south velocities at j=0 equal to the velocities at j=2 and the velocities at j=y_grid+1 equal to the velocities at j=y_grid-1. This makes a boundary at the northern and southern ends and doesn’t allow mass to leave the system. On the east/west velocities we included a boundary that any mass leaving one side would enter in on the other side. In the code when i was equal to either of the endpoints on the grid, i plus 1 was equal to the other endpoint. Also in the east/west velocities we made all velocities at j=0 equal to the velocities at j=2, and the velocities at j=y_grid+1 equal to velocities at j=y_grid-1. Our boundary condition at the north and south ends of the vortex allows us to maintain geostrophic balance in our simulation. Fixing this thickness would be too restrictive, so we allow it to change as the east-west velocity changes. This flexibility allows it to change while still allowing for conservation of momentum.
PROGRAMMING OF EQUATIONS
To begin our simulation, we have broken up the area around the GRS into little boxes, or grid cells. This allows us to assume the GRS is a two- dimensional feature (the GRS is nearly two- dimensional because it is relatively shallow as compared to its length and width). The grids were labeled with points in the form (i,j) where i is the position in the x- dimension and j is the position in the y- dimension. Each point (i,j) is dependent upon the points immediately adjacent to it.
The fluid flow must be advanced in time. To accomplish this, a time-stepping method was used. We define the superscript n as the current time period, n+1 as the subsequent time period, and n-1 as the previous time period. The values at each point (i,j) depend on the values at the previous time period from the same point (i,j).
We wrote a C++ program with a 100 by 225 grid that integrated for 200,000 time steps, which is about 3.6 years, real time. Each time step is approximately 571.5 seconds, or 9 minutes, 31 seconds. The limit of 200,000 time steps was chosen because the original version of the code became unstable shortly thereafter. We also felt this was a long enough time period to give a good representation of the stability of the vortex. For each variable in the code, we have three different, two-dimensional arrays representing the time levels (old, current, and new). In order to keep the program advancing in time, each array was moved down one time level (new becomes current, current becomes old). Data was outputted to files at 1, 5000, and every 20,000 time-steps from 20,000. Using our initial conditions we initialized our arrays. To give our equations stability, we calculated a time step that was small enough to eliminate noise.
This time stepping was accomplished by using a leapfrog method, in which our change in time was equal to the smaller of the two equations:
Dx/u and Dy/v
We multiplied the result by .01 to keep any instabilities from forming in the initial startup of the simulation. During the first time step, we had to divide this change in time by two, because we were not performing a full leapfrog. This is because we did not yet have the values at the n+1 time level for u, v, and H.
We used the Borland C++ Version 4.5 compiler to write most of the code. As the memory capacity needed for our program grew, we had to begin using the super-computer on Mode and Pi to compile our code. Our code is around 500 lines long. It took approximately three hours to write the first draft of the code. However, it took months to gain a working simulation. Some problems that we ran into included having a time step that was too large to be stable. This made the output numbers progress to infinity very quickly. We also discovered that initial conditions and boundary conditions had to be particular for the GRS to begin the spiral. In the beginning, we did not have the diffusion term or the density term, which ended up being very important to the final coding. Once everything was finally in place, it took about three hours to run our code on Pi. We did try running it on a personal computer, but after about 21 hours of processing, we decided to leave it up to Pi.
REFERENCES
Beebe, Reta. "Characteristic Zonal Winds and Long-lived Vortices in the Atmospheres of the
Outer Planets." Chaos Volume 4 Issue 2 pages 113-122 (1994)
Brown, Robert A., Fluid Mechanics of the Atmosphere. (Academic Press, Inc., San Diego,
1991). pages 1-11, 294-297
Gutfinger, C and Pnueli, D. Fluid Mechanics. (Cambridge University Press, Cambridge, 1992).
Pages 1-14, 162-165
Ingersoll, Andrew P. "Models of Jovian Vortices." Nature. Volume 331. February 25, 1988.
Pages 654-655
Marcus, Phillip S. "Jupiter’s Great Red Spot and Other Vortices." Annual Review of Astronomy
and Astrophysics. Volume 31, pages 523-73 (1993)
Marcus, Phillip S. "Numerical Simulation of Jupiter’s Great Red Spot." Nature. Volume 331.
February 25, 1988. Pages 691-696