H2 bulk to surface exceeds 2 Monolayers

[GijsVermarien]

Title: H2 surface bulk asymetry
Description:

If the bulk is shrinking, we can get more than 2 monolayers by H2 moving from bulk to surface (H2 can form in reactions in the bulk very slowly, even if it never can move into the bulk)

Steps to Reproduce:

Disable new

See: test_ice_dependent_desorption_changes_chemistry

Environment:
NA, develop branch

Additional Context:

_______________ test_ice_dependent_desorption_changes_chemistry ________________
[gw0] linux -- Python 3.12.13 /opt/hostedtoolcache/Python/3.12.13/x64/bin/python

temp_output_dir = '/tmp/uclchem_surface_test_ght28ybk'

    def test_ice_dependent_desorption_changes_chemistry(temp_output_dir):
        """Test that ice-coverage-dependent desorption actually affects chemistry.
    
        Verifies:
        - Model runs successfully with new features
        - Ice builds up over time
        - Surface chemistry evolves as ice grows
        - Chemical desorption efficiency varies with coverage
    
        """
        param_dict = {
            "endAtFinalDensity": False,
            "freefall": False,
            "initialDens": 1.0e4,
            "initialTemp": 10.0,
            "finalTime": 1.0e6,
            "points": 1,  # Explicitly set to 0D mode to avoid state pollution from 1D tests
            "enable_radiative_transfer": False,  # Explicitly disable to avoid state pollution from 1D tests
            "desorb": True,
            "chemdesorb": True,
            "outputFile": str(Path(temp_output_dir) / "surface_test.dat"),
        }
    
        result = uclchem.model.Cloud(param_dict=param_dict)
    
        # Basic checks
        assert result is not None, "Model failed to run"
        result.check_error()
    
        # Get dataframe
        df = result.get_joined_dataframes()
        assert len(df) > 0, "No output produced"
    
        # Verify ice buildup (SURFACE + BULK = total ice)
        early_ice = df["SURFACE"].iloc[2] + df["BULK"].iloc[2]
        late_ice = df["SURFACE"].iloc[-1] + df["BULK"].iloc[-1]
        assert late_ice > early_ice * 10, (
            f"Ice should build up significantly: {early_ice:.2e} → {late_ice:.2e}"
        )
    
        # Verify chemistry evolves
        early_co = df["#CO"].iloc[2]
        late_co = df["#CO"].iloc[-1]
        assert late_co != early_co, "Surface CO should evolve"
    
        # Verify complex molecules form (tests that chemistry is active)
        final_ch4 = df["CH4"].iloc[-1]
        assert final_ch4 > 1e-20, "Complex molecules should form"
    
        settings = uclchem.advanced.GeneralSettings()
        num_monolayers_is_surface = list(
            settings.search("num_monolayers_is_surface", True, True).values()
        )[0].get()
        gas_dust_density_ratio = list(
            settings.search("gas_dust_density_ratio", True, True).values()
        )[0].get()
        num_sites_per_grain = list(
            settings.search("num_sites_per_grain", True, True).values()
        )[0].get()
    
        num_monolayers_in_run = (
            df["SURFACE"].max() * gas_dust_density_ratio / num_sites_per_grain
        )
    
        tol = 5e-2
>       assert num_monolayers_in_run <= num_monolayers_is_surface + tol, (
            f"Number of monolayers of surface should be less than {num_monolayers_is_surface}, but was {num_monolayers_in_run:.2f} at most"
        )
E       AssertionError: Number of monolayers of surface should be less than 2.0, but was 2.82 at most
E       assert np.float64(2.8220666051813117) <= (array(2.) + 0.05)

[GijsVermarien]

@TobiasDijkhuis I'd like to dicuss what your intended behavior is here!