Abstract. There are many applications of the least-squares finite element method for the numerical solution of partial differential equations because of a number of benefits that the leastsquares method has. However, one of most well-known drawbacks of the least-squares finite element method is the lack of exact discrete mass conservation, in some contexts, due to the fact that leastsquares method minimizes the continuity equation in L 2 norm. In this paper, we explore the reason of the mass loss and provide new approaches to retain the mass even in severely under-resolved grid. . On the other hand, least-squares methods for Stokes and Navier-Stokes equations lack an exact sense of discrete mass conservation. More precisely, least-squares methods do not naturally constrain the velocities to exactly satisfy any form of the mass conservation equations on discrete volumes as Galerkin and finite volume methods typically do. Instead, these methods merely strike a balance between the L 2 errors in each equation, with no special attention to the one that expresses mass conservation. Conservation error generally tends to vanish as the resolution or approximation order increases, but some applications require a level of discrete conservation that substantially exceeds the overall accuracy obtained by the discretization. In particular, for problems in which the domain is long and thin and only velocity boundary conditions are prescribed on all but a relatively small portion of the boundary, the use of low-order finite element spaces with a relatively coarse grid without further mitigation can result in severe mass loss. This deficiency can be overcome in two-dimensional problems by using higher-order finite element spaces, as the results in [15] show. However, in three dimensions, further measures may be necessary, which is the focus of this paper.In [9], mass conservation was improved by introducing newly defined variables and corresponding boundary conditions. Here, we instead suggest more fundamental ways of achieving improved mass conservation. We first identify a cause of the loss of mass in the least-squares finite element method in three-dimensional problems. In section 3, we show that boundary conditions near edges parallel to the flow contribute to mass-loss, even when high-order finite elements are employed. In section 3.1, a