Skip to content

Add total alkalinity response to benthic-pelagic exchange of nutrients #183

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
yutiPML merged 4 commits into from
Aug 18, 2026
Merged

Add total alkalinity response to benthic-pelagic exchange of nutrients #183

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
eea8e69
Add total alkalinity response to benthic-pelagic exchange of nutrients
hpowley Aug 11, 2026
1f24dd5
Fix bug in transferring code
hpowley Aug 11, 2026
a39b115
addedd comment on new parameters
yutiPML Aug 11, 2026
3f61187
Updating expected results
wathen Aug 13, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion github-actions/fabm0d-gotm-ersem/expected.json
View file
Open in desktop

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion github-actions/gotm-fabm-ersem/expected.json
View file
Open in desktop

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion github-actions/gotm-fabm-ersem/expected_state.json
View file
Open in desktop

Large diffs are not rendered by default.

57 changes: 42 additions & 15 deletions github-actions/regen_expected_results.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,8 @@
import argparse
import json
import netCDF4 as nc
from numpy import ndarray, interp
from numpy import asarray, ndarray, interp


class NumpyEncoder(json.JSONEncoder):
def default(self, obj):
Expand All @@ -10,6 +11,38 @@ def default(self, obj):
return temp
return json.JSONEncoder.default(self, obj)


def extract_gotm_time_series(data, variable_name, depth=0.0):
var = data.variables[variable_name]
values = var[:].squeeze()

if values.ndim == 1:
return [float(v) for v in values]

zi = data.variables['zi'][:].squeeze()
z = data.variables['z'][:].squeeze()
var_time_series = []
for i in range(var.shape[0]):
depth_offset = depth + zi[i].reshape(-1)[-1]
z_row = z[i].reshape(-1)
var_row = values[i].reshape(-1)
var_time_series.append(float(interp(depth_offset, z_row, var_row)))
return var_time_series


def extract_final_state(data, variable_name):
values = data.variables[variable_name][:].squeeze()

if values.ndim == 1:
return float(values[-1])

if values.ndim == 2:
return asarray(values[-1, :])

raise RuntimeError(
f"Unsupported squeezed shape for {variable_name}: {values.shape}"
)

parser = argparse.ArgumentParser()
parser.add_argument('-p', '--data-path', type=str, required=True,
help='Path to output file from model run')
Expand Down Expand Up @@ -51,6 +84,9 @@ def default(self, obj):
expected_results = {}
for v in items:
try:
if v != "dates" and v not in data.variables:
print(f"Skipping missing variable: {v}")
continue
if v == "dates":
times = data.variables['time']
dates = nc.num2date(times[:],
Expand All @@ -59,23 +95,14 @@ def default(self, obj):
dates = [str(d).split(" ")[0] for d in dates]
expected_results[v] = dates
elif key == "expected" and model_run == "gotm":
depth = 0.0
var = data.variables[v]
zi = data.variables['zi'][:].squeeze()
z = data.variables['z'][:].squeeze()
var_time_series = []
for i in range(var.shape[0]):
depth_offset = depth + zi[i, -1]
var_time_series.append(interp(depth_offset, z[i, :], var[i, :].squeeze()))
expected_results[v] = var_time_series
expected_results[v] = extract_gotm_time_series(data, v)

elif model_run == "gotm":
expected_results[v] = extract_final_state(data, v)
elif data.variables[v].ndim == 4:
expected_results[v] = data.variables[v][:].squeeze() \
if model_run == "fabm0d" else expected_results[v][-1, :]
expected_results[v] = data.variables[v][:].squeeze()
elif data.variables[v].ndim == 3:
expected_results[v] = data.variables[v][:].squeeze() \
if model_run == "fabm0d" else \
float(data.variables[v][:].squeeze()[-1])
expected_results[v] = data.variables[v][:].squeeze()
else:
raise RuntimeError
except Exception as e:
Expand Down
14 changes: 13 additions & 1 deletion src/benthic_column_dissolved_matter.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,8 +34,10 @@ module ersem_benthic_column_dissolved_matter
type (type_bottom_state_variable_id) :: id_layer ! depth of bottom interface of own layer (where own concentration drops to zero)
type (type_horizontal_diagnostic_variable_id) :: id_conc_eq(nlayers) ! mean equilibrium pore water concentration in individual layers
type (type_horizontal_diagnostic_variable_id) :: id_conc_tot(nlayers) ! mean pore water concentration in individual layers
type (type_state_variable_id) :: id_TA

real(rk) :: ads(nlayers)
real(rk) :: relax, minD
real(rk) :: relax, minD, TA_factor
integer :: last_layer
logical :: correction
type (type_single_constituent),allocatable :: constituents(:)
Expand Down Expand Up @@ -102,7 +104,11 @@ subroutine benthic_dissolved_matter_initialize(self,configunit)
write (index,'(i0)') ilayer
call self%get_parameter(self%ads(ilayer),'ads'//trim(index),'-','adsorption in layer '//trim(index)//' (total:dissolved)',default=1.0_rk)
end do

call self%get_parameter(self%correction,'correction','','move losses in oxygenic layer to deeper layers if pelagic concentration is limiting',default=.false.)

call self%get_parameter(self%TA_factor,'TA_factor','','contribution to TA',default=0._rk)


! Create model that computes concentrations per benthic layer.
allocate(profile)
Expand Down Expand Up @@ -146,6 +152,7 @@ subroutine benthic_dissolved_matter_initialize(self,configunit)
end do
call self%request_coupling (self%id_Dm(nlayers), depth_of_sediment_column)
call profile%request_coupling(profile%id_Dm(nlayers),depth_of_sediment_column)
call self%register_state_dependency(self%id_TA,standard_variables%alkalinity_expressed_as_mole_equivalent)

end subroutine benthic_dissolved_matter_initialize

Expand Down Expand Up @@ -367,6 +374,9 @@ subroutine process_constituent(self,_ARGUMENTS_DO_BOTTOM_,info)
! Net change in column-integrated mass must equal column-integrated production - surface exchange.
! Thus, surface exchange = column-integrated production - net change (net change = relaxation)
_SET_BOTTOM_EXCHANGE_(info%id_pel,sms-(c_int_eq-(c_int+c_int_deep))/self%relax)

_SET_BOTTOM_EXCHANGE_(self%id_TA,(sms-(c_int_eq-(c_int+c_int_deep))/self%relax)*self%TA_factor)

_SET_HORIZONTAL_DIAGNOSTIC_(info%id_pbf,sms-(c_int_eq-(c_int+c_int_deep))/self%relax)
else
! Apply a "technical correction" in case flux from the oxygenated
Expand Down Expand Up @@ -426,6 +436,8 @@ subroutine process_constituent(self,_ARGUMENTS_DO_BOTTOM_,info)
norm_res_int = poro*sum(self%ads*residual_per_layer)
P_res_int = (c_int-c_int_eq)/norm_res_int*Dm(nlayers)
_SET_BOTTOM_EXCHANGE_(info%id_pel,sms+P_res_int) ! Equilibrium flux = depth-integrated production sms + residual flux P_res_int
_SET_BOTTOM_EXCHANGE_(self%id_TA,(sms+P_res_int)*self%TA_factor)

_SET_HORIZONTAL_DIAGNOSTIC_(info%id_pbf,sms+P_res_int)
_SET_BOTTOM_ODE_(info%id_int,-P_res_int-sms)

Expand Down
3 changes: 3 additions & 0 deletions testcases/ersem-15.06/fabm-ersem-15.06-L4-ben-docdyn-iop-denit.yaml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -752,6 +752,7 @@ instances:
ads1: 100.0 # adsorption in layer 1 (total:dissolved) (-), default = 1.0
ads2: 100.0 # adsorption in layer 2 (total:dissolved) (-), default = 1.0
ads3: 2.0 # adsorption in layer 3 (total:dissolved) (-), default = 1.0
TA_factor: -1. # scaling factor (sign) of the effect of benthic-pelagic fluxes on TA
initialization:
p: 6.8 # phosphorus (mmol/m^2)
coupling:
Expand All @@ -765,6 +766,7 @@ instances:
last_layer: 2 # sediment layer where concentration drops to zero, default = 3
relax: 5.0 # rate of relaxation towards equilibrium concentration profile (1/d)
minD: 0.0001 # minimum depth of zero-concentration isocline (m)
TA_factor: -1. # scaling factor (sign) of the effect of benthic-pelagic fluxes on TA
initialization:
n: 0.584 # nitrogen (mmol/m^2)
coupling:
Expand All @@ -779,6 +781,7 @@ instances:
ads2: 3.0 # adsorption in layer 2 (total:dissolved) (-), default = 1.0
ads3: 3.0 # adsorption in layer 3 (total:dissolved) (-), default = 1.0
correction: true # move losses in oxygenic layer to deeper layers if pelagic concentration is limiting, default = false
TA_factor: 1. # scaling factor (sign) of the effect of benthic-pelagic fluxes on TA
initialization:
n: 25.0 # nitrogen (mmol/m^2)
coupling:
Expand Down
3 changes: 3 additions & 0 deletions testcases/ersem-15.06/fabm-ersem-15.06-L4-ben-docdyn-iop-n2o.yaml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -745,6 +745,7 @@ instances:
ads1: 100.0 # adsorption in layer 1 (total:dissolved) (-), default = 1.0
ads2: 100.0 # adsorption in layer 2 (total:dissolved) (-), default = 1.0
ads3: 2.0 # adsorption in layer 3 (total:dissolved) (-), default = 1.0
TA_factor: -1. # scaling factor (sign) of the effect of benthic-pelagic fluxes on TA
initialization:
p: 6.8 # phosphorus (mmol/m^2)
coupling:
Expand All @@ -758,6 +759,7 @@ instances:
last_layer: 2 # sediment layer where concentration drops to zero, default = 3
relax: 5.0 # rate of relaxation towards equilibrium concentration profile (1/d)
minD: 0.0001 # minimum depth of zero-concentration isocline (m)
TA_factor: -1. # scaling factor (sign) of the effect of benthic-pelagic fluxes on TA
initialization:
n: 0.584 # nitrogen (mmol/m^2)
coupling:
Expand All @@ -772,6 +774,7 @@ instances:
ads2: 3.0 # adsorption in layer 2 (total:dissolved) (-), default = 1.0
ads3: 3.0 # adsorption in layer 3 (total:dissolved) (-), default = 1.0
correction: true # move losses in oxygenic layer to deeper layers if pelagic concentration is limiting, default = false
TA_factor: 1. # scaling factor (sign) of the effect of benthic-pelagic fluxes on TA
initialization:
n: 25.0 # nitrogen (mmol/m^2)
coupling:
Expand Down
3 changes: 3 additions & 0 deletions testcases/ersem-15.06/fabm-ersem-15.06-L4-ben-docdyn-iop.yaml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -735,6 +735,7 @@ instances:
ads1: 100.0 # adsorption in layer 1 (total:dissolved) (-), default = 1.0
ads2: 100.0 # adsorption in layer 2 (total:dissolved) (-), default = 1.0
ads3: 2.0 # adsorption in layer 3 (total:dissolved) (-), default = 1.0
TA_factor: -1. # scaling factor (sign) of the effect of benthic-pelagic fluxes on TA
initialization:
p: 6.8 # phosphorus (mmol/m^2)
coupling:
Expand All @@ -748,6 +749,7 @@ instances:
last_layer: 2 # sediment layer where concentration drops to zero, default = 3
relax: 5.0 # rate of relaxation towards equilibrium concentration profile (1/d)
minD: 0.0001 # minimum depth of zero-concentration isocline (m)
TA_factor: -1. # scaling factor (sign) of the effect of benthic-pelagic fluxes on TA
initialization:
n: 0.584 # nitrogen (mmol/m^2)
coupling:
Expand All @@ -762,6 +764,7 @@ instances:
ads2: 3.0 # adsorption in layer 2 (total:dissolved) (-), default = 1.0
ads3: 3.0 # adsorption in layer 3 (total:dissolved) (-), default = 1.0
correction: true # move losses in oxygenic layer to deeper layers if pelagic concentration is limiting, default = false
TA_factor: 1. # scaling factor (sign) of the effect of benthic-pelagic fluxes on TA
initialization:
n: 25.0 # nitrogen (mmol/m^2)
coupling:
Expand Down
3 changes: 3 additions & 0 deletions testcases/ersem-Blackford2024/fabm-ersem-Blackford2004-L4-ben-stdbac-stdlight.yaml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -718,6 +718,7 @@ instances:
ads1: 100.0 # adsorption in layer 1 (total:dissolved) (-), default = 1.0
ads2: 100.0 # adsorption in layer 2 (total:dissolved) (-), default = 1.0
ads3: 2.0 # adsorption in layer 3 (total:dissolved) (-), default = 1.0
TA_factor: -1. # scaling factor (sign) of the effect of benthic-pelagic fluxes on TA
initialization:
p: 6.8 # phosphorus (mmol/m^2)
coupling:
Expand All @@ -731,6 +732,7 @@ instances:
last_layer: 2 # sediment layer where concentration drops to zero, default = 3
relax: 5.0 # rate of relaxation towards equilibrium concentration profile (1/d)
minD: 0.0001 # minimum depth of zero-concentration isocline (m)
TA_factor: -1. # scaling factor (sign) of the effect of benthic-pelagic fluxes on TA
initialization:
n: 0.584 # nitrogen (mmol/m^2)
coupling:
Expand All @@ -745,6 +747,7 @@ instances:
ads2: 3.0 # adsorption in layer 2 (total:dissolved) (-), default = 1.0
ads3: 3.0 # adsorption in layer 3 (total:dissolved) (-), default = 1.0
correction: true # move losses in oxygenic layer to deeper layers if pelagic concentration is limiting, default = false
TA_factor: 1. # scaling factor (sign) of the effect of benthic-pelagic fluxes on TA
initialization:
n: 25.0 # nitrogen (mmol/m^2)
coupling:
Expand Down
3 changes: 3 additions & 0 deletions testcases/fabm-ersem-26.02-dvm-n2o.yaml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -751,6 +751,7 @@ instances:
ads1: 100.0 # adsorption in layer 1 (total:dissolved) (-), default = 1.0
ads2: 100.0 # adsorption in layer 2 (total:dissolved) (-), default = 1.0
ads3: 2.0 # adsorption in layer 3 (total:dissolved) (-), default = 1.0
TA_factor: -1. # scaling factor (sign) of the effect of benthic-pelagic fluxes on TA
initialization:
p: 6.8 # phosphorus (mmol/m^2)
coupling:
Expand All @@ -764,6 +765,7 @@ instances:
last_layer: 2 # sediment layer where concentration drops to zero, default = 3
relax: 5.0 # rate of relaxation towards equilibrium concentration profile (1/d)
minD: 0.0001 # minimum depth of zero-concentration isocline (m)
TA_factor: -1. # scaling factor (sign) of the effect of benthic-pelagic fluxes on TA
initialization:
n: 0.584 # nitrogen (mmol/m^2)
coupling:
Expand All @@ -778,6 +780,7 @@ instances:
ads2: 3.0 # adsorption in layer 2 (total:dissolved) (-), default = 1.0
ads3: 3.0 # adsorption in layer 3 (total:dissolved) (-), default = 1.0
correction: true # move losses in oxygenic layer to deeper layers if pelagic concentration is limiting, default = false
TA_factor: 1. # scaling factor (sign) of the effect of benthic-pelagic fluxes on TA
initialization:
n: 25.0 # nitrogen (mmol/m^2)
coupling:
Expand Down
3 changes: 3 additions & 0 deletions testcases/fabm-ersem-26.02-dvm.yaml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -741,6 +741,7 @@ instances:
ads1: 100.0 # adsorption in layer 1 (total:dissolved) (-), default = 1.0
ads2: 100.0 # adsorption in layer 2 (total:dissolved) (-), default = 1.0
ads3: 2.0 # adsorption in layer 3 (total:dissolved) (-), default = 1.0
TA_factor: -1. # scaling factor (sign) of the effect of benthic-pelagic fluxes on TA
initialization:
p: 6.8 # phosphorus (mmol/m^2)
coupling:
Expand All @@ -754,6 +755,7 @@ instances:
last_layer: 2 # sediment layer where concentration drops to zero, default = 3
relax: 5.0 # rate of relaxation towards equilibrium concentration profile (1/d)
minD: 0.0001 # minimum depth of zero-concentration isocline (m)
TA_factor: -1. # scaling factor (sign) of the effect of benthic-pelagic fluxes on TA
initialization:
n: 0.584 # nitrogen (mmol/m^2)
coupling:
Expand All @@ -768,6 +770,7 @@ instances:
ads2: 3.0 # adsorption in layer 2 (total:dissolved) (-), default = 1.0
ads3: 3.0 # adsorption in layer 3 (total:dissolved) (-), default = 1.0
correction: true # move losses in oxygenic layer to deeper layers if pelagic concentration is limiting, default = false
TA_factor: 1. # scaling factor (sign) of the effect of benthic-pelagic fluxes on TA
initialization:
n: 25.0 # nitrogen (mmol/m^2)
coupling:
Expand Down
Loading