adjust Number of Fractures to 65 due to better stats fit. Calculate stats in sensitivity analysis - WIP: TBD whether to include them in analysis doc

Author: softwareengineerprogrammer

docs/Fervo_Project_Red.md.jinja

@@ -172,7 +172,7 @@ The variance analysis (results displayed in legend captions) evaluates the predi
 Both models demonstrate high predictive fidelity, tracking steady-state flowing temperatures within 1.5°C of the empirical data.
 
 * **Overall Fit:** GEOPHIRES mathematically achieves a tighter overall fit, yielding a lower Root Mean Square Error (RMSE) and a higher coefficient of determination (R²).
-* **Systematic Bias:** The Fervo model exhibits slightly less systemic underestimation, with a cold bias of {{ fervo_bias_degc }}°C compared to the GEOPHIRES cold bias of {{ geophires_bias_degc }}°C.
+* **Systematic Bias:** The GEOPHIRES model exhibits slightly less systemic underestimation, with a cold bias of {{ geophires_bias_degc }}°C compared to the Fervo cold bias of {{ fervo_bias_degc }}°C.
 * **R² Context:** While the Fervo model yields a relatively low R² ({{ fervo_r2 }}), GEOPHIRES achieves a stronger R² of {{ geophires_r2 }}, indicating it more successfully captures the underlying physical trend of the data (the slight thermal drawdown) rather than simply averaging the noise. However, it is important to note that the absolute R² ceiling for both models is inherently suppressed by the dataset. Because the steady-state temperature profile is essentially a flat plateau, natural sensor variance and minor reservoir oscillations may account for a disproportionately large portion of the total variance, keeping the R² scores modest despite the favorable absolute error (RMSE).
 
 <div id="long-term-forecast-section"></div>

docs/_images/fervo_project_red-2026_production-temperature-data-vs-modeling-1.png

@@ -504,6 +504,7 @@ def _generate_fracture_sensitivity_graph(
     steady_state_start_years: float,
     sensitivity_graph_path: Path,
     show_excluded_measured_temperatures: bool = False,
+    calculate_stats: bool = True,
 ) -> None:
     _log.info('Running 8-year fracture sensitivity analysis...')
 
@@ -542,7 +543,13 @@ def _generate_fracture_sensitivity_graph(
     base_input_params: GeophiresInputParameters = get_project_red_input_params_and_result()[0]
     base_number_of_fractures = int(_get_input_parameters_dict(base_input_params)[number_of_fractures_param_name])
 
-    fracture_counts = [base_number_of_fractures, 57, 60, 66, 69]
+    fracture_counts = [
+        base_number_of_fractures,
+        base_number_of_fractures - 6,
+        base_number_of_fractures - 3,
+        base_number_of_fractures + 3,
+        base_number_of_fractures + 6,
+    ]
     client = GeophiresXClient()
 
     colors = {
@@ -556,6 +563,13 @@ def _generate_fracture_sensitivity_graph(
         base_number_of_fractures: '-.',
     }
 
+    if calculate_stats:
+        _log.info(
+            f'--- FRACTURE SENSITIVITY STATISTICAL ALIGNMENT (Steady-State > {steady_state_start_years} Years) ---'
+        )
+        y_true = df_included['Temperature_C'].values
+        ss_tot = float(np.sum((y_true - np.mean(y_true)) ** 2))
+
     for frac_count in fracture_counts:
         input_params: GeophiresInputParameters = ImmutableGeophiresInputParameters(
             from_file_path=base_input_params.as_file_path(),
@@ -574,6 +588,17 @@ def _generate_fracture_sensitivity_graph(
         geophires_x = [float(step) / float(time_steps_per_year) for step, _ in enumerate(profile)]
         geophires_y = [q.magnitude for q in profile]
 
+        if calculate_stats and not df_included.empty:
+            geo_interp = interp1d(geophires_x, geophires_y, kind='linear', fill_value='extrapolate')
+            y_geo = geo_interp(df_included['Time_Years'])
+
+            rmse_g = float(np.sqrt(((y_true - y_geo) ** 2).mean()))
+            bias_g = float((y_geo - y_true).mean())
+            ss_res_g = float(np.sum((y_true - y_geo) ** 2))
+            r2_g = 1.0 - (ss_res_g / ss_tot) if ss_tot != 0.0 else 0.0
+
+            _log.info(f'{frac_count} Fractures: RMSE={rmse_g:.2f}°C, R²={r2_g:.4f}, Bias={bias_g:.2f}°C')
+
         label_prefix = 'GEOPHIRES: ' if frac_count == base_number_of_fractures else ''
         label_suffix = ' (Baseline)' if frac_count == base_number_of_fractures else ''
 
@@ -598,7 +623,6 @@ def _generate_fracture_sensitivity_graph(
 
     def _savefig(version: int | str) -> None:
         fig.savefig(
-            # sensitivity_graph_path,
             sensitivity_graph_path.with_name(f'{sensitivity_graph_path.stem}-{version}').with_suffix(
                 sensitivity_graph_path.suffix
             ),

docs/_images/fervo_project_red-2026_production-temperature-data-vs-modeling-2.png

@@ -17,7 +17,7 @@ Reservoir Thermal Conductivity, 2.7
 Number of Segments, 1
 Gradient 1, 76.1
 
-Number of Fractures, 63, -- Fervo estimates between 75 and 100 fractures were created. This value is de-rated in the model to account for the physical reality of imperfect flow distribution and uneven utilization across the entire stimulated rock volume (Norbeck and Latimer, 2023).
+Number of Fractures, 65, -- Fervo estimates between 75 and 100 fractures were created. This value is de-rated in the model to account for the physical reality of imperfect flow distribution and uneven utilization across the entire stimulated rock volume (Norbeck and Latimer, 2023).
 Fracture Shape, 4, -- Rectangular geometry (Shape 4), representing standard transverse hydraulic fractures along a horizontal wellbore.
 Fracture Height, 300 foot, -- Estimated vertical propagation of the stimulated fracture network.
 Fracture Width, 365 foot, -- Set to match the distance between the injection and production wellbores, assuming a dipole flow field directly connecting the laterals.

docs/_images/fervo_project_red-2026_production-temperature-data-vs-modeling-fracture-sensitivity-1.png

@@ -45,4 +45,4 @@ def _vuq(v_u: dict[str, Any]) -> PlainQuantity:
 
         avg_production_temp_q = _vuq(r.result['RESERVOIR SIMULATION RESULTS']['Average Production Temperature'])
         self.assertGreater(avg_production_temp_q, _q(346, 'degF'))
-        self.assertLess(avg_production_temp_q, _q(356, 'degF'))
+        self.assertLess(avg_production_temp_q, _q(357, 'degF'))