Organic-matter effects on soil water holding

Organic matter and bulk density are the two soil properties practitioners can most directly influence through management. This notebook shows how increasing organic matter — and easing compaction — shifts water into plant-available storage (field capacity − wilting point) and drainable storage (saturation − field capacity), the fast-draining pore space that matters most for stormwater infiltration and detention. Gains are largest in coarse-textured soils and diminish in clays.

Two methods are compared across all 12 USDA texture classes:

• Section 2 — ROSETTA + Minasny & McBratney (2018): ROSETTA’s mineral baseline plus Minasny’s & McBratney’s empirical organic-carbon increments (the OC sensitivity preferred over Saxton–Rawls). Two sliders: mineral bulk density and organic matter (0–8 %).
• Section 3 — Comparison to Saxton & Rawls (2006): an independent, self-contained PTF taking sand/clay/OM directly, with a validation of its OC sensitivity against Minasny & McBratney (§3.1).

Section 1 (collapsible below) documents how we established ROSETTA’s mineral-baseline organic carbon — the anchor point for the Section 2 blend.

import numpy as np
import pandas as pd
import hvplot.pandas  # noqa: F401  (registers the .hvplot accessor)
import holoviews as hv
from scipy.stats import linregress

# Shared display helper (and the shared pandas float_format, set on import). See notebooks/_helpers.py.
from _helpers import show, soil_water_texture_band_diagram, soil_water_bd_om_blend_table, soil_water_table_html, VB, OC_BASELINE_PCT, oc_baseline_for_bd, OC_BD_SLOPE, OC_BD_INTERCEPT, MM_SLOPES, MM_GROUP

# ROSETTA texture × bulk-density baseline produced by rosetta_porosity_by_texture.ipynb
result = pd.read_csv("rosetta_porosity_by_texture.csv")

# Reconstruct the shared constants from the baseline table (canonical sand -> clay order)
TEXTURE_CLASSES = {
    row.texture_class: (row.sand_pct, row.silt_pct, row.clay_pct)
    for row in result[["texture_class", "sand_pct", "silt_pct", "clay_pct"]]
    .drop_duplicates()
    .itertuples(index=False)
}
HYDROLOGIC_SOIL_GROUP = dict(
    result[["texture_class", "hydrologic_soil_group"]].drop_duplicates().itertuples(index=False, name=None)
)
bulk_densities = np.sort(result["bulk_density_g_cm3"].unique())

texture_x = {cls: i for i, cls in enumerate(TEXTURE_CLASSES)}  # sand=0 ... clay=11
texture_ticks = [(i, f"{cls} ({HYDROLOGIC_SOIL_GROUP[cls]})") for cls, i in texture_x.items()]
INCHES_PER_FOOT = 12.0  # volumetric water content (cm3/cm3) -> inches per foot of soil depth

print(f"loaded {len(result)} ROSETTA rows: {len(TEXTURE_CLASSES)} textures x {len(bulk_densities)} bulk densities")
show(result)

loaded 144 ROSETTA rows: 12 textures x 12 bulk densities

	texture_class	hydrologic_soil_group	sand_pct	silt_pct	clay_pct	bulk_density_g_cm3	total_porosity	field_capacity_porosity	wilting_point_porosity	available_water_capacity	ksat_cm_day	ksat_in_hr	k_fc_cm_day	k_wp_cm_day	theta_r	vg_alpha_1cm	vg_n	k0_cm_day	mualem_L	implausible_bd
0	sand	A	92	5	3	0.800	0.560	0.141	0.072	0.069	992.244	16.277	0.008	0.000	0.070	0.029	1.848	53.148	-0.574	False
1	sand	A	92	5	3	0.900	0.531	0.129	0.068	0.062	874.440	14.344	0.006	0.000	0.066	0.026	1.930	36.508	-0.473	False
2	sand	A	92	5	3	1.000	0.503	0.113	0.064	0.049	803.505	13.181	0.004	0.000	0.063	0.025	2.033	29.219	-0.477	False
3	sand	A	92	5	3	1.100	0.475	0.096	0.061	0.036	763.127	12.518	0.003	0.000	0.060	0.025	2.154	24.929	-0.542	False
4	sand	A	92	5	3	1.200	0.450	0.082	0.058	0.024	731.398	11.998	0.002	0.000	0.058	0.026	2.289	22.146	-0.638	False
5	sand	A	92	5	3	1.300	0.425	0.072	0.056	0.016	690.465	11.327	0.002	0.000	0.056	0.027	2.430	20.779	-0.734	False
6	sand	A	92	5	3	1.400	0.401	0.064	0.054	0.010	630.566	10.344	0.001	0.000	0.054	0.029	2.563	20.654	-0.809	False
7	sand	A	92	5	3	1.500	0.376	0.059	0.052	0.007	544.378	8.930	0.001	0.000	0.052	0.031	2.653	21.556	-0.857	False
8	sand	A	92	5	3	1.600	0.352	0.057	0.050	0.006	424.038	6.956	0.001	0.000	0.050	0.032	2.638	23.625	-0.874	False
9	sand	A	92	5	3	1.700	0.327	0.056	0.049	0.007	301.562	4.947	0.001	0.000	0.049	0.033	2.532	27.468	-0.876	False
10	sand	A	92	5	3	1.800	0.303	0.055	0.048	0.007	218.069	3.577	0.002	0.000	0.048	0.035	2.449	36.333	-0.883	False
11	sand	A	92	5	3	1.900	0.281	0.054	0.047	0.007	173.663	2.849	0.003	0.000	0.047	0.038	2.415	83.843	-0.888	False
12	loamy sand	A	82	12	6	0.800	0.554	0.253	0.097	0.156	487.050	7.990	0.011	0.000	0.074	0.018	1.545	16.326	-0.346	False
13	loamy sand	A	82	12	6	0.900	0.525	0.228	0.088	0.140	421.489	6.914	0.009	0.000	0.071	0.018	1.589	14.527	-0.322	False
14	loamy sand	A	82	12	6	1.000	0.498	0.203	0.081	0.123	369.826	6.067	0.008	0.000	0.068	0.018	1.636	13.651	-0.361	False
15	loamy sand	A	82	12	6	1.100	0.472	0.180	0.074	0.106	321.878	5.280	0.008	0.000	0.066	0.019	1.684	13.228	-0.433	False
16	loamy sand	A	82	12	6	1.200	0.447	0.159	0.070	0.090	274.748	4.507	0.007	0.000	0.064	0.020	1.729	13.048	-0.520	False
17	loamy sand	A	82	12	6	1.300	0.423	0.141	0.066	0.076	229.151	3.759	0.006	0.000	0.061	0.021	1.769	13.056	-0.614	False
18	loamy sand	A	82	12	6	1.400	0.400	0.127	0.062	0.065	186.012	3.051	0.006	0.000	0.059	0.022	1.799	13.213	-0.706	False
19	loamy sand	A	82	12	6	1.500	0.377	0.116	0.060	0.057	145.188	2.382	0.006	0.000	0.057	0.024	1.810	13.454	-0.796	False
20	loamy sand	A	82	12	6	1.600	0.353	0.111	0.058	0.053	105.214	1.726	0.006	0.000	0.055	0.025	1.782	13.723	-0.890	False
21	loamy sand	A	82	12	6	1.700	0.328	0.112	0.058	0.054	69.610	1.142	0.007	0.000	0.054	0.027	1.716	14.156	-0.996	False
22	loamy sand	A	82	12	6	1.800	0.305	0.112	0.058	0.054	45.577	0.748	0.008	0.000	0.053	0.028	1.649	16.243	-1.096	False
23	loamy sand	A	82	12	6	1.900	0.282	0.109	0.058	0.051	32.713	0.537	0.015	0.000	0.052	0.030	1.609	31.984	-1.179	False
24	sandy loam	A	65	25	10	0.800	0.544	0.339	0.129	0.210	345.058	5.660	0.013	0.000	0.081	0.010	1.455	4.963	0.012	False
25	sandy loam	A	65	25	10	0.900	0.517	0.313	0.119	0.194	261.554	4.291	0.012	0.000	0.078	0.010	1.472	4.803	-0.047	False
26	sandy loam	A	65	25	10	1.000	0.490	0.289	0.110	0.178	196.156	3.218	0.011	0.000	0.075	0.011	1.486	4.834	-0.143	False
27	sandy loam	A	65	25	10	1.100	0.466	0.266	0.103	0.163	145.698	2.390	0.010	0.000	0.073	0.012	1.496	4.954	-0.257	False
28	sandy loam	A	65	25	10	1.200	0.442	0.246	0.097	0.148	107.304	1.760	0.009	0.000	0.071	0.013	1.501	5.111	-0.381	False
29	sandy loam	A	65	25	10	1.300	0.419	0.228	0.093	0.135	78.162	1.282	0.009	0.000	0.068	0.014	1.502	5.289	-0.509	False
30	sandy loam	A	65	25	10	1.400	0.396	0.213	0.089	0.124	56.362	0.925	0.008	0.000	0.066	0.015	1.496	5.509	-0.641	False
31	sandy loam	A	65	25	10	1.500	0.374	0.201	0.087	0.114	39.936	0.655	0.008	0.000	0.064	0.016	1.481	5.772	-0.776	False
32	sandy loam	A	65	25	10	1.600	0.351	0.192	0.086	0.106	27.110	0.445	0.008	0.000	0.063	0.017	1.451	6.082	-0.926	False
33	sandy loam	A	65	25	10	1.700	0.328	0.188	0.089	0.099	17.322	0.284	0.008	0.000	0.061	0.018	1.406	6.551	-1.113	False
34	sandy loam	A	65	25	10	1.800	0.306	0.184	0.092	0.092	10.893	0.179	0.008	0.000	0.060	0.019	1.360	7.607	-1.347	False
35	sandy loam	A	65	25	10	1.900	0.284	0.179	0.095	0.083	7.269	0.119	0.011	0.000	0.060	0.021	1.322	11.189	-1.619	True
36	loam	B	40	40	20	0.800	0.559	0.421	0.163	0.258	256.843	4.213	0.022	0.000	0.098	0.005	1.453	1.999	0.334	False
37	loam	B	40	40	20	0.900	0.530	0.395	0.154	0.241	165.697	2.718	0.018	0.000	0.095	0.005	1.456	1.745	0.281	False
38	loam	B	40	40	20	1.000	0.502	0.371	0.147	0.224	104.354	1.712	0.016	0.000	0.093	0.006	1.456	1.624	0.180	False
39	loam	B	40	40	20	1.100	0.476	0.348	0.141	0.207	65.280	1.071	0.014	0.000	0.091	0.006	1.452	1.585	0.054	False
40	loam	B	40	40	20	1.200	0.451	0.326	0.137	0.190	41.158	0.675	0.012	0.000	0.089	0.006	1.444	1.588	-0.083	False
41	loam	B	40	40	20	1.300	0.427	0.307	0.133	0.174	26.290	0.431	0.011	0.000	0.087	0.007	1.432	1.613	-0.227	False
42	loam	B	40	40	20	1.400	0.404	0.290	0.131	0.159	17.009	0.279	0.010	0.000	0.086	0.007	1.416	1.657	-0.377	False
43	loam	B	40	40	20	1.500	0.381	0.275	0.130	0.145	11.107	0.182	0.009	0.000	0.084	0.008	1.393	1.731	-0.538	False
44	loam	B	40	40	20	1.600	0.358	0.262	0.130	0.131	7.236	0.119	0.008	0.000	0.083	0.009	1.362	1.852	-0.732	False
45	loam	B	40	40	20	1.700	0.336	0.250	0.133	0.117	4.656	0.076	0.008	0.000	0.081	0.009	1.324	2.071	-0.994	False
46	loam	B	40	40	20	1.800	0.314	0.239	0.136	0.103	3.015	0.049	0.008	0.000	0.080	0.010	1.285	2.502	-1.350	False
47	loam	B	40	40	20	1.900	0.293	0.228	0.139	0.090	2.038	0.033	0.008	0.000	0.080	0.011	1.252	3.372	-1.803	True
48	silt loam	B	20	65	15	0.800	0.557	0.458	0.146	0.313	330.144	5.416	0.069	0.000	0.092	0.003	1.572	1.515	0.746	False
49	silt loam	B	20	65	15	0.900	0.527	0.431	0.138	0.293	203.429	3.337	0.056	0.000	0.088	0.003	1.574	1.286	0.741	False
50	silt loam	B	20	65	15	1.000	0.499	0.404	0.131	0.273	127.725	2.095	0.046	0.000	0.085	0.003	1.572	1.166	0.668	False
51	silt loam	B	20	65	15	1.100	0.473	0.378	0.125	0.253	81.420	1.336	0.038	0.000	0.083	0.003	1.565	1.109	0.553	False
52	silt loam	B	20	65	15	1.200	0.449	0.354	0.121	0.233	52.349	0.859	0.032	0.000	0.080	0.004	1.554	1.095	0.422	False
53	silt loam	B	20	65	15	1.300	0.427	0.331	0.117	0.214	33.895	0.556	0.027	0.000	0.078	0.004	1.538	1.110	0.284	False
54	silt loam	B	20	65	15	1.400	0.405	0.310	0.115	0.195	22.196	0.364	0.023	0.000	0.077	0.004	1.516	1.144	0.141	False
55	silt loam	B	20	65	15	1.500	0.384	0.291	0.113	0.178	14.735	0.242	0.019	0.000	0.075	0.005	1.489	1.195	-0.010	False
56	silt loam	B	20	65	15	1.600	0.364	0.274	0.114	0.161	9.942	0.163	0.016	0.000	0.074	0.005	1.455	1.269	-0.178	False
57	silt loam	B	20	65	15	1.700	0.343	0.259	0.115	0.144	6.902	0.113	0.013	0.000	0.073	0.006	1.416	1.390	-0.376	False
58	silt loam	B	20	65	15	1.800	0.323	0.244	0.117	0.127	5.086	0.083	0.012	0.000	0.073	0.007	1.376	1.603	-0.617	True
59	silt loam	B	20	65	15	1.900	0.303	0.230	0.118	0.111	4.470	0.073	0.011	0.000	0.073	0.008	1.340	2.003	-0.915	True
60	silt	B/D	7	88	5	0.800	0.562	0.351	0.103	0.249	1381.028	22.655	0.039	0.000	0.083	0.006	1.711	3.927	0.390	False
61	silt	B/D	7	88	5	0.900	0.533	0.380	0.102	0.278	690.764	11.331	0.064	0.000	0.078	0.004	1.707	2.692	0.642	False
62	silt	B/D	7	88	5	1.000	0.506	0.377	0.100	0.277	369.838	6.067	0.067	0.000	0.074	0.004	1.701	2.103	0.713	False
63	silt	B/D	7	88	5	1.100	0.482	0.361	0.096	0.265	215.558	3.536	0.059	0.000	0.070	0.004	1.693	1.801	0.686	False
64	silt	B/D	7	88	5	1.200	0.459	0.341	0.093	0.248	136.101	2.233	0.050	0.000	0.068	0.004	1.680	1.669	0.601	False
65	silt	B/D	7	88	5	1.300	0.438	0.319	0.090	0.229	92.347	1.515	0.041	0.000	0.066	0.004	1.664	1.632	0.482	False
66	silt	B/D	7	88	5	1.400	0.418	0.297	0.088	0.209	68.169	1.118	0.033	0.000	0.064	0.004	1.644	1.661	0.342	False
67	silt	B/D	7	88	5	1.500	0.398	0.274	0.086	0.189	58.825	0.965	0.026	0.000	0.063	0.005	1.620	1.760	0.182	False
68	silt	B/D	7	88	5	1.600	0.379	0.252	0.084	0.167	70.476	1.156	0.020	0.000	0.062	0.006	1.591	1.956	-0.007	False
69	silt	B/D	7	88	5	1.700	0.360	0.228	0.084	0.145	129.258	2.120	0.015	0.000	0.062	0.008	1.560	2.337	-0.241	True
70	silt	B/D	7	88	5	1.800	0.341	0.201	0.082	0.120	571.274	9.371	0.011	0.000	0.063	0.010	1.540	3.146	-0.535	True
71	silt	B/D	7	88	5	1.900	0.323	0.174	0.080	0.094	1463.279	24.004	0.009	0.000	0.064	0.015	1.527	4.956	-0.883	True
72	sandy clay loam	C	60	13	27	0.800	0.586	0.386	0.186	0.200	231.903	3.804	0.010	0.000	0.112	0.013	1.350	6.082	-0.785	False
73	sandy clay loam	C	60	13	27	0.900	0.560	0.367	0.176	0.191	174.015	2.855	0.009	0.000	0.108	0.013	1.359	4.962	-0.751	False
74	sandy clay loam	C	60	13	27	1.000	0.534	0.348	0.168	0.181	128.898	2.114	0.008	0.000	0.106	0.013	1.366	4.333	-0.755	False
75	sandy clay loam	C	60	13	27	1.100	0.509	0.330	0.160	0.170	93.746	1.538	0.007	0.000	0.103	0.013	1.370	3.974	-0.787	False
76	sandy clay loam	C	60	13	27	1.200	0.484	0.313	0.154	0.159	66.541	1.092	0.006	0.000	0.101	0.014	1.371	3.742	-0.843	False
77	sandy clay loam	C	60	13	27	1.300	0.459	0.298	0.149	0.149	45.805	0.751	0.006	0.000	0.099	0.014	1.368	3.551	-0.921	False
78	sandy clay loam	C	60	13	27	1.400	0.434	0.285	0.146	0.138	30.569	0.501	0.006	0.000	0.097	0.014	1.359	3.371	-1.023	False
79	sandy clay loam	C	60	13	27	1.500	0.408	0.273	0.145	0.128	19.687	0.323	0.005	0.000	0.096	0.014	1.342	3.212	-1.159	False
80	sandy clay loam	C	60	13	27	1.600	0.381	0.264	0.147	0.117	12.022	0.197	0.005	0.000	0.094	0.015	1.313	3.111	-1.363	False
81	sandy clay loam	C	60	13	27	1.700	0.354	0.257	0.152	0.105	6.918	0.113	0.005	0.000	0.094	0.015	1.275	3.166	-1.697	False
82	sandy clay loam	C	60	13	27	1.800	0.327	0.249	0.157	0.091	3.975	0.065	0.005	0.000	0.093	0.015	1.239	3.572	-2.195	True
83	sandy clay loam	C	60	13	27	1.900	0.301	0.237	0.159	0.078	2.444	0.040	0.006	0.000	0.093	0.016	1.211	4.640	-2.809	True
84	clay loam	D	30	35	35	0.800	0.599	0.448	0.205	0.244	215.577	3.536	0.016	0.000	0.121	0.007	1.379	2.611	-0.202	False
85	clay loam	D	30	35	35	0.900	0.570	0.427	0.196	0.231	137.810	2.261	0.013	0.000	0.118	0.007	1.380	2.049	-0.189	False
86	clay loam	D	30	35	35	1.000	0.542	0.406	0.189	0.216	85.477	1.402	0.010	0.000	0.115	0.007	1.378	1.713	-0.237	False
87	clay loam	D	30	35	35	1.100	0.514	0.385	0.184	0.202	51.959	0.852	0.009	0.000	0.113	0.007	1.374	1.511	-0.320	False
88	clay loam	D	30	35	35	1.200	0.487	0.366	0.179	0.187	31.207	0.512	0.008	0.000	0.112	0.007	1.366	1.386	-0.418	False
89	clay loam	D	30	35	35	1.300	0.461	0.349	0.176	0.172	18.636	0.306	0.007	0.000	0.111	0.007	1.355	1.299	-0.530	False
90	clay loam	D	30	35	35	1.400	0.434	0.332	0.175	0.157	11.137	0.183	0.006	0.000	0.110	0.008	1.340	1.232	-0.663	False
91	clay loam	D	30	35	35	1.500	0.408	0.317	0.174	0.142	6.693	0.110	0.006	0.000	0.109	0.008	1.320	1.186	-0.830	False
92	clay loam	D	30	35	35	1.600	0.381	0.302	0.176	0.127	4.037	0.066	0.005	0.000	0.109	0.008	1.294	1.174	-1.058	False
93	clay loam	D	30	35	35	1.700	0.355	0.289	0.178	0.110	2.447	0.040	0.005	0.000	0.109	0.008	1.263	1.228	-1.399	False
94	clay loam	D	30	35	35	1.800	0.330	0.274	0.181	0.094	1.524	0.025	0.005	0.000	0.109	0.009	1.231	1.412	-1.890	True
95	clay loam	D	30	35	35	1.900	0.305	0.260	0.181	0.078	1.009	0.017	0.006	0.000	0.110	0.009	1.204	1.848	-2.537	True
96	silty clay loam	D	10	56	34	0.800	0.609	0.479	0.199	0.280	229.209	3.760	0.033	0.000	0.122	0.005	1.436	2.158	0.129	False
97	silty clay loam	D	10	56	34	0.900	0.578	0.460	0.192	0.267	138.366	2.270	0.026	0.000	0.118	0.004	1.438	1.567	0.221	False
98	silty clay loam	D	10	56	34	1.000	0.549	0.438	0.186	0.252	81.976	1.345	0.021	0.000	0.115	0.004	1.437	1.219	0.207	False
99	silty clay loam	D	10	56	34	1.100	0.520	0.416	0.181	0.236	47.675	0.782	0.017	0.000	0.114	0.004	1.432	1.013	0.132	False
100	silty clay loam	D	10	56	34	1.200	0.492	0.395	0.176	0.218	27.436	0.450	0.015	0.000	0.112	0.004	1.424	0.895	0.030	False
101	silty clay loam	D	10	56	34	1.300	0.465	0.374	0.174	0.200	15.812	0.259	0.013	0.000	0.112	0.005	1.412	0.825	-0.087	False
102	silty clay loam	D	10	56	34	1.400	0.439	0.354	0.172	0.182	9.222	0.151	0.011	0.000	0.111	0.005	1.394	0.779	-0.221	False
103	silty clay loam	D	10	56	34	1.500	0.413	0.336	0.172	0.164	5.486	0.090	0.009	0.000	0.111	0.005	1.372	0.752	-0.379	False
104	silty clay loam	D	10	56	34	1.600	0.388	0.318	0.172	0.146	3.339	0.055	0.008	0.000	0.111	0.005	1.343	0.748	-0.573	False
105	silty clay loam	D	10	56	34	1.700	0.362	0.301	0.174	0.127	2.092	0.034	0.007	0.000	0.111	0.006	1.311	0.779	-0.829	True
106	silty clay loam	D	10	56	34	1.800	0.338	0.285	0.175	0.109	1.377	0.023	0.006	0.000	0.111	0.006	1.279	0.868	-1.167	True
107	silty clay loam	D	10	56	34	1.900	0.315	0.268	0.176	0.093	0.974	0.016	0.006	0.000	0.111	0.007	1.250	1.064	-1.599	True
108	sandy clay	D	50	7	43	0.800	0.616	0.408	0.222	0.186	240.026	3.937	0.009	0.000	0.133	0.018	1.302	8.437	-1.512	False
109	sandy clay	D	50	7	43	0.900	0.590	0.398	0.216	0.183	166.103	2.725	0.008	0.000	0.130	0.016	1.305	6.062	-1.395	False
110	sandy clay	D	50	7	43	1.000	0.564	0.386	0.209	0.176	115.164	1.889	0.006	0.000	0.126	0.016	1.305	4.659	-1.351	False
111	sandy clay	D	50	7	43	1.100	0.539	0.372	0.204	0.168	79.444	1.303	0.005	0.000	0.124	0.015	1.304	3.821	-1.345	False
112	sandy clay	D	50	7	43	1.200	0.514	0.358	0.199	0.159	54.335	0.891	0.005	0.000	0.122	0.015	1.300	3.312	-1.370	False
113	sandy clay	D	50	7	43	1.300	0.488	0.345	0.195	0.149	36.772	0.603	0.004	0.000	0.120	0.015	1.293	2.968	-1.434	False
114	sandy clay	D	50	7	43	1.400	0.463	0.332	0.194	0.139	24.361	0.400	0.004	0.000	0.119	0.015	1.282	2.687	-1.543	False
115	sandy clay	D	50	7	43	1.500	0.436	0.321	0.193	0.127	15.598	0.256	0.004	0.000	0.118	0.015	1.265	2.441	-1.723	True
116	sandy clay	D	50	7	43	1.600	0.408	0.310	0.196	0.115	9.416	0.154	0.003	0.000	0.117	0.015	1.242	2.248	-2.029	True
117	sandy clay	D	50	7	43	1.700	0.380	0.300	0.200	0.100	5.387	0.088	0.003	0.000	0.117	0.015	1.214	2.186	-2.542	True
118	sandy clay	D	50	7	43	1.800	0.350	0.288	0.203	0.084	3.138	0.051	0.003	0.000	0.118	0.015	1.186	2.423	-3.307	True
119	sandy clay	D	50	7	43	1.900	0.322	0.272	0.203	0.069	2.041	0.033	0.004	0.000	0.120	0.015	1.165	3.192	-4.266	True
120	silty clay	D	7	47	46	0.800	0.635	0.473	0.224	0.249	217.292	3.565	0.019	0.000	0.137	0.007	1.372	3.245	-0.458	False
121	silty clay	D	7	47	46	0.900	0.605	0.460	0.219	0.241	131.925	2.164	0.015	0.000	0.133	0.007	1.372	2.269	-0.334	False
122	silty clay	D	7	47	46	1.000	0.575	0.443	0.213	0.230	78.732	1.292	0.012	0.000	0.130	0.006	1.369	1.657	-0.298	False
123	silty clay	D	7	47	46	1.100	0.546	0.425	0.209	0.216	45.845	0.752	0.010	0.000	0.128	0.006	1.364	1.305	-0.330	False
124	silty clay	D	7	47	46	1.200	0.518	0.406	0.205	0.201	26.160	0.429	0.008	0.000	0.127	0.006	1.356	1.093	-0.405	False
125	silty clay	D	7	47	46	1.300	0.489	0.388	0.203	0.185	14.805	0.243	0.007	0.000	0.127	0.006	1.344	0.954	-0.503	False
126	silty clay	D	7	47	46	1.400	0.461	0.370	0.202	0.168	8.421	0.138	0.006	0.000	0.126	0.006	1.328	0.858	-0.624	False
127	silty clay	D	7	47	46	1.500	0.433	0.352	0.201	0.151	4.871	0.080	0.005	0.000	0.126	0.006	1.309	0.789	-0.786	False
128	silty clay	D	7	47	46	1.600	0.405	0.335	0.202	0.133	2.887	0.047	0.005	0.000	0.127	0.007	1.285	0.749	-1.013	True
129	silty clay	D	7	47	46	1.700	0.377	0.318	0.203	0.114	1.773	0.029	0.004	0.000	0.128	0.007	1.259	0.749	-1.338	True
130	silty clay	D	7	47	46	1.800	0.350	0.300	0.204	0.096	1.152	0.019	0.004	0.000	0.129	0.007	1.233	0.817	-1.791	True
131	silty clay	D	7	47	46	1.900	0.324	0.282	0.203	0.079	0.818	0.013	0.004	0.000	0.132	0.008	1.209	1.016	-2.381	True
132	clay	D	20	20	60	0.800	0.647	0.449	0.252	0.197	207.818	3.409	0.008	0.000	0.151	0.016	1.291	6.239	-1.596	False
133	clay	D	20	20	60	0.900	0.619	0.442	0.248	0.194	130.710	2.144	0.007	0.000	0.147	0.014	1.288	4.412	-1.453	False
134	clay	D	20	20	60	1.000	0.591	0.432	0.244	0.188	82.545	1.354	0.006	0.000	0.143	0.013	1.285	3.219	-1.393	False
135	clay	D	20	20	60	1.100	0.563	0.418	0.239	0.179	52.323	0.858	0.005	0.000	0.141	0.012	1.280	2.469	-1.391	False
136	clay	D	20	20	60	1.200	0.536	0.404	0.236	0.168	33.173	0.544	0.004	0.000	0.139	0.012	1.273	2.010	-1.427	False
137	clay	D	20	20	60	1.300	0.509	0.390	0.233	0.157	20.956	0.344	0.004	0.000	0.137	0.011	1.264	1.717	-1.501	False
138	clay	D	20	20	60	1.400	0.481	0.375	0.231	0.144	13.182	0.216	0.003	0.000	0.137	0.011	1.252	1.507	-1.628	True
139	clay	D	20	20	60	1.500	0.453	0.360	0.230	0.131	8.264	0.136	0.003	0.000	0.136	0.011	1.238	1.350	-1.831	True
140	clay	D	20	20	60	1.600	0.425	0.345	0.230	0.115	5.162	0.085	0.003	0.000	0.136	0.011	1.219	1.248	-2.163	True
141	clay	D	20	20	60	1.700	0.395	0.330	0.230	0.099	3.240	0.053	0.002	0.000	0.138	0.011	1.198	1.240	-2.681	True
142	clay	D	20	20	60	1.800	0.366	0.313	0.230	0.083	2.247	0.037	0.003	0.000	0.140	0.012	1.178	1.417	-3.430	True
143	clay	D	20	20	60	1.900	0.338	0.295	0.228	0.067	2.749	0.045	0.003	0.000	0.144	0.012	1.159	1.973	-4.424	True

1. Mineral-baseline organic carbon (from UNSODA 2.0)

NoteFor researchers: ROSETTA mineral-baseline organic carbon (UNSODA)

ROSETTA has no OC input, but its training samples (UNSODA + others) are mineral-dominated soils that do carry organic carbon — so the ROSETTA prediction is a nominal baseline at OC > 0, not a true organic-free (OC = 0) soil. Before blending in any organic-matter effect (Section 2), we estimate that baseline OC from the 367 UNSODA 2.0 samples that report both bulk density and organic-matter content (OC = 0.58·OM, van Bemmelen). OC falls clearly with bulk density (Pearson r ≈ −0.6); the scatter below overlays OLS regressions of OC on bulk density fit separately for topsoil, subsoil, and all mineral data, with the slopes/intercepts/R²/p reported in the accompanying table. For the mineral subset (OM ≤ 20 %): all-horizon mean ≈ 0.9 % OC; mineral topsoil (≤15 cm) median ≈ 1.1 % OC. Rather than a single anchor, we use the all-mineral fit as a BD-dependent baseline — OC_base(BD) = OC_BD_SLOPE·BD + OC_BD_INTERCEPT (floored at 0) — so ROSETTA’s prediction carries the OC mineral soils typically have at that bulk density (≈ 2.9 % at BD 0.8, ≈ 1.0 % at the mean BD ≈ 1.4, 0 by BD ≈ 1.72). Section 2 applies the Minasny & McBratney increments relative to that baseline. (OC_BASELINE_PCT ≈ 1 % is retained only as the reference value at the mean BD.)

data_temp/ is git-ignored; run pixi run python notebooks/fetch_unsoda.py to (re)create the UNSODA extract read below.

import os

# Reference mineral-baseline OC at the mean mineral BD ≈ 1.4 (UNSODA mineral subset: all-horizon
# mean ~0.9 %, topsoil median ~1.1 %). The blend (Section 2) uses the BD-DEPENDENT baseline
# oc_baseline_for_bd(BD) from the all-mineral OC~BD fit below, not this scalar.
OC_BASELINE_PCT = 1.0  # % organic carbon, reference only (value of the all-mineral fit at BD ≈ 1.4)

_unsoda_path = "data_temp/unsoda_bd_om.csv"
if os.path.exists(_unsoda_path):
    unsoda = pd.read_csv(_unsoda_path)
    mineral = unsoda[unsoda["is_mineral"]].copy()
    mineral["horizon_group"] = np.where(
        mineral["depth_upper"] <= 15, "topsoil (≤15 cm)", "subsoil (>15 cm)"
    )
    print("UNSODA mineral subset (OM ≤ 20%) — organic carbon %, OC = 0.58·OM:")
    print(mineral.groupby("horizon_group")["OC_pct"].agg(["count", "mean", "median"]).round(2))
    print(f"OC~BD Pearson r = {mineral['OC_pct'].corr(mineral['bulk_density_g_cm3']):.2f}")

    # OLS regressions of organic carbon on bulk density: topsoil, subsoil, and all mineral data.
    _reg_groups = {
        "topsoil (≤15 cm)": mineral[mineral["horizon_group"] == "topsoil (≤15 cm)"],
        "subsoil (>15 cm)": mineral[mineral["horizon_group"] == "subsoil (>15 cm)"],
        "all mineral": mineral,
    }
    _reg_colors = {"topsoil (≤15 cm)": "#1f77b4", "subsoil (>15 cm)": "#ff7f0e", "all mineral": "black"}
    _bd_line = np.array([mineral["bulk_density_g_cm3"].min(), mineral["bulk_density_g_cm3"].max()])
    _reg_rows, _reg_lines, _all_mineral_lr = [], [], None
    for _name, _g in _reg_groups.items():
        _lr = linregress(_g["bulk_density_g_cm3"], _g["OC_pct"])
        if _name == "all mineral":
            _all_mineral_lr = _lr
        _reg_rows.append({
            "regression": _name,
            "n": len(_g),
            "slope (%OC per g/cm³)": _lr.slope,
            "intercept (%OC)": _lr.intercept,
            "r": _lr.rvalue,
            "R²": _lr.rvalue ** 2,
            "p_value": _lr.pvalue,
        })
        _reg_lines.append(
            hv.Curve((_bd_line, _lr.intercept + _lr.slope * _bd_line), label=f"{_name} fit").opts(
                color=_reg_colors[_name], line_width=2,
                line_dash="solid" if _name == "all mineral" else "dashed",
            )
        )
    reg_df = pd.DataFrame(_reg_rows)

    _scatter = mineral.hvplot.scatter(
        x="bulk_density_g_cm3", y="OC_pct", by="horizon_group",
        xlabel="bulk density (g/cm³)", ylabel="organic carbon  OC = 0.58·OM  (%)",
        title="UNSODA 2.0: organic carbon vs. bulk density (mineral soils)",
        width=760, height=460, legend="top_right", alpha=0.6, size=25, ylim=(0, 10), grid=True,
    )
    # The all-mineral fit (black solid) IS the BD-dependent ROSETTA baseline used by the blend.
    # Verify the live UNSODA fit still matches the constants hard-coded in _helpers (the home page
    # relies on those, since it never runs this regression); fail loudly if UNSODA has drifted.
    assert np.isclose(_all_mineral_lr.slope, OC_BD_SLOPE, atol=1e-2) and \
        np.isclose(_all_mineral_lr.intercept, OC_BD_INTERCEPT, atol=1e-2), (
        f"all-mineral OC~BD fit ({_all_mineral_lr.slope:.4f}, {_all_mineral_lr.intercept:.4f}) "
        f"drifted from _helpers (OC_BD_SLOPE={OC_BD_SLOPE}, OC_BD_INTERCEPT={OC_BD_INTERCEPT}); "
        f"update the constants in notebooks/_helpers.py.")
    display(_scatter * hv.Overlay(_reg_lines))

    # Regression parameters table (p-values formatted in scientific notation to keep precision).
    _reg_show = reg_df.copy()
    _reg_show["p_value"] = _reg_show["p_value"].map(lambda p: f"{p:.2e}")
    print("OC = slope · bulk_density + intercept  (OLS):")
    display(show(_reg_show, height=160))
else:
    print(f"{_unsoda_path} not found — run `pixi run python notebooks/fetch_unsoda.py` to regenerate it.")
    print(f"Proceeding with the documented default OC_BASELINE_PCT = {OC_BASELINE_PCT} % OC.")

UNSODA mineral subset (OM ≤ 20%) — organic carbon %, OC = 0.58·OM:
                  count  mean  median
horizon_group                        
subsoil (>15 cm)    240 0.660   0.310
topsoil (≤15 cm)    119 1.460   1.130
OC~BD Pearson r = -0.63

OC = slope · bulk_density + intercept  (OLS):

	regression	n	slope (%OC per g/cm³)	intercept (%OC)	r	R²	p_value
0	topsoil (≤15 cm)	119	-3.337	5.959	-0.703	0.494	5.06e-19
1	subsoil (>15 cm)	240	-2.668	4.561	-0.548	0.300	3.53e-20
2	all mineral	359	-3.160	5.427	-0.634	0.402	9.28e-42

2. ROSETTA + organic-matter modifier (Minasny & McBratney 2018)

The interactive diagram and line plots below show how available and drainable water change as you add organic matter and adjust compaction across all 12 USDA texture classes. Use the bulk density slider to represent compaction (higher = more compacted) and the organic matter slider to explore realistic management scenarios. The low-BD + high-OM corner represents a healthy, well-structured soil; the high-BD + low-OM corner a compacted, depleted one. Greyed texture columns mark physically implausible BD × OM combinations.

NoteFor researchers: Minasny & McBratney slopes, blend method, and caveats

The approach keeps ROSETTA’s texture + bulk-density skill for the mineral soil baseline (Section 1), then adds the empirical organic-carbon increments from Minasny & McBratney (2018) Table 2 — an OC sensitivity derived from >50,000 measurements and preferred here over Saxton–Rawls (see §3.1). Per +1 % organic carbon (= +10 g C kg⁻¹), by USDA texture group:

Minasny & McBratney group	ΔWP	ΔAWC	ΔSAT (mm 100 mm⁻¹ per 1 % OC)
Coarse	0.86	1.94	4.59
Medium	0.68	1.79	3.59
Fine	0.54	1.41	3.23

Blend (volumetric, cm³/cm³); ROSETTA gives the baseline at the mineral bulk density set by the BD slider, anchored at the BD-dependent baseline OC, OC_base(BD), from the UNSODA all-mineral OC~BD fit (Section 1; ≈ 2.9 % at BD 0.8, ≈ 1.0 % at BD 1.4, 0 by BD ≈ 1.72):

• WP(OC) = WP_ROSETTA + (ΔWP/100)·(OC − OC_base)
• AWC(OC) = AWC_ROSETTA + (ΔAWC/100)·(OC − OC_base) (AWC_ROSETTA = FC − WP from ROSETTA)
• SAT(OC) = θₛ_ROSETTA + (ΔSAT/100)·(OC − OC_base)
• FC(OC) = WP + AWC ; drainable = SAT − FC

We apply Minasny & McBratney’s WP, AWC and SAT slopes — the three quantities that define the unavailable / available / drainable bands — and derive FC = WP + AWC. This reproduces Minasny & McBratney’s headline AWC sensitivity exactly; the drainable response follows from ΔSAT − ΔFC. (Because Minasny & McBratney regressed each property independently, ΔAWC ≠ ΔFC − ΔWP; anchoring on ΔFC instead would understate the AWC response by ~25–50 %, especially in fine soils.)

Caveats. (1) ROSETTA’s prediction is a nominal baseline at the BD-dependent baseline OC, OC_base(BD) (all-mineral fit, Section 1 — ≈ 1 % near BD 1.4, higher at low BD, 0 by BD ≈ 1.72), not OC = 0; the Minasny & McBratney increments are applied relative to it, so the slider’s 0 % end is a truly organic-free mineral soil (drier than ROSETTA), clamped at ≥ 0. The firebrick reference line in the line plots marks OC_base(BD) and moves with the BD slider. (2) The BD and OC sliders are independent “what-if” axes; in reality organic matter lowers bulk density (the low-BD ↔︎ high-OC diagonal is the realistic region), and a low-BD + high-OC corner double-counts porosity, so don’t read the extreme corners as coupled predictions. (3) The modifier is linear, whereas Minasny & McBratney found diminishing returns (largest gains 0→1 % OC), so it may overstate gains at high OC; their data span OC < 10 %. (4) OM ≈ OC / 0.58 (van Bemmelen); the line-plot OC axis is capped at 5 % (≈ 8.6 % OM) and the diagram’s OM slider spans 0–8 %. (5) In the diagram, texture columns are greyed where the blended saturation exceeds the BD-implied pore space (1 − BD/2.65) — physically impossible, i.e. extrapolation at that BD × OM (mirrors Notebook 1’s implausible_bd flag).

# ROSETTA mineral baseline (per bulk density) + additive Minasny & McBratney (2018) OC
# increments, applied RELATIVE to the BD-dependent baseline oc_baseline_for_bd(BD), across the full
# bulk-density range. VB, MM_SLOPES, MM_GROUP, oc_baseline_for_bd are imported from _helpers.

oc_values = np.round(np.arange(0.0, 8.0 + 1e-9, 0.5), 2)  # organic carbon %, 0–8 (≈ 0–14% OM)

base = result.set_index(["texture_class", "bulk_density_g_cm3"])
blend_rows = []
for bd in bulk_densities:
    oc_base_bd = float(oc_baseline_for_bd(bd))   # BD-dependent mineral-baseline OC (all-mineral fit)
    for oc in oc_values:
        for cls in TEXTURE_CLASSES:
            s = MM_SLOPES[MM_GROUP[cls]]
            sat0 = base.loc[(cls, bd), "total_porosity"]
            fc0 = base.loc[(cls, bd), "field_capacity_porosity"]
            wp0 = base.loc[(cls, bd), "wilting_point_porosity"]
            awc0 = fc0 - wp0
            d_oc = oc - oc_base_bd               # increments relative to the BD-dependent baseline OC
            wp = max(wp0 + s["WP"] / 100 * d_oc, 0.0)
            awc = max(awc0 + s["AWC"] / 100 * d_oc, 0.0)   # Minasny & McBratney AWC slope; floor at 0
            sat = max(sat0 + s["SAT"] / 100 * d_oc, wp + awc)  # keep SAT >= FC
            fc = wp + awc                       # derive FC so AWC matches Minasny & McBratney exactly
            blend_rows.append(
                {
                    "texture_class": cls,
                    "hydrologic_soil_group": HYDROLOGIC_SOIL_GROUP[cls],
                    "mm_group": MM_GROUP[cls],
                    "bulk_density_g_cm3": bd,
                    "oc_pct": oc,
                    "om_pct_approx": round(oc / VB, 2),
                    "wilting_point_porosity": wp,
                    "field_capacity_porosity": fc,
                    "total_porosity": sat,
                    "available_water_capacity": awc,
                    "drainable_water": sat - fc,
                }
            )
blend_df = pd.DataFrame(blend_rows)
blend_df["texture_class"] = pd.Categorical(blend_df["texture_class"], categories=list(TEXTURE_CLASSES), ordered=True)
# order each (BD, texture) series ascending in OC, then flag extrapolation: a blended saturation
# exceeding the BD-implied pore space (1 - BD/2.65) is physically impossible (greyed in the plots).
blend_df = blend_df.sort_values(["bulk_density_g_cm3", "texture_class", "oc_pct"]).reset_index(drop=True)
blend_df["implausible"] = blend_df["total_porosity"] > (1.0 - blend_df["bulk_density_g_cm3"] / 2.65)
_blend_grp = blend_df.groupby(["bulk_density_g_cm3", "texture_class"], observed=True, sort=False)["implausible"]
_blend_extrap_mask = blend_df["implausible"] | _blend_grp.shift(-1).fillna(False)  # keep boundary so segments join

print(f"{len(blend_df)} rows  ({len(TEXTURE_CLASSES)} textures x {len(bulk_densities)} BD x {len(oc_values)} OC)")
show(blend_df[(blend_df["bulk_density_g_cm3"] == 1.5) & (blend_df["oc_pct"].isin([0.0, 2.0, 4.0]))].round(3))

2448 rows  (12 textures x 12 BD x 17 OC)

	texture_class	hydrologic_soil_group	mm_group	bulk_density_g_cm3	oc_pct	om_pct_approx	wilting_point_porosity	field_capacity_porosity	total_porosity	available_water_capacity	drainable_water	implausible
1428	sand	A	coarse	1.500	0.000	0.000	0.046	0.046	0.345	0.000	0.299	False
1432	sand	A	coarse	1.500	2.000	3.450	0.063	0.096	0.437	0.033	0.341	True
1436	sand	A	coarse	1.500	4.000	6.900	0.080	0.152	0.529	0.071	0.377	True
1445	loamy sand	A	coarse	1.500	0.000	0.000	0.054	0.097	0.345	0.043	0.248	False
1449	loamy sand	A	coarse	1.500	2.000	3.450	0.071	0.153	0.437	0.082	0.284	True
1453	loamy sand	A	coarse	1.500	4.000	6.900	0.088	0.209	0.529	0.121	0.319	True
1462	sandy loam	A	coarse	1.500	0.000	0.000	0.081	0.181	0.342	0.101	0.161	False
1466	sandy loam	A	coarse	1.500	2.000	3.450	0.098	0.237	0.434	0.139	0.197	True
1470	sandy loam	A	coarse	1.500	4.000	6.900	0.115	0.293	0.526	0.178	0.233	True
1479	loam	B	medium	1.500	0.000	0.000	0.125	0.258	0.356	0.133	0.099	False
1483	loam	B	medium	1.500	2.000	3.450	0.138	0.307	0.428	0.169	0.121	False
1487	loam	B	medium	1.500	4.000	6.900	0.152	0.356	0.500	0.204	0.143	True
1496	silt loam	B	medium	1.500	0.000	0.000	0.109	0.274	0.360	0.166	0.085	False
1500	silt loam	B	medium	1.500	2.000	3.450	0.122	0.324	0.431	0.201	0.108	False
1504	silt loam	B	medium	1.500	4.000	6.900	0.136	0.373	0.503	0.237	0.130	True
1513	silt	B/D	medium	1.500	0.000	0.000	0.081	0.257	0.373	0.176	0.116	False
1517	silt	B/D	medium	1.500	2.000	3.450	0.095	0.307	0.445	0.212	0.138	True
1521	silt	B/D	medium	1.500	4.000	6.900	0.108	0.356	0.517	0.248	0.161	True
1530	sandy clay loam	C	coarse	1.500	0.000	0.000	0.139	0.254	0.376	0.115	0.122	False
1534	sandy clay loam	C	coarse	1.500	2.000	3.450	0.157	0.310	0.468	0.153	0.158	True
1538	sandy clay loam	C	coarse	1.500	4.000	6.900	0.174	0.366	0.560	0.192	0.194	True
1547	clay loam	D	medium	1.500	0.000	0.000	0.170	0.300	0.383	0.130	0.083	False
1551	clay loam	D	medium	1.500	2.000	3.450	0.183	0.349	0.455	0.166	0.106	True
1555	clay loam	D	medium	1.500	4.000	6.900	0.197	0.399	0.527	0.202	0.128	True
1564	silty clay loam	D	medium	1.500	0.000	0.000	0.167	0.319	0.389	0.152	0.070	False
1568	silty clay loam	D	medium	1.500	2.000	3.450	0.180	0.368	0.460	0.188	0.092	True
1572	silty clay loam	D	medium	1.500	4.000	6.900	0.194	0.418	0.532	0.223	0.115	True
1581	sandy clay	D	fine	1.500	0.000	0.000	0.190	0.307	0.414	0.118	0.106	False
1585	sandy clay	D	fine	1.500	2.000	3.450	0.201	0.346	0.478	0.146	0.132	True
1589	sandy clay	D	fine	1.500	4.000	6.900	0.211	0.385	0.543	0.174	0.158	True
1598	silty clay	D	fine	1.500	0.000	0.000	0.198	0.339	0.411	0.141	0.072	False
1602	silty clay	D	fine	1.500	2.000	3.450	0.208	0.378	0.475	0.169	0.098	True
1606	silty clay	D	fine	1.500	4.000	6.900	0.219	0.417	0.540	0.198	0.123	True
1615	clay	D	fine	1.500	0.000	0.000	0.226	0.347	0.431	0.121	0.084	False
1619	clay	D	fine	1.500	2.000	3.450	0.237	0.386	0.496	0.149	0.110	True
1623	clay	D	fine	1.500	4.000	6.900	0.248	0.425	0.560	0.177	0.136	True

# AVAILABLE water capacity vs. organic carbon, one line per texture class, with a BD slider.
# Each line is solid where plausible and grey-dashed where the (BD, OC) state is an extrapolation
# (blended saturation > 1 - BD/2.65); the greyed region grows as the BD slider increases.
hv.output(widget_location="bottom")

om8_oc = 8.0 * VB  # 8% OM ≈ 4.64% OC marks the top of the primary range


def blend_line(ycol, ylabel, title, ylim):
    common = dict(
        x="oc_pct", y=ycol, by="texture_class", groupby="bulk_density_g_cm3",
        dynamic=False, width=820, height=500, xlim=(0, 5), ylim=ylim,
    )
    solid = blend_df.assign(**{ycol: blend_df[ycol].where(~blend_df["implausible"])}).hvplot.line(
        xlabel="soil organic carbon (% by weight)   [OM ≈ OC / 0.58]",
        ylabel=ylabel, title=title, legend="right", grid=True, **common,
    )
    dashed = (
        blend_df.assign(**{ycol: blend_df[ycol].where(_blend_extrap_mask)})
        .hvplot.line(**common)
        .opts(hv.opts.Curve(color="lightgray", line_dash="dashed", alpha=0.9))
        .opts(show_legend=False)
    )
    y_lab = ylim[1]  # labels just below the top of the plot

    # The baseline-OC reference is now BD-dependent (all-mineral fit), so it must move with the BD
    # slider: build one ref frame per bulk density as a HoloMap keyed on the same dimension as the
    # solid/dashed line HoloMaps, so they overlay frame-by-frame.
    def _refs_for_bd(bd):
        ocb = float(oc_baseline_for_bd(bd))
        return (
            hv.VLine(ocb).opts(color="firebrick", line_dash="dashed", line_width=1)   # BD-dependent baseline OC
            * hv.VLine(om8_oc).opts(color="black", line_dash="dotted", line_width=1)   # 8% OM
            * hv.Text(ocb, y_lab, f" ROSETTA baseline OC ≈ {ocb:.1f}%", halign="left", valign="top").opts(
                text_color="firebrick", text_font_size="8pt")
            * hv.Text(om8_oc, y_lab, "8% OM ", halign="right", valign="top").opts(
                text_color="black", text_font_size="8pt")
        )

    refs = hv.HoloMap({bd: _refs_for_bd(bd) for bd in bulk_densities}, kdims="bulk_density_g_cm3")
    return (solid * dashed * refs).redim(
        bulk_density_g_cm3=hv.Dimension("Bulk density, g/cm³ (higher is more compacted)", default=1.4, value_format=lambda v: f"{v:.1f}")
    )


blend_line(
    "available_water_capacity",
    "available water capacity (cm³/cm³)",
    "ROSETTA + Minasny & McBratney blend: AVAILABLE water vs. organic carbon — {dimensions}",
    (0, 0.42),
)

Takeaway: In well-managed, low-compaction soils, even modest organic matter gains (1–2 % OC) measurably increase plant-available water — most in sandy and loamy soils, least in heavy clays.

# DRAINABLE water (saturation − field capacity) vs. organic carbon — the rapidly draining pore
# space that matters for stormwater storage / infiltration. BD slider; extrapolated (BD, OC)
# states are grey-dashed (same flag as the AVAILABLE-water plot above).
hv.output(widget_location="bottom")

blend_line(
    "drainable_water",
    "drainable water  SAT − FC  (cm³/cm³)",
    "ROSETTA + Minasny & McBratney blend: DRAINABLE water vs. organic carbon — {dimensions}",
    (0, 0.60),
)

Takeaway: Organic matter increases the fast-draining (macro)pore space most strongly in coarse soils — the same soils that most benefit for stormwater infiltration — while heavy clays see smaller and less consistent gains.

# FAO-style transposed diagram for the ROSETTA + Minasny & McBratney blend, with TWO sliders: mineral bulk
# density and organic MATTER (the way the audience thinks about it; OM ≈ OC / 0.58). Computed
# directly from the ROSETTA baseline + Minasny & McBratney increments so the slider can carry round OM values.
# dynamic=False embeds every (BD, OM) frame so it works without a live kernel. (OM capped at
# 8% on a 1% grid to bound the frame count / latency.)
import panel as pn

pn.extension()

om_grid = [round(float(o), 1) for o in np.arange(0.0, 8.0 + 1e-9, 1.0)]  # organic matter %, 0–8

# One shared pair of sliders drives BOTH the figure and the table; placed BETWEEN them. The Panel
# layout is embedded (embed=True) so every (BD, OM) state works in static HTML without a kernel.
# The styled HTML table (soil_water_table_html) gives bold/wrapped headers, a hidden index, full
# texture names on one line, and all rows — formatting a Bokeh DataTable cannot express.
_SLIDER_CSS = [":host, :host * { font-size: 0.85rem; }"]  # sized so the BD slider title stays on one line
bd_slider = pn.widgets.DiscreteSlider(
    name="Bulk density, g/cm³ (higher is more compacted)",
    options=[round(float(b), 1) for b in bulk_densities], value=1.4, width=380, stylesheets=_SLIDER_CSS)
om_slider = pn.widgets.DiscreteSlider(
    name="Organic matter (% by weight)", options=om_grid, value=2.0, width=380, stylesheets=_SLIDER_CSS)


def _blend_figure(bd, om):
    t = soil_water_bd_om_blend_table(result, bd, om)
    x = np.array([texture_x[cls] for cls in TEXTURE_CLASSES])
    pwp = t["wilting_point_porosity"].to_numpy() * INCHES_PER_FOOT
    fc = t["field_capacity_porosity"].to_numpy() * INCHES_PER_FOOT
    por = t["total_porosity"].to_numpy() * INCHES_PER_FOOT
    ov = soil_water_texture_band_diagram(
        x, pwp, fc, por,
        texture_labels=t["texture_class"].astype(str).tolist(),
        implausible=t["implausible"].to_numpy(),
    )
    return ov.opts(
        hv.opts.Overlay(width=720, height=520, toolbar="right", legend_position="top_left",
                        xlabel="Texture class (Hydrologic Soil Group); coarse → fine",
                        ylabel="Water Storage Capacity (inches per foot of soil depth)",
                        title=f"Soil water vs. texture for bulk density = {bd:.1f} g/cm³ & organic matter = {om:g}%"),
        hv.opts.Curve(xticks=texture_ticks, xrotation=45, ylim=(0, 10)),
        hv.opts.Area(xticks=texture_ticks, xrotation=45, ylim=(0, 10)),
    )


def _blend_table(bd, om):
    t = soil_water_bd_om_blend_table(result, bd, om)
    out = pd.DataFrame({
        "texture (HSG)": [f"{c} ({HYDROLOGIC_SOIL_GROUP[c]})" for c in t["texture_class"]],
        "wilting point": t["wilting_point_porosity"].to_numpy(),
        "field capacity": t["field_capacity_porosity"].to_numpy(),
        "saturation": t["total_porosity"].to_numpy(),
        "available water": t["available_water_capacity"].to_numpy(),
        "drainable water": t["drainable_water"].to_numpy(),
    }).round(3)
    return pn.pane.HTML(soil_water_table_html(out), width=720)


blend_layout = pn.Column(
    pn.panel(pn.bind(_blend_figure, bd_slider, om_slider)),
    pn.Row(bd_slider, om_slider),
    pn.pane.HTML(
        "<b>Soil water retention by texture class (porosity volume fraction, cm³/cm³)</b>",
        styles={"font-size": "1rem", "margin": "8px 0 2px"},
    ),
    pn.panel(pn.bind(_blend_table, bd_slider, om_slider)),
    width=790,
)

# embed() bakes every (BD, OM) state into static HTML so the sliders work with no kernel, but it
# SILENTLY truncates if its limits are exceeded (→ dead/stale frames). Assert the bounds so a
# future finer grid fails loudly here instead (max_opts > options per slider; max_states > product).
assert len(bulk_densities) <= 40 and len(om_grid) <= 40, "max_opts=40 too low for slider options"
assert len(bulk_densities) * len(om_grid) <= 2000, "max_states=2000 too low for BD×OM combos"
blend_layout.embed(max_states=2000, max_opts=40, progress=False)

Takeaway: Adding organic matter and easing compaction together shift water into both plant-available and drainable storage — gains are strongest in coarse-textured soils and diminish toward clays. The realistic management path runs along the low-BD + high-OM diagonal.

3. Organic-matter sensitivity comparison to Saxton & Rawls (2006)

The “Soil Water Characteristic Estimates by Texture and Organic Matter for Hydrologic Solutions” publication by Saxton & Rawls (2006) is an earlier study that has been commonly used to estimate how changes in soil organic matter affect soil water storage. Here we compare findings from this older study with the newer results from Minasny & McBratney (2018) that we use in section 2 above.

NoteFor researchers: Saxton–Rawls comparison

The Saxton & Rawls (2006) pedotransfer functions take sand, clay, and organic-matter % directly and were developed from USDA/NRCS data for the continental USA. Self-contained (no ROSETTA baseline) — though (see §3.1) it gives a smaller, and for clays negative, OC effect than Minasny & McBratney.

Restricted to its calibrated range, OM ≤ 8 % by weight (≈ 4.6 % organic carbon).

def saxton_rawls(sand_frac, clay_frac, om_pct):
    """Saxton & Rawls (2006), SSSAJ 70:1569-1578 — soil-water characteristics from texture and
    organic matter. Inputs: sand & clay as decimal mass fractions (0-1), organic matter in % by
    weight. Returns volumetric water contents (cm³/cm³): theta_1500 (permanent wilting point),
    theta_33 (field capacity), theta_S (total porosity / saturation). Calibrated for OM up to
    ~8 %; higher OM is extrapolation. Scalars or numpy arrays (broadcast together).
    """
    S, C, OM = sand_frac, clay_frac, om_pct
    t1500t = -0.024 * S + 0.487 * C + 0.006 * OM + 0.005 * (S * OM) - 0.013 * (C * OM) + 0.068 * (S * C) + 0.031
    t1500 = t1500t + (0.14 * t1500t - 0.02)
    t33t = -0.251 * S + 0.195 * C + 0.011 * OM + 0.006 * (S * OM) - 0.027 * (C * OM) + 0.452 * (S * C) + 0.299
    t33 = t33t + (1.283 * t33t**2 - 0.374 * t33t - 0.015)
    tS33t = 0.278 * S + 0.034 * C + 0.022 * OM - 0.018 * (S * OM) - 0.027 * (C * OM) - 0.584 * (S * C) + 0.078
    tS33 = tS33t + (0.636 * tS33t - 0.107)
    tS = t33 + tS33 - 0.097 * S + 0.043
    return t1500, t33, tS


OM_VALID_MAX = 8.0  # Saxton–Rawls organic-matter calibration limit (% by weight)
om_values = np.round(np.arange(0.0, OM_VALID_MAX + 1e-9, 1.0), 1)  # 0–8 %, 1 % steps (calibrated range)

sr_rows = []
for om in om_values:
    for cls, (sand, silt, clay) in TEXTURE_CLASSES.items():
        pwp, fc, por = saxton_rawls(sand / 100, clay / 100, om)
        sr_rows.append(
            {
                "texture_class": cls,
                "hydrologic_soil_group": HYDROLOGIC_SOIL_GROUP[cls],
                "om_pct": om,
                "wilting_point_porosity": pwp,
                "field_capacity_porosity": fc,
                "total_porosity": por,
                "available_water_capacity": fc - pwp,
                "drainable_water": por - fc,
            }
        )
sr_df = pd.DataFrame(sr_rows)
sr_df["texture_class"] = pd.Categorical(sr_df["texture_class"], categories=list(TEXTURE_CLASSES), ordered=True)

print(f"{len(sr_df)} rows  ({len(TEXTURE_CLASSES)} textures x {len(om_values)} OM levels)")
show(sr_df[sr_df["om_pct"].isin([0.0, 2.0, 4.0, 8.0])])

108 rows  (12 textures x 9 OM levels)

	texture_class	hydrologic_soil_group	om_pct	wilting_point_porosity	field_capacity_porosity	total_porosity	available_water_capacity	drainable_water
0	sand	A	0.000	0.009	0.049	0.417	0.040	0.368
1	loamy sand	A	0.000	0.030	0.085	0.399	0.055	0.313
2	sandy loam	A	0.000	0.058	0.144	0.384	0.086	0.240
3	loam	B	0.000	0.122	0.253	0.394	0.131	0.141
4	silt loam	B	0.000	0.095	0.277	0.392	0.181	0.115
5	silt	B/D	0.000	0.041	0.278	0.366	0.237	0.088
6	sandy clay loam	C	0.000	0.161	0.253	0.392	0.092	0.139
7	clay loam	D	0.000	0.210	0.345	0.435	0.136	0.090
8	silty clay loam	D	0.000	0.204	0.370	0.456	0.166	0.086
9	sandy clay	D	0.000	0.257	0.368	0.429	0.111	0.061
10	silty clay	D	0.000	0.271	0.417	0.501	0.146	0.083
11	clay	D	0.000	0.352	0.474	0.528	0.122	0.054
24	sand	A	2.000	0.032	0.077	0.460	0.044	0.383
25	loamy sand	A	2.000	0.051	0.114	0.445	0.062	0.332
26	sandy loam	A	2.000	0.076	0.172	0.437	0.096	0.265
27	loam	B	2.000	0.134	0.274	0.446	0.140	0.172
28	silt loam	B	2.000	0.107	0.299	0.461	0.192	0.162
29	silt	B/D	2.000	0.054	0.306	0.458	0.252	0.152
30	sandy clay loam	C	2.000	0.174	0.273	0.424	0.099	0.151
31	clay loam	D	2.000	0.216	0.355	0.469	0.139	0.113
32	silty clay loam	D	2.000	0.209	0.377	0.499	0.169	0.122
33	sandy clay	D	2.000	0.264	0.376	0.441	0.112	0.066
34	silty clay	D	2.000	0.272	0.414	0.525	0.142	0.111
35	clay	D	2.000	0.350	0.461	0.522	0.110	0.061
48	sand	A	4.000	0.056	0.107	0.505	0.051	0.398
49	loamy sand	A	4.000	0.073	0.144	0.494	0.071	0.350
50	sandy loam	A	4.000	0.094	0.201	0.491	0.107	0.289
51	loam	B	4.000	0.146	0.296	0.499	0.150	0.203
52	silt loam	B	4.000	0.118	0.323	0.532	0.204	0.209
53	silt	B/D	4.000	0.067	0.336	0.551	0.268	0.215
54	sandy clay loam	C	4.000	0.186	0.293	0.457	0.107	0.164
55	clay loam	D	4.000	0.223	0.366	0.502	0.143	0.137
56	silty clay loam	D	4.000	0.213	0.385	0.543	0.171	0.158
57	sandy clay	D	4.000	0.270	0.383	0.453	0.113	0.070
58	silty clay	D	4.000	0.273	0.411	0.549	0.138	0.138
59	clay	D	4.000	0.349	0.447	0.516	0.099	0.068
96	sand	A	8.000	0.102	0.175	0.604	0.073	0.429
97	loamy sand	A	8.000	0.115	0.211	0.598	0.096	0.387
98	sandy loam	A	8.000	0.131	0.264	0.603	0.133	0.339
99	loam	B	8.000	0.171	0.343	0.607	0.172	0.264
100	silt loam	B	8.000	0.142	0.372	0.674	0.230	0.303
101	silt	B/D	8.000	0.093	0.398	0.739	0.304	0.342
102	sandy clay loam	C	8.000	0.211	0.335	0.525	0.124	0.190
103	clay loam	D	8.000	0.236	0.386	0.570	0.150	0.184
104	silty clay loam	D	8.000	0.223	0.400	0.630	0.177	0.230
105	sandy clay	D	8.000	0.284	0.398	0.477	0.114	0.079
106	silty clay	D	8.000	0.275	0.404	0.597	0.130	0.192
107	clay	D	8.000	0.345	0.421	0.504	0.076	0.083

# Saxton–Rawls AVAILABLE and DRAINABLE water vs. organic matter, one line per texture class,
# over the calibrated 0–8 % OM range.
hv.output(widget_location="right")

sr_awc = sr_df.hvplot.line(
    x="om_pct", y="available_water_capacity", by="texture_class",
    xlabel="soil organic matter (% by weight)",
    ylabel="available water capacity (cm³/cm³)",
    title="Saxton–Rawls: AVAILABLE water vs. organic matter",
    width=820, height=460, legend="right", grid=True, ylim=(0, 0.22),
)
sr_drain = sr_df.hvplot.line(
    x="om_pct", y="drainable_water", by="texture_class",
    xlabel="soil organic matter (% by weight)",
    ylabel="drainable water  SAT − FC  (cm³/cm³)",
    title="Saxton–Rawls: DRAINABLE water vs. organic matter",
    width=820, height=460, legend="right", grid=True, ylim=(0, 0.45),
)
(sr_awc + sr_drain).cols(1)

# FAO-style transposed diagram, Saxton–Rawls, with an organic-matter slider (calibrated 0–8 % OM).
hv.output(widget_location="bottom")


def _sr_profile(om):
    d = (
        sr_df[sr_df["om_pct"] == om]
        .set_index("texture_class")
        .reindex(list(TEXTURE_CLASSES))
        .reset_index()
    )
    x = d["texture_class"].map(texture_x).to_numpy()
    pwp = d["wilting_point_porosity"].to_numpy() * INCHES_PER_FOOT
    fc = d["field_capacity_porosity"].to_numpy() * INCHES_PER_FOOT
    por = d["total_porosity"].to_numpy() * INCHES_PER_FOOT
    return soil_water_texture_band_diagram(
        x, pwp, fc, por, texture_labels=d["texture_class"].astype(str).tolist()
    )


sr_profiles = hv.HoloMap(
    {om: _sr_profile(om) for om in om_values},
    kdims=[hv.Dimension("Soil organic matter (% by weight)", default=2.0, value_format=lambda v: f"{v:.1f}")],
)

sr_profiles.opts(
    hv.opts.Overlay(
        width=820,
        height=520,
        legend_position="top_left",
        xlabel="texture class (hydrologic soil group); coarse → fine",
        ylabel="Water Storage Capacity (inches per foot of soil depth)",
        title="Soil water vs. texture (Saxton–Rawls) — {dimensions}",
    ),
    hv.opts.Curve(xticks=texture_ticks, xrotation=45, ylim=(0, 10)),
    hv.opts.Area(xticks=texture_ticks, xrotation=45, ylim=(0, 10)),
)

3.1 Validation: ΔAWC/ΔOC vs. Minasny & McBratney (2018)

NoteFor researchers: Saxton–Rawls vs. Minasny & McBratney AWC sensitivity comparison

Do the two PTF families agree on how much available water organic carbon adds? We compute the Saxton–Rawls ΔAWC per +1 % organic carbon for each texture class — over the same OC 0.5 % → 1.5 % interval Minasny & McBratney used for PTF-derived slopes (OC = 0.58·OM) — then average by their coarse/medium/fine groups and compare against Table 2.

Expect only order-of-magnitude agreement: both say the effect is small and decreases from coarse to fine textures, but Saxton–Rawls is systematically lower and turns negative for clays — a known feature of the Rawls/Saxton–Rawls lineage (Minasny & McBratney note their neural net “did not show a negative effect with an increase in OC for clay content larger than 60 %”). This is why Section 2 builds the blend on the Minasny & McBratney increments rather than Saxton–Rawls.

# Saxton–Rawls ΔAWC/ΔOC vs. Minasny & McBratney (2018) Table 2.
MM_AWC_SLOPE = {"general": 1.16, "coarse": 1.94, "medium": 1.79, "fine": 1.41}


def sr_awc_slope_per_pct_oc(sand_frac, clay_frac):
    """Saxton–Rawls ΔAWC over OC 0.5%->1.5% (Minasny & McBratney's interval), in mm H2O 100 mm-1 (= vol%)."""
    om_lo, om_hi = 0.5 / VB, 1.5 / VB  # OC% -> OM%
    p_lo, f_lo, _ = saxton_rawls(sand_frac, clay_frac, om_lo)
    p_hi, f_hi, _ = saxton_rawls(sand_frac, clay_frac, om_hi)
    return ((f_hi - p_hi) - (f_lo - p_lo)) * 100.0  # cm3/cm3 over 1% OC -> mm/100mm


val_df = pd.DataFrame(
    [
        {
            "texture_class": cls,
            "mm_group": MM_GROUP[cls],
            "saxton_rawls_dAWC": sr_awc_slope_per_pct_oc(sand / 100, clay / 100),
        }
        for cls, (sand, silt, clay) in TEXTURE_CLASSES.items()
    ]
)

grp_cmp = val_df.groupby("mm_group", sort=False)["saxton_rawls_dAWC"].mean().reset_index()
grp_cmp["mm_group"] = pd.Categorical(grp_cmp["mm_group"], categories=["coarse", "medium", "fine"], ordered=True)
grp_cmp = grp_cmp.sort_values("mm_group")
grp_cmp["minasny_mcbratney"] = grp_cmp["mm_group"].map(MM_AWC_SLOPE)

print("Mean ΔAWC/ΔOC by texture group (mm H₂O 100 mm⁻¹ per +1% OC):")
print(grp_cmp.round(2).to_string(index=False))
val_df.round(2)

Mean ΔAWC/ΔOC by texture group (mm H₂O 100 mm⁻¹ per +1% OC):
mm_group  saxton_rawls_dAWC minasny_mcbratney
  coarse              0.660             1.940
  medium              0.740             1.790
    fine             -0.430             1.410

	texture_class	mm_group	saxton_rawls_dAWC
0	sand	coarse	0.480
1	loamy sand	coarse	0.670
2	sandy loam	coarse	0.860
3	loam	medium	0.820
4	silt loam	medium	0.990
5	silt	medium	1.350
6	sandy clay loam	coarse	0.640
7	clay loam	medium	0.290
8	silty clay loam	medium	0.240
9	sandy clay	fine	0.070
10	silty clay	fine	-0.350
11	clay	fine	-1.010

# Grouped bars: Saxton–Rawls vs Minasny & McBratney mean ΔAWC/ΔOC, by texture group.
hv.output(widget_location="right")

grp_long = grp_cmp.melt(
    id_vars="mm_group",
    value_vars=["saxton_rawls_dAWC", "minasny_mcbratney"],
    var_name="method",
    value_name="dAWC_mm",
)
grp_long["method"] = grp_long["method"].map(
    {"saxton_rawls_dAWC": "Saxton–Rawls (2006)", "minasny_mcbratney": "Minasny & McBratney (2018)"}
).astype(str)

bars = grp_long.hvplot.bar(
    x="mm_group",
    y="dAWC_mm",
    by="method",
    color=["#4c78a8", "#f58518"],
    xlabel="texture group (Minasny & McBratney classes)",
    ylabel="ΔAWC / ΔOC  (mm H₂O 100 mm⁻¹ per +1% OC)",
    title="AWC sensitivity to organic carbon: Saxton–Rawls vs. Minasny & McBratney (2018)",
    width=760,
    height=470,
    legend="top_right",
    ylim=(-1.3, 2.3),
)
refs = (
    hv.HLine(0).opts(color="black", line_width=1)
    * hv.HLine(MM_AWC_SLOPE["general"]).opts(color="gray", line_dash="dotted", line_width=1)
)
bars * refs