A combination of Landsat 8 and Sentinel-2 offers a high frequency of observations (3-5 days) at moderate spatial resolution (10-30 m), which is essential for crop yield studies. Existing methods traditionally apply vegetation indices (VIs) that incorporate surface reflectances (SRs) in two or more spectral bands into a single variable, and rarely address the incorporation of SRs into empirical regression models of crop yield. In this work, we address these issues by normalizing satellite data (both VIs and SRs) derived from NASA's Harmonized Landsat Sentinel-2 (HLS) product, through a phenological fitting. We apply a quadratic function to fit VIs or SRs against accumulated growing degree days (AGDDs), which affects the rate of crop development. The derived phenological metrics for VIs and SRs, namely peak, area under curve (AUC), and fitting coefficients from a quadratic function, were used to build empirical regression winter wheat models at a regional scale in Ukraine for three years, 2016-2018. The best results were achieved for the model with near infrared (NIR) and red spectral bands and derived AUC, constant, linear, and quadratic coefficients of the quadratic model. The best model yielded a root mean square error (RMSE) of 0.201 t/ha (5.4%) and coefficient of determination R 2 = 0.73 on cross-validation.is that methods for generating products from coarse spatial resolution sensors can be ported to moderate (Landsat 8/OLI, Sentinel-2/MSI) or high (Planet/PlanetScope) spatial resolution sensors. However, the practice shows that such a transition is not always straightforward due to larger data gaps because of clouds and uneven coverage, sensor characteristics, and increased spatial resolution (at least at the order of 100, when going from 250 to 30 m).Consider, for example, a crop yield assessment/forecasting application [4]. The hypothesis is that satellite-based features, such as vegetation indices (VIs) or biophysical parameters derived at a single date or accumulated over some time period, can be correlated to crop yields [5,6]. Since the reference data on crop yields are mainly available at the regional scale, the corresponding empirical models are built by averaging satellite-based features over those regions and correlating these derived variables with crop yields [7][8][9]. It is assumed that there is a homogeneity within the region in terms of crops grown and agricultural practices applied and, therefore, the averaging should be performed for satellite data acquired at the same (or approximately the same) stages of crop growth, meaning that the data are normalized. This is usually the case for coarse spatial resolution remote sensing sensors, which enable a higher likelihood of obtaining a high temporal frequency of cloud-free data over the Earth's surface [10,11]. This is also evidenced by multiple successful applications of coarse spatial resolution satellite data to crop yield assessment and forecasting [5,[7][8][9][12][13][14][15].However, this is not the case for moderate spatial resolution sat...