Laminar separation bubbles (LSB's) over the suction surface of a wing at low Reynolds number ( O (10 4 ) - O (10 6 ) based on the airfoil chord length) can significantly affect the aerodynamic performance of the wing, and pose a unique challenge for the predictive capabilities of simulation tools due to their high sensitivity to flow environments and wing surface conditions. In this work a series of two-dimensional (2D) and threedimensional (3D) low -order, and high -order accurate unstructured -grid -based numerical methods with varying model fidelity levels were used to study LSB physics over a NACA 0012 airfoil both in a clean freestream and in a turbulent freestream at a chord -based Reynolds number of 12,000. Lift production and time -averaged flow fields were compared with available experimental results. A major discovery is that in clean freestream flow a 3D high -order numerical scheme is necessary to capture LSB physics. This is due to the sensitivity of LSB-induced laminar -turbulent transition to flow conditions and boundary geometry at low Reynolds number. In freestream flows with moderate background turbulence (similar to 5%), 2D simulations failed to capture subtle 3D flow physics due to their intrinsic limitation, but can reasonably predict time -averaged airfoil performance. Similarity and distinction between freestream vortex-LSB interaction in 2D and eddy-LSB interaction in 3D were explained. The role of the Kelvin -Helmholtz instability and Klebanoff modes in the transition of 3D airfoils were shown to be critical for understanding laminar -turbulent transition and LSB formation on airfoils in clean and turbulent freestreams.