The function is_mortar_conductive should not change the values of lower-dimensional interfaces, since they are assumed to always be conductive.

Typically, this is done through the specific volume, which can be obtained from the typical length of the subdomain (i.e., the aperture). Ask Ivar for help on this one.

In the volume conservation equation for lower-dimensional subdomains: Should we use the integrated mortar fluxes (the one given by default by PorePy) or should we scale by the mortar cell volume and effectively use the normal velocities?

Specifically, we would like to avoid init_res_norm to be potentially non-defined:

# Compute 'error' as norm of the residual
res_norm = np.linalg.norm(b, 2)
if itr == 0:
    init_res_norm = res_norm
else:
    inti_res_norm = max(res_norm, init_res_norm)
rel_res = res_norm / init_res_norm

Ideally, one should be able to pass only pressure_threshold, proj_h_high and proj_h_low.

The trace operator currently takes into account only the projection of the traces from the internal boundary of the bulk onto the highest-dimensional interfaces. If one wants to project from 1D to 0D, these are filled with zeros. Naturally, this creates a difference in hydraulic head and consequently an artificial mortar flux.

A correct definition of the projection operators should solve this issue.

Relative permeability has to be included in the definition of the mortar fluxes

To initialize it, use is_mortar_conductive = np.zeros(gb.num_mortar_cells(), dtype=bool).

The value corresponding to dry cell within a fracture domain must be set such that:

h_frac = pressure_threshold + datum

The function is_water_volume_negative can be deprecated and replaced by np.any(vol(h_frac.evaluate(dof_manager).val) < 0) in all scripts.

The pressure threshold in a fracture is not well defined when the fracture neighbors two different types of soils. This introduces a problem since the water volume correction depends on this value.

A possible solution is to introduce a representative fracture pressure threshold. This will be a single number representing a sort of equivalent pressure threshold for the fracture.

The natural way to define it will be as a weighted average, where the weights include the area occupied by the specific type of soil.

Currently, the conductivity state (1 = conductive, 0 = blocking) in lower-dimensional mortar cells is not computed correctly. This depends on wheter the projected hydraulic heads from the lower- and higher-dimensional neighboring subdomains are either dry or wet. However, there is currently an inconsistency.

Indeed, the problem seems to be deeper and related to the values of "dry" hydraulic head being different when using ghost cells.

Plausible solutions:

1. Use the value of pressure head to determine wheter a mortar cell is conductive or blocking. We might even consider using the pressure head in the interface flux calculation (think a bit more about this, does it make sense?)

This can be done by setting gb.pressure_threshold = val. It might not be the best option, but it will do the work.

Currently, if t = sim_time and the solution did not convergence (this can be because we have not satisfied the tolerance, or we found negative volume), the simulation ends, rather than be recomputed.

Check if this is related with the time loop (user-side) or a bug in the time-stepping scheme from the core of PorePy.

A new class DataExporter must be created, for each numerical example. This class has to be able to effectively store primary and secondary variables, i.e., hydraulic_head, pressure_head, water_volume, water_table_height, water_content, saturation, etc.