Vertical axis turbine (VAT) arrays can achieve larger power generation per land area than their horizontal axis counterparts, due to the positive synergy from clustering VATs in close proximity. The VATs generate a three-dimensional wake that evolves unevenly over the vertical and transverse directions according to two governing length scales, namely the rotor's diameter and height. Theoretical wake models need to capture such a complex wake dynamics to enable reliable array design that maximises energy output. This paper presents two new theoretical VAT wake models based on super-Gaussian and Gaussian shape functions, which account for the three-dimensional velocity deficit distribution in the wake. The super-Gaussian model represents the initial elliptical shape with the superposition of vertical and lateral shape functions that progressively converge into an axisymmetric circular-shaped wake at a downstream distance that depends on the rotor's height-to-diameter aspect ratio. Our Gaussian model improves the initial wake width prediction taking into account the rectangular rotor's cross-section. Our models were well validated with large-eddy simulations (LES) of single VATs with varying aspect ratios and thrust coefficients operating in an atmospheric boundary layer. The super-Gaussian model attained a good agreement with LES in both near and far wake, whilst the Gaussian model represented well the far-wake region.