## 2 Methods and Process Procedures
### Brief Description of PSRFM  
The reconstruction quality of a cloud-free Sentinel-2 image time-series using our proposed approach essentially depends on two key processes. The first is to generate high quality synthetic Sentinel-2 images for cloudy and cloud shadow pixel replacement and the second is to identify cloudy and cloud shadow pixels that replace the corresponding pixels of the synthetic images. Herein, PSRFM [(24; 26)] was used to generate cloud-free synthetic Sentinel-2 image time-series, i.e., the first key process of the reconstruction procedure. We selected PSRFM for the image fusion processing because it had been evaluated [(24; 26)] against different landscape environments and dynamics, as well as being compared to published and well-known image fusion models, such as STARFM [(18)] and ESTARFM [(19)]. The evaluations and comparisons indicated that PSRFM outperformed STARFM and ESTARFM visually and quantitively.  
When interpolating estimates of surface reflectance, PSRFM requires two co-located image pairs of clear-sky fine and coarse resolution images acquired at start date \(t_{0}\) and end date \(t_{n}\) and a time-series of cloud-free coarse spatial resolution images acquired between the start and end dates donated as \(t_{1}\), \(t_{1}\),..., \(t_{n-1}\) for synthetic fine resolution image prediction. PSRFM can also be applied for extrapolation when only one pair of clear-sky fine and coarse spatial resolution images are available. For example, for predicting images earlier than the date of the first clear-sky image pair at \(t_{0}\) or images later than the date of the last clear-sky image pair at \(t_{n}\), a synthetic image can be generated by the forward or backward prediction.  
Figure 1: Reconstruction process flowchart. Date \(t_{k}\) indicates all acquisition dates when clear-sky images were not available.  
The modelling process involved a spatial partition of the fine resolution pixels sharing spectral similarity and fitting a linear transformation relating fine and coarse resolution spectra by unmixing each coarse resolution pixel based on land cover clusters. Then, the spectral change between the reference (\(t_{0}\) and \(t_{n}\)) and prediction dates (\(t_{k}\)) was predicted forward and backward, respectively. Following that, a smoothing process combined the forward and backward predictions through an average filter weighted by either the time elapsed between reference and predicted date or the uncertainties of the two predicted images [24].  
For the reconstruction process, we downloaded all Sentinel-2 images for which a cloud-free image time-series was required and then grouped every image into one of three classes: clear-sky images (cloud coverage \(<5\%\)), partially cloud-covered images (cloud coverage \(<50\%\)), and fully cloud-covered images (cloud coverage \(\geq 50\%\)). Here, an image with greater than \(50\%\) of cloud coverage was considered as full cloud coverage. This was based on the requirement of following cloud and cloud shadow identification methods for a reliable result. After the images were grouped, the clear-sky images were used as reference images to predict synthetic Sentinel-2 images at those dates when the images were partially or fully cloud-covered. For an image with full cloud coverage, it was entirely replaced by its corresponding synthetic image and for a partially cloud-covered image, its cloudy and cloud shadow pixels were replaced by the corresponding pixels of the synthetic image. Consequently, the final reconstructed cloud-free image time series consisted of three types of images: the clear-sky observations, the blended images that replaced the fully cloud-covered images, and the observation images with cloudy and cloud shadow pixels replacement from the synthetic images.