  Gaussian process is an indispensable tool in clustering functional data,
owing to it's flexibility and inherent uncertainty quantification. However,
when the functional data is observed over a large grid (say, of length $p$),
Gaussian process clustering quickly renders itself infeasible, incurring
$O(p^2)$ space complexity and $O(p^3)$ time complexity per iteration; and thus
prohibiting it's natural adaptation to large environmental applications. To
ensure scalability of Gaussian process clustering in such applications, we
propose to embed the popular Vecchia approximation for Gaussian processes at
the heart of the clustering task, provide crucial theoretical insights towards
algorithmic design, and finally develop a computationally efficient expectation
maximization (EM) algorithm. Empirical evidence of the utility of our proposal
is provided via simulations and analysis of polar temperature anomaly
(\href{https://www.ncei.noaa.gov/access/monitoring/climate-at-a-glance/global/time-series}{noaa.gov})
data-sets.
