12.747 Lectures 19&20: Section 2c:
Upper Ocean 1-D Seasonal Models
File last modified 9 May 2002
19.2c The Physical Model: Doing It
Well, enough fooling around, let's do it! First, let's look at the code. It consists of the main m-file, pwp.m plus a bunch of ancillary m-files which do a number of tasks. They are
To run the model, you need to download all these little files (remember to shift-click) and save them in your current directory. Then run MATLAB. Running this model in MATLAB on a 100 MHz pentium takes about 15 minutes per model year, or about 45 minutes for a complete, three year run. This is unfortunately slow, but because the code is not easily vectorized. You can achieve about a factor of 5-10 speedup by converting the MATLAB code into FORTRAN or C. That's up to you. This is what pwp.m looks like:
inifctr; % Initialize useful factors inihydro; % initialize water column profile % % Main time step loop % for it=1:nt T(1)=T(1)+thf(it); % add sensible + latent heat flux S(1)=S(1)-FWFlux(it); % flux fresh water (latent/precip) T=T+rhf(it)*dRdz; % add radiant heat to profile dostins; % do static instability adjustment addmom; % add wind stress induced momentum dobrino; % do bulk Ri No Adjustment dogrino; % do gradient Ri No Adjustment advdif; % advect and diffuse dooutput; % if time, save data end save run1 Ta Sa zmld nnstin nnbri nngri % save run info.
When the code runs, there is a significant period of time while it is "initializing" (in inifctr), where it is building the forcing vectors, etc. but then it spends most of its time in the time iterations. Once all the forcing vectors are built, and the initial vertical profile is specified (in inihydro), then the main time loop takes place. The surface-most cell is first cooled by the heat flux (sensible + latent + outgoing longwave radiation), then fresh water flux applied. Next, radiant heating is to the whole water column. Once this forcing is applied, then the water column is adjusted for static stability (dostins). The reason for doing that here is that the cooling process tends to produce density instabilities. Next, wind stress is applied to the mixed layer (addmom). After that, the bulk Richardson Number stability is applied, and then the gradient Richardson number mixing. Finally, vertical advection-diffusion is applied.
The routine "dooutput.m" contains code which prints to the screen, but it is up to you so save the model output from the run in a form which is useable for analysis of the run in the end. The single line at the end of the program simply saves some of the data for that purpose. What is saved is the accumulated temperature and salinity profiles (Ta and Sa) the mixed layer depth history (zmld) and counters indicating water column adjustment activity (see previous section). You can build in other things that you think might be informative or interesting, and save that as well. For example, the horizontal velocities are kept in a two column vector "UV" which you might save to see what the velocity shear is doing during the run.
Well, let's look at the model behavior. Below is a plot for a three year run with Kz = 1 x 10-4 m2/s.
You can see that the model seems to get into some kind of stable "loop" like a seasonal cycle, with the formation of a summer warm layer, and a winter deep mixing event. However, there is that kind of long term drift that we mentioned: note how the seasonal penetration of the 20 degree isotherm (at about 120 m depth) seems to shoal during the run. Realistically, we cannot expect the model to balance over long times, because we have ignored non-vertical, non-local processes such as horizontal advection of heat. Thus we have to cut our losses and limit our simulations to a few years. The curves above look qualitatively like what is observed:
Comparison of the model results with observations is a critical phase in the process of modeling. The tricky part is finding the appropriate data to make the comparison, and then to choose the right aspects of the model output to compare with the data. What are the critical modes of model behavior? What is important to match with "reality" and what can we safely ignore? The answers are complex, and rooted in what your goals are. The only tunable parameter in the model, as formulated is the vertical mixing rate "Kz", although certainly some of our choices of climatological forcing need to be looked at more carefully. For the sake of instruction, however, we'll stick to the one parameter. Now let's look at some comparisons.
We compare the model temperatures at three fixed depths (the surface, 50 meters and 100 meters depths) as solid lines vs data taken from time series stations near Bermuda between 1980 and 1990 (the red dots)..
Three values of vertical diffusivity are shown (the color coded lines). First, the surface temperature seems to be extremely sensitive to our choice of Kz: the low Kz values tend to trap heat in the upper layer, producing a very hot, thin mixed layer. Thus the first order observation is that the optimal value of Kz for the model is somewhere around 1.5 x 10-4 m2/s. Getting the model surface temperatures right, aside from being an indicator of subtler aspects of model performance, is obviously critical if we are to do the gas budgets right, because the gas fluxes are driven by solubility effects, which are in turn dependent on temperature. The thermal history seems to lock us into a relatively narrow range of acceptable Kz.
This seems to be the case for the 50 m depth as well, although data scatter makes it difficult to make a definitive comparison. What seems to happen, however, is that the low Kz cases generate a thermal overshoot at that depth in the autumn. The 50 m comparison is also critical, because the trapping and transport of heat is a useful diagnostic of the trapping and transport of gases (as we will see in the next section) so we definitely need to get that right. Finally, all the 100 m model curves appear warm compared to the data, although the higher Kz case seems to be better than the others. This offset should raise little alarm bells in your head, but may also be due to subtle problems with our choice of initial hydrography and bottom boundary conditions. What may be more important at this level is the amplitude of variation.
The data scatter is a little problematic in comparing model results. We can also compare to hydrographic climatologies. Below is a comparison to the Levitus hydrographic atlas data:
which shows more-or-less the same thing that we discussed before.
Another way is to look at the mean difference between the climatological data and the model results over the entire model domain. This can be a little misleading, since a residual offset may be introduced by imperfect choice of initial conditions (if the initial heat content of the model is not quite right, and the model is thermally closed, i.e. the heat fluxes average out to zero over the annual cycle, the initial heat content will always be wrong). However a better way to compare is to look at the mean square difference. We show two view of the mean square difference as a function of model diffusivity (the X-axis) and the stochastic forcing amplitude (Y-axis). The latter is simply a multiplier applied to random winds added to the climatology. A value of 0 means that the climatological mean winds only are used (no random wind forcing is added). A value of 1 means that the average amplitude of the random winds is equal to the Trenberth climatological wind stress standard deviation.
There are two aspects of this comparison which are interesting. The first is that the comparison is not perfect: the minimum mean squared difference is about 0.15 degrees-C squared, corresponding to an RMS offset of about 0.2 degrees C. This is symptomatic of the effects mentioned above. Second, there appears to be a compensatory trade-off between vertical diffusion and random wind forcing. The optimum fit for purely smooth, climatological mean forcing is about 10-4 m2/s, while for the more realistic case with climatological stochastic forcing, it is more like 0.9 x 10-4 m2/s. Thus the diffusivity is masquerading as random wind events, which serve to episodically deepen and shoal the mixed layer. This is characteristic of this kind of problem.
Now the final concern is how the mixed layer is simulated. This, of course, will affect the ability of the model to transfer gases to and from the atmosphere. One must be careful, however, since one's choice of criterion for what exactly you call the mixed layer depth may influence your "success". What we mean is that actual ocean observations do not usually show a clean point at which the mixed layer ends, and the stratified ocean below begins. This is result of a complex meteorological history, and just plain orneriness. Here we've used a compromise criterion that the the mixed layer depth is half way between the depth at which the temperature is more than 0.15 degrees lower than the surface, and the depth above. So, for example, if the first depth at which the temperature drops below the critical value is 150m, and if the depth sampled immediately above that is 120m, then the mixed layer depth will be 135m. Here the model doesn't do too badly, but it still has its problems. Below is a comparison with the Levitus climatology:
The comparison again tells us that the very high or very low diffusivity cases do not do a good job in emulating the system, since the winter mixed layers are not deep enough. The Kz around 1.0 x 10-4 m2/s seems to do best. There is, however, a problem with the model forming a shallow spring mixed layer too soon (compared to the climatology), and not eroding fast enough in the autumn. The climatology, actually, may be wrong in the spring, since it may "smear out" in time an annually variable onset of spring warming. Thus we aren't too concerned with the spring disparity. However, the autumnal deepening curve does not match well, although the model eventually "catches up".
In summary, it looks like we have a physical model which does a half-way decent job of pretending it's an ocean. Now we'll go ahead and try to put gases into it.
GoTo Next Section