DA with Random Initial Conditions on Soil Layers

This notebook demonstrates how to create an ensemble of models with varying initial conditions for different soil layers. By introducing randomness into the initial conditions, we can better simulate the natural variability in soil properties and improve the robustness of data assimilation (DA) processes.

Estimated time to run the notebook = 2min

import os
import numpy as np
import pyvista as pv
from pyCATHY.DA.perturbate import perturbate_parm
from pyCATHY.DA import perturbate
from pyCATHY.DA.cathy_DA import DA
from pyCATHY.plotters import cathy_plots as CTp

---

nlay = 6
scenario = {
            'per_name':['ic'],
            'per_type': [[None]*nlay],
            'per_nom':[[-15]*nlay],
            'per_mean':[[-15]*nlay],
            'per_sigma': [[3.75]*nlay],
            'per_bounds': [[None]*nlay],
            'sampling_type': [['normal']*nlay],
            'transf_type':[[None]*nlay],
            'listUpdateParm': ['St. var.']
            }

---

simuWithDA = DA(
                        dirName='./',
                        prj_name= 'DA_with_non_uniform_ic',
                        notebook=True,
                    )

# linear z depth
# ---------------
zb = np.linspace(0, 2, nlay)
nstr = len(zb)
zr = list(np.ones(len(zb))/nstr)

simuWithDA.update_prepo_inputs(
                                nstr=nstr,
                                zratio=zr,
                                base=max(zb),
                                )

🏁 Initiate CATHY object
🔄 Update hap.in file
🔄 Update dem_parameters file
🔄 Update dem_parameters file

---

simuWithDA.update_dem_parameters()
simuWithDA.update_veg_map()
simuWithDA.update_soil()

Traceback (most recent call last):
  File "/home/runner/work/pycathy_wrapper/pycathy_wrapper/examples/DA/plot_1_DA_NoUni_ic (conflicted copy 2025-05-21 093615).py", line 57, in <module>
    simuWithDA.update_veg_map()
  File "/home/runner/work/pycathy_wrapper/pycathy_wrapper/pyCATHY/cathy_tools.py", line 3353, in update_veg_map
    indice_veg, str_hd_rootmap = in_CT.read_root_map(os.path.join(
  File "/home/runner/work/pycathy_wrapper/pycathy_wrapper/pyCATHY/importers/cathy_inputs.py", line 529, in read_root_map
    with open(os.path.join(rootmapfile), "r") as f:  # open the file for reading
FileNotFoundError: [Errno 2] No such file or directory: '/home/runner/work/pycathy_wrapper/pycathy_wrapper/examples/DA/./DA_with_non_uniform_ic/input/veg_map'

---

simuWithDA.NENS = 3
list_pert = perturbate.perturbate(
                                    simuWithDA,
                                    scenario,
                                    simuWithDA.NENS,
                                 )

var_per_dict_stacked = {}
for dp in list_pert:
    var_per_dict_stacked = perturbate_parm(
                                var_per_dict_stacked,
                                parm=dp,
                                type_parm = dp['type_parm'], # can also be VAN GENUCHTEN PARAMETERS
                                mean =  dp['mean'],
                                sd =  dp['sd'],
                                sampling_type =  dp['sampling_type'],
                                ensemble_size =  dp['ensemble_size'], # size of the ensemble
                                per_type= dp['per_type'],
                                nlayers = nlay,
                                savefig= os.path.join(simuWithDA.workdir,
                                                      simuWithDA.project_name,
                                                      simuWithDA.project_name + dp['savefig'])
                                )

This in normally directly called when using run_DA_sequential()

simuWithDA._create_subfolders_ensemble()
simuWithDA.update_ENS_files(var_per_dict_stacked,
                            var_per_dict_stacked.keys(),
                            cycle_nb=0
                            )

---

pl = pv.Plotter(shape=(1,2))
for i, ensi in enumerate([1,3]):
    DApath = f'DA_Ensemble/cathy_{ensi}/vtk/'
    path = os.path.join(simuWithDA.workdir,
                        simuWithDA.project_name,
                        DApath,
                        simuWithDA.project_name + '.vtk'
                        )

    pl.subplot(0,i)
    CTp.show_vtk(path,
                 'ic_nodes',
                 ax=pl,
                 clim = [-25,-5],
                 #show_scalar_bar=True,
                 )
    _ = pl.add_legend('')
    pl.add_title(f'Ensemble nb:{ensi}')

pl.show()