We consider boundary element method (BEM) formulation of three-dimensional transient Stokes flow. We derive a representation of the Stokeslet and stresslet in terms of incomplete gamma functions and investigate the nature of the singularity of the single- and double layer potentials. Further, we derive the Green's representation formula for changing surfaces and investigate a problem with two moving cubes. Nystr\"om discretization does not work well for Stokes flow therefore we develop Galerkin discretization methods. We give analytical formulas for the time integration which results in a sequence of space-dependent integral operators and develop Galerkin schemes with tensor product piecewise polynomial ansatz functions. The time integrals of the matrix coefficients of the discretized linear system are expressed in terms of incomplete gamma functions and numerical results show right convergence rates. Furthermore, we implement the message passing interface (MPI) to solve a large size problem and see the efficiency of the parallel algorithm.