The wrinkling phenomenon occurring in nonlinear elastic membrane structures has been studied by numerous methods for many decades. However, it remains a challenging issue to model wrinkled membranes appropriately because they exhibit minimal resistance to compression and bending. As a result, most compressive stresses within membranes are released through localized instability. It leads to out-of-plane displacements, which are typically referred to as wrinkles. In general, the previous numerical approaches dealing with wrinkling can be divided into two categories. The first relies on the buckling theory to resolve the wrinkle details explicitly with a highly dense mesh of shell elements, causing expensive computation costs. An alternative is to use a coarse mesh of membrane elements enhanced with proper wrinkling models grounded in the so-called tension field theory [1]. It is primarily concerned with predicting average stress and displacement fields rather than providing adequate information regarding wrinkle geometry. Moreover, this treatment may lead to convergence problems due to the absence of bending stiffness and sudden changes in the tangent stiffness matrix. Hence, it is necessary to develop a robust and suitable wrinkling model. We propose a novel variationally consistent membrane wrinkling model for analyzing the mechanical responses of wrinkled thin membranes. The elastic strain energy density is split into tensile and compressive terms via a spectral decomposition of the strain tensor. Tensile and compressive parts of the stress and tangent material matrices are then obtained via consistent variation from the respective strain energies. Considering only the positive part of the strain energy in the variational formulation, we obtain a membrane with zero compressive stiffness. By adding the negative strain energy multiplied with a reduction factor, we further obtain residual compressive stiffness, which improves stability and allows handling also states of slackening. Comparison with results from analytical, numerical, and experimental examples from the literature on membrane wrinkling problems demonstrate the great performance and capability of the proposed approach, which is also compatible with commercial finite element software.