A novel deep learning approach for data assimilation of complex hydrological systems is introduced, addressing the limitations of traditional methods that rely on Gaussian assumptions or suffer from low computational efficiency. The proposed method, called DA(DL), leverages deep learning (DL) to model non-linear relationships and recognize complex patterns, making it suitable for problems involving non-linearity, high-dimensionality, and non-Gaussianity. DA(DL) generates training data from a prior ensemble and trains a DL model to update system knowledge, such as model parameters. For highly non-linear models, an iterative form of DA(DL) is implemented. Strategies like data augmentation and local updating are introduced to enhance DA(DL) for small ensemble sizes and equifinality issues. In two hydrological DA cases involving Gaussian and non-Gaussian distributions, DA(DL) shows promising performance compared to ensemble smoother methods like ES(K) and ES(DL). Potential improvements to DA(DL) include designing better DL model architectures, imposing physical constraints, and further updating important variables like model states, forcings, and error terms. The method is tested in a non-Gaussian groundwater flow case, where DA(DL) outperforms ES(K) and ES(DL) in estimating the spatial distribution of hydraulic conductivity. In a Gaussian case, DA(DL) also performs well, demonstrating its versatility and practicality as a reliable DA method. The results highlight the effectiveness of DA(DL) in handling complex hydrological systems with non-linear and non-Gaussian characteristics.A novel deep learning approach for data assimilation of complex hydrological systems is introduced, addressing the limitations of traditional methods that rely on Gaussian assumptions or suffer from low computational efficiency. The proposed method, called DA(DL), leverages deep learning (DL) to model non-linear relationships and recognize complex patterns, making it suitable for problems involving non-linearity, high-dimensionality, and non-Gaussianity. DA(DL) generates training data from a prior ensemble and trains a DL model to update system knowledge, such as model parameters. For highly non-linear models, an iterative form of DA(DL) is implemented. Strategies like data augmentation and local updating are introduced to enhance DA(DL) for small ensemble sizes and equifinality issues. In two hydrological DA cases involving Gaussian and non-Gaussian distributions, DA(DL) shows promising performance compared to ensemble smoother methods like ES(K) and ES(DL). Potential improvements to DA(DL) include designing better DL model architectures, imposing physical constraints, and further updating important variables like model states, forcings, and error terms. The method is tested in a non-Gaussian groundwater flow case, where DA(DL) outperforms ES(K) and ES(DL) in estimating the spatial distribution of hydraulic conductivity. In a Gaussian case, DA(DL) also performs well, demonstrating its versatility and practicality as a reliable DA method. The results highlight the effectiveness of DA(DL) in handling complex hydrological systems with non-linear and non-Gaussian characteristics.