The aim of this study was to propose a methodology that provides a detailed description of the argillic zone of a hydrothermal field, based on satellite multispectral data. More specifically, we developed a method based on spectral unmixing where hydroxyl-bearing alteration is represented by a single endmember (representing clays) and the three (nearly) non-altered primary volcanic lithologies, namely, two types of lava flows (basic and acidic compositions) and the loose materials (alluvial/beach deposits, scree, pyroclastic deposits, etc.), are represented by three endmembers. We also used one endmember representing elemental sulfur that is present in fumarolic vents hosted by active hydrothermal craters. The methodology was applied in the south part of Lakki plain inside the Nisyros volcano caldera (Greece), using Sentinel-2, Landsat-8/OLI, and ASTER satellite multispectral datasets. Specifically, it was applied separately to each one of the three datasets. The spectral unmixing results, combined with the relative geological map, provide quantitative estimations of the primary volcanic and loose material areas affected by alteration. In addition, pixels with high abundance values of hydroxyl-bearing alteration corresponded to mapped areas with strong hydrothermal alteration. The developed methodology is superior to conventional approaches (e.g., alteration spectral index) in terms of its ability to describe the overall pattern of the hydrothermal field. The most accurate results were taken when applied to ASTER or Sentinel-2 MSI data.