The existence of two regimes for El Niño (EN) events, moderate and strong, has been previously shown in the GFDL CM2.1 climate model and also suggested in observations. The two regimes have been proposed to originate from the nonlinearity in the Bjerknes feedback, associated with a threshold in sea surface temperature (Tc) that needs to be exceeded for deep atmospheric convection to occur in the eastern Pacific. However, although the recent 2015–16 EN event provides a new data point consistent with the sparse strong EN regime, it is not enough to statistically reject the null hypothesis of a unimodal distribution based on observations alone. Nevertheless, we consider the possibility suggestive enough to explore it with a simple theoretical model based on the nonlinear Bjerknes feedback. In this study, we implemented this nonlinear mechanism in the recharge-discharge (RD) ENSO model and show that it is sufficient to produce the two EN regimes, i.e. a bimodal distribution in peak surface temperature (T) during EN events. The only modification introduced to the original RD model is that the net damping is suppressed when T exceeds Tc, resulting in a weak nonlinearity in the system. Due to the damping, the model is globally stable and it requires stochastic forcing to maintain the variability. The sustained low-frequency component of the stochastic forcing plays a key role for the onset of strong EN events (i.e. for T> Tc), at least as important as the precursor positive heat content anomaly (h). High-frequency forcing helps some EN events to exceed Tc, increasing the number of strong events, but the rectification effect is small and the overall number of EN events is little affected by this forcing. Using the Fokker–Planck equation, we show how the bimodal probability distribution of EN events arises from the nonlinear Bjerknes feedback and also propose that the increase in the net feedback with increasing T is a necessary condition for bimodality in the RD model. We also show that the damping strength determines both the adjustment time-scale and equilibrium value of the ensemble spread associated with the stochastic forcing.