% vertical_pred_prey is a Matlab script for calculating % the interactions among particle flux, flux feeder, and prey % % % Time-stamp: <2004-05-17 16:53:03 adrian> % % % USAGE: % Type vertical_pred_prey at the Matlab prompt. % % HISTORY: % **-**-01: Original code by GAJ. % 25-04-03: Code re-written & re-organized by ABB % % NOTES: % For theory, see % Jackson, G.A. & A.B. Burd (2002) Deep-Sea Res. II 49:193-218 % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % George A. Jackson % Department of Oceanography, Texas A&M University. % % Adrian B. Burd % Department of Marine Sciences, University of Georgia % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Clear workspace close all clear all %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % MODEL PARAMETERS % % Define dimensional (and some dimensionless) parameters and from % them, calculate required dimensionless parameter values. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Define dimensional parameters a = 0.3; % efficiency of P eating particles b = 123; % clearance rate for P [m^3/gC_p/d] c = 0.3; % efficiency of Q eating P d = 0.01; % Predator specific mortality [1/d] F_surf = 0.1; % Surface particle flux [gC/m^2/d] sigma_p = 3.95; % Interaction cross section [m^2/gC_p] lambda = 6.25e-04; % Bacterial degradation rate [1/d] v_f = 30; % Particle settling vely [m/d] omega_f = 2*pi/365; % Flux forcing frequency [1/d] % Define some non-dimensional parameters gamma = 0.25; % relative amplitude of seasonal flux oscillations % Calculate some dimensional constants beta_0 = sqrt(a*sigma_p*d*F_surf); % Lotka Voleterra frequency p_ss = d/(b*c); % Steady-state P [gC/m^3] % Calculate non-dimensional parameters psi = sigma_p*v_f/(b*c); nu = b*c*lambda/(d*v_f*sigma_p); beta_0 = beta_0/d; omega_f = omega_f/d; beta_0 = omega_f/0.5; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % MODEL GRID % % Set up a grid from 100 -> 1000 m depth with 100 layers and 0 -> 4 % years in time with the time step determined by the CFL % condition. Note convention that dimensional depths and time have % "depth" and "t" in their names, dimensionless version have "z" % and "tau". % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% t_start = 0; % Starting day [d] t_final = 12*365; % End day [d] top_depth = 100; % Depth of top layer [m] bottom_depth = 1000; % Depth of bottom layer [m] n_layers = 100; % Number of layers tau_start = t_start*d; % Dimensionless start time tau_final = t_final*d; % Dimensionless final time z_top = top_depth/(v_f/d); % Dimensionless top layer z_bot = bottom_depth/(v_f/d); % Dimensionless bottom layer z = linspace(z_top, z_bot, n_layers); % Dimensionless depths z = z'; delta_z = (z_bot - z_top)/(n_layers-1); % Dimensionless delta z delta_tau = delta_z; n_tau = floor((tau_final - tau_start)/delta_tau); tau = (0:(n_tau-1))*delta_tau; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % THE INTEGRATION % % Following G, use a method of lines by discretizing space (use % backward difference for flux and central difference for % diffusion) and an Euler scheme for the time derivative. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Some functions we need - these are the time variations of the % particle flux at the surface and the time variations of the % vertical migration: NOTE, the definition of ekz is different from % that in G's code, my version is 1/his. f0 = 1 + gamma*sin(omega_f*tau); % Flux forcing at surface ekz = exp(-psi*(1 + nu)*z); % Dimensionless exp(-k*z) beta_sq = beta_0^2.*ekz; % Define a matrix of coefficients for the central differences used % in the second order derivatives. ckz = 1e-4; Akz = -2*diag(ones(n_layers,1)) + diag(ones((n_layers-1),1),1) + ... diag(ones(n_layers-1,1),-1); Akz = sparse((delta_z^(-2)*ckz)*Akz); % The dimensionless initial values - i.e. the initial depth % distributions. Here we choose the steady state sans migration % case for f and q and a Gaussian for p. This translates into a % constant for the dimensionless f f_init = ones(n_layers,1); q_init = f_init.*ekz; p_init = f_init; % Set up indexing vectors for the finite difference z_range_1 = 1:n_layers-1; z_range_2 = 2:n_layers; % Set up storage for the solutions and put the initial conditions % in the storage ft = zeros(n_layers, n_tau); pt = ft; qt = ft; ft(:,1) = f_init; pt(:,1) = p_init; qt(:,1) = q_init; % Now loop through each time step for it = 2 : n_tau dfdt = zeros(n_layers,1); dpdt = dfdt; dqdt = dfdt; % The equation for the derivative of the flux dfdt = psi*ft(:,it-1).*(1 - pt(:,it-1))*delta_tau; % The prey equation. This includes a diffusive term. Note the % difference in the formulation of the (f-q) term when we don't % have the diffusion in. dpdt = beta_sq.*pt(:,it-1).*(ft(:,it-1) - qt(:,it-1)./ekz) ... + Akz*pt(:,it-1); % The predator equations dqdt = qt(:,it-1).*(pt(:,it-1) - 1) + Akz*qt(:,it-1); % Update the flux, prey and predator depth profiles ft(1,it) = f0(it); ft(z_range_2,it) = ft(z_range_1,it-1) + dfdt(z_range_1); pt(:,it) = pt(:,it-1) + delta_tau*dpdt; qt(:,it) = qt(:,it-1) + delta_tau*dqdt; % Impose no-flux conditions at the boundaries pt([1, n_layers],it) = pt([2, n_layers-1],it); qt([1, n_layers],it) = qt([2, n_layers-1],it); end % Plot some things figure(1) subplot(3,1,1) plot(tau/d/365, ft(1,:), 'b', tau/d/365, ft(n_layers,:), 'r--') ylabel('particle flux') title(['\psi = ' num2str(psi) ' \nu = ', num2str(nu), ' \gamma =' ... ' ' num2str(gamma)]) subplot(3,1,2) plot(tau/d/365, pt(1,:), 'b', tau/d/365, pt(n_layers,:), 'r--') ylabel('particle feeder p') subplot(3,1,3) plot(tau/d/365, qt(1,:), 'b', tau/d/365, qt(n_layers,:), 'r--') ylabel('predator q'); xlabel('time (yr)') figure(2) plot(pt(1,:), qt(1,:), 'b') plot(pt(n_layers,:), qt(n_layers,:), 'r') xlabel('particle feeder p') ylabel('predator q') title(['\psi = ' num2str(psi) ' \nu = ', num2str(nu), ' \gamma =' ... ' ' num2str(gamma)]) %------------ Last line of vertical_pred_prey.m