Four-dimensional variational (4D-Var) data assimilation is an effective method for obtaining physically consistent time-varying states. In this study, a method using a neural network surrogate model obtained by machine learning is proposed to solve one of the most serious challenges in 4D-Var: to construct an adjoint model. The feasibility of the proposed method was demonstrated by a 4D-Var experiment using a surrogate model for the Lorenz 96 model. In the method, several effective procedures have been proposed to obtain an accurate surrogate model and the assimilated initial conditions, including two-stage learning (i.e., single- and multi-step learning) of neural networks, limiting the target states of the surrogate model to a small subspace of the state phase space, and updating the surrogate model during 4D-Var iterations.