Stommel's Gulf Stream Model

This worksheet is a MAPLE translation of Reza Malek-Madani's Mathematica worksheet.

In this computer assignment you will be reproducing a fundamental result of ocean circulation, the
wind-driven circulation of Henry Stommel (1948). In Stommel's model we encounter a system of
nonlinear differential equations of the form

(*) x'=f(x, y), y'= g(x, y),

where x and y are the east-west and north-south coordinates of a typical fluid particle in the Atlantic

ocean. The functions f and g take into account the balance of forces due to the prevailing trade winds

in the Atlantic as well as the force due to Coriolis.

At this stage of the course we are not ready to derive the equations that I will cite shortly; We are
able, however, to solve the initial-value problem that results once f and g are given to us.

The Atlantic ocean of Stommel's model consists of a unit square. The functions f and g have two
parameters a and b in them that are affected by the strength of the prevailing winds and the Coriolis
parameter. We start our program by defining a constant A by

A = (exp(a)-1)/(1-exp(b))

in terms of which we define the stream function P by

P = sin(pi y) (exp(ax) + A exp(bx) - 1 -A).

Finally, the f and g of (*) are defined by

f = dP/dy , g = - dP/dx.

At this stage we have not studied what a stream function is so you are not expected to understand
the logic behind the previous three mathematical statements. Suffice it to say that by the end of this
semester the mathematical process we just went through to get at f and g will be second nature to
you.

Now that f and g are known explicitly, we are ready to plot the phase portrait of (*). The following
program does just that for the initial data

((0.5,0.3), (0.5, 0.2), (0.5, 0.1), (0.5, 0.01));

I have arbitrarily chosen a=-10, b=1, and t_final = 2. Here is the program.

sm311o_6.mws,wdj,1-20-98

> restart;

> with(plots):

> a:=-10; b:=1; t_final:=2;

> A:=(exp(a)-1)/(1-exp(b));

> P:= (u,v)->sin(Pi*v)*(exp(a*u)+A*exp(b*u) -1 -A);

> f:=(u,v)->subs({x=u,y=v},diff(P(x,y),y));
g:=(u,v)->-subs({x=u,y=v},diff(P(x,y),x));

> eqns:=
diff(x(t),t)=f(x(t),y(t)),diff(y(t),t)=g(x(t),y(t));
ics:=[[0.5,0.4], [0.5,0.3], [0.5, 0.2], [0.5, 0.1], [0.5, 0.01]];

> for i from 1 to nops(ics) do nsoln[i]:=dsolve({eqns,x(0)=ics[i,1],y(0)=ics[i,2]},{x(t),y(t)}, type=numeric);
od;

> for i from 1 to nops(ics) do
de_plot[i]:=odeplot(nsoln[i],[x(t),y(t)],0..2,numpoints=100):
od:

1. Here is the output of the program. Which way is the fluid circulating?

> display([seq(de_plot[i],i=1..nops(ics))]);

2. Change the parameter values of a and b to -1 and 1, respectively. Plot the phase portrait. How
does it differ from the one above?

3. Let's go back to a=-10 and b=1. Change t_final =1. What is the effect of this parameter on the
orbits of the fluid particles?


4. Let b=1 and t_final = 2. Let a vary from -1 to -10 at increments of -1. Plot the graph of the orbit
of (0.5, 0.1) for each a. Combine all orbits on the same screen. What is the effect of varying a on
the orbit of this particle?


5. Replace sin(pi y) with sin(2 pi y) in the definition of P. Modify the above program to allow for this

change. Execute it with parameter values t_final = 2, a=-10, and b=1. Report on how the output

differs from the original case. You may need to alter the content of data to get a representative picture

of what is going on in this case.

6. Replace sin(pi y) with sin(3 pi y) in the definition of P. Modify the program to allow for this

change. Execute it with parameter values t_final = 2, a=-10, and b=1. Report on how the output

differs from the previous two cases. You may need to alter the content of data to get a representative

picture of what is going on in this case.