Four-dimensional variational (4D-Var) data assimilation is an effective assimilation method for obtaining physically consistent time-varying states. In this study, I propose a method using a neural network surrogate model obtained by machine learning to solve one of the most serious challenges in 4D-Var, which is to construct an adjoint model. The feasibility of the method was demonstrated by a 4D-Var experiment using a surrogate model for the Lorenz 96 model. Several effective procedures have been proposed to obtain an accurate surrogate model and the assimilated initial conditions: 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.