The current study seeks a fundamental explanation to the development of twodimensional coherent structures (2DCSs) in shallow mixing layers. A nonlinear numerical model based on the depth-averaged shallow water equations is used to investigate the temporal evolution of shallow mixing layers, where the mapping from temporal to spatial results is made using the velocity at the center of the mixing layers. The flow is periodic in the streamwise direction. Transmissive boundary conditions are used in the cross-stream boundaries to prevent reflections. Numerical results are compared to linear stability analysis, mean-field theory, and secondary stability analysis. Results suggest that the onset and development of 2DCS in shallow mixing layers are the result of a sequence of instabilities governed by linear theory, meanfield theory, and secondary stability theory. The linear instability of the shearing velocity gradient gives the onset of 2DCS. When the perturbations reach a certain amplitude, the flow field of the perturbations changes from a wavy shape to a vortical (2DCS) structure because of nonlinearity. The development of the vertical 2DCS does not appear to follow weakly nonlinear theory; instead, it follows mean-field theory. After the formation of 2DCS, separate 2DCSs merge to form larger 2DCS. In this way, 2DCSs grow and shallow mixing layers develop and grow in scale. The merging of 2DCS in shallow mixing layers is shown to be caused by the secondary instability of the 2DCS. Eventually 2DCSs are dissipated by bed friction. The sequence of instabilities can cause the upscaling of the turbulent kinetic energy in shallow mixing layers.