Geochemical Perspectives Intercept Therapeutic Properties of Saturnia Thermal Waters: A Retrospective Cohort Study

1. Introduction

The ancient practice of balneotherapy, etymologically derived from the Latin balneum (bath) and the Greek therapeia (healing), has undergone a notable renaissance within contemporary medicine. This resurgence represents a compelling intersection of time-honored empirical knowledge and modern molecular science. Balneotherapy exerts its therapeutic effects via two synergistic pathways: first, the inherent physical properties of water, including hydrostatic pressure and thermal modulation; and second, the unique chemical compositions of individual thermal springs, particularly those enriched with hydrogen sulfide (H_₂S), which has emerged as a key molecular mediator in human physiology [1].

Recent molecular studies have shed light on the immunomodulatory potential of sulfur-rich thermal waters. These waters appear to fine-tune inflammatory pathways, bolster innate immune responses, exhibit antimicrobial activity, and reinforce epidermal barrier integrity through multiple, interdependent mechanisms [2]. Among such natural phenomena, the Saturnia thermal springs in Tuscany are especially noteworthy due to their singular hydrogeochemical profile. Emerging from a subterranean course through ancient limestone formations, the Saturnia waters carry a stable and distinct mineral composition that sets them apart from conventional thermal sources. Their calcium-sulfate dominance, precisely balanced H_₂S concentrations, and trace levels of boron, lithium, and strontium contribute to a uniquely bioactive environment. The waters ‘surface consistently at 37.5°C—a physiologically favorable temperature—and their composition has been shaped over millennia by geological filtration and equilibrium processes [3]. The mineral profile of Saturnia waters is characterized by calcium concentrations of 522–600 mg/L, sulfate levels between 1470–1480 mg/L, and H_₂S concentrations ranging from 7.56 to 11.10 mg/L. This biochemical constellation supports a therapeutic landscape that has demonstrated noteworthy clinical outcomes [4]. Recent advances in molecular biology have elucidated the mechanisms through which H_₂S exerts its therapeutic functions. Recognized as a gas transmitter, H_₂S influences critical physiological processes including inflammation, oxidative stress regulation, and tissue regeneration [5]. Its role in immune modulation is particularly well-documented, with evidence pointing to its capacity to modulate T-cell differentiation and cytokine expression [6]. Despite the increasing body of literature supporting the clinical utility of Saturnia waters and the mechanistic insights surrounding H_₂S, significant questions remain. Notably, the translation of specific chemical constituents into reproducible clinical outcomes, the influence of exposure duration on therapeutic efficacy, and the relative contribution of individual mineral elements remain insufficiently characterized [6]. The emerging evidence of H_₂S’s signaling and immunological roles offers a compelling rationale for investigating the relationships between defined geochemical parameters and measurable health benefits. A better understanding of these associations could facilitate the optimization of therapeutic protocols tailored to specific patient needs [7]. To address this knowledge gap, our retrospective study investigated whether the unique geochemical profile of Saturnia waters correlates with improvements in patient-reported quality-of-life (QoL) metrics. By comparing short-term versus prolonged exposure protocols, we aimed not only to validate the efficacy of these waters but also to elucidate the underlying mechanisms through which they exert their benefits. Specifically, our study addressed three key research questions:

• What is the correlation between individual geochemical parameters—especially H₂S concentration—and measurable therapeutic outcomes?
• Is there a dose–response relationship between exposure duration and clinical efficacy?
• How do specific mineral constituents differentially contribute to improvements in defined health domains?

2. Material and Methods

2.1. Study Design

This retrospective analysis examined clinical records from the Saturnia Thermal Springs Centre database covering the period from January 2016 to March 2024. The investigation utilized pre-existing clinical data reflecting two naturally occurring treatment modalities: acute exposure (2–3 hours, n = 114) and sustained exposure (3–14 days, n = 76). Data extraction and analysis were conducted between June and August 2024 in accordance with established standards for retrospective clinical research.

2.2. Ethics and Data Protection

The retrospective nature of this analysis, conducted within the framework of standard clinical practice at our institution, precluded the requirement for formal ethics committee approval. Nevertheless, the investigation rigorously adhered to the Declaration of Helsinki principles while implementing comprehensive data protection measures. The methodological framework incorporated a multi-tiered approach to data security and participant confidentiality. Patient records underwent systematic anonymization and de-identification prior to analysis, utilizing standardized protocols developed for retrospective clinical investigations. Data extraction procedures followed empirically validated methodologies to ensure consistency across the longitudinal dataset. The secure research database architecture implemented restricted access protocols, limiting data availability exclusively to authorized research personnel. All analyses were conducted exclusively on anonymized datasets, with results presented in aggregate form to maintain the highest standards of participant confidentiality while preserving scientific validity.

2.3. Patient Selection and Data Stratification

Records from the Saturnia Thermal Springs Centre database were screened to identify eligible cases between January 2016 and March 2024. Inclusion criteria required patients to be ≥18 years of age, to have undergone balneotherapy, and to have completed standardized QoL assessments. Records were eligible only if they included comprehensive medical documentation, complete QoL questionnaires (in Italian or English), and clear records of treatment duration and frequency. The Saturnia center follows standardized clinical screening protocols that exclude individuals with severe cardiovascular disease, acute infections, active malignancies, or pregnancy. Therefore, the dataset inherently reflects a medically screened population who completed prescribed thermal therapy. Records were categorized into two cohorts based on treatment duration: an acute exposure group (single session lasting 2–3 hours) and a sustained exposure group (multiple sessions over 3–14 days). This stratification mirrored typical usage patterns at the facility, allowing examination of both immediate and cumulative therapeutic effects (Figure 1).

2.4. Data Collection and Sample Characteristics

Using a systematic retrospective sampling strategy, clinical records were extracted from the Saturnia Thermal Springs Centre database (2016–2024). The final sample comprised 190 patients selected based on completeness and documentation quality. Post-hoc power analysis using G*Power 3.1 confirmed adequacy for multiple regression analyses (f² = 0.15, α = 0.05, 1-β = 0.80, predictors = 10). The sample (N = 190) was sufficient to detect clinically meaningful correlations while supporting subgroup comparisons. Of these, 114 (60%) represented acute exposure and 76 (40%) sustained exposure. This distribution is consistent with regional patterns of thermal facility utilization and ensures balanced cohorts for comparative evaluation.

2.5. Data Extraction and Management

Our analysis utilized standardized assessment data routinely collected at the Saturnia Thermal Springs Centre as part of standard clinical practice. The facility’s established protocol includes systematic QoL assessments at two time points: baseline (T0) and post-treatment evaluation (T1). For acute exposure cases (2-3 hours), T1 assessments were conducted immediately following treatment completion, while for sustained exposure cases (3-14 days), final evaluations were performed on the last treatment day. The facility’s standard assessment protocol involves qualified healthcare professionals administering validated questionnaires in a dedicated clinical area, ensuring consistent evaluation conditions across all cases. QoL assessments are integrated into routine clinical documentation, with standardized procedures for data collection and storage.

2.6. Outcome Assessment

2.6.1. Primary Outcomes

QoL was assessed using validated Italian versions of the SF-36 and WHOQOL-BREF instruments at T0 and T1. The SF-36, with strong internal consistency (α = 0.82–0.94) and test-retest reliability (ICC = 0.85–0.92), provided Physical and Mental Component Summary scores. A 5-point change constituted the Minimal Clinically Important Difference (MCID), with scores normalized on a 0–100 scale. The WHOQOL-BREF demonstrated reliability across its physical (α = 0.89, ICC = 0.88), psychological (α = 0.87, ICC = 0.83), social (α = 0.85, ICC = 0.76), and environmental (α = 0.93, ICC = 0.87) domains. MCID was defined as a 10% improvement from baseline.

2.6.2. Secondary Outcomes

Pain assessment utilized the Visual Analogue Scale (VAS), employing a 100mm horizontal line with anchors denoting “no pain” to “worst possible pain.” This instrument demonstrated excellent test-retest reliability (ICC = 0.90), with an MCID established at 15mm reduction. Dermatological quality of life assessment employed the dermatological quality of life (DLQI), comprising a total score range of 0-30, with demonstrated internal consistency (α = 0.88-0.92) and test-retest reliability (ICC = 0.84). The MCID was established at a 4-point reduction. Global therapeutic response was evaluated using the Patient Global Impression of Change (PGIC) scale at T1, utilizing a 7-point ordinal scale referenced to pre-treatment status. Treatment success was defined as scores ≤2, reflecting significant improvement from baseline.

2.6.3. Psychometric Properties of Assessment Tools

The SF-36 assessed eight domains: Physical Functioning (10 items, α = 0.93), Role-Physical (4 items, α = 0.89), Bodily Pain (2 items, α = 0.90), General Health (5 items, α = 0.82), Vitality (4 items, α = 0.86), Social Functioning (2 items, α = 0.85), Role-Emotional (3 items, α = 0.88), and Mental Health (5 items, α = 0.94).

The DLQI provided clinically validated cut-offs: 0–1 (no effect), 2–5 (small effect), 6–10 (moderate effect), 11–20 (very large effect), and 21–30 (extremely large effect).

Table 1. SF-36 Health Survey Structure and Psychometric Properties.

Domain	Items (n)	Score Range	Reliability (α)	Description
Physical Functioning	10	0-100	0.93	Extent of physical limitations in daily activities
Role-Physical	4	0-100	0.89	Impact of physical health on work/daily role performance
Bodily Pain	2	0-100	0.90	Pain intensity and interference with normal activities
General Health	5	0-100	0.82	Overall perception of health status and well-being
Vitality	4	0-100	0.86	Energy levels, fatigue, and vigor
Social Functioning	2	0-100	0.85	Impact of health on social activities and relationships
Role-Emotional	3	0-100	0.88	Effect of emotional health on work/daily activities
Mental Health	5	0-100	0.94	Psychological distress and well-being assessment

Note: All domains normalized to 0-100 scale; higher scores indicate better health status.

Table 1. SF-36 Health Survey Structure and Psychometric Properties.

Domain	Items (n)	Score Range	Reliability (α)	Description
Physical Functioning	10	0-100	0.93	Extent of physical limitations in daily activities
Role-Physical	4	0-100	0.89	Impact of physical health on work/daily role performance
Bodily Pain	2	0-100	0.90	Pain intensity and interference with normal activities
General Health	5	0-100	0.82	Overall perception of health status and well-being
Vitality	4	0-100	0.86	Energy levels, fatigue, and vigor
Social Functioning	2	0-100	0.85	Impact of health on social activities and relationships
Role-Emotional	3	0-100	0.88	Effect of emotional health on work/daily activities
Mental Health	5	0-100	0.94	Psychological distress and well-being assessment

Note: All domains normalized to 0-100 scale; higher scores indicate better health status.

Table 2. Dermatological quality of life (DLQI) Comprehensive Structure and Scoring.

Domain	Items (n)	Score Range	Key Questions	Clinical Interpretation
Symptoms and feelings	2	0-6	Itching, self-consciousness	0-1: No effect
Daily activities	2	0-6	Shopping, home management	2-5: Small effect
Leisure	2	0-6	Social, recreational activities	6-10: Moderate effect
Work and school	1	0-3	Work/study impact	11-20: Very large effect
Personal relationships	2	0-6	Intimate relationships, social interactions	21-30: Extremely large effect
Treatment	1	0-3	Treatment burden

Note: Each item scored 0-3 (0=not at all, 3=very much); total score interpretation based on validated thresholds.

Table 2. Dermatological quality of life (DLQI) Comprehensive Structure and Scoring.

Domain	Items (n)	Score Range	Key Questions	Clinical Interpretation
Symptoms and feelings	2	0-6	Itching, self-consciousness	0-1: No effect
Daily activities	2	0-6	Shopping, home management	2-5: Small effect
Leisure	2	0-6	Social, recreational activities	6-10: Moderate effect
Work and school	1	0-3	Work/study impact	11-20: Very large effect
Personal relationships	2	0-6	Intimate relationships, social interactions	21-30: Extremely large effect
Treatment	1	0-3	Treatment burden

Note: Each item scored 0-3 (0=not at all, 3=very much); total score interpretation based on validated thresholds.

Table 3. Patient Global Impression of Change (PGIC) Scale with Clinical Interpretations and Response Criteria.

Score	Category	Clinical Interpretation	Treatment Response
1	Very much improved	Complete or near-complete resolution	Success
2	Much improved	Substantial improvement with minimal symptoms	Success
3	Minimally improved	Slight but noticeable improvement	Partial Response
4	No change	No detectable change from baseline	Non-response
5	Minimally worse	Slight deterioration	Negative Response
6	Much worse	Substantial deterioration	Negative Response
7	Very much worse	Severe deterioration	Negative Response

Note: Assessment period: comparing current status to pre-treatment baseline; Success defined as scores ≤2.

Table 3. Patient Global Impression of Change (PGIC) Scale with Clinical Interpretations and Response Criteria.

Score	Category	Clinical Interpretation	Treatment Response
1	Very much improved	Complete or near-complete resolution	Success
2	Much improved	Substantial improvement with minimal symptoms	Success
3	Minimally improved	Slight but noticeable improvement	Partial Response
4	No change	No detectable change from baseline	Non-response
5	Minimally worse	Slight deterioration	Negative Response
6	Much worse	Substantial deterioration	Negative Response
7	Very much worse	Severe deterioration	Negative Response

Note: Assessment period: comparing current status to pre-treatment baseline; Success defined as scores ≤2.

2.7. Standard Treatment Protocols and Bias Consideration

The therapeutic protocols implemented at the Saturnia Thermal Springs Centre reflect decades of refined clinical experience, offering a naturalistic framework for evaluating treatment outcomes. Our analysis focused on two distinct therapeutic patterns that emerged from standard clinical practice: acute exposure protocols, characterized by single-session immersions lasting 2-3 hours, and sustained exposure protocols, involving daily sessions over periods ranging from 3 to 14 days. This natural stratification provides valuable insight into real-world treatment outcomes while maintaining rigorous safety standards.

2.8. Geochemical Data

The geochemical analysis of Saturnia thermal waters was based on comprehensive historical data [14]. This approach was justified by the documented long-term stability of the thermal system’s composition [5]. The waters’ composition is characterized by a distinctive calcium-sulfate dominance, with significant H₂S (Table 4). Of therapeutic relevance are the elevated concentrations of SO₄²⁻ (1475.00 ± 5.00 mg/L) and Ca²⁺ (561.00 ± 39.00 mg/L), alongside therapeutically significant levels of H₂S (14.00 ± 1.20 mg/L). The mineral profile is further enriched by notable concentrations of trace elements, including boron (19315.0 ± 872.0 μg/L), lithium (633.1 ± 28.4 μg/L), and strontium (11384.0 ± 334.0 μg/L), each potentially contributing to the observed therapeutic effects. Physical parameters remained remarkably stable throughout the study period, with temperature maintained at 37.5 ± 0.2°C and pH at 6.85 ± 0.05, creating optimal conditions for therapeutic applications.

2.9. Statistical Analysis

Initial data preprocessing encompassed a thorough quality assessment, including evaluation of outliers through standardized scores and box plot analysis, coupled with normality testing via Q-Q plots and the Shapiro-Wilk test. QoL scores underwent standardization following validated guidelines, with appropriate transformations applied to non-normally distributed variables to ensure robust statistical analysis. The management of missing data followed a systematic approach, beginning with Little’s MCAR test for pattern evaluation. We implemented a tiered imputation strategy, employing single imputation through expectation-maximization for variables with minimal missingness (<5%) and multiple imputation generating twenty datasets for variables with moderate missingness (5-20%). Records exceeding 20% missing data were excluded to maintain analytical integrity. Our primary analytical framework examined the relationships between standardized geochemical parameters and QoL outcomes through multiple complementary approaches. Descriptive statistics provided initial characterization of demographic patterns and QoL measures, while changes in scores underwent analysis through paired t-tests or Wilcoxon signed-rank tests, depending on distribution characteristics. Comparisons between acute and sustained exposure cohorts employed ANCOVA methodology, with careful adjustment for baseline scores to account for initial differences between groups. To investigate the predictive power of geochemical parameters while controlling for demographic factors and exposure duration, we developed hierarchical multiple regression models. These underwent comprehensive diagnostic evaluation, including assessments of linearity, residual normality, homoscedasticity, and error independence. The models incorporated carefully selected interaction terms based on theoretical considerations and statistical significance, with final model selection guided by improvements in adjusted R² values and the maintenance of theoretical coherence. To address multiple comparison concerns, we implemented the Benjamini-Hochberg procedure with a pre-specified false discovery rate of 0.05, applying this correction to all primary statistical tests conducted across QoL domains and geochemical parameters. This approach resulted in 18 of 24 tests maintaining significance after correction, supporting the robustness of our primary findings. Effect sizes were calculated and reported with standardized 95% confidence intervals throughout, utilizing Cohen’s d with bias correction for between-group comparisons and partial eta-squared for repeated measures analyses. Model validation followed a rigorous protocol, incorporating assessments of multicollinearity through variance inflation factors, residual normality via the Shapiro-Wilk test, and homoscedasticity through the Breusch-Pagan test. All analyses were conducted using R version 4.1.2 (R Core Team, 2021), with statistical significance consistently set at α = 0.05 (two-tailed), providing a comprehensive framework for evaluating the relationship between thermal water composition and therapeutic outcomes.

3. Results

3.1. Characteristics of the Study Population

Our retrospective analysis encompassed 190 cases, naturally stratified into two therapeutic modalities: acute exposure (n=114) and sustained exposure (n=76). The study population demonstrated a mean age of 52.3 years [95% CI: 50.4, 54.2] with a moderate female predominance (58.0%). Baseline characteristics between exposure groups exhibited comparable distributions, validating subsequent comparative analyses (Table 5).

3.2. Quality of Life Outcomes

Analysis of longitudinal data revealed significant improvements across QoL domains following thermal water exposure. The WHOQOL-BREF physical health domain exhibited the most pronounced enhancement (Δ = 18.4 points [95% CI: 15.2, 21.6], d = 0.85 [95% CI: 0.72, 0.98], p < .001). This substantial effect size was corroborated by parallel improvements in SF-36 Physical Component Summary scores (Δ = 11.7 points [95% CI: 9.8, 13.6], d = 0.67 [95% CI: 0.54, 0.80], p < .001). Duration-dependent analysis revealed distinct therapeutic patterns. The sustained exposure cohort consistently demonstrated superior outcomes, particularly in physical health domains (between-group difference = 6.4 points [95% CI: 3.8, 9.0], η² = 0.15 [95% CI: 0.08, 0.22]). This dose-response relationship manifested across multiple QoL dimensions, though with varying effect magnitudes (Table 6).

3.3. Geochemical Parameters and Clinical Response

H_₂S concentration emerged as the strongest predictor of therapeutic response (β = 0.28 [95% CI: 0.15, 0.41], p < .001), explaining 45.8% of outcome variance in the adjusted model. Multiple regression analysis, controlling for demographic factors, demonstrated significant predictive power for QoL improvements (adjusted R² = 0.435 [95% CI: 0.338, 0.532], F (5,184) = 31.12, p < .001) (Table 7). Tables S1 and S2 in Supplementary Materials presents the results of the regression analysis examining predictors of QoL improvement.

3.4. Correlation Analysis and Effect Modifiers

Examination of the relationship between specific geochemical parameters and therapeutic outcomes revealed a complex pattern of associations. H_₂S concentration demonstrated the strongest correlation with physical well-being (r = 0.45 [95% CI: 0.36, 0.54], p < .001), followed by SO₄²⁻ levels (r = 0.38 [95% CI: 0.29, 0.47], p < .001). These associations remained robust after adjustment for potential confounders through partial correlation analysis.

The trace element profile exhibited distinct patterns of therapeutic association. Boron demonstrated consistent moderate correlations across both physical (r = 0.31 [95% CI: 0.21, 0.41], p < .01) and psychological domains (r = 0.33 [95% CI: 0.23, 0.43], p < .01), suggesting a potential role in integrated therapeutic response.

3.5. Exposure-Dependent Response Patterns

Multiple regression analysis revealed significant interaction effects between exposure duration and geochemical parameters. The sustained exposure cohort demonstrated enhanced therapeutic responses particularly evident in:

Physical Component Summary: Δ = 11.7 [95% CI: 9.8, 13.6], d = 0.67 [95% CI: 0.54, 0.80] Mental Component Summary: Δ = 10.8 [95% CI: 8.9, 12.7], d = 0.62 [95% CI: 0.49, 0.75]

The magnitude of these effects remained significant after Benjamini-Hochberg correction for multiple comparisons (all p < .001).

3.6. Subgroup Analysis and Therapeutic Specificity

In the dermatological subgroup (n=73), calcium concentration emerged as a significant predictor (β = 0.31 [95% CI: 0.05, 0.57], p = .022), suggesting mineral-specific therapeutic mechanisms. This relationship-maintained significance after controlling for treatment duration and baseline severity scores (adjusted R² = 0.392 [95% CI: 0.307, 0.477]).

3.7. Temporal Analysis of Therapeutic Response

Time-series analysis of QoL metrics revealed distinct response patterns between acute and sustained exposure cohorts. The sustained exposure group demonstrated progressive improvement across all domains, with effect sizes increasing over the treatment period. The acute exposure group showed immediate improvements that, while significant, were of lesser magnitude:

WHOQOL-BREF Physical Domain:

• Acute: Δ = 15.2 [95% CI: 12.4, 18.0], d = 0.71
• Sustained: Δ = 21.6 [95% CI: 18.4, 24.8], d = 0.94

These temporal patterns suggest cumulative therapeutic effects that vary by exposure duration and specific outcome domain (Table 8 and Table 9 and Table S3 in Supplemetary Materials).

3.8. Multivariate Analysis

3.8.1. Hierarchical Regression Models

Hierarchical multiple regression analysis revealed complex interactions between geochemical parameters and therapeutic outcomes. The final model, incorporating both main effects and interaction terms, demonstrated substantial explanatory power (R² = 0.458 [95% CI: 0.384, 0.532], adjusted R² = 0.443, F (5,184) = 31.12, p < .001). H_₂S concentration emerged as the primary predictor (β = 0.28 [95% CI: 0.15, 0.41], p < .001), followed by exposure duration (β = 0.25 [95% CI: 0.12, 0.38], p < .001) (Table S4 in Supplementary Materials.

3.8.2. Component Analysis and Variable Relationships

Principal component analysis of the geochemical parameters identified three distinct therapeutic clusters explaining 78.6% of total variance (Table 10):

• Sulfur-dominant effects (eigenvalue = 3.24, variance explained = 41.8%)
• Mineral-trace element interactions (eigenvalue = 2.17, variance explained = 24.3%)
• Temporal exposure patterns (eigenvalue = 1.86, variance explained = 12.5%)

3.8.3. Interaction Effects and Effect Modifiers

Significant interaction effects emerged between H₂S concentration and exposure duration (β = 0.19 [95% CI: 0.06, 0.32], p = .004), suggesting synergistic therapeutic mechanisms. Age demonstrated modest negative moderation of treatment response (β = -0.17 [95% CI: -0.30, -0.04], p = .015), while gender showed no significant effect modification (β = -0.08 [95% CI: -0.21, 0.05], p = .232) (Table S5 and Table S6, Supplementary Materials).

3.8.4. Model Diagnostics and Validation

Comprehensive diagnostic procedures confirmed model robustness:

• Variance Inflation Factors remained below 2.5, indicating absence of problematic multicollinearity
• Residual analysis demonstrated normal distribution (Shapiro-Wilk, p = .412)
• Homoscedasticity was confirmed through Breusch-Pagan testing (χ² = 3.24, p = .198)
• Cook’s distance values identified no influential outliers (maximum D = 0.034)

The model’s predictive validity was assessed through k-fold cross-validation (k = 10), yielding stable performance metrics (mean RMSE = 0.142 [95% CI: 0.128, 0.156]).

3.9. Sensitivity Analyses

3.9.1. Parameter Stability Analysis

Sensitivity analyses assessed the stability of primary findings under varying analytical conditions. Sequential removal of influential cases (Cook’s distance > 0.029) demonstrated robust parameter estimates, with H₂S coefficients maintaining significance (βrange = 0.26-0.31, all p < .001). The relationship between exposure duration and therapeutic response remained stable across multiple analytical frameworks (η²range = 0.14-0.17).

3.9.2. Threshold Effects and Non-Linear Relationships

Examination of potential non-linear relationships through restricted cubic splines revealed threshold effects in H₂S therapeutic response. The inflection point occurred at 8.45 mg/L [95% CI: 7.89, 9.01], above which therapeutic effects demonstrated logarithmic progression (F = 24.31, p < .001). This finding suggests optimal therapeutic ranges for geochemical parameters (Tables S7–S9 in Supplementary Materials).

3.9.3. Temporal Stability and Seasonal Effects

Analysis of temporal stability across the eight-year study period revealed consistent effect sizes despite seasonal variations in thermal water parameters. The coefficient of variation for primary outcomes remained stable (CVrange = 0.11-0.14), indicating robust therapeutic effects independent of temporal fluctuations (Table S10 and Table S11 in Supplementary Materials).

3.9.4. Subgroup Consistency

Stratified analyses across demographic and clinical subgroups demonstrated consistent effect patterns:

Age-stratified analysis (quartiles):

• Q1 (≤45 years): β = 0.29 [95% CI: 0.15, 0.43]
• Q2 (46-52 years): β = 0.27 [95% CI: 0.13, 0.41]
• Q3 (53-59 years): β = 0.28 [95% CI: 0.14, 0.42]
• Q4 (≥60 years): β = 0.26 [95% CI: 0.12, 0.40]

The consistency of effect sizes across subgroups supports the generalizability of primary findings while identifying potential effect modifiers (Table S12, Supplementary Materials).

3.9.5. Methodological Robustness

Alternative analytical approaches, including weighted regression analyses and robust estimation methods, corroborated primary findings:

• Weighted least squares: β = 0.27 [95% CI: 0.14, 0.40]
• Robust regression (M-estimation): β = 0.29 [95% CI: 0.16, 0.42]
• Mixed-effects modeling: β = 0.28 [95% CI: 0.15, 0.41]

The convergence of results across methodological frameworks supports the robustness of observed therapeutic relationships (Table S13 and Table S14, Supplementary Materials).

4. Discussion

Our comprehensive analysis of Saturnia thermal waters reveals a significant association between their distinctive geochemical profile and measurable improvements in QoL parameters. The observed correlations between specific mineral constituents, particularly H₂S, and enhanced well-being metrics suggest underlying mechanisms of considerable complexity. These findings enrich the current understanding of the mineral composition–therapeutic outcome relationship by offering empirical evidence that extends beyond traditional observational research [1,2,3,4]. The biological activity of sulfurous waters, especially via H₂S as a molecular mediator, highlights the intersection of geochemistry and human physiology.

While the therapeutic use of thermal waters is well-documented historically, our study offers quantitative evidence linking molecular pathways to clinical outcomes. The therapeutic environment at Saturnia is characterized by exceptional geochemical stability, validated through long-term monitoring. Waters consistently emerge at 37.5°C (±0.5°C), establishing physiologically favorable conditions. The calcium-sulfate dominant profile (Ca: 560 ±20 mg/L; SO4: 1475 ±5 mg/L) is complemented by therapeutically relevant H₂S concentrations (9.33 ±1.77 mg/L). Additional trace elements, including boron (4.5 ±0.3 mg/L), lithium (0.18 ±0.02 mg/L), and strontium (8.8 ±0.4 mg/L), further contribute to the therapeutic milieu. The geochemical stability that characterizes the Saturnia thermal system—supported by ISO 17025-certified monthly analyses and coefficients of variation under 5%—provides a reliable basis for clinical application.

All treatments were conducted under standardized safety protocols, beginning with comprehensive medical screening and continuing through regular monitoring of physiological responses. The acute exposure protocol emphasizes careful observation of immediate therapeutic responses, while the sustained exposure protocol incorporates progressive adaptation periods and regular medical assessments to optimize therapeutic benefits while ensuring patient safety. This combination of stable geochemical parameters and standardized therapeutic protocols creates a well-controlled environment for evaluating treatment outcomes, while still reflecting the real-world conditions of clinical practice.

The natural variation in treatment patterns provides valuable insights into the relationship between exposure duration and therapeutic efficacy, while the consistent geochemical profile allows for reliable assessment of specific mineral-dependent effects. Monitoring of physical-chemical parameters, such as pH (7.0 ±0.2) and electrical conductivity (2450 μS/cm ±50 μS/cm), providing additional context for understanding the therapeutic environment. These measurements, combined with regular chemical analyses, ensure that all treatments occur within optimal therapeutic parameters while maintaining the highest standards of safety and efficacy. The dose-dependent nature of these effects, particularly evident in our extended exposure cohort, demonstrates biological responses that are precisely regulated by varying concentrations of active compounds. Recent molecular investigations by Andrés et al. [15] characterized how H₂S transforms into colloidal sulfur within the cutaneous matrix, thereby initiating specific cellular responses. Our findings extend this research, demonstrating significant correlations between sulfide concentrations and QoL improvements that correspond to these molecular events.

The relationship between exposure duration and therapeutic benefit suggests a temporally sensitive mechanism, potentially reflecting cumulative physiological adaptations [16]. H₂S appears to interact with MAP kinase pathways [17,18], a hypothesis supported by Mirandola et al. [19], and contributes to keratinocyte modulation and vascular responses [19,20]. These biochemical cascades provide plausible explanations for the observed improvements in both physical and psychological domains.

The magnitude of improvement in physical health domains (Cohen’s d = 0.85) observed in our study differs notably from prior reports on sulfurous thermal waters, highlighting the distinctiveness of the Saturnia composition. This robust effect size warrants careful contextualization within established clinical significance thresholds. While the aggregated Cohen’s d indicates a substantial therapeutic effect, the observed variability across different QoL dimensions necessitates nuanced interpretation. Notably, the WHOQOL-BREF physical domain scores improved by 18.4 points, far exceeding the MCID of 10%, thereby suggesting a clinically meaningful benefit.

Differences in response patterns between acute and sustained exposure cohorts further underscore the potential dose-response relationship. A mean difference of 6.4 points (95% CI: 3.8, 9.0) supports the hypothesis of duration-dependent therapeutic efficacy, although this finding should be interpreted with caution due to the potential for selection bias. While studies by Fioravanti et al. and Kovács et al. [21,22] reported comparable directional trends, discrepancies in effect sizes may be attributed to the unique mineral composition of Saturnia waters.

The role of trace elements in shaping therapeutic outcomes is particularly noteworthy. Correlations between boron concentrations and psychological well-being indicators introduce a novel dimension to thermal water research, potentially implicating trace element biochemistry in neuromodulatory processes. Similarly, the association between lithium levels and mental health domains presents intriguing avenues for future research [23].

Improvements in dermatological parameters also appear to be influenced by specific geochemical components—notably calcium and sulfur. Li et al. [24] describe the formation of pentathionic acid in the epidermis, offering a plausible explanation for the antimicrobial and barrier-enhancing properties of sulfurous waters. When coupled with calcium’s role in keratinocyte differentiation, a mechanistic framework begins to emerge for interpreting the dermatological benefits observed.

The dose-response patterns revealed in this study contribute meaningfully to the broader discourse on optimal thermal water exposure protocols. Although Karagülle and Karagülle [25] previously emphasized the benefits of extended treatment durations, our comparative analysis between acute and sustained exposure offers empirical support for cumulative therapeutic effects. Nevertheless, the underlying mechanisms remain to be elucidated through prospective molecular studies.

Particularly compelling are the relationships identified between specific trace element concentrations and distinct QoL domains. The elevated boron and lithium levels found in Saturnia waters, compared to other thermal sources, may underlie the psychological well-being improvements documented here. These findings build upon those of Gálvez et al. [26], who posited neuromodulatory effects of mineral-rich waters, although causality remains to be demonstrated via controlled trials.

The observed associations likely reflect the synergistic action of multiple geochemical constituents, rather than isolated mineral effects. This holistic perspective aligns with emerging themes in balneotherapy research, which emphasize the therapeutic relevance of comprehensive mineral profiles [27].

However, several methodological considerations must temper the interpretation of our findings. The retrospective nature of the dataset, while rich in ecological validity, limits causal inference. Real-world stratification of exposure groups introduces potential confounders, and unmeasured psychosocial or lifestyle variables may also contribute to the observed outcomes. For example, the enhanced results associated with sustained exposure may partially reflect environmental or psychological benefits unrelated to mineral composition per se.

While validated QoL instruments underpin our assessments, reliance on self-reported data suggests a need for future investigations incorporating objective biochemical and physiological measures. The durability of the observed benefits remains unknown, underscoring the importance of longitudinal follow-up studies. These limitations highlight key avenues for future research.

Prospective studies employing controlled exposure protocols will be essential to explore the temporal dynamics underlying mineral-outcome associations. The development of biomarker panels specific to thermal water exposure could serve to validate subjective improvements with objective physiological metrics. Additionally, comparative analyses across diverse thermal water compositions may clarify the relative contributions of specific mineral elements.

The therapeutic potential of thermal water exposure to complement conventional interventions warrants systematic investigation. Our findings support the hypothesis that therapeutic outcomes are best understood as the product of integrated mineral actions. The unique geochemical profile of Saturnia waters raises critical questions regarding generalizability. Establishing minimal threshold concentrations—particularly for hydrogen sulfide (9.33 ±1.77 mg/L)—across different thermal sources will be essential to define evidence-based treatment standards. The exceptional geochemical stability observed may represent an optimal therapeutic condition, although acceptable variability ranges for clinical efficacy remain to be determined.

Developing standardized, evidence-informed protocols for thermal therapy will require integrated approaches that combine clinical, geochemical, and molecular data. Our study contributes to this objective by offering detailed insight into the mineral-outcome relationships that define the therapeutic landscape of Saturnia thermal waters.

4.1. Limitations and Future Directions

While our findings underscore robust associations between geochemical composition and clinical outcomes, several methodological constraints merit discussion. The retrospective design limits causal interpretation, and natural group stratification introduces potential confounding. The retrospective design, while provides substantial ecological validity, inherently constrains causal inference. The natural stratification of exposure groups, though reflecting authentic therapeutic utilization patterns, introduces potential selection factors that warrant consideration in interpretation.

Absence of molecular or biochemical markers precludes direct mechanistic insight. Though our dataset offers strong ecological validity, self-report instruments—despite their validation—limit biological specificity.

These methodological considerations point toward several promising research directions. Prospective studies could employ controlled exposure protocols while incorporating biological sampling to examine the molecular basis of the mineral-outcome associations our findings suggest. The development of biomarker panels specific to thermal water exposure could provide objective correlates to subjective well-being measures. Additionally, Comparative studies across thermal sources may clarify the contribution of distinct mineral elements to therapeutic effects.

5. Conclusions

Our results advance the scientific framework of balneotherapy by demonstrating dose-response patterns and mineral-outcome associations specific to Saturnia thermal waters. The distinct geochemical composition, including high and stable H₂S concentrations, appears to underpin meaningful therapeutic benefits. Nevertheless, controlled trials are required to establish causality, define minimal effective concentrations, and assess the generalizability of these findings across diverse thermal environments. Integrated research approaches combining clinical, geochemical, and molecular perspectives are essential to develop evidence-based thermal therapy protocols.