The transient seepage in anisotropic media can be moddled by an elliptical partial differential equation with kinematic boundary conditions. Introducing an eigen-relation, the general solution and the boundary conditions of the plane problem can be expressed by terms of complex analytical function vector as $f_k(z_k)$ with complex eigenvalues $z_k=x+p_k y$. A new conjugate boundary element method for holomorphic function (CBEM) has been proposed by author recently (see DOI:10.21656/1000-0887.370001[in Chinese]). Employing CBEM and transforming the original plane defined in into several isotropic planes defined in with , a number of weighted residual equations are established as the equivalent weak form to the original governing equations. Following the way similar to the conventional BEM of potential problem, a series of linear equation system with constraned equations can be obtained, in which multiple BEM equations of potential problems are employed and Cauchy-Riemann realtions are used as constrained linear equations in discrete space. Overcoming the mathematical difficulty arising from domain decomposition, domain transformationt and solution of large linear system with weak ill-conditioned matrices, a complete MCBEM for seepage in anisotropic media will be finally established. Three transient free-surface seepage in phreatic surface studied by Rafiezadeh and Ataie-Ashtiani (ref. to EABE 46(2014):51-66) are used as numerical examples, and the results show the reliability and flexibility of the proposed method. Comparing with the other numerical methods to seepage analysis, the most remarkable advantage of the proposed method is suitable for parallel computing, which is discussed at the end of this paper.