For the first time three dimensional simulation is conducted for vertical zone-melting crystal growth, simultaneously considering the time-dependent fluid flow, heat transfer, and moving interfaces. The numerical method is based on an efficient finite-volume method with front tracking. The growth of 2-cm diameter GaAs in a quartz ampoule is used as an example. Under an axisymmetric operation, the two-dimensional solution can bifurcate subcritically to three-dimensional ones as the convection or superheating increases. The three-dimensional modes, which can have either one- or two-fold symmetry, depend on the equilibrium aspect ratio of the molten zone. Under normal gravity, the axisymmeric mode exists only for low superheating. Furthermore, the multiple slates, including other multi-fold modes, are path dependent, which are further illustrated though the dynamics of heating and cooling operations.