Let $${{\mathcal {Q}}}\subset {{\mathbb {C}}}\left $$ be an arbitrary set of polynomials in noncommutative indeterminates such that $$q(0)=0$$ for all $$q\in {{\mathcal {Q}}}$$ . The noncommutative variety $$\begin{aligned} {{\mathcal {V}}}_{f,{{\mathcal {Q}}}}^m({{\mathcal {H}}}):=\left\{ { X}=(X_1,\ldots , X_n)\in \mathbf{D}_f^m({{\mathcal {H}}}): \ q({ X})=0 \ \text { for all } q\in {{\mathcal {Q}}}\right\} , \end{aligned}$$ where $$\mathbf{D}_f^m({{\mathcal {H}}})$$ is a noncommutative regular domain in $$B({{\mathcal {H}}})^n$$ and $$B({{\mathcal {H}}})$$ is the algebra of bounded linear operators on a Hilbert space $${{\mathcal {H}}}$$ , admits a universal model $$B^{(m)}=(B_1^{(m)},\ldots , B_n^{(m)})$$ such that $$q(B^{(m)})=0$$ , $$q\in {{\mathcal {Q}}}$$ , acting on a model space which is a subspace of the full Fock space with n generators. In this paper, we obtain a Beurling type characterization of the joint invariant subspaces under the operators $$B_1^{(m)},\ldots , B_n^{(m)}$$ , in terms of partially isometric multi-analytic operators acting on model spaces. More generaly, a Beurling-Lax-Halmos type representation is obtained and used to parameterize the wandering subspaces of the joint invariant subspaces under $$B_1^{(m)}\otimes I_{{\mathcal {E}}},\ldots , B_n^{(m)}\otimes I_{{\mathcal {E}}}$$ , and to characterize when they are generating for the corresponding invariant subspaces. Similar results are obtained for any pure n-tuple $$(X_1,\ldots , X_n)$$ in the noncommutative variety $${{\mathcal {V}}}_{f,{{\mathcal {Q}}}}^m({{\mathcal {H}}})$$ . We characterize the elements in the noncommutative variety $${{\mathcal {V}}}_{f,{{\mathcal {Q}}}}^m({{\mathcal {H}}})$$ which admit characteristic functions, develop an operator model theory for the completely non-coisometric elements, and prove that the characteristic function is a complete unitary invariant for this class of elements. This extends the classical Sz.-Nagy–Foias functional model for completely non-unitary contractions, based on characteristic functions. Our results apply, in particular, when $${{\mathcal {Q}}}$$ consists of the noncommutative polynomials $$Z_iZ_j-Z_jZ_i$$ , $$i,j=1,\ldots , n$$ . In this case, the model space is a symmetric weighted Fock space, which is identified with a reproducing kernel Hilbert space of holomorphic functions on a Reinhardt domain in $${{\mathbb {C}}}^n$$ , and the universal model is the n-tuple $$(M_{\lambda _1},\ldots , M_{\lambda _n})$$ of multipliers by the coordinate functions.