A helpful reference for this assignment is given below:
Roache, P.J. (1976) Computational Fluid Dynamics, Chap. 3 is particularly useful.
1.) (a)Consider a well mixed water mass (read that as a "box model") which is ventilated at a rate of 0.1 y-1. That is, it has a "true" age of 10 years. Now run the model using the tritium source function natrsf.dat, from 1950 to 1996 in one year steps, setting it up as we have done in the lecture notes. Don't forget to start at the "steady state" values for the two isotopes. For simplicity, use a mean decay life for tritium of 18 years. This time, however, we want you to plot as a function of time, the tritium, the 3He and the tritium-helium age, defined by
where T is the tritium concentration, and the log is the natural logarithm (not base 10). Combine all of these in a single plot (use the subplot(3,1,1), subplot(3,1,2) and subplot(3,1,3) commands). Please explain in general terms what you are seeing.
(b) Now, do the same calculation again, but instead of using a constant renewal rate, use a climatically varying renewal time given by
where Tau is the renewal time in years and t is the year
(e.g.1981). Do the same plots as above, and comment on what
parameter has the greatest sensitivity to the climatic variations.
2.) Using the von Neumann stability analysis
technique, derive the appropriate stability criteria for the first
upwind differencing method (FUDM). Since this derivation would be
difficult to do in MATLAB or in a standard editor, please do this on
real paper (or e-mail us a PDF file, please, no
Word files) and hand it in on Thursday.
3.) Using the FUDM, perform an impulse experiment to demonstrate empirically what the numerical diffusivity is for this technique. Do this in the following way: set up a 1 dimensional system which is 25,000 km long. Choose a reasonable delta-X (say 50 km) and stable delta-T. Do not use a time increment greater than a month. Use the purely advective model (no explicit diffusivity). Start with zero concentration everywhere, then set the concentration of the grid point in the middle of the domain to be 1. Run the model with a velocity of .01 m/s for 10 years. As you run the model, generate annual plots of the following:
The last two plots should be cumulative (i.e. should retain the previous years' points, do this with "hold on" for them) and the first one should just be the most recent year (don't use "hold on" for it). In principle, the initial spike should evolve as a gaussian with time according to the relation
where K is the effective (numerical) diffusivity, t is the time, and x0 is the (moving) center of the peak. Now the area under the peak is given by (yes, this time we've got it right!)
now in a perfect world (model), this area, which is the total amount of stuff in the domain, should remain constant, and since we started with a control volume of size delta-X and a concentration of 1, then we have
or, squaring this, we have
Now if you plot the LHS vs. time, the slope will be what is in the parentheses on the RHS. So if you can estimate the slope, you have
Now do the above and calculate from the observed slope the apparent numeric diffusivity. Compare it to the "theoretical" numeric diffusivity. Repeat the calculation after first doubling the velocity and once again after doubling the space step. Interpret your results.