Numerical studies of decaying two-dimensional (2D) magnetohydrodynamic (MHD) turbulence are presented. Results concern mainly the decay laws of integral quantities and the scaling properties of spectra and structure functions. The self-similarity of the decay is characterized by the invariance of the ratios of resistive and viscous dissipation rates, kinetic and magnetic energies, and mean parallel and perpendicular velocities. The energy decay rate is primarily determined by the conservation of the mean-square magnetic potential, but the process is somewhat retarded at intermediate times by current density condensation in the magnetic eddies. In contrast to the 3D case, the energy spectrum agrees with the Iroshnikov–Kraichnan (IK) law Ek∼k−3/2, but is modified by a strong anomalous bottleneck effect. Structure functions of the Elsässer fields satisfy the property of extended self-similarity, which provides accurate values of the relative scaling exponents. Also the structure functions themselves show clear scaling properties, particularly at the highest Reynolds numbers studied. The scaling exponent ζ2 is consistent with the IK spectrum, but ζ3 and ζ4 remain distinctly below unity, which is surprising in view of an exact relation suggesting ζ3=1. Finally, the role of the tearing mode in 2D MHD turbulence is revisited.