Generally, the sea ice prediction skills can be improved by assimilating the observed sea ice data into a numerical forecast model to update the initial fields of the model. Meanwhile, it is necessary to assimilate sea ice thickness (SIT) while assimilating sea ice concentration (SIC) to keep the two harmony in assimilation. However, due to the lack of the SIT from satellite remote sensing observation, it cannot meet the bivariate assimilation requirement in Arctic melting season. In order to solve this problem, an easy-to-implement bivariate regression mode of SIT is tentatively established based on the grid reanalysis data of SIC and SIT, through which the SIT field is statistically constructed. Then, the ice-ocean coupled numerical forecast experiment is carried out in which both the observed SIC and the constructed SIT are jointly assimilated using the spatial multi-scale recursive filter (SMRF) method. Results show the joint assimilation of SIC and the constructed SIT can greatly improve the forecast accuracy of sea ice elements especially in the multi-year ice region of Arctic center, where the average absolute error between the SIT forecast and in situ observations is about 0.14 m. Further, effects of the bivariate assimilation on the ocean elements are also deeply investigated in melting season. The higher forecast skill of sea surface temperature and drift flow can be obtained via the bivariate assimilation scheme considering the ice-ocean coupled dynamics and the feedback process between them.