How to: Apply Observed Water Level as Initial Condition in a 1D D-HYDRO (Delft3D FM) Model

Background

To generate accurate forecasts, it's crucial to represent the current system state as accurately as possible. In operational systems, the system state is often prepared by running a hindcast simulation using past data and generating a state file at the end. This “warm state” is then used as the initial condition for the forecast simulation.

However, when the quality and availability of the model input data are insufficient, it is possible that the uncertainty propagates into the forecast. In this cases, it is beneficial to correct the state with data assimilation techniques, or may be more effective to use direct real-time water level observations as the initial condition for your forecast simulation.

In this How-to, we demonstrate how to setup a 1D D-HYDRO(Delft3DFM) model to run with observed water level as initial conditions, using a python script.

When is this approach suitable?

This method is particularly helpful when:

The model domain is relatively flat and controlled (e.g., canal systems).
There are sufficient water level observations.
Water level gradients are minimal, so velocity initialization can be approximated or ignored.

⚠️ Caution: If the system is highly dynamic (e.g., during or right after a flood event), using only observed water levels may cause instabilities at the model start and increase the required warm-up time.

Step-by-Step Guide

1. Configure the D-HYDRO (Delft3DFM) model to use initial water level file

Set the model to read an initial water level file (InitialWaterLevel.ini) in the Initial Fields file (initialFields.ini). Example:

[General]
fileVersion = 2.00
fileType = iniField
[Initial]
quantity = waterlevel
dataFile = InitialWaterLevel.ini
dataFileType = 1dField
Setup the InitialWaterLevel.ini file by associating the water level gauges to the branch and chainage in the model. Example:

[General]
fileVersion = 2.00
fileType = 1dField
[Global]
quantity = WaterLevel
unit = m
value = 3.000
[Branch]
branchId = Canal 01
values = $Canal_01_upstream$

[Branch]
branchId = Canal 02
numLocations = 2
chainage = 0.000 75.000
values = $Canal_02_upstream$ $Canal_02_downstream$

2. In Delft-FEWS, export the observed water level data as a netcdf file

Export a NetCDF file containing recent observations of the water level gauges that are needed by the model.

3. Replace placeholders using a Python script

Since FEWS currently doesn’t support direct value replacement in NetCDF initial condition files, you can use a python script to:
- Extract the latest observed water levels.
- Replace the corresponding values in the initialWatrerlevel.ini file.

Example:

Generate_initialcondition.py

# Configure logging
log_file = 'Generate_initialcondition.log'
pi_output_file = "Generate_initialcondition.xml"
logging.basicConfig(level=logging.DEBUG,
                    format='%(asctime)s - %(levelname)s - %(message)s',
                    handlers=[
                        logging.FileHandler(log_file),
                        logging.StreamHandler()
                    ])


def convert_to_pi_diagnostics(log_file, pi_output_file):
    if not os.path.exists(log_file):
        logging.error(f"Log file not found: {log_file}")
        return

    namespaces = {
        '': 'http://www.wldelft.nl/fews/PI',
        'xsi': 'http://www.w3.org/2001/XMLSchema-instance'
    }
    root = ET.Element("Diag", version="1.2", xmlns="http://www.wldelft.nl/fews/PI", 
                      xmlns_xsi="http://www.w3.org/2001/XMLSchema-instance",
                      xsi_schemaLocation="http://www.wldelft.nl/fews/PI https://fewsdocs.deltares.nl/schemas/version1.0/pi-schemas/pi_diag.xsd")

    with open(log_file, "r") as file:
        for line in file:
            match = re.match(r"(?P<timestamp>.*?) - (?P<level>\w+) - (?P<message>.*)", line)
            if match:
                timestamp = match.group("timestamp")
                level = match.group("level")
                message = match.group("message")

                level_mapping = {
                    "ERROR": 0,
                    "WARNING": 2,
                    "INFO": 3,
                    "DEBUG": 4
                }
                level_num = level_mapping.get(level, 4)
                event_code = f"EventCode.{level.upper()}"

                ET.SubElement(root, "line", level=str(level_num), 
                              description=message, eventCode=event_code)

    tree = ET.ElementTree(root)
    pi_output_file = Path(pi_output_file).with_suffix('.xml')  # Ensure proper extension
    tree.write(pi_output_file, encoding="utf-8", xml_declaration=True)
    logging.info(f"Pi Diagnostics file created: {pi_output_file}")



def fill_ini_with_nc_data(nc_file_path, ini_file_path, output_file_path="output.ini"):
    # Step 1: Read netCDF file
    ds = nc.Dataset(nc_file_path)

    # Extract station IDs and waterlevel data
    station_ids = ds.variables['station_id'][:]
    waterlevels = ds.variables['waterlevel'][:]

    # Decode bytes if needed
    if station_ids.dtype.kind in {'S', 'O'}:
        station_ids = [b''.join(station_ids_i).decode('utf-8').strip('\x00') for station_ids_i in station_ids]
    
    # Decode masked array
    waterlevels = np.ma.filled(waterlevels, fill_value=np.nan)


    # Take the last time step
    last_values = waterlevels[-1, :]

    # Build dictionary
    station_value_dict = dict(zip(station_ids, last_values))
    logging.debug("NC file parsed. ")

    # Step 2: Read the ini file
    with open(ini_file_path, 'r', encoding='utf-8') as f:
        ini_text = f.read()
    logging.debug("INI file parsed. ")

    # Step 3: Replace placeholders
    def replacement(match):
        key = match.group(1)
        value = station_value_dict.get(key)
        if value is None or np.ma.is_masked(value) or np.isnan(value):
            logging.warning(f"Warning: No value found for {key}")
            return match.group(0)  # Keep original placeholder if missing or invalid
        return f"{value:.3f}"

    filled_ini_text = re.sub(r'\$(.*?)\$', replacement, ini_text)
    logging.debug("INI file filled. ")

    # Step 4: Save to output file
    with open(output_file_path, 'w', encoding='utf-8') as f:
        f.write(filled_ini_text)

    logging.info(f"Filled INI file written to {output_file_path}")


if __name__ == "__main__":
    parser = argparse.ArgumentParser(description='Generate initial condition required by Delft3DFM based on water level in the netcdf file.')
    parser.add_argument('nc_file_path', type=Path, help='Path to the waterlevel.nc. Must contain station_id that will be used by InitialWaterLevel.ini. ')
    parser.add_argument('ini_file_path', type=Path, help='Path to the InitialWaterLevel.ini. Must contain station_id in the netcdf file as property keywords in between ##.')
    args = parser.parse_args()
    
    nc_file_path = args.nc_file_path
    ini_file_path = args.ini_file_path
    output_file_path = ini_file_path
    
    try:
        fill_ini_with_nc_data(nc_file_path, ini_file_path, output_file_path)
    except Exception as e:
        logging.error("An error occurred: %s", e)

    convert_to_pi_diagnostics(log_file, pi_output_file)

4. Integrate in the general adaptor

Add the Python script to the executeActivities in the General Adapter configuration.