  In fluid flow simulation, the multi-continuum model is a useful strategy.
When the heterogeneity and contrast of coefficients are high, the system
becomes multiscale, and some kinds of reduced-order methods are demanded.
Combining these techniques with nonlinearity, we will consider in this paper a
dual-continuum model which is generalized as a multi-continuum model for a
coupled system of nonlinear Richards equations as unsaturated flows, in complex
heterogeneous fractured porous media; and we will solve it by a novel
multiscale approach utilizing the constraint energy minimizing generalized
multiscale finite element method (CEM-GMsFEM). In particular, such a nonlinear
system will be discretized in time and then linearized by Picard iteration
(whose global convergence is proved theoretically). Subsequently, we tackle the
resulting linearized equations by the CEM-GMsFEM and obtain proper offline
multiscale basis functions to span the multiscale space (which contains the
pressure solution). More specifically, we first introduce two new sources of
samples, and the GMsFEM is used over each coarse block to build local auxiliary
multiscale basis functions via solving local spectral problems, that are
crucial for detecting high-contrast channels. Second, per oversampled coarse
region, local multiscale basis functions are created through the CEM as
constrainedly minimizing an energy functional. Various numerical tests for our
approach reveal that the error converges with the coarse-grid size alone and
that only a few oversampling layers, as well as basis functions, are needed.
