Freak/rogue waves are considered to be the causes of marine accidents and their generation mechanism is closely related to the formation of wave groups. However, observations that capture the spatiotemporal evolution of coherent wave groups in directional windsea are rather limited. The paper presents a new technique known as the surface wave reconstruction by ensemble adjoint-free data assimilation (SWEAD) method that enables reconstruction of a spatiotemporal wave field covering a large area from wave records limited in observational density and spatial extent. We reconstructed spatiotemporal profiles of nonlinear surface gravity waves from virtual observational data using the adjoint-free four-dimensional variational data assimilation (a4DVar) scheme. The higher-order spectral method (HOSM) is used as a forward deep-water nonlinear wave model in a realistic sea state. The a4DVar scheme uses perturbed ensemble simulations to calculate the cost function gradient and Hessian; thus, construction of an adjoint model is not needed. A few extensions of the a4DVar scheme are proposed in this study. For efficient wave reconstruction, perturbed ensemble simulation results are reused by increasing the searching direction dimension at each iteration while assuring conformity to the perturbed model’s linearity. For regularization, Fourier coefficient magnitudes are constrained by a known power spectrum from the phase-averaged wave model. Twin experiments were conducted for a unidirectional wave with virtual wave gauge data and a multidirectional wave with virtual stereo camera imaging data. For both unidirectional and multidirectional cases, nonlinear freak wave–related wave groups were well reproduced, which is impossible using a linear model.