We present a study of a large-scale energy spectrum and integral invariants in temporally developing grid turbulence at mesh Reynolds numbers of ReM = 10 000 and 20 000 by employing direct numerical simulations (DNSs) in a periodic box. The simulations are initialized with a velocity field that approximates the wakes induced by the bars of conventional square grids. The turbulence statistics obtained in the temporal DNS agree well with those of the previous experiments in both the production and decay regions. The temporally developing grid turbulence also has a so-called non-equilibrium region, which is consistent with its spatially developing counterpart, where the normalized dissipation rate of turbulence kinetic energy (TKE), Cε, increases as the turbulence decays. The decay exponent n of TKE is n = 1.22 at ReM = 20 000 and n = 1.35 at ReM = 10 000, which are close to the values for the Saffman turbulence [i.e., 6/5 for ReM = 20 000 and 6(1 + p)/5 ≈ 1.36 for ReM = 10 000 with p ≈ 0.13 obtained by Cε ∼ tp at large t]. The longitudinal integral length scale and the TKE dissipation rate also exhibit temporal evolutions consistent with the Saffman turbulence for both ReM. The Saffman integral directly evaluated in the grid turbulence tends to be time-independent after the turbulence evolves for about 200 times of characteristic time scale defined by mesh size divided by the mean velocity of a fluid passing the grid. A direct evaluation of the TKE spectrum E(k) shows that E(k) ≈ Lk2/4π2 is valid for a finite range of low wavenumbers.