In this paper we addressed that the estimation of ocean surface heat flux from ocean temperature data is an ill-posed inverse problem, just like a well-known ill-posed problem: differentiation of noisy data. We reviewed engineering literature on such problems. Using a Mellor-Yamada level 2.5 closure model and simulated temperature data, we conducted numerical experiments of retrieving heat fluxes using variational data assimilation. This study shows that estimated nonsolar heat fluxes with high time resolution have large errors even if the observation errors are small. We also discussed the ill-posedness based on sensitivity coefficients. However, by applying some regularization methods and a longer assimilation window, it is feasible to estimate heat fluxes with reasonable accuracies, at least at some lower temporal resolution. On the basis of our model and experiment configurations, the high nonlinearity of the ocean mixed layer model also can cause difficulty in optimization when the assimilation window is long (e.g., 7 days). We proposed a modified an adjoint method approach and yielded good results.