% Stommel's 1961 model of convection in coupled boxes. % % Compute the temperature and salinity difference, here represented % by y and x, between two well-stirred boxes that are both % in contact with a reservoir that has (y, x) = (1, 1). Both % boxes conduct y and x at rates 1 and delta (delta < 1), and the % density difference between the boxes is given by d = -y + R*x, where % we are assuming R > 1, typically. There can be advection (or % flushing) between the boxes at a rate d*lambda, where the % flushing does not depend upon the sign of the density anomaly. % % This code also has 1) flushing with inertia (not especially % interesting), 2) a flickering or random temperature % perturbation (slightly interesting), and 3) an oscillating % reservoir temperature (hoping for but not yet finding chaos; % may be present in parameter regimes not checked). % % This code has not been fully tested, but seems to reproduce % S61's two cases fairly well when all three of the 'new' things % are turned off, of course. % % Written by Jim Price, April 28, 1999. Public domain. % set(0,'DefaultLineLineWidth',1.2) set(0,'DefaultTextFontSize',14) set(0,'DefaultAxesLineWidth',1.6) set(0,'DefaultAxesFontSize',14) clear; format compact; nn = 0; % set the model parameters (follows S61, aside from 'new') R = 2.0; % abs of the ratio of the expansion coefficients, x/y delta = 1/6; % conduction rate of salinity wrt temperature lambda = 0.2; % inverse non-d flushing rate; = inf for no flushing q = 0.; % initial flushing rate (0 to 1) 'new' qdelta = 100.; % time constant (inertia) for flushing; 'new' % set = 1/dtau for equilibrium flushing as in S61 % set = 0.2 for slowly responding flushing yres = 1.; % steady reservoir y, = 1 for S61 case 'new' resosc = 0.; % amplitude of reservoir y oscillation 'new' dtau = 0.01; % the time step of non-d time; 0.01 seems OK nstep = 1500; % number of time steps; 1500 is usually % enough to insure convergence to a steady state % % This model version is set up to do integration over a range % of initial T,S or y,x from 0 to 1. The increment of T and S % are set by delT and delS = 1/ni, where ni is the number of % integrations (n1 = 1 to 20 is reasonable). yres0 = yres; ni = 6; delT = 1/ni; % make ni = 1 or 2 to reduce the number of % integrations to be done delS = delT; for n1=0:delT:1 % n1 and n2 are used to set the initial T,S for n2=0:delS:1 if n1==0 | n1==1 | n2==0 | n2==1 % skip all non-boundary initial points x(1) = n1; % set the initial temperature y(1) = n2; % set the initial salinity for m = 2:nstep tau(m) = m*dtau; % the non-d time % evaluate the reservoir temperature (y); note that % this temperature is steady if resosc = 0. (the S61 case) yres = yres0 + resosc*sin(tau(m)*pi); % the first part of a second order R-K method; time step forward % by half a time step to make a first guess at the new times dr = abs(R*x(m-1) - y(m-1)); % the density anomaly qequil = dr/lambda; % the equilibrium flushing yh = y(m-1) + dtau*(yres - y(m-1))/2 - ... dtau*y(m-1)*q/2; % time step the temperature xh = x(m-1) + dtau*delta*(1 - x(m-1))/2 - ... dtau*x(m-1)*q/2; % time step the salinity qh = q + dtau*qdelta*(qequil - q)/2; % time step the flushing % the second part; use the half time step values to make a full step dr = abs(R*xh - yh); qequil = dr/lambda; y(m) = y(m-1) + dtau*(yres - yh) - ... dtau*qh*yh; x(m) = x(m-1) + dtau*delta*(1 - xh) - ... dtau*qh*xh; q = q + dtau*qdelta*(qequil - qh); % now add on a flickering temperature if you want to (or comment out) % tflickamp = 0.01; % set the amplitude here % tflick = tflickamp*unifrnd(-1., 1.); % y(m) = y(m) + tflick; end % end of time step loop d = R*x - y; % evaluate the density if nn == 0; nn = 1; % make a time series plot of the first case only figure(1) clf reset subplot(2,1,1) plot(tau, x,'-', tau, y,'--') legend('salinity', 'temperature') ylabel('T, S diff, non-d') title('Experiment 1,1') subplot(2,1,2) plot(tau, d) xlabel('time, non-d') ylabel('density diff') % contour the density (or flushing rate), and add the (T,S) % trajectories on top of the contours ym = 0:0.1:1; xm = 0:0.1:1; for k1 = 1:11 for k2 = 1:11 dm(k1,k2) = (1/lambda)*(R*xm(k2) - ym(k1)); end end figure(2) clf reset dc = -10:2:20; c = contour(xm, ym, dm, dc,'k'); clabel(c); xlabel('salinity diff, non-d') ylabel('temp diff, non-d') hold on end % the if on nn = 0 [m1 m2] = size(x); % plot the individual trajectories. % color code according to which equilibrium point % the trajectory ends up on. this will likely have be % reset if the model parameters (R, delta, lambda) are changed. if d(m2) >= 0 plot(x, y, 'r--') plot(x(m2), y(m2), '*r') else plot(x, y, 'g') plot(x(m2), y(m2), '*g') end clear x y d end % or to gate out non-boundary values of n1, n2 end % loop on n1 end % loop on n2 % make some plots to show where roots (equil. points) are for k=1:60 f(k) = (k-30)*0.1; lhs(k) = lambda*f(k); rhs(k) = (R/(1 + abs(f(k))/delta)) - 1/(1 + abs(f(k))); end figure(3) clf reset plot(f, rhs, f, lhs) xlabel('f, flow rate') ylabel('lhs(f), rhs(f)') title('roots of S61 model') grid % end of the script