Single-phase natural circulation thermosiphon loops have been attracting increased interest as they represent the prototype of passive safety systems. However, the stability properties of thermosiphon loops, which can affect and compromise their functionality, are still actively investigated. Traditionally, the stability analysis of thermosiphon loops has been simplified to one-dimensional (1D) calculations, on the argument that the flow would be monodimensional when the diameter of the pipe D is orders of magnitude smaller than the length of the loop Lt. However, at lower Lt/D ratios, rectangular thermosiphon loops show that the flow presents three-dimensional (3D) effect, which also has been confirmed by stability analyses in toroidal loops. In this paper, we performed a series of high-fidelity simulations using the spectral-element code nek5000 to investigate the stability behavior of the flow in rectangular thermosiphon loops. A wide range of Lt/D ratio from 10 to 200 has been considered, and the results show many different outcomes compared to previous 1D analytical calculations or stability theory. Moreover, we analyzed the flow in rectangular thermosiphon loops using proper orthogonal decomposition (POD), and we observed that the cases without flow reversal are characterized by swirl modes typical of bent pipes and high-frequency oscillation of the related time coefficients obtained by Galerkin projection. However, the swirl mode was not observed in cases with flow reversals, and these cases are characterized by symmetric flow field at second POD mode and the similarity of low-frequency oscillation in the projection of POD modes.