Quantifying $pCO_2$ Bias from Mean-State Errors¶

1. Decomposition of $pCO_2$ Bias¶

First, we quantify how biases in the modeled mean state (DIC, alkalinity, and temperature) propagate into errors in simulated $pCO_2$ .

Using a first-order Taylor expansion, the total bias in $pCO_2$ ( $\Delta pCO_2^{total}$ ) can be approximated as the sum of the linear contributions from each state variable:

\Delta pCO_2^{total} \approx \left( \frac{\partial pCO_2}{\partial DIC} \right) \Delta DIC + \left( \frac{\partial pCO_2}{\partial Alk} \right) \Delta Alk + \left( \frac{\partial pCO_2}{\partial T} \right) \Delta T

(1)

Where:

$\Delta DIC, \Delta Alk, \Delta T$ are the mean-state biases (Model $-$ Observed).
$\frac{\partial pCO_2}{\partial X}$ are the partial derivatives (sensitivities) evaluated at the observed reference state.

We find that the alkalinity and DIC biases in the model largely cancel out, but any temperature bias in the model will propagate into $pCO_2$ space

import numpy as np
import PyCO2SYS as pyco2
import pandas as pd

# -------------------------------
# helper: compute pCO2
# -------------------------------
def compute_pCO2(alk_m3, dic_m3, sal, temp, pres):
    rho = 1026.0 
    out = pyco2.sys(
        par1=(alk_m3 / rho) * 1000,
        par2=(dic_m3 / rho) * 1000,
        par1_type=1, par2_type=2,
        salinity=sal, temperature=temp, pressure=pres,
    )
    return out["pCO2"]

# -------------------------------
# derivatives
# -------------------------------
def get_derivatives(state, eps=30, eps_t=0.1):
    alk, dic = state["alk"], state["dic"]
    s, t, p = state["sal"], state["temp"], state["pres"]
    
    # Chemistry First derivatives
    dpc_dDIC = (compute_pCO2(alk, dic + eps, s, t, p) - 
                compute_pCO2(alk, dic - eps, s, t, p)) / (2 * eps)
    
    dpc_dALK = (compute_pCO2(alk + eps, dic, s, t, p) - 
                compute_pCO2(alk - eps, dic, s, t, p)) / (2 * eps)
    
    # Temperature First derivative
    dpc_dT = (compute_pCO2(alk, dic, s, t + eps_t, p) - 
              compute_pCO2(alk, dic, s, t - eps_t, p)) / (2 * eps_t)
    
    return dpc_dDIC, dpc_dALK, dpc_dT

# -------------------------------
# main calculation
# -------------------------------
def dpCO2_bias_analysis(obs, model):
    # 1. Calculate Sensitivities at Observed State
    dDIC_sens, dALK_sens, dTemp_sens = get_derivatives(obs)

    # 2. Mean state biases
    delta_DIC = model["dic"] - obs["dic"]
    delta_ALK = model["alk"] - obs["alk"]
    delta_Temp = model["temp"] - obs["temp"]

    # 3. Individual components
    error_DIC = dDIC_sens * delta_DIC
    error_ALK = dALK_sens * delta_ALK
    error_Temp = dTemp_sens * delta_Temp
    
    # 4. Actual total for reference
    actual_total = compute_pCO2(model['alk'], model['dic'], model['sal'], model['temp'], model['pres']) - \
                   compute_pCO2(obs['alk'], obs['dic'], obs['sal'], obs['temp'], obs['pres'])

    # Summary table
    data = {
        "Component": [
            "DIC First-Order", 
            "ALK First-Order", 
            "Temp First-Order", 
            "Total Approximation", 
            "Actual Bias"
        ],
        "Sensitivity": [
            dDIC_sens, 
            dALK_sens, 
            dTemp_sens, 
            np.nan, 
            np.nan
        ],
        "Bias (ΔX)": [
            delta_DIC, 
            delta_ALK, 
            delta_Temp, 
            np.nan, 
            np.nan
        ],
        "pCO2 Error (uatm)": [
            error_DIC, 
            error_ALK, 
            error_Temp,
            error_DIC + error_ALK + error_Temp,
            actual_total
        ]
    }
    
    return pd.DataFrame(data)

# -------------------------------
# example usage
# -------------------------------
model = dict(alk=2320, dic=2110.0, sal=33.75, temp=11.5, pres=0)
obs = dict(alk=2380, dic=2160.0, sal=33.75, temp=10, pres=0)

df_results = dpCO2_bias_analysis(obs, model)
df_results

Full Multi-Parameter Analysis of $pCO_2$ Sensitivity¶

We now consider $pCO_2$ as a function of Dissolved Inorganic Carbon ( $DIC$ ), Total Alkalinity ( $Alk$ ), Temperature ( $T$ ), and Salinity ( $S$ ):

f = f(DIC, Alk, T, S)

(2)

The change in $pCO_2$ resulting from the alkalinity intervention $I$ is:

\delta f = f(DIC, Alk + I, T, S) - f(DIC, Alk, T, S)

(3)

We approximate $\delta f$ using the sensitivity of $pCO_2$ to alkalinity, $\theta$ :

\theta = \frac{\partial pCO_2}{\partial Alk}

(4)

\delta f \approx \theta \cdot I

(5)

To determine how the contribution of biases in all four background parameters ( $\Delta DIC, \Delta Alk, \Delta T, \Delta S$ ) affects the predicted change $\delta f$ , we calculate the bias in sensitivity $\Delta \theta$ and multiply by the intervention:

\Delta(\delta f) = \Delta \theta \cdot I = \left( \frac{\partial \theta}{\partial DIC} \Delta DIC + \frac{\partial \theta}{\partial Alk} \Delta Alk + \frac{\partial \theta}{\partial T} \Delta T + \frac{\partial \theta}{\partial S} \Delta S \right) \cdot I

(6)

import numpy as np
import pandas as pd
import PyCO2SYS as pyco2

# ---------------------------------------------------------
# 1. Setup states and intervention
# ---------------------------------------------------------
obs = dict(alk=2380.0, dic=2160.0, sal=33.75, temp=10.0, pres=0)
model = dict(alk=2320.0, dic=2110.0, sal=33.75, temp=11.5, pres=0)
I = 30.0  # Intervention magnitude (umol/kg)

# ---------------------------------------------------------
# helper: compute pCO2
# ---------------------------------------------------------
def compute_pCO2(alk, dic, sal, temp, pres):
    out = pyco2.sys(
        par1=alk,
        par2=dic,
        par1_type=1, 
        par2_type=2,
        salinity=sal, 
        temperature=temp, 
        pressure=pres,
    )
    return out["pCO2"]

def get_theta(alk, dic, temp, sal, pres=0):
    """Calculates theta = dpCO2/dAlk using a central difference."""
    d_alk = 0.01
    pCO2_high = compute_pCO2(alk + d_alk, dic, sal, temp, pres)
    pCO2_low  = compute_pCO2(alk - d_alk, dic, sal, temp, pres)
    return (pCO2_high - pCO2_low) / (2 * d_alk)

def get_isoQ(alk, dic, sal, temp, pres):
    out = pyco2.sys(
        par1=alk,
        par2=dic,
        par1_type=1,
        par2_type=2,
        salinity=sal,
        temperature=temp,
        pressure=pres,
    )
    return 1/out["isocapnic_quotient"]   # key name in PyCO2SYS
    
isocapnic_quotient = get_isoQ(
    obs['alk'], obs['dic'], obs['sal'], obs['temp'], obs['pres']
)
# ---------------------------------------------------------
# 2. Calculate Sensitivities (dTheta/dX) at OBS state
# ---------------------------------------------------------
eps = 1e-2 

dTheta_dAlk = (get_theta(obs['alk']+eps, obs['dic'], obs['temp'], obs['sal']) - 
               get_theta(obs['alk']-eps, obs['dic'], obs['temp'], obs['sal'])) / (2*eps)

dTheta_dDic = (get_theta(obs['alk'], obs['dic']+eps, obs['temp'], obs['sal']) - 
               get_theta(obs['alk'], obs['dic']-eps, obs['temp'], obs['sal'])) / (2*eps)

dTheta_dTemp = (get_theta(obs['alk'], obs['dic'], obs['temp']+eps, obs['sal']) - 
                get_theta(obs['alk'], obs['dic'], obs['temp']-eps, obs['sal'])) / (2*eps)

dTheta_dSal = (get_theta(obs['alk'], obs['dic'], obs['temp'], obs['sal']+eps) - 
               get_theta(obs['alk'], obs['dic'], obs['temp'], obs['sal']-eps)) / (2*eps)


# 3. Compute Biases and Expected Errors (Linearized)
# ---------------------------------------------------------
variables = ['alk', 'dic', 'temp', 'sal']
sensitivities = [dTheta_dAlk, dTheta_dDic, dTheta_dTemp, dTheta_dSal]

results = []
for var, sens in zip(variables, sensitivities):
    bias = model[var] - obs[var]
    error_pco2 = sens * bias * I
    results.append({
        "Variable": var,
        "Sensitivity (∂θ/∂X)": sens,
        "Total Bias (ΔX)": bias,
        "Linearized Error [μatm]": error_pco2,
        "Isocapnic Quotient": isocapnic_quotient
    })

df = pd.DataFrame(results)
sum_linear_error = df["Linearized Error [μatm]"].sum()

# ---------------------------------------------------------
# 4. Calculate ACTUAL Error (Direct Computation)
# ---------------------------------------------------------
# Delta f in OBS state
df_obs = (compute_pCO2(obs['alk'] + I, obs['dic'], obs['sal'], obs['temp'], obs['pres']) - 
          compute_pCO2(obs['alk'], obs['dic'], obs['sal'], obs['temp'], obs['pres']))

# Delta f in MODEL state
df_model = (compute_pCO2(model['alk'] + I, model['dic'], model['sal'], model['temp'], model['pres']) - 
            compute_pCO2(model['alk'], model['dic'], model['sal'], model['temp'], model['pres']))

# Change this line in your script:
actual_error = df_model - df_obs 

# And the attribution:
error_pco2 = sens * bias * I  # where bias = model - obs

# ---------------------------------------------------------
# 5. Output Results
# ---------------------------------------------------------
print(f"Intervention Intensity (I): {I} umol/kg\n")
print(df.to_string(index=False))
print("-" * 60)
print(f"Expected Error (Sum of terms):   {sum_linear_error:10.4f} μatm")
print(f"Actual Error (Direct f_obs-f_mod): {actual_error:10.4f} μatm")
print(f"Residual (Higher-order terms):    {actual_error - sum_linear_error:10.4f} μatm")

Intervention Intensity (I): 30.0 umol/kg

Variable  Sensitivity (∂θ/∂X)  Total Bias (ΔX)  Linearized Error [μatm]  Isocapnic Quotient
     alk             0.011500            -60.0               -20.700032            0.858088
     dic            -0.012470            -50.0                18.705461            0.858088
    temp            -0.065056              1.5                -2.927531            0.858088
     sal            -0.035555              0.0                -0.000000            0.858088
------------------------------------------------------------
Expected Error (Sum of terms):      -4.9221 μatm
Actual Error (Direct f_obs-f_mod):    -4.5408 μatm
Residual (Higher-order terms):        0.3813 μatm

import numpy as np
import pandas as pd
import PyCO2SYS as pyco2

# ---------------------------------------------------------
# 1. Setup states and intervention
# ---------------------------------------------------------
obs = dict(alk=2380.0, dic=2160.0, sal=33.75, temp=10.0, pres=0)
model = dict(alk=2320.0, dic=2110.0, sal=33.75, temp=11.5, pres=0)
I = 30.0  # Intervention magnitude (umol/kg)

# ---------------------------------------------------------
# Helpers: compute carbonate system properties
# ---------------------------------------------------------
def compute_pCO2(alk, dic, sal, temp, pres):
    out = pyco2.sys(
        par1=alk,
        par2=dic,
        par1_type=1, 
        par2_type=2,
        salinity=sal, 
        temperature=temp, 
        pressure=pres,
    )
    return out["pCO2"]

def get_theta(alk, dic, temp, sal, pres=0):
    """Calculates theta = dpCO2/dAlk using a central difference."""
    d_alk = 0.01
    pCO2_high = compute_pCO2(alk + d_alk, dic, sal, temp, pres)
    pCO2_low  = compute_pCO2(alk - d_alk, dic, sal, temp, pres)
    return (pCO2_high - pCO2_low) / (2 * d_alk)

def get_isoQ(alk, dic, temp, sal, pres=0):
    """Extracts the isocapnic quotient from PyCO2SYS."""
    out = pyco2.sys(
        par1=alk,
        par2=dic,
        par1_type=1,
        par2_type=2,
        salinity=sal,
        temperature=temp,
        pressure=pres,
    )
    # Using the standard PyCO2SYS key name
    return out["isocapnic_quotient"] 
    
# ---------------------------------------------------------
# 2. Calculate Sensitivities (dIsoQ/dX) at OBS state
# ---------------------------------------------------------
eps = 1e-2 

# Central difference derivatives for the Isocapnic Quotient (Q)
dQ_dAlk = (get_isoQ(obs['alk']+eps, obs['dic'], obs['temp'], obs['sal']) - 
           get_isoQ(obs['alk']-eps, obs['dic'], obs['temp'], obs['sal'])) / (2*eps)

dQ_dDic = (get_isoQ(obs['alk'], obs['dic']+eps, obs['temp'], obs['sal']) - 
           get_isoQ(obs['alk'], obs['dic']-eps, obs['temp'], obs['sal'])) / (2*eps)

dQ_dTemp = (get_isoQ(obs['alk'], obs['dic'], obs['temp']+eps, obs['sal']) - 
            get_isoQ(obs['alk'], obs['dic'], obs['temp']-eps, obs['sal'])) / (2*eps)

dQ_dSal = (get_isoQ(obs['alk'], obs['dic'], obs['temp'], obs['sal']+eps) - 
           get_isoQ(obs['alk'], obs['dic'], obs['temp'], obs['sal']-eps)) / (2*eps)

# ---------------------------------------------------------
# 3. Compute Biases and Expected Errors for IsoQ (Linearized)
# ---------------------------------------------------------
variables = ['alk', 'dic', 'temp', 'sal']
sensitivities_Q = [dQ_dAlk, dQ_dDic, dQ_dTemp, dQ_dSal]

results = []
for var, sens in zip(variables, sensitivities_Q):
    bias = model[var] - obs[var]
    
    # Linearized error framework applied to the Quotient parameter
    error_Q = sens * bias 
    
    results.append({
        "Variable": var,
        "Sensitivity (∂Q/∂X)": sens,
        "Total Bias (ΔX)": bias,
        "Linearized Error [Quotient Units]": error_Q
    })

df = pd.DataFrame(results)
sum_linear_error_Q = df["Linearized Error [Quotient Units]"].sum()

# ---------------------------------------------------------
# 4. Calculate ACTUAL Error in Quotient (Direct Computation)
# ---------------------------------------------------------
# Value of Q in the baseline OBS state
Q_obs = get_isoQ(obs['alk'], obs['dic'], obs['temp'], obs['sal'])

# Value of Q in the biased MODEL state
Q_model = get_isoQ(model['alk'], model['dic'], model['temp'], model['sal'])

# Direct baseline difference mismatch
actual_error_Q = Q_model - Q_obs 

# ---------------------------------------------------------
# 5. Output Results
# ---------------------------------------------------------
print(f"Intervention Intensity (I): {I} umol/kg\n")
print(df.to_string(index=False))
print("-" * 60)
print(f"Expected Error in Q (Sum of terms):   {sum_linear_error_Q:10.6f}")
print(f"Actual Error in Q (Direct mod-obs):    {actual_error_Q:10.6f}")
print(f"Residual (Higher-order terms):         {actual_error_Q - sum_linear_error_Q:10.6f}")

Intervention Intensity (I): 30.0 umol/kg

Variable  Sensitivity (∂Q/∂X)  Total Bias (ΔX)  Linearized Error [Quotient Units]
     alk             0.000596            -60.0                          -0.035771
     dic            -0.000649            -50.0                           0.032432
    temp             0.000300              1.5                           0.000450
     sal            -0.000168              0.0                          -0.000000
------------------------------------------------------------
Expected Error in Q (Sum of terms):    -0.002889
Actual Error in Q (Direct mod-obs):     -0.002949
Residual (Higher-order terms):          -0.000061

import numpy as np
import pandas as pd
import PyCO2SYS as pyco2

# ---------------------------------------------------------
# 1. Setup states and intervention
# ---------------------------------------------------------
obs = dict(alk=2380.0, dic=2160.0, sal=33.75, temp=10.0, pres=0)
model = dict(alk=2320.0, dic=2110.0, sal=33.75, temp=11.5, pres=0)
I = 30.0  # Intervention magnitude (umol/kg)

# ---------------------------------------------------------
# Helpers: compute carbonate system properties
# ---------------------------------------------------------
def compute_pCO2(alk, dic, sal, temp, pres):
    out = pyco2.sys(
        par1=alk,
        par2=dic,
        par1_type=1, 
        par2_type=2,
        salinity=sal, 
        temperature=temp, 
        pressure=pres,
    )
    return out["pCO2"]

def get_theta(alk, dic, temp, sal, pres=0):
    """Calculates theta = dpCO2/dAlk using a central difference."""
    d_alk = 0.01
    pCO2_high = compute_pCO2(alk + d_alk, dic, sal, temp, pres)
    pCO2_low  = compute_pCO2(alk - d_alk, dic, sal, temp, pres)
    return (pCO2_high - pCO2_low) / (2 * d_alk)

def get_isoQ(alk, dic, temp, sal, pres=0):
    """Extracts the isocapnic quotient from PyCO2SYS."""
    out = pyco2.sys(
        par1=alk,
        par2=dic,
        par1_type=1,
        par2_type=2,
        salinity=sal,
        temperature=temp,
        pressure=pres,
    )
    return out["isocapnic_quotient"] 
    
# ---------------------------------------------------------
# 2. Compute baseline eta_max (Q_obs)
# ---------------------------------------------------------
eta_max_obs = get_isoQ(obs['alk'], obs['dic'], obs['temp'], obs['sal'])
eta_max_model = get_isoQ(model['alk'], model['dic'], model['temp'], model['sal'])

# ---------------------------------------------------------
# 3. Calculate Sensitivities (dIsoQ/dX) at OBS state
# ---------------------------------------------------------
eps = 1e-2 

dQ_dAlk = (get_isoQ(obs['alk']+eps, obs['dic'], obs['temp'], obs['sal']) - 
           get_isoQ(obs['alk']-eps, obs['dic'], obs['temp'], obs['sal'])) / (2*eps)

dQ_dDic = (get_isoQ(obs['alk'], obs['dic']+eps, obs['temp'], obs['sal']) - 
           get_isoQ(obs['alk'], obs['dic']-eps, obs['temp'], obs['sal'])) / (2*eps)

dQ_dTemp = (get_isoQ(obs['alk'], obs['dic'], obs['temp']+eps, obs['sal']) - 
            get_isoQ(obs['alk'], obs['dic'], obs['temp']-eps, obs['sal'])) / (2*eps)

dQ_dSal = (get_isoQ(obs['alk'], obs['dic'], obs['temp'], obs['sal']+eps) - 
           get_isoQ(obs['alk'], obs['dic'], obs['temp'], obs['sal']-eps)) / (2*eps)

# ---------------------------------------------------------
# 4. Compute Biases and Percent Errors (Linearized)
# ---------------------------------------------------------
variables = ['alk', 'dic', 'temp', 'sal']
sensitivities_Q = [dQ_dAlk, dQ_dDic, dQ_dTemp, dQ_dSal]

results = []
for var, sens in zip(variables, sensitivities_Q):
    bias = model[var] - obs[var]
    error_Q = sens * bias 
    
    # Calculate error as a percentage of baseline eta_max
    pct_error_Q = (error_Q / eta_max_obs) * 100
    
    results.append({
        "Variable": var,
        "Sensitivity (∂Q/∂X)": sens,
        "Total Bias (ΔX)": bias,
        "Linearized Error [Units]": error_Q,
        "Error (% of eta_max)": pct_error_Q
    })

df = pd.DataFrame(results)
sum_linear_error_Q = df["Linearized Error [Units]"].sum()
sum_pct_error_Q = df["Error (% of eta_max)"].sum()

# ---------------------------------------------------------
# 5. Calculate ACTUAL Percent Error (Direct Computation)
# ---------------------------------------------------------
actual_error_Q = eta_max_model - eta_max_obs
actual_pct_error_Q = (actual_error_Q / eta_max_obs) * 100

# ---------------------------------------------------------
# 6. Output Results
# ---------------------------------------------------------
print(f"Observed eta_max (Baseline Q): {eta_max_obs:.6f}")
print(f"Model eta_max (Biased Q):     {eta_max_model:.6f}\n")

print(df.to_string(index=False, formatters={
    "Sensitivity (∂Q/∂X)": "{:.6f}".format,
    "Total Bias (ΔX)": "{:.2f}".format,
    "Linearized Error [Units]": "{:.6f}".format,
    "Error (% of eta_max)": "{:+.2f}%".format
}))
print("-" * 75)
print(f"Expected Bias Summary (Sum of terms): {sum_linear_error_Q:10.6f} ({sum_pct_error_Q:+.2f}%)")
print(f"Actual Bias Summary (Direct mod-obs):  {actual_error_Q:10.6f} ({actual_pct_error_Q:+.2f}%)")
print(f"Residual (Higher-order terms):         {actual_error_Q - sum_linear_error_Q:10.6f} ({(actual_pct_error_Q - sum_pct_error_Q):+.2f}%)")

Observed eta_max (Baseline Q): 1.165382
Model eta_max (Biased Q):     1.162433

Variable Sensitivity (∂Q/∂X) Total Bias (ΔX) Linearized Error [Units] Error (% of eta_max)
     alk            0.000596          -60.00                -0.035771               -3.07%
     dic           -0.000649          -50.00                 0.032432               +2.78%
    temp            0.000300            1.50                 0.000450               +0.04%
     sal           -0.000168            0.00                -0.000000               -0.00%
---------------------------------------------------------------------------
Expected Bias Summary (Sum of terms):  -0.002889 (-0.25%)
Actual Bias Summary (Direct mod-obs):   -0.002949 (-0.25%)
Residual (Higher-order terms):          -0.000061 (-0.01%)

print(get_isoQ(2380.0, 2160.0, 10.0, 33.75))

1.1653819662966605

import numpy as np
import pandas as pd
import PyCO2SYS as pyco2

# ---------------------------------------------------------
# 1. Setup states and intervention
# ---------------------------------------------------------
obs = dict(alk=2380.0, dic=2160.0, sal=33.75, temp=10.0, pres=0)
model = dict(alk=2320.0, dic=2110.0, sal=33.75, temp=11.5, pres=0)
I = 30.0  # Intervention magnitude (umol/kg)

# ---------------------------------------------------------
# Helper: Compute pCO2 to find baseline target
# ---------------------------------------------------------
def get_baseline_pCO2(alk, dic, sal, temp, pres=0):
    out = pyco2.sys(
        par1=alk, par2=dic, par1_type=1, par2_type=2,
        salinity=sal, temperature=temp, pressure=pres
    )
    return out["pCO2"]

# ---------------------------------------------------------
# Core Helper: Compute eta_max explicitly holding pCO2 constant
# ---------------------------------------------------------
def compute_eta_max(alk, dic, sal, temp, pres=0, d_alk=1.0):
    """
    Computes eta_max = dDIC / dAlk holding pCO2 constant.
    1. Finds target baseline pCO2.
    2. Adds d_alk.
    3. Uses PyCO2SYS to calculate the new equilibrium DIC at that target pCO2.
    """
    # Step 1: Get target pCO2
    target_pCO2 = get_baseline_pCO2(alk, dic, sal, temp, pres)
    
    # Step 2 & 3: Compute state after Alk addition holding pCO2 constant
    # par1 = Alk (type 1), par2 = pCO2 (type 4)
    out_equilibrium = pyco2.sys(
        par1=alk + d_alk,
        par2=target_pCO2,
        par1_type=1,
        par2_type=4,
        salinity=sal,
        temperature=temp,
        pressure=pres
    )
    
    new_dic = out_equilibrium["dic"]
    delta_dic = new_dic - dic
    
    # eta_max = delta_DIC / delta_Alk
    return delta_dic / d_alk

# ---------------------------------------------------------
# 2. Compute Baseline eta_max values for Obs and Model
# ---------------------------------------------------------
eta_max_obs = compute_eta_max(obs['alk'], obs['dic'], obs['sal'], obs['temp'])
eta_max_model = compute_eta_max(model['alk'], model['dic'], model['sal'], model['temp'])

# ---------------------------------------------------------
# 3. Calculate Sensitivities (dEta_max / dX) at OBS state
# ---------------------------------------------------------
eps = 1e-2 

dEta_dAlk = (compute_eta_max(obs['alk']+eps, obs['dic'], obs['sal'], obs['temp']) - 
             compute_eta_max(obs['alk']-eps, obs['dic'], obs['sal'], obs['temp'])) / (2*eps)

dEta_dDic = (compute_eta_max(obs['alk'], obs['dic']+eps, obs['sal'], obs['temp']) - 
             compute_eta_max(obs['alk'], obs['dic']-eps, obs['sal'], obs['temp'])) / (2*eps)

dEta_dTemp = (compute_eta_max(obs['alk'], obs['dic'], obs['sal'], obs['temp']+eps) - 
              compute_eta_max(obs['alk'], obs['dic'], obs['sal'], obs['temp']-eps)) / (2*eps)

dEta_dSal = (compute_eta_max(obs['alk'], obs['dic'], obs['sal']+eps, obs['temp']) - 
             compute_eta_max(obs['alk'], obs['dic'], obs['sal']-eps, obs['temp'])) / (2*eps)

# ---------------------------------------------------------
# 4. Compute Biases and Percent Errors (Linearized)
# ---------------------------------------------------------
variables = ['alk', 'dic', 'temp', 'sal']
sensitivities_eta = [dEta_dAlk, dEta_dDic, dEta_dTemp, dEta_dSal]

results = []
for var, sens in zip(variables, sensitivities_eta):
    bias = model[var] - obs[var]
    error_eta = sens * bias 
    
    # Express the error contribution as a percentage of the true baseline eta_max
    pct_error_eta = (error_eta / eta_max_obs) * 100
    
    results.append({
        "Variable": var,
        "Sensitivity (∂η/∂X)": sens,
        "Total Bias (ΔX)": bias,
        "Linearized Error [Units]": error_eta,
        "Error (% of eta_max)": pct_error_eta
    })

df = pd.DataFrame(results)
sum_linear_error_eta = df["Linearized Error [Units]"].sum()
sum_pct_error_eta = df["Error (% of eta_max)"].sum()

# ---------------------------------------------------------
# 5. Calculate ACTUAL Percent Error
# ---------------------------------------------------------
actual_error_eta = eta_max_model - eta_max_obs
actual_pct_error_eta = (actual_error_eta / eta_max_obs) * 100

# ---------------------------------------------------------
# 6. Output Results
# ---------------------------------------------------------
print(f"Observed eta_max (Baseline dDIC/dAlk): {eta_max_obs:.6f}")
print(f"Model eta_max (Biased dDIC/dAlk):     {eta_max_model:.6f}\n")

print(df.to_string(index=False, formatters={
    "Sensitivity (∂η/∂X)": "{:.6f}".format,
    "Total Bias (ΔX)": "{:.2f}".format,
    "Linearized Error [Units]": "{:.6f}".format,
    "Error (% of eta_max)": "{:+.2f}%".format
}))
print("-" * 75)
print(f"Expected Bias Summary (Sum of terms): {sum_linear_error_eta:10.6f} ({sum_pct_error_eta:+.2f}%)")
print(f"Actual Bias Summary (Direct mod-obs):  {actual_error_eta:10.6f} ({actual_pct_error_eta:+.2f}%)")
print(f"Residual (Higher-order terms):         {actual_error_eta - sum_linear_error_eta:10.6f} ({(actual_pct_error_eta - sum_pct_error_eta):+.2f}%)")

Observed eta_max (Baseline dDIC/dAlk): 0.858073
Model eta_max (Biased dDIC/dAlk):     0.860250

Variable Sensitivity (∂η/∂X) Total Bias (ΔX) Linearized Error [Units] Error (% of eta_max)
     alk           -0.000439          -60.00                 0.026339               +3.07%
     dic            0.000478          -50.00                -0.023881               -2.78%
    temp           -0.000221            1.50                -0.000332               -0.04%
     sal            0.000124            0.00                 0.000000               +0.00%
---------------------------------------------------------------------------
Expected Bias Summary (Sum of terms):   0.002127 (+0.25%)
Actual Bias Summary (Direct mod-obs):    0.002177 (+0.25%)
Residual (Higher-order terms):           0.000050 (+0.01%)

Quantifying pCO_2 Bias from Mean-State Errors

Quantifying pCO2pCO_2pCO2​ Bias from Mean-State Errors¶

1. Decomposition of pCO2pCO_2pCO2​ Bias¶

Full Multi-Parameter Analysis of pCO2pCO_2pCO2​ Sensitivity¶

Quantifying $pCO_2$ Bias from Mean-State Errors¶

1. Decomposition of $pCO_2$ Bias¶

Full Multi-Parameter Analysis of $pCO_2$ Sensitivity¶