Many models have been used to examine subduction dip angle but none of which is both, dynamic and time dependent. Numerical models are developed to follow the thermal evolution of the subduction zone. The energy equation includes advection, the adiabatic gradient, viscous dissipation, radiogenic sources and variable diffusivity. The pressure is eliminated by combining the two momentum equations into a fourth, order stream-function equation. Both the energy and stream-function equations are solved numerically on a special non-uniform offset grid using the Altemating-Directon Implicit method. In a convection cell with Herring-Nabarro viscosity, circulation will be restricted to the left side of the cell while the descending limh is on the right. Altering the aspect ratio by changing the depth, has a much more pronounced effect on the total circulation and velocities than by changing the width. A convection cell will reorganize its circulation pattern to achieve steady-state when a highly viscous block, is "inserted” into the cell to disrupt normal circulation. Finally, the convective pattern in a system in which the boundary conditions are perpetually changing is determined by both its present boundary conditions and also by those preceding it by several millions of years.