Skip to content

add support for sectional dry deposition - #230

Open
zdaq12 wants to merge 6 commits into
compdyn:masterfrom
zdaq12:sect_drydep
Open

zdaq12 wants to merge 6 commits into
compdyn:masterfrom
zdaq12:sect_drydep

Conversation

@zdaq12

@zdaq12 zdaq12 commented Sep 9, 2026

Copy link
Copy Markdown
Contributor

Summary

Adds dry deposition support to the sectional representation. Also corrects the
deposition velocity formula shared by all representations to match
Zhang et al. (2001) and Emerson et al. (2020) as published.

Sectional dry deposition

scenario_binned_loss() applies per-bin loss each timestep via the same
scenario_loss_rate() dispatcher already used by the particle-resolved
and modal paths, so the sectional model supports every
loss_function_type (constant, volume, drydep, chamber), not just dry
deposition.

Each bin decays as m(t+dt) = m(t) * exp(-rate*dt) — the exact solution
of dm/dt = -rate*m. Mass (vol_conc) is
prognostic and num_conc is diagnosed, matching the convention already used by
sectional coagulation.

output_sectional()/input_sectional() gain a scenario argument so
dry deposition parameters are written to and read from NetCDF output,
mirroring output_modal()/input_modal().

Deposition velocity formula change

before: V_d = V_s + 1/(R_a + R_s + R_a*R_s*V_s)
after:  V_d = V_s + 1/(R_a + R_s)

The R_a*R_s*V_s term comes from the two-layer flux balance
(Seinfeld and Pandis, 1998); however, it is not a physically representative term, but
rather a direct artifact of the mathematical manipulation of the multilayer resistance analogy
system of equations with settling acting through all layers. Strictly speaking, neither
Zhang et al. (2001) nor Emerson et al. (2020) use this term. Both follow Slinn (1982)'s
derivation where settling is treated purely in parallel. Since drydep_params_t is explicitly
parameterized to switch between Z01-and E20 coefficient sets, retaining the cross
term would implement neither parameterization as published.

Applied in all three places: scenario_loss_rate_drydep
(particle/sectional), scenario_integrated_loss_rate_drydep (modal,
analytic), and dep_vel_integrand (modal, QUADPACK).

New spec-file option

do_aero_dilution (sectional runs) skips both background-exchange
dilution and the mixing-layer-growth correction when set to no. The
modal representation has no dilution treatment at all, so this is needed
for the modal/sectional comparison and lets dry deposition be
studied in isolation.

Tests

  • test_loss_16 — sectional dry deposition regression test
  • test_loss_17 — checks sectional against particle-resolved

Reference files for test_loss_13, test_loss_14 (QUADPACK), and
test_loss_16 are regenerated due to formula update.

Example scenario

scenarios/7_drydep now demonstrates both modal and sectional representations. The
sectional case is added as drydep_sect.spec with 1_run_sect.sh and
2_process_sect.sh, mirroring the existing modal scripts. It uses the
standard extract_sectional_* tools rather than a dedicated
post-processor. A README documents both cases.

@codecov

codecov Bot commented Sep 9, 2026

Copy link
Copy Markdown

Codecov Report

❌ Patch coverage is 70.31250% with 19 lines in your changes missing coverage. Please review.
✅ Project coverage is 78.47%. Comparing base (7984020) to head (75a7cae).

Files with missing lines Patch % Lines
src/scenario.F90 68.75% 10 Missing ⚠️
src/run_sect.F90 70.83% 7 Missing ⚠️
src/output.F90 71.42% 2 Missing ⚠️
Additional details and impacted files
@@            Coverage Diff             @@
##           master     #230      +/-   ##
==========================================
- Coverage   78.49%   78.47%   -0.03%     
==========================================
  Files          56       56              
  Lines        9755     9805      +50     
==========================================
+ Hits         7657     7694      +37     
- Misses       2098     2111      +13     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@zdaq12
zdaq12 marked this pull request as ready for review September 9, 2026 16:46
@jcurtis2
jcurtis2 self-requested a review September 9, 2026 20:40

@jcurtis2 jcurtis2 left a comment

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Mostly looks straightforward to get this merged. Here are my comments:

  • As I can't offer suggestions on files you haven't touched: quickly just add the output directory scenarios/7_drydep/out to the .gitignore. Looks like we've gotten lazy about this for new scenarios but lets at least fix this one here.
  • You've added this do_aero_dilution to the sectional spec file to disable dilution because the modal lacks dilution entirely so that you can compare the two fairly. In general, this is likely unnecessary as we typically disable dilution by having our input files set in a particular way that causes this process to be zero. Dilution depends on (1) background mixing rate and (2) mixing height changes So for (1) we've controlled it by having rate = 0 in aero_back.dat so that there is no mixing, but we've keep mixing from changes in box height in height.dat. Mixing height does impact the dry deposition rates and the results will change if we change height.dat to lose our boundary layer evolution. I don't know if we consider that important for this scenario at this time. So this feels like a workaround for a deficiency in modal model by introducing something permanent in the sectional model. But this all isn't 100% without impact, its just whether or not that impact matters.
  • Possibly worth a sentence in the PR description that you change the particle size distribution (larger size, broader distribution) where you discuss the example scenario.
  • Not directly related to your changes but since you are making edits in that scenario directory and comparing the modal to sectional, you have the units of tot_mass_conc in the modal processing code to be micrograms per cubic meter.

Comment thread src/run_sect.F90
Comment on lines 149 to 163
call bin_kernel(bin_grid_size(bin_grid), bin_grid%centers, aero_data, &
run_sect_opt%coag_kernel_type, env_state, k_bin)
call smooth_bin_kernel(bin_grid_size(bin_grid), k_bin, ck)
do i = 1,bin_grid_size(bin_grid)
do j = 1,bin_grid_size(bin_grid)
ck(i,j) = ck(i,j) * 1d6 ! m^3/s to cm^3/s
end do
end do

! multiply kernel with constant timestep and logarithmic grid distance
do i = 1,bin_grid_size(bin_grid)
do j = 1,bin_grid_size(bin_grid)
ck(i,j) = ck(i,j) * run_sect_opt%del_t * bin_grid%widths(i)
end do
end do

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not your mistake, but I think we should take this time to wrap this code in a if (run_sect_opt%do_coagulation) to ensure that this is only called when coagulation is on.

Doing this change, you can fix your sectional spec files that required you to have do_coagulation yes and setting the kernel to zero.

Comment thread src/scenario.F90

end do
else
return

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Change this to a die_msg to handle the other loss processes that aren't currently implemented so that if someone turns them on, they fail. Right now, this reads like someone could have these loss processes on and have the simulation run successfully.

The modal model (scenario_update_aero_modes) has a similar issue that I didn't catch in the previous review. It felt more complete when looking at it as it added a check for SCENARIO_LOSS_FUNCTION_CONSTANT (although it only returns) but it doesn't handle the others (there's no final else in that function).

Comment thread src/run_sect.F90 Outdated
Comment thread src/run_sect.F90 Outdated
Comment thread src/run_sect.F90
Comment on lines +316 to +320
call spec_file_read_logical(file, 'do_optical', run_sect_opt%do_optical)
if (run_sect_opt%do_optical) then
call spec_file_die_msg(527436819, file, &
"sectional run does not support optical properties calculation")
end if

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

While neither the optical properties or MOSAIC are supported, move this block inside the logic for checking MOSAIC. The particle model (in run_part.F90) checks for do_mosaic and then if do_mosaic is true, it will check if optical properties are true. It will just make it easier in the future when MOSAIC is supported.

partmc/src/run_part.F90

Lines 204 to 212 in 7984020

if (run_part_opt%do_mosaic) then
call mosaic_init(env_state, aero_data, run_part_opt%del_t, &
run_part_opt%do_optical)
if (run_part_opt%do_optical) then
call mosaic_aero_optical_init(env_state, aero_data, &
aero_state, gas_data, gas_state)
call mosaic_optical_wavelengths(aero_data)
end if
end if

Comment thread src/scenario.F90
Comment on lines +1181 to +1217
!> Updates an array (of size equal to the number of sections) containing
!> the dry deposition velocity for each bin in a sectional simulation.
!!
!! This is a diagnostic helper: it returns deposition velocities (m s^{-1})
!! rather than loss rates, by multiplying the rate from
!! scenario_loss_rate_drydep() back by the mixing layer height. It is not
!! used by the sectional time-stepping itself, which calls
!! scenario_binned_loss().
subroutine scenario_section_drydep_rates(scenario, bin_grid, aero_data, &
env_state, rates)

!> Scenario data.
type(scenario_t), intent(in) :: scenario
!> Bin grid.
type(bin_grid_t), intent(in) :: bin_grid
!> Aerosol data.
type(aero_data_t), intent(in) :: aero_data
!> Environmental state.
type(env_state_t), intent(in) :: env_state
!> Deposition velocities for each section/bin (m s^{-1}).
real(kind=dp), intent(inout) :: rates(:)

integer :: i_bin
real(kind=dp) :: density, vol

call assert_msg(516274839, size(rates) == bin_grid_size(bin_grid), &
"rates array size must match the number of bins")

density = aero_data%density(1)

do i_bin = 1,bin_grid_size(bin_grid)
vol = aero_data_rad2vol(aero_data, bin_grid%centers(i_bin))
rates(i_bin) = scenario_loss_rate_drydep(vol, density, &
aero_data, env_state, scenario) * env_state%height
end do

end subroutine scenario_section_drydep_rates

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is this worth keeping/maintaining?

Comment thread src/output.F90

if (present(scenario)) then
call drydep_params_input_netcdf(scenario%drydep, ncid)
end if

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is called unconditionally but should be guarded. The drydep_params_output_netcdf is only called when the loss process is dry deposition

    if (scenario%loss_function_type == SCENARIO_LOSS_FUNCTION_DRYDEP) then
       call drydep_params_output_netcdf(scenario%drydep, ncid)
    end if

but drydep_params_input_netcdf would always expect it provided scenario is passed to input_sectional.

Comment thread scenarios/7_drydep/README
Comment on lines +40 to +42
only active process. The sectional case runs with the zero
coagulation kernel. It also sets do_aero_dilution to no; the modal
representation has no dilution treatment.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just a note to edit this if changes are made regarding coagulation or dilution.

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You can remove this spec file entirely - this directory isn't doing any dilution.

Comment thread test/loss/test_loss_17.sh

Copy link
Copy Markdown
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Change the spec file and related output processing for the sectional run as the particle run has no dilution.

zdaq12 and others added 2 commits September 16, 2026 18:53
Co-authored-by: Jeffrey Curtis <jcurtis2@illinois.edu>
Co-authored-by: Jeffrey Curtis <jcurtis2@illinois.edu>

This branch has not been deployed

No deployments
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants