The unsteady Reynolds-averaged Navier–Stokes (URANS) equations provide a computationally efficient tool to simulate unsteady turbulent flows for a wide range of applications. To account for the errors introduced by the turbulence closure model, recent works have adopted data assimilation (DA) to enhance their predictive capabilities. Recognizing the challenges posed by the computational cost of four-dimensional variational DA for unsteady flows, we propose a three-dimensional DA framework that incorporates a time-discrete Fourier transform of the URANS equations, facilitating the use of the stationary discrete adjoint method in Fourier space. Central to our methodology is the introduction of a corrective, divergence-free, and unsteady forcing term, derived from a Fourier series expansion, into the URANS equations. This term aims at mitigating discrepancies in the modeled divergence of Reynolds stresses, allowing for the tuning of stationary parameters across different Fourier modes. While designed to accommodate multiple modes in general, the basic capabilities of our framework are demonstrated for a setup that is truncated after the first Fourier mode. The effectiveness of our approach is demonstrated through its application to turbulent flow around a two-dimensional circular cylinder at a Reynolds number of 3900. Our results highlight the method's ability to reconstruct mean flow accurately and improve the vortex shedding frequency (Strouhal number) through the assimilation of zeroth mode data. Additionally, the assimilation of first mode data further enhances the simulation's capability to capture low-frequency dynamics of the flow, and finally, it runs efficiently by leveraging a coarse mesh.