History matching is crucial for reliable numerical simulation of geological carbon storage (GCS) in deep subsurface aquifers. This study focuses on inferring highly complex aquifer permeability fields with multi- and intra-facies heterogeneity to improve the characterization of CO2 plume migration. We propose a deep learning (DL)-based parameterization strategy combined with the ensemble smoother with multiple data assimilation (ESMDA) algorithm to formulate an integrated inverse framework. The DL model is employed to parameterize non-Gaussian permeability fields using low-dimensional latent variables in a Gaussian distribution, thereby mitigating the non-Gaussianity issue faced by the ensemble-based ESMDA inverse method and simultaneously alleviating the computational burden of high-dimensional inversion. The efficacy of the integrated DL-ESMDA inverse framework is demonstrated using a 3-D GCS model, where it estimates the non-Gaussian permeability field characterized by multi- and intra-facies heterogeneity. Results show that the DL model is able to represent the highly complex and high-dimensional permeability fields using low-dimensional latent vectors. The DL-ESMDA framework sequentially updates these low-dimensional latent vectors instead of the original high-dimensional permeability field to obtain posterior estimations of the permeability field. The resulting CO2 plume migration closely matches historical measurements, suggesting a significantly improved model reliability after history matching. Additionally, a substantial reduction in uncertainty for future plume migration predictions beyond the history matching period is observed. The proposed framework provides an effective approach for reliable characterization of CO2 plume migration in highly heterogeneous aquifers, enhancing GCS project operation and risk analysis.