Given the equations derived last class, students worked on their MATLAB encodings in order to create a functional model for chemical release in the west passage.

The lateral flow in x₁ can be represented by:

U=(0.5*B/Mu).*Z.^2+((Uo/h)-(0.5*B*h/Mu)).*Z;

where Z is an array of depths.

A plot of the temperature distribution at any timestep can be plotted by:

figure
pcolor(squeeze(T(:,:,1)))
shading flat
colormap jet  %popular rainbow color gradient
xlabel('N-S Distance')
ylabel('Height above bottom, m')
text(140,32,'Time is 0 hours');
title('Temperature','FontSize', 18);
pause(.2)

where the third parameter in the T matrix (second line above) suggests which timestep to plot.

A chemical density matrix can be set up so that the full chemical release is at a specific x₁ pixel (icx) at depth (icz):

C(icz,icx,:)=Cmax;

The timestep looping code updates neighboring parcels for both temperature and chemical density in time via:

for k = 1:(ntimes-1)
    for i = 2:(nz-1)
        for j = 2:(ny-1)
            if (U(i)>0)
                T(i,j,k+1)=T(i,j,k)-U(i)*(dt/dy)*(T(i,j+1,k)-
                T(i,j,k))+Kh(j)*(dt/dy^2)*(T(i,j+1,k)-2*T(i,j,k)+T(i,j-1,k))+
                Kv(j)*(dt/dz^2)*(T(i+1,j,k)-2*T(i,j,k)+T(i-1,j,k));

                C(i,j,k+1)=C(i,j,k)-U(i)*(dt/dy)*(C(i,j+1,k)-
                C(i,j,k))+Kh(j)*(dt/dy^2)*(C(i,j+1,k)-2*C(i,j,k)+C(i,j-1,k))+
                Kv(j)*(dt/dz^2)*(C(i+1,j,k)-2*C(i,j,k)+C(i-1,j,k));
            else
                T(i,j,k+1)=T(i,j,k)-U(i)*(dt/dy)*(T(i,j,k)-
                T(i,j-1,k))+Kh(j)*(dt/dy^2)*(T(i,j+1,k)-2*T(i,j,k)+
                T(i,j-1,k))+Kv(j)*(dt/dz^2)*(T(i+1,j,k)-2*T(i,j,k)+T(i-1,j,k));

                C(i,j,k+1)=C(i,j,k)-U(i)*(dt/dy)*(C(i,j,k)-
                C(i,j-1,k))+Kh(j)*(dt/dy^2)*(C(i,j+1,k)-2*C(i,j,k)+
                C(i,j-1,k))+Kv(j)*(dt/dz^2)*(C(i+1,j,k)-2*C(i,j,k)+C(i-1,j,k));
            end
        end
    end
    
    T(1,:,k+1)= T(2,:,k+1);  %depth, dist, time
    T(31,:,k+1)=T(30,:,k+1); %assign top row of Ts

    for i=2:nz-1
       T(i,1,k+1)=  T(i,2,k+1);
       T(i,201,k+1)=T(i,200,k+1);
    end
end

And we can create interesting movies showing different plots by timestep by coding:

Before the time stepping loop:

v = VideoWriter('myvideo.avi');
open(v);

In the time stepping code, anywhere after the figure has the contents you want:

   frame = getframe(gcf);
   writeVideo(v,frame);

After the time stepping:

close(v);

A simple example of a video is a capture of the thermal and chemical gradients for the subject area by time (or grab the AVI format).

A Python-based solution is available in a format useful for a Jupyter Notebook.