Geophysical flows are often two-dimensional (2D) in nature. The decay dynamics of the canonical, symmetric initial conditions defined by the Taylor–Green (TG) flow is simulated in a 2D domain. The standard lattice Boltzmann equation method with the Bhatnagar–Gross–Krook collision model is employed. The vortex dynamics is investigated for bounded and unbounded flows with various sets of Gaussian vortices and varying Reynolds numbers. The formation of thin filaments and the generation of small-scale structures near the no-slip walls are found. To characterize the decay properties, the integral quantities, namely, kinetic energy, E(t), enstrophy, Ω(t) , and palinstrophy, P(t), are tracked as a function of the eddy-turnover-time. The initial vortex population in the domain dictates the decay rate of these integral quantities. The energy spectra, E(k), are computed to investigate the effect of domain boundaries and vortex populations on energy distribution among different length scales for Re=50000 . The evolution of TG flow in a bounded domain shows an exponent close to −3, which may be attributed to the generation of thin filaments and small-scale structures in vortex-wall interactions. Exponent changes with an increase in the vortex populations.