remain separated for longer distances (long bubble globally modifying wall-pressure distribution [1]). In both cases, airfoil lift, drag and moment are influenced by the presence of the bubble. RANS (Reynolds-Averaged Navier-Stokes) turbulence models require specific transition information and/or modification to handle the separation-induced transition which dominates the flow [7, 10]. Numerous early [1, 6, 9, 11, 12] and more recent [8, 13-15] studies, both experimental and computational, highlight the complex physics of the transitional flow inside LSBs, such as pronounced three-dimensionality, high unsteadiness, turbulent breakdown driven by multiple coherent structures and complex oscillator and amplifier instability mechanisms that may lead to notable changes in the flow topology. Several approaches of varying degree of empiricism have been used to predict this difficult flow. Ultimately, improved transport-equation based transition models [16] are expected to accurately predict such flows. Correlation-based models [17, 18] reformulate widely used transition correlations [19] into coefficients of transport equations for the intermittency [20, 21] and other transition parameters, to obtain local, and recently [22] Galilean invariant, transition models. Physics-based models [23-25] generally rely on the concept [26] of laminar (pre-transitional) kinetic energy k L ≈ 1 2 u 2 (essentially streamwise, in agreement with transition physics [27] and with the observed rapid increase of the Reynolds-stress tensor anisotropy in very low Reynolds number channel flows [28]), which is computed by a specific transport equation, to trigger and control transition in the turbulence model. Transport-equation based transition models are quite successful in mimicking transitional flow effects and in detecting transition [16]. Given the extreme difficulty of accurately predicting transition, algebraic nonlocal transition models are also developed [29], including specific adaptations of additional transport equations [30-32]. Numerous authors have contributed to RANS modeling of LSBs on airfoils, especially at low chord-based Reynolds number Re c. Yuan et al. [33] and Windte et al. [34] performed extensive RANS modeling of LSBs for the flow around a SD7003 airfoil at Re c =60,000, with various models, and concluded that Menter's BSL [35] k-ω performed quite well. Source terms were disabled in the laminar regions upstream of transition onset which was determined using an approximate envelope method [34]. Lian and Shyy [36] proposed a nonlocal intermittency model for the activation of Wilcox's 1994 k-ω [37], through direct weighting of the eddy viscosity ν effective