add print cell/s

        start = time.perf_counter()

        aquifer_ds = root[aquifer].to_dataset()

        updated_ds = compute_results_for_aquifer(aquifer_ds, aquifer, config, temperature_model)

        updated_ds, n_cells = compute_results_for_aquifer(aquifer_ds, aquifer, config, temperature_model)

        root[aquifer] = xr.DataTree(dataset=updated_ds, name=aquifer)

        end = time.perf_counter()

        elapsed = end - start

        print(f"Time required for aquifer {aquifer}: {elapsed:.2f} seconds")

        cells_per_sec = n_cells / elapsed

        print(f"Time required for aquifer {aquifer}: {elapsed:.2f} seconds. {cells_per_sec} cells per second")

    if config.results_dir != "":

        zarr_path = f"{config.results_dir}/ROSL_ROSLU_test.zarr"

                f"ntg={ntg}, depth={depth:.3f}, porosity={porosity:.3f}, "

                f"utc:{result.utc}"

            )

            result.utc = 1000

            result = result._replace(utc=1000)

        return tuple(result) if result is not None else tuple(np.nan for _ in output_names)

    aquifer_ds["temperature"] = temperature

    y_coord, x_coord = xr.broadcast(aquifer_ds["y"], aquifer_ds["x"])

    valid_mask = (

            aquifer_ds["thickness"].notnull()

            & aquifer_ds["transmissivity"].notnull()

            & aquifer_ds["transmissivity_with_ntg"].notnull()

            & aquifer_ds["ntg"].notnull()

            & aquifer_ds["depth"].notnull()

            & aquifer_ds["porosity"].notnull()

            & aquifer_ds["temperature"].notnull()

            & (aquifer_ds["transmissivity_with_ntg"] >= config.kh_cutoff)

    )

    n_cells = int(valid_mask.sum().compute())

    # use hc_accum mask

    hc_mask = aquifer_ds["hc_accum"] != 0

    vars_to_mask = [

    print(f"Computed results for {aquifer_name}")

    return aquifer_ds

    return aquifer_ds, n_cells

def compute_hydro_properties_mc(