The determination of the temperature in and above the slab in subduction
zones, using models where the top of the slab is precisely known, is important
to test hypotheses regarding the causes of arc volcanism and intermediate-depth
seismicity. While 2D and 3D models can predict the thermal structure with high
precision for fixed slab geometries, a number of regions are characterized by
relatively large geometrical changes. Examples include the flat slab segments
in South America that evolved from more steeply dipping geometries to the
present day flat slab geometry. We devise, implement, and test a numerical
approach to model the thermal evolution of a subduction zone with prescribed
changes in slab geometry over time. Our numerical model approximates the
subduction zone geometry by employing time dependent deformation of a B\'ezier
spline which is used as the slab interface in a finite element discretization
of the Stokes and heat equations. We implement the numerical model using the
FEniCS open source finite element suite and describe the means by which we
compute approximations of the subduction zone velocity, temperature, and
pressure fields. We compute and compare the 3D time evolving numerical model
with its 2D analogy at cross-sections for slabs that evolve to the present-day
structure of a flat segment of the subducting Nazca plate