A number of studies have found that the formation of double twins in low symmetry metals can lead to the onset of strain localizations, leading further to void nucleation and ultimately fracture. This work extends a recently developed three-dimensional crystal plasticity finite element framework [1] to explicitly model kinematics and kinetics of nucleation/formation, propagation, and thickening/growth of a discrete double twin lamella ({101¯1}−{101¯2}) in a magnesium alloy AZ31. With this approach, morphological and crystallographic reorientations as well as the shear transformation strains associated with strain accommodation by the double twinning sequence are modeled during simple compression and tension. The simulations predict that the distribution of local stress-strain fields during formation and growth of primary contraction twin creates the driving force for the formation of a secondary extension twin variant, which is consistent with experimental observations in both compression and tension. In particular, the contraction twin variant (01¯11)[01¯12¯] is predicted to form an internal extension twin variant (011¯2)[01¯11]. Furthermore, the prediction of the underlining crystallographic slip deformation mechanism reveals a substantial activity of basal slip within the contraction portion of the double twin, causing strain localization in its vicinity. Finally, the simulations reveal a gradient in the traction force field across twin-parent interface, suggesting that contraction twin-parent boundaries are weak links in the microstructure where voids can nucleate.