A Cluster Approach to Cloud Cover Classification over South America and Adjacent Oceans Using a k-means/k-means++ Unsupervised Algorithm on GOES IR Imagery

: An unsupervised k-means/k-means++ clustering algorithm was implemented on daily images of standardized anomalies of brightness temperature ( T b ) derived from the Geostationary Operational Environmental Satellite (GOES)-13 infrared data for the period 1 December 2010 to 30 November 2016. The goal was to decompose each individual T b image into four clusters that captures the characteristics of different cloud regimes. The extracted clusters were ordered by their mean value in an ascending fashion so that the lower the cluster order, the higher the clouds they represent. A linear regression between temperature and height with temperature used as the predictor was conducted to estimate cloud top heights (CTHs) from the T b values. The analysis of the results was performed in two different ways: sample dates and seasonal features. Cluster 1 is the less dominant one, representing clouds with the highest tops and variabilities. Cluster 4 is the most dominant one and represents a cloud regime that spans the lowest 2 km of the troposphere. Clusters 2 and 3 are entangled in the sense that both have their CTHs spanning the middle troposphere. Correlations between the monthly time series of the number of pixels in each cluster and of the entropy with several circulation indices are also introduced. Additionally, a fractal-related analysis was carried out on cluster 1 in order to resolve cirrus and cumulonimbus.


Introduction
In the recent years the scientific community has focused on cloud cover (CC) because of its role in the Earth's radiative balance (ERB). In 2013, the Intergovernmental Panel on Climate Change (IPCC) included for the first time a dedicated chapter to this topic alongside aerosols [1]. Clouds critically influence the radiation budget since they reflect a substantial portion of the incoming solar radiation back to space, and partly trap the outgoing terrestrial radiation and re-emit it to space at Northern Hemisphere, partly because of the number of validation opportunities the research community is presented with when compared with the Southern Hemisphere (SH). With a limited number of exceptions (e.g., [42]) no algorithms have hitherto been applied to, or particularly developed for the analysis of cloud systems in SA and adjacent oceans using GOES retrievals. The main reasons for undertaking an automated classification method for meteorological applications is thus linked to the huge volume of data that is available every day from the satellite and to the lack of studies in the region. In certain cases the analysis of several consecutive images is crucial to forecasting and nowcasting (e.g., severe weather events). As pointed out in [43] there is yet another motivation for this paper and it relies upon the identification of synoptic and subsynoptic cloud systems that follow the analysis of individual images. In order to provide further insights on CC classification in general, and for the SH in particular, this study is based upon an unsupervised k-means clustering algorithm with a k-means++ initializing algorithm that is individually applied on daily IR GOES-13 brightness temperature (T b ) data over Southern SA and the contiguous oceans for the 2010-2016 time span. The k-means algorithm has been extensively used to classify cloud top pressure/optical thickness joint histograms in order to obtain cloud regimes [44][45][46][47]. This k-means classification was extended with the incorporation of IR information [48]. All these studies were applied to gridded data. The current paper paves the way for the analysis of sequential images in the study region on a pixel-by-pixel basis, as well as for the characterization of synoptic and subsynoptic processes from a CC perspective.

Data
The GOES-13 dataset used in this study was obtained from the NOAA's Comprehensive Large Array-Data Stewardship System (CLASS, available at http://www.class.noaa.gov). Spanning the period 1 December 2010 to 30 November 2016 (2192 days), daily NetCDF-formatted channel 4 IR images taken at 11:45 UTC were selected for each of the days of the study period. Each pixel in the selected imagery represents an area of 4 km by 4 km; the central IR wavelength corresponds to 10.70 µm with a bandwidth of 1 µm [32] (pp. 11-12). Availability was over 97%, totaling 2136 images.
The basic IR quantities measured by the satellite sensors are raw digital counts, which are routinely converted to scaled radiances (SRs) packaged in 10-bit words in a format termed GVAR (for GOES Variable Format) [49]. Matlab's ncread function was used to read the GVAR values from the NetCDF files. A transformation of the GVAR values, i.e., calibration, is mandatory so as to express the SRs into a physical quantity. This calibration was carried out following NOAA's procedure to transform SR into T b [49]. It represents the intensity of the radiation at a certain frequency, and in this paper is the working variable. The calibration process requires three intermediate steps. The first one calculates the scene radiance (R). Expressed in mW m −2 sr −1 cm, R is calculated from SR using the linear transformation The relationship in (1) is extensive to all GOES IR channels; for the particular case of channel 4, the values of the scaling slope m and the scaling intercept b are 5.2285 and 15.6854, respectively [49]. The second calibration step transforms the R values into an effective temperature T eff using the inverse of the Planck function as follows In (2) v stands for the central wavenumber of the IR channel 4 (v=937.23449 cm -1 ) and c 1 and c 2 are equal to 1.191066 10 −5 mW m −2 sr −1 cm 4 and 1.438833 K cm, respectively [49]. The final transformation is the one from T eff to T b , which is made through a quadratic function T b =α+βT eff +γT eff 2 (3) For the IR channel 4 α, β and γ in (3) are equal to −0.52227011, 1.0023802 and −2.0798856 10 −6 , respectively [49].
It should be noted that the retrieval of some cloud properties from IR measurements at selected channels is sensitive not only to these desired properties but also to a number of atmospheric and surface parameters that may vary both spatially and temporally [50]. As a consequence, the T b values obtained through the aforementioned set of equations are subject to certain limitations due to the inherent sensitivity of the observed IR radiances. Measurements in the VIS and in the IR spectra can be overwhelmed by thicker/lower water clouds and so thin cirrus clouds have faint signals that are difficult to detect in both channels [51]. Following this, cirrus cloud top temperatures can be biased towards higher values and in consequence their top altitudes can be underestimated [52]. Moreover, subtle differences can emerge from the detection related to cloudy or cloud-free pixels depending on the band [53]. Further discussions regarding the use of a single IR channel are given below.
Each of the original 2136 individual images was transformed into T b images following the calibration process. Each resulting image was checked for consistency, and some of them were discarded for having either one or both of the following issues: (a) incomplete raw data information related to noise (49 images, 2.29%) and (b) images that have an anomalous scale of T b values (51 images, 2.39%). Examples of images with these two issues are shown in Figure 1. Excluding the images having issues (a) and (b), which accounted for 4.68% of the total, the number of analyzable images was 2036. Pixels that populate the lower blank portions at both sides of the original T b images have no valid latitudes or longitudes as they are beyond the edge marked by the Earth's curvature (see, e.g., Figure 1c). These pixels were removed from each original image prior to the analysis, so the input data consisted of non-zero values only. Prior to the clustering analysis, and in order to homogenize the set of usable images, each one of them was applied a month-by-month and pixel-by-pixel subtraction of the corresponding monthly mean and divided by the corresponding monthly standard deviation (SD). Figure S1 shows the monthly mean values and SDs. The homogenization process makes the input images comparable to each other. In particular, it removed the meridional gradient that is present in the mean values (cf. Figure S1).

The k-means/k-means++ Clustering Algorithm
Clustering is a widespread technique that is used in many applications, such as in data mining ([54] and references therein), in pattern recognition [55,56] and in machine learning problems ([57] and references therein). In particular, the k-means algorithm [58] is an exploratory procedure that is extensively used in pattern recognition and clustering ([59] and references therein). Given a field with data points p(i,j), i≤I, j≤J as input and an integer value K known a priori, the k-means method consists in finding K cluster centers (or centroids) so that the points assigned to a centroid are closer to it than they are to any other. The assignment is unique and the centroid values are usually obtained by weight-averaging all the points identified with this centroid [54]. Despite a number of shortcomings, simplicity and speed of convergence towards an optimal solution make this method one of the most useful across all disciplines [54]. Nevertheless, the performance of k-means was improved by combining it with the k-means++ initialization algorithm, which is a subprocess that seeds the centroids [60]. To this particular aim, k-means++ outperforms k-means both by achieving the classification task quicker and by a faster convergence to a minimal intra-class (intra-cluster) variance [60,61].
The k-means++ algorithm operates as follows: 1. A pixel from p(i,j) is randomly selected as the first centroid C 1 ; 2. All individual Euclidean distances from each pixel to C 1 , denoted as d(p,C 1 ), are computed; 3. The second centroid C 2 is randomly selected with probability.
Equation (4) shows that the farther the location of C 2 from C 1 the higher the probability of selecting it as the second centroid.
4. The process is repeated until all K centroids are obtained. Similar to step 3), the k-th centroid C k (3 ≤ k ≤ K) is selected from p(i,j) with probability.
A successful determination of the K centroids carries with the condition of the overall distance between each other being a maximum. Once the K centroids have been determined with the aid of k-means++, the process continues as with the standard k-means algorithm: 5. The distances d(p,C k ) for k ≤ K are computed; 6. Each grid point is assigned to the cluster with the closest centroid; 7. For k ≤ K, the average distance of all the pixels belonging to the k-th cluster is calculated so as to reassign this value to the corresponding centroid. The k-means process was repeated until the maximum number of predefined iterations was reached or when cluster assignments no longer changed. The kmeans/kmeans++ technique was individually implemented on each of the 2036 homogenized images using Matlab's kmeans and kmeans++ functions. Several replicates with random starting seeds were generated in order to avoid local minima and overfitting problems.
The establishment of the optimal value for the unknown parameter K that minimizes misclassification is a frequent issue. To this respect, the Caliński-Harabasz criterion (CHC) [62] was used here in order to estimate it. CHC establishes a variance ratio criterion (VRC) from which the optimal number of clusters K is estimated using the following quotient In (6) K represents the number of clusters, p is the number of pixels in the image, BGSS is the between-group (i.e., inter-cluster) sum of squares and WGSS is the within-group (i.e., intra-cluster) sum of squares. BGSS is made up from the summation of all the possible combinations of the squared distances calculated between two points in different clusters. The computation of WGSS involves the summation of all the possible squared distances in each individual cluster and then adding them up for all the K clusters. Since K is an unknown parameter the VRC was calculated using nine different values from K=2 to K=10 and the one that lead to the maximization of the VRC value was adopted. Maximization requires BGSS (WGSS) to be as greater (smaller) as possible in order to get a better separation of the groups. Since the optimal value of K is image-specific it was found to fluctuate in the 3-7 range for the 2036 input images. Figure 2 shows the way K distributes across the image dataset. There is a peak at K=4 (36% of the cases) followed by K=3 (28% of the cases). In order to use a single value of K the nearest integer from the average over the entire set of input images was established; this value turned out to be K=4 . A validation method using bootstrapping was carried out by randomly choosing pixels from the different input images in order to determine a single value of K for the resulting subsample image. This process was replicated several times and the outcome was also K=4 on average. Examples on two different dates whose input images have K=3 (4 July 2011) and K=7 (24 April 2013), i.e., the lowest and highest optimal values, respectively, are presented in order to confront their cluster sets with those for K=4. Individual cluster labeling was carried out over the mean value in an ascending fashion so that cluster 1 has the lowest mean in all sets. Figure 3 shows the input image and clusters 1-3 on 4 July 2011; Figure 4 shows the input image and clusters 1-4 on the same date. Similarly, Figure 5 shows the input image and clusters 1-7 on 24 April 2013 and Figure 6 shows the input image and clusters 1-4 on the same date. Table 1 shows the number of pixels in either set along with the percentage of pixels over the total on these two dates. The only commonality on either date is that the lowest (highest) order cluster includes the coldest (warmest) pixels. The change in the value of K yields a reassignment of pixels between the original and the new clusters. A split will occur when the value of K increases whereas a consolidation will take place when K decreases. On 4 July 2011, the rise in the value of K from 3 to 4 leads to a mild decrease in the number of pixels that populated the original cluster 1 that nets in a gain of pixels for the new cluster 2, the original cluster 2 is partly split between the new clusters 2 and 3, and the original cluster 3 is partly split between the new clusters 3 and 4. Initially, cluster 3 was the most populated cluster, including more than a half of the total number of pixels and representing the largest positive standardized anomalies. Following the splitting the most populated of the new clusters is cluster 3, and the largest positive standardized anomalies in cluster 4 account for a fifth of the total number of pixels. On the other hand, cluster 1 remains the less populated cluster after splitting. On 24 April 2013, the reduction in the value of K from 7 to 4 leads to the transfer of some of the pixels in the original cluster 2 to the new cluster 1, the rest of the pixels in the original cluster 2, the original cluster 3 and a small portion of the original cluster 4 are amalgamated into the new cluster 2, most of the contributions to the new cluster 3 comes from the amalgamation of the original clusters 4 and 5, and the original clusters 6 and 7 are amalgamated into the new cluster 4 with contributions from the original cluster 5 as well. After consolidation cluster 1 remains the less populated cluster and the new cluster 3 is the more populated cluster with half the total number of pixels. The new cluster 4, which includes the largest positive anomalies, represent more than a third of the total number of pixels. Table 1. Number of pixels included in each cluster for two selected dates that have the optimal number of clusters not coinciding with the fixed value of 4 used in the paper. Similar figures for the number of clusters set to 4 are also shown. The representativeness of each individual cluster within its group is presented in round brackets and it is estimated through the percentage over the total number of pixels in each image.

Cluster
In (7) p i stands for the "probability" of each cluster-estimated as the ratio between the number of pixels in the cluster and the total number of pixels in the input image. Given that the lower the value of the entropy the higher the information provided relation (7) can be used as an absolute value for comparisons between cluster sets regardless of the value of K. With the aid of Table 1 the calculation of the information entropy on 4 July 2011 gave 0.92 and 1.09 for K = 3 and K = 4, respectively. In other words, a fewer number of clusters as estimated by the VRC gave more information. On the other hand, the value of the entropy on 24 April 2013 was 1.68 and 1.09 for K = 7 and K = 4, respectively. This particular case stresses that an optimal number of clusters does not necessarily imply the best distribution of information across them.

Evaluation of Cloud Top Heights
The pixel value in each cluster was multiplied by the corresponding monthly SD and added the corresponding monthly mean in order to present the clusters in the original variable T b . A linear regression between temperature and height, using temperature expressed in Kelvin as the predictor, was conducted over reanalysis data in order to estimate cloud top heights (CTHs) from T b values. The reanalysis database used for these calculations consisted of air temperatures and geopotential heights retrieved from the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) reanalysis webpage (https://www.esrl.noaa.gov/psd/data/gridded/data.ncep.reanalysis.pressure.html) [64]. Each GOES pixel was associated to the closest grid point in the reanalysis. The regression was carried out for the later dataset between 1000 and 70 hPa so that information from both the troposphere and the lowermost stratosphere was considered. This was done for each of the 2036 individual dates and for the total number of 5,045,535 pixels in each image. The regression coefficients were used to associate T b values with CTHs. Given that the instrument on board the satellite did not collect information from below the Earth surface negative CTHs that resulted from the regression were not included in the analysis.

Results
The proper way to describe the outcome of the clustering process is to analyze the clusters of each image in the dataset individually. For the sake of brevity, only three sample dates were selected for this purpose. On the other hand, a description of the major features for the full analysis will be made from a seasonal perspective.

Weather Analysis
In order to provide a detailed description of the cluster classification results, three different sample dates were randomly selected. These dates are 1 August 2011 (austral winter), 4 January 2014 (austral summer) and 12 October 2016 (austral spring). Prior to analyzing the configuration of each cluster it is appropriate to describe at least the more prominent features of the CC that can be seen in the original T b images for these dates. An exhaustive analysis is however beyond the scope of the present work. Figure 7 shows the original T b monochromatic images on these three sample dates along with the corresponding colored T b values. It also shows the CTHs on 1 August 2011. . Mounted on a mid-latitude wavetrain, a frontal zone can be easily identified over the southern tip of SA, and there are visible signs of a cut-off circulation north of this front over central Argentina. Middle-to-high clouds seem to dominate in both of these two cloud systems. Even though Southern SA is under the influence of frontal activity all year round [66], generally speaking winter fronts are typically stronger due to an equatorward advance of baroclinicity, whose main driver is related to a strengthening of UTJs during the season across the entire SH [67]. Frequently associated to cyclogenesis [68], UTJs are more generally connected to the maintenance of the storm tracks due to their baroclinicity [69]. Figure S2 shows the geopotential height at ten different levels from 850 (lower troposphere) to 100 hPa (lower stratosphere) along with the wind field and the kinetic energy on 1 August 2011. It shows the presence of an UTJ, with maximum strength at 200 hPa in the leading part of an upper-level trough that is in a "tear-off" stage over central Argentina.
In lighter grey tones in Figure 7d, another frontal structure in a northwest-by-southeast position can be seen as a cloud band stretching from the east of the Ande-where several convective systems can be seen in the foreground-towards the South Atlantic. The T b values of the convective cells located in North-Western Bolivia and Eastern Brazil were around or below 210 K; convective-related clouds were also visible over Paraguay (Figure 7e). This cloud arrangement is characteristic of summer and is linked to an area of moisture flux convergence [70]. In fact, the region is associated to the local minima of outgoing longwave radiation where the well-known South Atlantic convergence zone (SACZ) develops as a result of the cold fronts reaching the region becoming stationary there ([71] and references therein). The south-easternmost tip of the cloud structure that sits along the SACZ spirals cyclonically (i.e., clockwise in the SH) well over the South Atlantic and bears a strong resemblance with an occlusion in line with the circulation features of the region [66].
The set of Figures 7f,g seems to have the largest number of distinct clouds between the three presented dates. Among these clouds, there were two contiguous mesoscale convective complexes (MCCs) situated in the eastern portion of central SA, approximately over Northern Uruguay and Southern Brazil, a region where MCCs tend to develop [72]. With an even texture, these very high clouds are characterized by lighter grey tones in both systems in Figure 7f, while Figure 7g shows that the T b values found at the top of these two MCCs were inclined towards the lowest values of the scale. The western, round-shaped MCC was smaller than the eastern one, from which cirrus plumes can be seen streaming off in its north-eastern portion following the westerly flow aloft. Figure 8 shows the configuration of clusters 1-4 on the three selected dates. Even though pixels represented by the lowest portion of the scale are found in each of the clusters their number generally decreased in favor of hotter pixels as the cluster order increased. Figure 9 shows the silhouette diagrams for these case studies. In general, each pixel was well matched to their own cluster. The information entropy for these dates was 0.96 (01/08/2011), 1.10 (04/01/2014) and 1.14 (12/10/2016; results not shown). These results mean that between the palettes of clouds present in these three days (cf. Figure 7) a better representation by means of a four-cluster set was achieved for the first of the dates. Another interpretation is that in terms of cloudiness the last date is the one that is more disordered.  According to the number of pixels in Figure 8a-c, cluster 1 accounts for approximately 8%, 4% and 7%, respectively, of the corresponding original T b images (not shown). Unlike with the rest of the clusters, pixels in cluster 1 on the selected dates did not extend beyond approximately 60 °S. Table 2 shows the mean values and the SDs for both T b and CTH for all the clusters in Figure 8. The mean T b value for Figure 8a-c was 232 K, 230 K and 228 K, respectively; SD equaled 11 K, 12 K and 11 K, respectively. In terms of CTH, the mean values were 10.71 km, 11.56 km and 10.26 km, respectively, and the SDs were 1.79 km, 2.01 km and 2.75 km, respectively. The distributions of T b and CTH for the clusters in Figure 8 are shown in Figure 10. Regarding cluster 1, Figure 10a shows that the T b pixels spanned the support's lower half, with values below approximately 260 K. The distributions in Figure 10b were more scattered than those in Figure 10a, highlighting the fact that the same temperature may be related to distinct heights at different regions. The Cis associated to the UTJ's curvature and the upper-level trough in Figure 7b or the convective cells over Bolivia, Brazil and Paraguay in Figure 7d were just examples of the CC represented by this cluster.  Cluster 2's number of pixels in Figure 8d-f represents 14%, 11% and 12% of the original T b images, respectively (not shown). The mean T b values for Figures 8d-f were 256 K, 253 K and 253 K, respectively, while the SDs were 15 K, 10 K and 13 K, respectively ( Table 2). As to CTH, the mean values were 6.04 km, 7.34 km and 6.06 km, respectively, and the SD values were 1.91 km, 1.63 km and 2.67 km, respectively. In line with the aforementioned results Figure 10c (10d) shows that the T b (CTH) distribution was shifted towards greater (lower) values of the corresponding scales. Figure  8d-f shows that pixels in cluster 2 spanned the entire domain. Nevertheless, the abundances for the colder pixels seemed to be greater for latitudes beyond 30 °S. Cloud decks that surround cluster 1's higher clouds in the region of the UTJ-associated Cis (cf. Figure 7a,b) were well represented by this cluster; most of the pixels along the SACZ and portions of the occlusion in the Atlantic (cf. Figure  7c,d) also were.

Cluster Analysis
Cluster 3's number of pixels in Figure 8g-i dramatically increased with respect to the former two clusters to account for 58%, 43% and 52%, respectively, of the total number of pixels in the original T b images (not shown). Pixels in these images spanned the entire domain. The mean T b values for Figure 8g-i were 278 K, 280 K and 279 K, while the SD was 13 K, 10 K and 11 K, respectively ( Table 2). In terms of CTH, the mean values were 2.16 km, 2.94 km and 2.98 km and the SD equaled 1.15 km, 1.39 km and 1.28 km, respectively. The T b and the CTH distributions ( Figure  10e,f) were in agreement with these results. Moreover, these two histograms were sharper than the previous ones. As with cluster 2, the colder pixels in Figure 8g-i spanned the domain's lower portion. A feature shared by these three figures was the presence of yellow-red pixels both off the Chilean and the Peruvian coasts and off the Angolan and the Namibian coasts (the African continent was out of sight in all the panels of Figure 8). These pixels were located in regions where persistent stratocumulus (Scs) occurs [73,74]. Scs are usually found below 2000 m independent of the Earth's region (WMO 2019). Some of the pixels that wrap around the occlusion in the Atlantic (cf. Figure 7d) were included in this cluster (Figure 8h). Figure 8j-l shows the configuration of cluster 4 on the selected dates. The number of pixels represented in them accounted for 20%, 42% and 29%, respectively, of the total number of pixels in the input images (not shown). Aside from Figure 8k, these percentages marked a difference in that they included fewer pixels than cluster 3. As with cluster 2 and 3, all latitudes were spanned by cluster 4's pixels on the selected dates. Reddish pixels were by far the most abundant in all three cases. The T b mean value for Figures 8j-l was 278 K, 285 K and 284 K, respectively, while the SD values were 12 K, 9 K and 11 K, respectively. As to CTH, the mean values were 1.22 km, 1.61 km and 1.44 km, respectively; the SD values equaled 0.92 km, 1 km and 0.84 km, respectively ( Table 2). The histograms of T b for clusters 3 and 4 bore a resemblance both in the shape and in the support. By contrast, the CTH cluster 4's distributions shifted towards the lowest values, with the class intervals that had a noticeable contribution extending through 5 km (Figure 10h).

Seasonal Features
The first part of this subsection is dedicated to the analysis of the mean values, SDs and frequency distributions from a seasonal perspective. Table 3 shows the seasonal grand-time mean values and SDs for both T b and CTH. The term grand-time mean was used here to reflect that the calculations were carried out within the entire period of analysis irrespective of the year, incorporating information from December, January and February (summer), March, April and May (autumn), June, July and August (winter), and September, October and November (spring). A general result from Table 3 is that the seasonal grand-time mean values for T b (CTH) increased (decreased) as the cluster order increased, and that the intra-seasonal variability decreased as the cluster order increased for both T b and CTH. With the exception of cluster 4, in terms of CTH winter presents the seasonal lowest mean values and autumn had the highest ones; this made the autumn-to-winter transition generally steeper than in any other season-to-season case. As to variability, it was greater in spring for clusters 1 and 2, in autumn for cluster 3 and in summer for cluster 4. The grand-time histograms for the four seasons are presented in Figure 11. They did not differ in a remarkable fashion from the corresponding ones shown in Figure 10.  Figure 11. As in Figure 10 but for the seasonal frequency distributions. These distributions were prepared including the entire set of 2036 analyzable images.
The seasonal character of the clusters represented in the columns of Table 3 warrants them to be termed "cloud regimes" to parallel similar plots that are present in the literature [75,76]. Clouds in cluster 1 represent high clouds including, but not restricted to, cirriform clouds, as this cluster is likely capturing the features of clouds that have a vertical extent as well, such as the cumulonimbus (Cbs). This is the less dominant cluster, representing, on average, 6% of the pixels. On the opposite side of the classification clouds in cluster 4 have their tops below 2 km. According to [77] they are classified as stratiform clouds. This is the second most dominant cluster as it represents, on average, 32% of the pixels. Clusters 2 (49%) and 3 (13%) spanned the middle troposphere. Altostratus and nimbostratus are examples of clouds that may populate them [77].
Given that seasonality implies that a number of different physical processes conspire to give the figures shown in Table 3 and Figure 11, it is difficult to draw conclusions about individual processes from these aggregated values. Monthly means were calculated for both the number of pixels that populate each cluster and the entropy. These values are shown in Figure 12. A noticeable feature in this figure is that the entropy closely followed the evolution of cluster 4. It can also be seen that the entropy increased when the number of pixels in clusters 3 and 4 were close to one another, i.e., when these two clusters became less distinguishable from each other. This tended to occur in fall and winter. In an endeavor to provide further insights on the intraseasonal variability and to shed light on the possible physical processes involved the mean values in Figure 12 were correlated with several circulation indices both on a grand-time seasonal basis and for the entire period, i.e., disregarding season or year. The calculations were carried out using the Spearman's correlation coefficient (ρ), which acts over the ranks of the pair of correlated variables (rather than over the raw values) and permits establishing whether there exists a relationship between the variables by means of a monotonic increasing or decreasing function [78]. It is worth noting that the method did not disclose the specific functional relationship between the variables, it just assessed the degree of monotonicity. The circulation indices that were chosen for this specific purpose were the Antarctic  Table 4, which only included significant values. The DMI was excluded from the table as the ρ values resulted as not significant in all cases.       with the cloudiness present in synoptic-scale processes (cf. Figure 12) and their pixels spreading 70 across the entire study region (cf. In the above equation C is a constant. The linearization of (9) leads to ln P = ln C + D(P) 2 ln A (10) from which both the intercept ln C and the slope D(P) 2 ⁄ can be easily obtained from a linear fit ones in [95], were used to classify the pixels in cluster 1.
106 Figure 13a shows the breakdown of Figure 8a into groups of pixels of Cis and Cbs. In this 107 particular case Cis and Cbs represent approximately 89% and 9% of cluster 1's total number of 108 pixels, respectively. The remaining 2% of the pixels correspond to the clouds that were unaccounted 109 for by the condition A > 200 km 2 . The largest cluster seen in Figure 13a was the one associated to Upper-level frontogenesis, which is most frequent in winter due to an increased baroclinicity (cf.
into the troposphere destabilize this latter layer and may trigger convection [100-102]. The detection of a greater number of Ci and Cb pixels is therefore predicted to occur in the winter months, as 123 shown in Figure 13b.

162
In principle, the four clusters can be associated to different types of clouds. For instance, cluster 4 163 seems to be associated to shallow clouds that span the lowest 2 km of the troposphere. On the other 164 hand, cluster 1 tends to represent clouds whose tops populate the upper troposphere. Clusters 2 and

171
Another point that is worth mentioning is the use of vertical profiles to establish, through linear 172 regressions, a relationship between temperature and height. Non-linear fits to the temperature 173 profiles could also be conducted. Even though this alternative is expected to improve the outcome 174 that relates the two variables, it is worth emphasizing that the results described in this paper are 175 qualitatively in correspondence with the existing literature.

176
Regarding the partition of cluster 1 into Ci and Cb pixels the results is encouraging considering