This study aims to identify the subsurface structure of the geothermal system in the Pacitan area, East Java, using gravity methods based on GGMPlus data. The analysis procedures include terrain correction to obtain the Complete Bouguer Anomaly (CBA), separation of regional and local anomalies through upward continuation, as well as two-dimensional (2D) modeling and Second Vertical Derivative (SVD) analysis to identify fault structure. The CBA values in the study area range from 107 to 115.5 mGal, where high anomalies correlate with volcanic breccia rocks and highlands, while low anomalies indicate the presence of sandstone in lowland areas. Geothermal manifestations in the form of hot springs exhibit gravity anomalies of approximately 110.5–111.5 mGal and are situated in the north of the fault in a northwest–southeast direction. The 2D modeling results reveal the presence of reservoir rocks consisting of sandstone and clay with densities ranging from 1.79 to 2.20 g/cm³, and lava rocks with densities of 2.80–2.90 g/cm³ acting as cap rocks. SVD analysis indicates the existence of a normal fault directly associated with geothermal fluid pathways. This study demonstrates that the GGMPlus gravity method is effective for preliminary exploration of geothermal systems in potential areas such as Pacitan.