Phytoplankton blooms in polar regions with seasonal sea ice cover show a unimodal seasonality. However, the bloom processes are controlled by multiple physical and biogeochemical factors, including sea ice, light availability, mixed layer depth, and nutrients; those may result in great uncertainties in simulating phytoplankton bloom by the Earth System Models (ESMs). In this study, the results of 11 Coupled Model Intercomparison Phase-6 (CMIP6) ESMs were analyzed and evaluated with various types of observational products in order to determine whether those ESMs can correctly model the phytoplankton blooms in three Arctic shelf seas, Barents Sea, Chukchi Sea, and Bering Sea. By calculating multiple indices that represent light and nutrient limitations, the error sources of simulated surface chlorophyll a concentrations were comprehensively analyzed. Our results show that the 11 ESMs can be divided into three groups based on ice-adjusted photoperiod, rate of change of mixed layer depth, and surface nitrate concentration. Some groups are characterized by the smallest bias between modeled indices and observation-based reference, and those ESMs perform best in simulating phytoplankton bloom characteristics. The other groups of ESMs differ significantly from the reference values in terms of surface nitrate and/or rate of change of mixed layer depth, resulting in delayed occurrences of annual chlorophyll a peak concentration and greater differences in corresponding peak values. In general, in addition to the two primary constraints of light and nutrients, the ESMs should also well represent the upper mixed layer controlled by temperature and salinity distributions, so as to accurately simulate the seasonal variation of surface chlorophyll a concentration. The above analyses indicate ESMs can be used in assessing polar planktonic ecosystems, and there is room for improving ecosystem-related parametrization in future ESM development.