Updating snow in CLM

[timhoar]

Use case

Assimilation experiments that update the snow in CLM require a partitioning function that redistributes the
posterior (diagnostic) SWE into the (prognostic) snow variables - which have an arbitrary number of layers.
Most snow observations are of snowcover fraction or snow water equivalent (SWE). It is useful to provide these
diagnostic quantities (frac_sno and H2OSNO, respectively) from the CLM history files as part of the DART vector.
However, updating each layer of the prognostic variables (H2OSOI_LIQ, H2OSOI_ICE) directly through the ensemble
correlation may not result in a posterior SWE that matches the posterior SWE in H2OSNO.

Note that we have never attempted to create snow from scratch. Nor have we ever attempted to change things like grain size or the amount of black carbon. Also, and perhaps most importantly, if one ensemble member did not have (unobserved) snow-related quantities, the (unobserved) snow-related values were left unchanged in ALL ensemble members. This is different than a failed forward operator. The forward operator can calculate a meaningful value if there is snow in the gridcell, but when the state is being updated - anything in the cutoff radius is considered. (It may be clearer to state that if there were only 79 values on the unobserved axis (think snow ice in layer 5) and 80 estimates of the observed SWE, the regression is ill-defined.)

Is your feature request related to a problem?

This is really a feature request, but if there is a problem, it is that updating the prognostic variables directly does not result in a SWE consistent with the observations.

Describe the your prefered solution

Maintain the original ratios of snow in each active layer, but modify the amount of snow in each layer such that the
resulting SWE matches the desired posterior SWE. The most logical place to do this is in dart_to_clm since it
reads both the posterior variables from DART and the original variables from CLM (which are needed to calculate the original ratios of snow in the layers). This removes the consideration that each ensemble member may have a different number of active snow layers in each column.

Eroding the snow is easy. Increasing snow is harder. Creating snow from scratch is impossible.

Describe any alternatives you have considered

Ideally, I suspect it is most appropriate to modify the atmospheric forcing to provide new snow.

[timhoar]

This is awful, but I found a prototype update_snow() routine in an ancient version of the CLM model_mod (circa 2015). This was committed to the classic_Tb branch (brightness temperature) and has extra routines in it for that. The git revision is commit 78007e2b86f25fc1d3bf423d46fc4abcca5cb94d

clm_model_mod_update_snow.f90.txt

I have no idea how useful it may be, and there is certainly lots to be removed from it and added directly to the dart_to_clm program.

[braczka]

Just for clarity -- the proposed partitioning function solution here which is to 'modify the amount of snow in each layer such that the resulting SWE matches the desired posterior SWE' can solve 2 potential problems: 1) the mismatch between the posterior SWE in H2OSNO (which we want) vs. the SWE constructed from the ensemble correlation of each layer of the prognostic variables (H2OSOI_LIQ, HSOSOI_ICE) and 2) remove ill-defined regressions in which each snow layer is not defined for all 80 ensemble members (therefore all snow-related values are left unchanged, no update applied) .

[braczka]

Update on the repartitioning of SWE:

1. I have migrated the subroutine update_snow from the clm_model_mod_ update_snow.f90.txt provided by Tim (from above) and inserted it into the current version of clm_main on one of my local feature branches. The functionality is currently being tested and I will create a PR when ready. I have preserved the repartitioning scheme from this older version, and made it compatible with the current clm_main code version in dart_to_clm.f90. The updated 'dart_to_clm.f90' repartitions snow layers of H2OSOI_ICE, H2OSOI_LIQ, DZSNO, and SNOW_DEPTH. The scheme is invoked as a namelist option from dart_to_clm_nml. For clm variables which have both above and below surface ice/liquid state variables (H2OSOI_LIQ, H2OSOI_ICE) the below surface layers are updated per the traditional dart posterior update. For the above surface layers the new repartitioning scheme is used.

2. In the current implementation of the repartitioning scheme It is unclear if other snow depth variables need to be repartitioned along with DZSNO (snow layer depth). For example repartitioning DZSNO may lead to inconsistencies with ZSNO and ZISNO. I need to look at the updated behavior and decide whether to add ZSNO (middle of snow layer depth) and ZISNO (top interface of snow layer) to the repartitioning update.

3. As currently implemented the repartitioning is a scientific advance -- but does not remove the need for missing values (See issue can we avoid missing values in CLM by changing what fields are in the state? #280). To do that the requirement of Fill_Vals/Missing_vals would have to be removed in clm_to_dart.f90 prior to when the repartitioning is implemented in dart_to_clm.f90. Also, the need for missing values in clm goes well beyond the snow layers and is needed for subsurface layers without initialized land properties (bedrock, ocean etc).

[braczka]

In order to repartition the snowcover SWE, dart_to_clm.f90 requires the posterior of the clm variable H2OSNO. This variable is a bit more challenging to deal with because depending upon the versions of CLM it can appear within the restart file (domain 1) or the history/vector_history files (domains 2,3). For consistency it is more straightforward to generate this variable within a vector_history file which is part of domain 3. However, the current functionality of dart_to_clm.f90 is to locate posterior values from DART within dart_posterior.nc which is linked to clm restart files (domain 1 only). Diagnostic variables within the history/vector_history files are not overwritten with the posterior values during the execution of filter and are not passed on to the output_state_file_list within filter_nml. However, the posterior value of H2OSNO (and other non-restart variables) can still be accessed through the filter stage analysis.nc file. This requires that the analysis stage be output within stages_to_write, but I think this is the preferred route to go in order to retrieve the posterior of H2OSNO. The alternative would be to try to change the filter settings such that other domains besides the restart variables pass the posterior values onto the output_state_file_list. Not sure that this is possible or a preferred approach...?

[braczka]

Update: Repartitioning of snow algorithm has been implemented on my local branch and tested for 7 day period with CLM across global test run. It functions successfully and as expected. Based on discussions with CGD folks the code should update not only 'DZSNO' (snow layer thickness) but also 'ZSNO' (middle layer depth) and 'ZISNO' (interface depth) for consistency. This is because CLM snow algorithm depends upon layer thickness/depth for the compaction/disaggregation algorithm to function correctly within CLM.

The above approach which repartitions all available snow layers for a given column, may inadvertently influence albedo state variables (e.g. black, organic carbon, dust) which could have unintended consequences for surface energy balance. To account for this, it was recommended by CGD to apply DART adjustment to the bottom snow layer only for a given column. This accounts for adjustments in ice/liquid water status, but also maintains radiative characteristics of top snow layers. Will add this approach as a separate partitioning scheme in DART code, and make partitioning scheme (all layers vs. bottom layer only) a user option. Plan to test impact of two schemes as part of Western US domain as part of Utah collaboration.

[hkershaw-brown]

fixed in #313