In this paper, an optimized complex nonnegative tensor factor 2D deconvolution (CNTF2D) is proposed to separate the sources that have been mixed in an underdetermined reverberant environment. Unlike conventional methods, the proposed model decomposition is performed directly on the statistics in the form of spectral covariance matrix instead of the data itself (i.e. the mixed signal). For faster convergence the model is adapted under the hybrid framework of the generalized expectation maximization and multiplicative update algorithms. This paper also proposes a solution to the issue of optimizing the model order i.e., number of components and convolutive parameters in the CNTF2D model. To this end, a latent-observation model based on Gamma-Exponential process is developed. In addition, the proposed Gamma-Exponential process can be used to initialize the parameterization of the CNTF2D. The proposed algorithm encodes a set of variable sparsity parameters derived from the Gibbs distribution. This permits a stable update and optimizes the CNTF2D with the correct degree of sparseness in the time-frequency domain. Experimental results on the underdetermined reverberant mixing environment have shown that the proposed algorithm is effective at separating the mixture with an average signal-to-distortion ratio of 2.5 dB.