Skip to content

API Reference

BRW09Forcing

Surface and ocean temperature data loaded from thermistor temperature record during the Barrow 2009 field study.

Source code in seaice3p/params/forcing.py 
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
@serde(type_check=coerce)
class BRW09Forcing:
    """Surface and ocean temperature data loaded from thermistor temperature record
    during the Barrow 2009 field study.
    """

    Barrow_top_temperature_data_choice: str = "air"

    def __post_init__(self):
        """populate class attributes with barrow dimensional air temperature
        and time in days (with missing values filtered out).

        Note the metadata explaining how to use the barrow temperature data is also
        in seaice3p/forcing_data. The indices corresponding to days and air temp are
        hard coded in as class variables.
        """
        DATA_INDICES = {
            "time": 0,
            "air": 8,
            "bottom_snow": 18,
            "top_ice": 19,
        }
        data = np.genfromtxt(
            Path(__file__).parent.parent / "forcing_data/BRW09.txt", delimiter="\t"
        )
        top_temp_index = DATA_INDICES[self.Barrow_top_temperature_data_choice]
        time_index = DATA_INDICES["time"]

        barrow_top_temp = data[:, top_temp_index]
        barrow_days = data[:, time_index] - data[0, time_index]
        barrow_top_temp, barrow_days = _filter_missing_values(
            barrow_top_temp, barrow_days
        )

        self.barrow_top_temp = barrow_top_temp
        self.barrow_days = barrow_days

__post_init__()

populate class attributes with barrow dimensional air temperature and time in days (with missing values filtered out).

Note the metadata explaining how to use the barrow temperature data is also in seaice3p/forcing_data. The indices corresponding to days and air temp are hard coded in as class variables.

Source code in seaice3p/params/forcing.py 
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
def __post_init__(self):
    """populate class attributes with barrow dimensional air temperature
    and time in days (with missing values filtered out).

    Note the metadata explaining how to use the barrow temperature data is also
    in seaice3p/forcing_data. The indices corresponding to days and air temp are
    hard coded in as class variables.
    """
    DATA_INDICES = {
        "time": 0,
        "air": 8,
        "bottom_snow": 18,
        "top_ice": 19,
    }
    data = np.genfromtxt(
        Path(__file__).parent.parent / "forcing_data/BRW09.txt", delimiter="\t"
    )
    top_temp_index = DATA_INDICES[self.Barrow_top_temperature_data_choice]
    time_index = DATA_INDICES["time"]

    barrow_top_temp = data[:, top_temp_index]
    barrow_days = data[:, time_index] - data[0, time_index]
    barrow_top_temp, barrow_days = _filter_missing_values(
        barrow_top_temp, barrow_days
    )

    self.barrow_top_temp = barrow_top_temp
    self.barrow_days = barrow_days

BRW09InitialConditions dataclass

values for bottom (ocean) boundary

Source code in seaice3p/params/dimensional/initial_conditions.py 
12
13
14
15
16
17
@serde(type_check=coerce)
@dataclass(frozen=True)
class BRW09InitialConditions:
    """values for bottom (ocean) boundary"""

    Barrow_initial_bulk_gas_in_ice: float = 1 / 5

BRW09OceanForcing

Ocean temperature provided by Barrow 2009 data at 2.4m and specify ocean fixed gas saturation state

Source code in seaice3p/params/ocean_forcing.py 
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
@serde(type_check=coerce)
class BRW09OceanForcing:
    """Ocean temperature provided by Barrow 2009 data at 2.4m and specify ocean
    fixed gas saturation state"""

    ocean_gas_sat: float = 1.0

    def __post_init__(self):
        """populate class attributes with barrow dimensional ocean temperature
        and time in days (with missing values filtered out).

        Note the metadata explaining how to use the barrow temperature data is also
        in seaice3p/forcing_data.
        """
        data = np.genfromtxt(
            Path(__file__).parent.parent / "forcing_data/BRW09.txt", delimiter="\t"
        )
        ocean_temp_index = 43
        time_index = 0

        barrow_bottom_temp = data[:, ocean_temp_index]
        barrow_ocean_days = data[:, time_index] - data[0, time_index]
        barrow_bottom_temp, barrow_ocean_days = _filter_missing_values(
            barrow_bottom_temp, barrow_ocean_days
        )

        self.barrow_bottom_temp = barrow_bottom_temp
        self.barrow_ocean_days = barrow_ocean_days

__post_init__()

populate class attributes with barrow dimensional ocean temperature and time in days (with missing values filtered out).

Note the metadata explaining how to use the barrow temperature data is also in seaice3p/forcing_data.

Source code in seaice3p/params/ocean_forcing.py 
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
def __post_init__(self):
    """populate class attributes with barrow dimensional ocean temperature
    and time in days (with missing values filtered out).

    Note the metadata explaining how to use the barrow temperature data is also
    in seaice3p/forcing_data.
    """
    data = np.genfromtxt(
        Path(__file__).parent.parent / "forcing_data/BRW09.txt", delimiter="\t"
    )
    ocean_temp_index = 43
    time_index = 0

    barrow_bottom_temp = data[:, ocean_temp_index]
    barrow_ocean_days = data[:, time_index] - data[0, time_index]
    barrow_bottom_temp, barrow_ocean_days = _filter_missing_values(
        barrow_bottom_temp, barrow_ocean_days
    )

    self.barrow_bottom_temp = barrow_bottom_temp
    self.barrow_ocean_days = barrow_ocean_days

Config dataclass

contains all information needed to run a simulation and save output

this config object can be saved and loaded to a yaml file.

Source code in seaice3p/params/params.py 
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
@serde(type_check=coerce)
@dataclass(frozen=True)
class Config:
    """contains all information needed to run a simulation and save output

    this config object can be saved and loaded to a yaml file."""

    name: str
    total_time: float
    savefreq: float

    physical_params: PhysicalParams
    bubble_params: BubbleParams
    brine_convection_params: BrineConvectionParams
    forcing_config: ForcingConfig
    ocean_forcing_config: OceanForcingConfig
    initial_conditions_config: InitialConditionsConfig
    numerical_params: NumericalParams = NumericalParams()
    scales: Scales | None = None

    def save(self, directory: Path):
        with open(directory / f"{self.name}.yml", "w") as outfile:
            outfile.write(to_yaml(self))

    @classmethod
    def load(cls, path):
        with open(path, "r") as infile:
            yaml = infile.read()
        return from_yaml(cls, yaml)

ConstantForcing dataclass

Constant temperature forcing

Source code in seaice3p/params/forcing.py 
37
38
39
40
41
42
@serde(type_check=coerce)
@dataclass(frozen=True)
class ConstantForcing:
    """Constant temperature forcing"""

    constant_top_temperature: float = -1.5

CubicLiquidus dataclass

Cubic fit to liquidus to give liquidus salinity in terms of temperature

S = a0 + a1 T + a2 T^2 + a3 T^3

defaults are taken from Notz PhD thesis for fit to Assur seawater data

Source code in seaice3p/params/dimensional/water.py 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
@serde(type_check=coerce)
@dataclass(frozen=True)
class CubicLiquidus:
    """Cubic fit to liquidus to give liquidus salinity in terms of temperature

    S = a0 + a1 T + a2 T^2 + a3 T^3

    defaults are taken from Notz PhD thesis for fit to Assur seawater data
    """

    eutectic_temperature: float = -21.1  # deg Celsius
    a0: float = -1.2
    a1: float = -21.8
    a2: float = -0.919
    a3: float = -0.0178

    def get_liquidus_salinity(self, temperature):
        return (
            self.a0
            + self.a1 * temperature
            + self.a2 * temperature**2
            + self.a3 * temperature**3
        )

    def get_liquidus_temperature(self, salinity):
        temperature = fsolve(
            lambda x: salinity - self.get_liquidus_salinity(x),
            np.full_like(salinity, -2),
        )
        if temperature.size == 1:
            return temperature[0]
        else:
            return temperature

DISEQPhysicalParams dataclass

Bases: BasePhysicalParams

non dimensional numbers for the mushy layer

Source code in seaice3p/params/physical.py 
50
51
52
53
54
55
56
@serde(type_check=coerce)
@dataclass(frozen=True)
class DISEQPhysicalParams(BasePhysicalParams):
    """non dimensional numbers for the mushy layer"""

    # only used in DISEQ model
    damkohler_number: float = 1

DimensionalBRW09OceanForcing dataclass

Ocean temperature provided by Barrow 2009 data at 2.4m and specify ocean fixed gas saturation state

Source code in seaice3p/params/dimensional/ocean_forcing.py 
56
57
58
59
60
61
62
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalBRW09OceanForcing:
    """Ocean temperature provided by Barrow 2009 data at 2.4m and specify ocean
    fixed gas saturation state"""

    pass

DimensionalConstantTurbulentFlux dataclass

Parameters for calculating the turbulent surface sensible and latent heat fluxes

NOTE: If you are running a simulation with ERA5 reanalysis forcing you must set the ref_height=2m as this is the appropriate value for the atmospheric reanalysis quantities

The windspeed given here will only be used with ERA5 forcing if the windspeed key is set to None in the forcing_data_file_keys dictionary.

Source code in seaice3p/params/dimensional/forcing.py 
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalConstantTurbulentFlux:
    """Parameters for calculating the turbulent surface sensible and latent heat
    fluxes

    NOTE: If you are running a simulation with ERA5 reanalysis forcing you must set
    the ref_height=2m as this is the appropriate value for the atmospheric reanalysis
    quantities

    The windspeed given here will only be used with ERA5 forcing if the windspeed key
    is set to None in the forcing_data_file_keys dictionary.
    """

    ref_height: float = 10  # m
    windspeed: float = 5  # m/s
    air_temp: float = 0  # deg C
    specific_humidity: float = 3.6e-3  # kg water / kg air
    atm_pressure: float = 101.325  # KPa

    air_density: float = 1.275  # kg/m3
    air_heat_capacity: float = 1005  # J/kg K
    air_latent_heat_of_vaporisation: float = 2.501e6  # J/kg

DimensionalERA5Forcing

read ERA5 data from netCDF file located at data_path.

Simulation will take atmospheric forcings from the start date specified in the string format YYYY-MM-DD

forcing_data_file_keys is a mapping of the descriptive names of the forcing data to be provided to the simulationa and the values are the corresponding strings giving the name of that variable in the netCDF file. The default values are the ERA5 variable names and the SnowModel-LG snow depth name.

Note that if you pass "sd" for the snow depth the simulation will assume you have provided snow depth in m of water equivalent and you must provide a snow density for the conversion.

If you pass None for the snow depth the simulation will procede with no snow layer. If you pass None for the windspeed the simulation will use the constant windspeed defined in the turbulent flux forcing parameters.

Source code in seaice3p/params/dimensional/forcing.py 
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
@serde(type_check=coerce)
class DimensionalERA5Forcing:
    """read ERA5 data from netCDF file located at data_path.

    Simulation will take atmospheric forcings from the start date specified in the
    string format YYYY-MM-DD

    forcing_data_file_keys is a mapping of the descriptive names of the
    forcing data to be provided to the simulationa and the values are the corresponding
    strings giving the name of that variable in the netCDF file.
    The default values are the ERA5 variable names and the SnowModel-LG snow depth name.

    Note that if you pass "sd" for the snow depth the simulation will assume you have
    provided snow depth in m of water equivalent and you must provide a snow density for
    the conversion.

    If you pass None for the snow depth the simulation will procede with no snow layer.
    If you pass None for the windspeed the simulation will use the constant windspeed
    defined in the turbulent flux forcing parameters.
    """

    data_path: Path
    start_date: str  # YYYY-MM-DD
    forcing_data_file_keys: ERA5FileKeys = ERA5FileKeys()
    SW_forcing: DimensionalSWForcing = DimensionalConstantSWForcing()
    LW_forcing: DimensionalLWForcing = DimensionalConstantLWForcing()
    turbulent_flux: DimensionalTurbulentFlux = DimensionalConstantTurbulentFlux()
    oil_heating: DimensionalOilHeating = DimensionalBackgroundOilHeating()

DimensionalFixedHeatFluxOceanForcing dataclass

Provides constant ocean heat flux at the bottom of the domain

Parameters:

Name Type Description Default
ocean_heat_flux float

The constant heat flux at the bottom of the domain in W/m2

 1
Source code in seaice3p/params/dimensional/ocean_forcing.py 
14
15
16
17
18
19
20
21
22
23
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalFixedHeatFluxOceanForcing:
    """Provides constant ocean heat flux at the bottom of the domain

    Args:
        ocean_heat_flux: The constant heat flux at the bottom of the domain in W/m2
    """

    ocean_heat_flux: float = 1

DimensionalFixedTempOceanForcing dataclass

Fixed temperature and gas saturation ocean boundary condition

Source code in seaice3p/params/dimensional/ocean_forcing.py 
 6
 7
 8
 9
10
11
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalFixedTempOceanForcing:
    """Fixed temperature and gas saturation ocean boundary condition"""

    ocean_temp: float = -1

DimensionalMonoBubbleParams dataclass

Bases: DimensionalBaseBubbleParams

Source code in seaice3p/params/dimensional/bubble.py 
15
16
17
18
19
20
21
22
23
24
25
26
27
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalMonoBubbleParams(DimensionalBaseBubbleParams):
    bubble_radius: float = 1e-3  # bubble radius in m

    @property
    def bubble_radius_scaled(self):
        r"""calculate the bubble radius divided by the pore scale

        .. math:: \Lambda = R_B / R_0

        """
        return self.bubble_radius / self.pore_radius

bubble_radius_scaled property

calculate the bubble radius divided by the pore scale

.. math:: \Lambda = R_B / R_0

DimensionalMonthlyHeatFluxOceanForcing dataclass

Provides constant ocean heat flux at the bottom of the domain in each month

Proivde an average monthly ocean heat flux with the entries i=0, 1, 2, 3, ...., 11 in the tuple corresponding to the months January, February, March, April, ...., December

Parameters:

Name Type Description Default
monthly_ocean_heat_flux Tuple[float, ...]

Tuple of ocean heat flux values in each month in W/m2

 (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0)
Source code in seaice3p/params/dimensional/ocean_forcing.py 
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalMonthlyHeatFluxOceanForcing:
    """Provides constant ocean heat flux at the bottom of the domain in each month

    Proivde an average monthly ocean heat flux with the entries
    i=0, 1, 2, 3, ...., 11
    in the tuple corresponding to the months
    January, February, March, April, ...., December

    Args:
        monthly_ocean_heat_flux: Tuple of ocean heat flux values in each month in W/m2
    """

    monthly_ocean_heat_flux: Tuple[float, ...] = (
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
    )

DimensionalParams dataclass

Contains all dimensional parameters needed to calculate non dimensional numbers.

To see the units each input should have look at the comment next to the default value.

Source code in seaice3p/params/dimensional/dimensional.py 
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalParams:
    """Contains all dimensional parameters needed to calculate non dimensional numbers.

    To see the units each input should have look at the comment next to the default
    value.
    """

    name: str
    total_time_in_days: float
    savefreq_in_days: float
    lengthscale: float

    gas_params: DimensionalEQMGasParams | DimensionalDISEQGasParams
    bubble_params: DimensionalMonoBubbleParams | DimensionalPowerLawBubbleParams
    brine_convection_params: DimensionalRJW14Params | NoBrineConvection
    forcing_config: DimensionalRadForcing | DimensionalBRW09Forcing | DimensionalConstantForcing | DimensionalYearlyForcing | DimensionalRobinForcing | DimensionalERA5Forcing
    ocean_forcing_config: DimensionalBRW09OceanForcing | DimensionalFixedTempOceanForcing | DimensionalFixedHeatFluxOceanForcing | DimensionalMonthlyHeatFluxOceanForcing
    initial_conditions_config: DimensionalOilInitialConditions | UniformInitialConditions | BRW09InitialConditions | PreviousSimulation

    water_params: DimensionalWaterParams = DimensionalWaterParams()
    numerical_params: NumericalParams = NumericalParams()
    frame_velocity_dimensional: float = 0  # velocity of frame in m/day
    gravity: float = 9.81  # m/s2

    @property
    def damkohler_number(self):
        r"""Return damkohler number as ratio of thermal timescale to nucleation
        timescale
        """
        if isinstance(self.gas_params, DimensionalEQMGasParams):
            return None

        return (
            (self.lengthscale**2) / self.water_params.thermal_diffusivity
        ) / self.gas_params.nucleation_timescale

    @property
    def total_time(self):
        """calculate the total time in non dimensional units for the simulation"""
        return self.total_time_in_days / self.scales.time_scale

    @property
    def savefreq(self):
        """calculate the save frequency in non dimensional time"""
        return self.savefreq_in_days / self.scales.time_scale

    @property
    def frame_velocity(self):
        """calculate the frame velocity in non dimensional units"""
        return self.frame_velocity_dimensional / self.scales.velocity_scale

    @property
    def B(self):
        r"""calculate the non dimensional scale for buoyant rise of gas bubbles as

        .. math:: \mathcal{B} = \frac{\rho_l g R_0^2 h}{3 \mu \kappa}

        """
        stokes_velocity = (
            (self.water_params.liquid_density - self.gas_params.gas_density)
            * self.gravity
            * self.bubble_params.pore_radius**2
            / (3 * self.water_params.liquid_viscosity)
        )
        velocity_scale_in_m_per_second = (
            self.water_params.thermal_diffusivity / self.lengthscale
        )
        return stokes_velocity / velocity_scale_in_m_per_second

    @property
    def Rayleigh_salt(self):
        r"""Calculate the haline Rayleigh number as

        .. math:: \text{Ra}_S = \frac{\rho_l g \beta \Delta S H K_0}{\kappa \mu}

        """
        match self.brine_convection_params:
            case DimensionalRJW14Params():
                return (
                    self.water_params.liquid_density
                    * self.gravity
                    * self.water_params.haline_contraction_coefficient
                    * self.water_params.salinity_difference
                    * self.lengthscale
                    * self.brine_convection_params.reference_permeability
                    / (
                        self.water_params.thermal_diffusivity
                        * self.water_params.liquid_viscosity
                    )
                )
            case NoBrineConvection():
                return None

    @property
    def expansion_coefficient(self):
        r"""calculate

        .. math:: \chi = \rho_l \xi_{\text{sat}} / \rho_g

        """
        return (
            self.water_params.liquid_density
            * self.gas_params.saturation_concentration
            / self.gas_params.gas_density
        )

    @property
    def lewis_gas(self):
        r"""Calculate the lewis number for dissolved gas, return np.inf if there is no
        dissolved gas diffusion.

        .. math:: \text{Le}_\xi = \kappa / D_\xi

        """
        if self.gas_params.gas_diffusivity == 0:
            return np.inf

        return self.water_params.thermal_diffusivity / self.gas_params.gas_diffusivity

    @property
    def scales(self):
        """return a Scales object used for converting between dimensional and non
        dimensional variables."""
        return Scales(
            self.lengthscale,
            self.water_params.thermal_diffusivity,
            self.water_params.liquid_thermal_conductivity,
            self.water_params.ocean_salinity,
            self.water_params.salinity_difference,
            self.water_params.ocean_freezing_temperature,
            self.water_params.temperature_difference,
            self.gas_params.gas_density,
            self.water_params.liquid_density,
            self.water_params.ice_density,
            self.gas_params.saturation_concentration,
            self.bubble_params.pore_radius,
            self.water_params.haline_contraction_coefficient,
        )

    def save(self, directory: Path):
        """save this object to a yaml file in the specified directory.

        The name will be the name given with _dimensional appended to distinguish it
        from a saved non-dimensional configuration."""
        with open(directory / f"{self.name}_dimensional.yml", "w") as outfile:
            outfile.write(to_yaml(self))

    @classmethod
    def load(cls, path):
        """load this object from a yaml configuration file."""
        with open(path, "r") as infile:
            yaml = infile.read()
        return from_yaml(cls, yaml)

B property

calculate the non dimensional scale for buoyant rise of gas bubbles as

.. math:: \mathcal{B} = \frac{\rho_l g R_0^2 h}{3 \mu \kappa}

Rayleigh_salt property

Calculate the haline Rayleigh number as

.. math:: \text{Ra}_S = \frac{\rho_l g \beta \Delta S H K_0}{\kappa \mu}

damkohler_number property

Return damkohler number as ratio of thermal timescale to nucleation timescale

expansion_coefficient property

calculate

.. math:: \chi = \rho_l \xi_{\text{sat}} / \rho_g

frame_velocity property

calculate the frame velocity in non dimensional units

lewis_gas property

Calculate the lewis number for dissolved gas, return np.inf if there is no dissolved gas diffusion.

.. math:: \text{Le}\xi = \kappa / D\xi

savefreq property

calculate the save frequency in non dimensional time

scales property

return a Scales object used for converting between dimensional and non dimensional variables.

total_time property

calculate the total time in non dimensional units for the simulation

load(path) classmethod

load this object from a yaml configuration file.

Source code in seaice3p/params/dimensional/dimensional.py 
198
199
200
201
202
203
@classmethod
def load(cls, path):
    """load this object from a yaml configuration file."""
    with open(path, "r") as infile:
        yaml = infile.read()
    return from_yaml(cls, yaml)

save(directory)

save this object to a yaml file in the specified directory.

The name will be the name given with _dimensional appended to distinguish it from a saved non-dimensional configuration.

Source code in seaice3p/params/dimensional/dimensional.py 
190
191
192
193
194
195
196
def save(self, directory: Path):
    """save this object to a yaml file in the specified directory.

    The name will be the name given with _dimensional appended to distinguish it
    from a saved non-dimensional configuration."""
    with open(directory / f"{self.name}_dimensional.yml", "w") as outfile:
        outfile.write(to_yaml(self))

DimensionalPowerLawBubbleParams dataclass

Bases: DimensionalBaseBubbleParams

Source code in seaice3p/params/dimensional/bubble.py 
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalPowerLawBubbleParams(DimensionalBaseBubbleParams):
    bubble_distribution_power: float = 1.5
    minimum_bubble_radius: float = 1e-6
    maximum_bubble_radius: float = 1e-3

    @property
    def minimum_bubble_radius_scaled(self):
        r"""calculate the bubble radius divided by the pore scale

        .. math:: \Lambda = R_B / R_0

        """
        return self.minimum_bubble_radius / self.pore_radius

    @property
    def maximum_bubble_radius_scaled(self):
        r"""calculate the bubble radius divided by the pore scale

        .. math:: \Lambda = R_B / R_0

        """
        return self.maximum_bubble_radius / self.pore_radius

maximum_bubble_radius_scaled property

calculate the bubble radius divided by the pore scale

.. math:: \Lambda = R_B / R_0

minimum_bubble_radius_scaled property

calculate the bubble radius divided by the pore scale

.. math:: \Lambda = R_B / R_0

DimensionalRobinForcing dataclass

This forcing imposes a Robin boundary condition of the form surface_heat_flux=heat_transfer_coefficient * (restoring_temp - surface_temp)

Source code in seaice3p/params/dimensional/forcing.py 
164
165
166
167
168
169
170
171
172
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalRobinForcing:
    """This forcing imposes a Robin boundary condition of the form
    surface_heat_flux=heat_transfer_coefficient * (restoring_temp - surface_temp)
    """

    heat_transfer_coefficient: float = 6.3  # W/m2K
    restoring_temperature: float = -30  # deg C

DimensionalWaterParams dataclass

Source code in seaice3p/params/dimensional/water.py 
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalWaterParams:
    liquid_density: float = 1028  # kg/m3
    ice_density: float = 916  # kg/m3
    ocean_salinity: float = 34  # g/kg
    liquidus: LinearLiquidus | CubicLiquidus = LinearLiquidus()
    latent_heat: float = 334e3  # latent heat of fusion for ice in J/kg
    liquid_specific_heat_capacity: float = 4184  # J/kg degC
    solid_specific_heat_capacity: float = 2009  # J/kg degC
    liquid_thermal_conductivity: float = 0.54  # water thermal conductivity in W/m deg C
    solid_thermal_conductivity: float = 2.22  # ice thermal conductivity in W/m deg C
    snow_thermal_conductivity: float = 0.31  # snow thermal conductivity in W/m deg C
    snow_density: float = 150  # snow density kg/m3

    eddy_diffusivity: float = 0

    salt_diffusivity: float = 0  # molecular diffusivity of salt in water in m2/s
    # used to calculate Rayleigh number for convection and density contraction in liquid equation of state
    haline_contraction_coefficient: float = 7.5e-4  # 1/ppt

    # calculated from moreau et al 2014 value of kinematic viscosity for sewater 2.7e-6
    # dynamic liquid_viscosity = 2.7e-6 * liquid_density
    liquid_viscosity: float = 2.78e-3  # dynamic liquid viscosity in Pa.s

    @property
    def eutectic_salinity(self):
        if isinstance(self.liquidus, LinearLiquidus):
            return self.liquidus.eutectic_salinity
        if isinstance(self.liquidus, CubicLiquidus):
            return self.liquidus.get_liquidus_salinity(
                self.liquidus.eutectic_temperature
            )

        raise NotImplementedError

    @property
    def eutectic_temperature(self):
        if isinstance(self.liquidus, LinearLiquidus) or isinstance(
            self.liquidus, CubicLiquidus
        ):
            return self.liquidus.eutectic_temperature

        raise NotImplementedError

    @property
    def salinity_difference(self):
        r"""calculate difference between eutectic salinity and typical ocean salinity

        .. math:: \Delta S = S_E - S_i

        """
        return self.eutectic_salinity - self.ocean_salinity

    @property
    def ocean_freezing_temperature(self):
        """calculate salinity dependent freezing temperature using linear liquidus with
        ocean salinity

        .. math:: T_i = T_L(S_i) = T_E S_i / S_E

        or using a cubic fit for the liquidus curve

        """
        if isinstance(self.liquidus, LinearLiquidus):
            return (
                self.eutectic_temperature * self.ocean_salinity / self.eutectic_salinity
            )
        if isinstance(self.liquidus, CubicLiquidus):
            return self.liquidus.get_liquidus_temperature(self.ocean_salinity)

        raise NotImplementedError

    @property
    def temperature_difference(self):
        r"""calculate

        .. math:: \Delta T = T_i - T_E

        """
        return self.ocean_freezing_temperature - self.eutectic_temperature

    @property
    def concentration_ratio(self):
        r"""Calculate concentration ratio as

        .. math:: \mathcal{C} = S_i / \Delta S

        """
        return self.ocean_salinity / self.salinity_difference

    @property
    def stefan_number(self):
        r"""calculate Stefan number

        .. math:: \text{St} = L / c_p \Delta T

        """
        return self.latent_heat / (
            self.temperature_difference * self.liquid_specific_heat_capacity
        )

    @property
    def thermal_diffusivity(self):
        r"""Return thermal diffusivity in m2/s

        .. math:: \kappa = \frac{k}{\rho_l c_p}

        """
        return self.liquid_thermal_conductivity / (
            self.liquid_density * self.liquid_specific_heat_capacity
        )

    @property
    def conductivity_ratio(self):
        r"""Calculate the ratio of solid to liquid thermal conductivity

        .. math:: \lambda = \frac{k_s}{k_l}

        """
        return self.solid_thermal_conductivity / self.liquid_thermal_conductivity

    @property
    def specific_heat_ratio(self):
        r"""Calculate the ratio of solid to liquid specific heat capacities

        .. math:: \lambda = \frac{c_{p,s}}{c_{p,l}}

        """
        return self.solid_specific_heat_capacity / self.liquid_specific_heat_capacity

    @property
    def eddy_diffusivity_ratio(self):
        r"""Calculate the ratio of eddy diffusivity to thermal diffusivity in
        the liquid phase

        .. math:: \lambda = \frac{\kappa_\text{turbulent}}{\kappa_l}

        """
        return self.eddy_diffusivity / self.thermal_diffusivity

    @property
    def snow_conductivity_ratio(self):
        r"""Calculate the ratio of snow to liquid thermal conductivity

        .. math:: \lambda = \frac{k_{sn}}{k_l}

        """
        return self.snow_thermal_conductivity / self.liquid_thermal_conductivity

    @property
    def lewis_salt(self):
        r"""Calculate the lewis number for salt, return np.inf if there is no salt
        diffusion.

        .. math:: \text{Le}_S = \kappa / D_s

        """
        if self.salt_diffusivity == 0:
            return np.inf

        return self.thermal_diffusivity / self.salt_diffusivity

concentration_ratio property

Calculate concentration ratio as

.. math:: \mathcal{C} = S_i / \Delta S

conductivity_ratio property

Calculate the ratio of solid to liquid thermal conductivity

.. math:: \lambda = \frac{k_s}{k_l}

eddy_diffusivity_ratio property

Calculate the ratio of eddy diffusivity to thermal diffusivity in the liquid phase

.. math:: \lambda = \frac{\kappa_\text{turbulent}}{\kappa_l}

lewis_salt property

Calculate the lewis number for salt, return np.inf if there is no salt diffusion.

.. math:: \text{Le}_S = \kappa / D_s

ocean_freezing_temperature property

calculate salinity dependent freezing temperature using linear liquidus with ocean salinity

.. math:: T_i = T_L(S_i) = T_E S_i / S_E

or using a cubic fit for the liquidus curve

salinity_difference property

calculate difference between eutectic salinity and typical ocean salinity

.. math:: \Delta S = S_E - S_i

snow_conductivity_ratio property

Calculate the ratio of snow to liquid thermal conductivity

.. math:: \lambda = \frac{k_{sn}}{k_l}

specific_heat_ratio property

Calculate the ratio of solid to liquid specific heat capacities

.. math:: \lambda = \frac{c_{p,s}}{c_{p,l}}

stefan_number property

calculate Stefan number

.. math:: \text{St} = L / c_p \Delta T

temperature_difference property

calculate

.. math:: \Delta T = T_i - T_E

thermal_diffusivity property

Return thermal diffusivity in m2/s

.. math:: \kappa = \frac{k}{\rho_l c_p}

EQMPhysicalParams dataclass

Bases: BasePhysicalParams

non dimensional numbers for the mushy layer

Source code in seaice3p/params/physical.py 
44
45
46
47
@serde(type_check=coerce)
@dataclass(frozen=True)
class EQMPhysicalParams(BasePhysicalParams):
    """non dimensional numbers for the mushy layer"""

ERA5Forcing

Forcing parameters for simulation forced with atmospheric variables from reanalysis data in netCDF file located at data_path.

Never create this object directly but instead initialise from a dimensional simulation configuration as we must pass it the simulation timescale to correctly read the atmospheric variables from the netCDF file.

Source code in seaice3p/params/forcing.py 
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
@serde(type_check=coerce)
class ERA5Forcing:
    """Forcing parameters for simulation forced with atmospheric variables
    from reanalysis data in netCDF file located at data_path.

    Never create this object directly but instead initialise from a dimensional
    simulation configuration as we must pass it the simulation timescale to correctly
    read the atmospheric variables from the netCDF file.
    """

    data_path: Path
    start_date: str
    timescale_in_days: float
    forcing_data_file_keys: ERA5FileKeys = ERA5FileKeys()
    snow_density: Optional[float] = None
    SW_forcing: DimensionalSWForcing = DimensionalConstantSWForcing()
    LW_forcing: DimensionalLWForcing = DimensionalConstantLWForcing()
    turbulent_flux: DimensionalTurbulentFlux = DimensionalConstantTurbulentFlux()
    oil_heating: DimensionalOilHeating = DimensionalBackgroundOilHeating()

    def __post_init__(self):
        data = xr.open_dataset(self.data_path)
        DATES = getattr(data, self.forcing_data_file_keys.time).to_numpy()
        DIMLESS_TIMES = (1 / self.timescale_in_days) * np.array(
            [
                (date - np.datetime64(self.start_date)) / np.timedelta64(1, "D")
                for date in DATES
            ]
        )

        # convert to deg C
        T2M = (
            getattr(data, self.forcing_data_file_keys.temperature_at_2m_in_K).to_numpy()
            - 273.15
        )
        D2M = (
            getattr(data, self.forcing_data_file_keys.dewpoint_at_2m_in_K).to_numpy()
            - 273.15
        )

        LW = getattr(
            data, self.forcing_data_file_keys.longwave_radiation_in_W_m2
        ).to_numpy()
        SW = getattr(
            data, self.forcing_data_file_keys.shortwave_radiation_in_W_m2
        ).to_numpy()

        # convert to KPa
        ATM = (
            getattr(data, self.forcing_data_file_keys.surface_pressure_in_Pa).to_numpy()
            / 1e3
        )

        wind_key = self.forcing_data_file_keys.windspeed_at_2m_in_m_s
        if wind_key is None:
            WIND = np.full_like(DIMLESS_TIMES, self.turbulent_flux.windspeed)
        else:
            WIND = getattr(data, wind_key).to_numpy()

        # Calculate specific humidity in kg/kg from dewpoint temperature
        SPEC_HUM = _calculate_specific_humidity(ATM, D2M)

        snow_key = self.forcing_data_file_keys.snow_depth_in_m
        # if ERA5 standard short name for snow depth in m of water equivalent use snow
        # density to convert to m of snow
        if snow_key == "sd":
            if self.snow_density is None:
                raise ValueError("No snow density provided")
            SNOW_DEPTH = getattr(data, "sd").to_numpy() * (1000 / self.snow_density)

        # If snow key is another name assume snow depth is just in m of snow
        elif snow_key is not None:
            SNOW_DEPTH = getattr(data, snow_key).to_numpy()

        # If snow key is None assume no snow
        else:
            SNOW_DEPTH = np.zeros_like(DIMLESS_TIMES)

        # Provide functions to interpolate forcing data at non-dimensional times
        # during simulation
        self.get_2m_temp = partial(
            np.interp, xp=DIMLESS_TIMES, fp=T2M, left=np.nan, right=np.nan
        )
        self.get_LW = partial(
            np.interp, xp=DIMLESS_TIMES, fp=LW, left=np.nan, right=np.nan
        )
        self.get_SW = partial(
            np.interp, xp=DIMLESS_TIMES, fp=SW, left=np.nan, right=np.nan
        )
        self.get_ATM = partial(
            np.interp, xp=DIMLESS_TIMES, fp=ATM, left=np.nan, right=np.nan
        )
        self.get_spec_hum = partial(
            np.interp, xp=DIMLESS_TIMES, fp=SPEC_HUM, left=np.nan, right=np.nan
        )
        self.get_snow_depth = partial(
            np.interp, xp=DIMLESS_TIMES, fp=SNOW_DEPTH, left=np.nan, right=np.nan
        )
        self.get_windspeed = partial(
            np.interp, xp=DIMLESS_TIMES, fp=WIND, left=np.nan, right=np.nan
        )

FixedHeatFluxOceanForcing dataclass

Provides constant dimensionless ocean heat flux at the bottom of the domain and fixed gas saturation state.

Source code in seaice3p/params/ocean_forcing.py 
27
28
29
30
31
32
33
34
@serde(type_check=coerce)
@dataclass(frozen=True)
class FixedHeatFluxOceanForcing:
    """Provides constant dimensionless ocean heat flux at the bottom of the domain and fixed gas
    saturation state."""

    ocean_heat_flux: float = 1
    ocean_gas_sat: float = 1.0

FixedTempOceanForcing dataclass

Fixed temperature and gas saturation ocean boundary condition

Source code in seaice3p/params/ocean_forcing.py 
18
19
20
21
22
23
24
@serde(type_check=coerce)
@dataclass(frozen=True)
class FixedTempOceanForcing:
    """Fixed temperature and gas saturation ocean boundary condition"""

    ocean_temp: float = 0.1
    ocean_gas_sat: float = 1.0

MonoBubbleParams dataclass

Bases: BaseBubbleParams

Parameters for population of identical spherical bubbles.

Source code in seaice3p/params/bubble.py 
24
25
26
27
28
29
@serde(type_check=coerce)
@dataclass(frozen=True)
class MonoBubbleParams(BaseBubbleParams):
    """Parameters for population of identical spherical bubbles."""

    bubble_radius_scaled: float = 1.0

MonthlyHeatFluxOceanForcing

Provides constant dimensionless ocean heat flux at the bottom of the domain in each month

and ocean gas saturation state.

Proivde an average monthly ocean heat flux with the entries i=0, 1, 2, 3, ...., 11 in the tuple corresponding to the months January, February, March, April, ...., December

Parameters:

Name Type Description Default
monthly_ocean_heat_flux

Tuple of dimensionless ocean heat flux values in

 required
Source code in seaice3p/params/ocean_forcing.py 
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
@serde(type_check=coerce)
class MonthlyHeatFluxOceanForcing:
    """Provides constant dimensionless ocean heat flux at the bottom of the domain in
    each month

    and ocean gas saturation state.

    Proivde an average monthly ocean heat flux with the entries
    i=0, 1, 2, 3, ...., 11
    in the tuple corresponding to the months
    January, February, March, April, ...., December

    Args:
        monthly_ocean_heat_flux: Tuple of dimensionless ocean heat flux values in
        each month
    """

    start_date: str
    timescale_in_days: float
    monthly_ocean_heat_flux: Tuple[float, ...] = (
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
    )
    ocean_gas_sat: float = 1.0

    def __post_init__(self):
        # Provide functions to interpolate day of year to monthly ocean heat flux
        self._interpolate_heat_flux = partial(
            np.interp,
            xp=np.array([15, 46, 74, 105, 135, 166, 196, 227, 258, 288, 319, 349]),
            fp=np.array(self.monthly_ocean_heat_flux),
            period=365,
        )

    def get_ocean_heat_flux(self, simulation_time: float) -> float:
        start_datetime = datetime.strptime(self.start_date, "%Y-%m-%d")
        current_datetime = start_datetime + timedelta(
            days=self.timescale_in_days * simulation_time
        )
        current_day = (
            current_datetime - datetime(start_datetime.year, 1, 1)
        ).total_seconds() / 86400
        return self._interpolate_heat_flux(current_day)

NoBrineConvection dataclass

No brine convection

Source code in seaice3p/params/dimensional/convection.py 
5
6
7
8
@serde(type_check=coerce)
@dataclass(frozen=True)
class NoBrineConvection:
    """No brine convection"""

NumericalParams dataclass

parameters needed for discretisation and choice of numerical method

Source code in seaice3p/params/dimensional/numerical.py 
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
@serde(type_check=coerce)
@dataclass(frozen=True)
class NumericalParams:
    """parameters needed for discretisation and choice of numerical method"""

    I: int = 50
    regularisation: float = 1e-6
    solver_choice: str = "RK23"  # scipy.integrate.solve_IVP solver choice

    @property
    def step(self):
        return 1 / self.I

OilInitialConditions dataclass

values for bottom (ocean) boundary

Source code in seaice3p/params/initial_conditions.py 
12
13
14
15
16
17
18
19
20
21
22
23
@serde(type_check=coerce)
@dataclass(frozen=True)
class OilInitialConditions:
    """values for bottom (ocean) boundary"""

    # Non dimensional parameters for summer initial conditions
    initial_ice_depth: float = 0.5
    initial_ocean_temperature: float = -0.05
    initial_ice_temperature: float = -0.1
    initial_oil_volume_fraction: float = 1e-7
    initial_ice_bulk_salinity: float = -0.1
    initial_oil_free_depth: float = 0

PowerLawBubbleParams dataclass

Bases: BaseBubbleParams

Parameters for population of bubbles following a power law size distribution between a minimum and maximum radius.

Source code in seaice3p/params/bubble.py 
32
33
34
35
36
37
38
39
40
41
@serde(type_check=coerce)
@dataclass(frozen=True)
class PowerLawBubbleParams(BaseBubbleParams):
    """Parameters for population of bubbles following a power law size distribution
    between a minimum and maximum radius.
    """

    bubble_distribution_power: float = 1.5
    minimum_bubble_radius_scaled: float = 1e-3
    maximum_bubble_radius_scaled: float = 1

RJW14Params dataclass

Parameters for the RJW14 parameterisation of brine convection

Source code in seaice3p/params/convection.py 
 6
 7
 8
 9
10
11
12
13
14
15
16
@serde(type_check=coerce)
@dataclass(frozen=True)
class RJW14Params:
    """Parameters for the RJW14 parameterisation of brine convection"""

    Rayleigh_salt: float = 44105
    Rayleigh_critical: float = 2.9
    convection_strength: float = 0.13
    couple_bubble_to_horizontal_flow: bool = False
    couple_bubble_to_vertical_flow: bool = False
    advective_heat_flux_in_ocean: bool = True

RadForcing dataclass

Forcing parameters for radiative transfer simulation with oil drops

we have not implemented the non-dimensionalisation for these parameters yet and so we just pass the dimensional values directly to the simulation

Source code in seaice3p/params/forcing.py 
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
@serde(type_check=coerce)
@dataclass(frozen=True)
class RadForcing:
    """Forcing parameters for radiative transfer simulation with oil drops

    we have not implemented the non-dimensionalisation for these parameters yet
    and so we just pass the dimensional values directly to the simulation"""

    SW_forcing: DimensionalSWForcing = DimensionalConstantSWForcing()
    LW_forcing: DimensionalLWForcing = DimensionalConstantLWForcing()
    turbulent_flux: DimensionalTurbulentFlux = DimensionalConstantTurbulentFlux()
    oil_heating: DimensionalOilHeating = DimensionalBackgroundOilHeating()

RobinForcing dataclass

Dimensionless forcing parameters for Robin boundary condition

Source code in seaice3p/params/forcing.py 
221
222
223
224
225
226
227
@serde(type_check=coerce)
@dataclass(frozen=True)
class RobinForcing:
    """Dimensionless forcing parameters for Robin boundary condition"""

    biot: float = 12
    restoring_temperature: float = -1.3

Scales dataclass

Source code in seaice3p/params/convert.py 
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
@serde(type_check=coerce)
@dataclass(frozen=True)
class Scales:
    lengthscale: float  # domain height in m
    thermal_diffusivity: float  # m2/s
    liquid_thermal_conductivity: float  # W/m deg C
    ocean_salinity: float  # g/kg
    salinity_difference: float  # g/kg
    ocean_freezing_temperature: float  # deg C
    temperature_difference: float  # deg C
    gas_density: float  # kg/m3
    liquid_density: float  # kg/m3
    ice_density: float  # kg/m3
    saturation_concentration: float  # kg(gas)/kg(liquid)
    pore_radius: float  # m
    haline_contraction_coefficient: float  # 1/ppt

    @property
    def time_scale(self):
        """in days"""
        return SECONDS_TO_DAYS * self.lengthscale**2 / self.thermal_diffusivity

    @property
    def velocity_scale(self):
        """in m /day"""
        return self.lengthscale / self.time_scale

    def convert_from_dimensional_temperature(self, dimensional_temperature):
        """Non dimensionalise temperature in deg C"""
        return (
            dimensional_temperature - self.ocean_freezing_temperature
        ) / self.temperature_difference

    def convert_to_dimensional_temperature(self, temperature):
        """get temperature in deg C from non dimensional temperature"""
        return (
            self.temperature_difference * temperature + self.ocean_freezing_temperature
        )

    def convert_from_dimensional_grid(self, dimensional_grid):
        """Non dimensionalise domain depths in meters"""
        return dimensional_grid / self.lengthscale

    def convert_to_dimensional_grid(self, grid):
        """Get domain depths in meters from non dimensional values"""
        return self.lengthscale * grid

    def convert_from_dimensional_time(self, dimensional_time):
        """Non dimensionalise time in days"""
        return dimensional_time / self.time_scale

    def convert_to_dimensional_time(self, time):
        """Convert non dimensional time into time in days since start of simulation"""
        return self.time_scale * time

    def convert_from_dimensional_bulk_salinity(self, dimensional_bulk_salinity):
        """Non dimensionalise bulk salinity in g/kg"""
        return (
            dimensional_bulk_salinity - self.ocean_salinity
        ) / self.salinity_difference

    def convert_to_dimensional_bulk_salinity(self, bulk_salinity):
        """Convert non dimensional bulk salinity to g/kg"""
        return self.salinity_difference * bulk_salinity + self.ocean_salinity

    def convert_from_dimensional_bulk_gas(self, dimensional_bulk_gas):
        """Non dimensionalise bulk gas content in kg/m3"""
        return dimensional_bulk_gas / self.gas_density

    def convert_to_dimensional_bulk_gas(self, bulk_gas):
        """Convert dimensionless bulk gas content to kg/m3"""
        return self.gas_density * bulk_gas

    def convert_dimensional_bulk_air_to_argon_content(self, dimensional_bulk_gas):
        """Convert kg/m3 of air to micromole of Argon per Liter of ice"""
        mass_ratio_of_argon_in_air = 0.01288
        micromoles_of_argon_in_a_kilogram_of_argon = 1 / (3.9948e-8)
        liters_in_a_meter_cubed = 1e3
        return (
            dimensional_bulk_gas
            * mass_ratio_of_argon_in_air
            * micromoles_of_argon_in_a_kilogram_of_argon
            / liters_in_a_meter_cubed
        )

    def convert_from_dimensional_dissolved_gas(self, dimensional_dissolved_gas):
        """convert from dissolved gas in kg(gas)/kg(liquid) to dimensionless"""
        return dimensional_dissolved_gas / self.saturation_concentration

    def convert_to_dimensional_dissolved_gas(self, dissolved_gas):
        """convert from non dimensional dissolved gas to dimensional dissolved gas in
        kg(gas)/kg(liquid)"""
        return self.saturation_concentration * dissolved_gas

    def convert_from_dimensional_heating(self, dimensional_heating):
        """convert from heating rate in W/m3 to dimensionless units"""
        return (
            dimensional_heating
            * self.lengthscale**2
            / (self.liquid_thermal_conductivity * self.temperature_difference)
        )

    def convert_from_dimensional_heat_flux(self, dimensional_heat_flux):
        """convert from heat flux in W/m2 to dimensionless units"""
        return (
            dimensional_heat_flux
            * self.lengthscale
            / (self.liquid_thermal_conductivity * self.temperature_difference)
        )

    def convert_to_dimensional_heat_flux(self, heat_flux):
        """convert from dimensionless heat flux to heat flux in W/m2"""
        return (
            heat_flux
            * (self.liquid_thermal_conductivity * self.temperature_difference)
            / self.lengthscale
        )

time_scale property

in days

velocity_scale property

in m /day

convert_dimensional_bulk_air_to_argon_content(dimensional_bulk_gas)

Convert kg/m3 of air to micromole of Argon per Liter of ice

Source code in seaice3p/params/convert.py 
81
82
83
84
85
86
87
88
89
90
91
def convert_dimensional_bulk_air_to_argon_content(self, dimensional_bulk_gas):
    """Convert kg/m3 of air to micromole of Argon per Liter of ice"""
    mass_ratio_of_argon_in_air = 0.01288
    micromoles_of_argon_in_a_kilogram_of_argon = 1 / (3.9948e-8)
    liters_in_a_meter_cubed = 1e3
    return (
        dimensional_bulk_gas
        * mass_ratio_of_argon_in_air
        * micromoles_of_argon_in_a_kilogram_of_argon
        / liters_in_a_meter_cubed
    )

convert_from_dimensional_bulk_gas(dimensional_bulk_gas)

Non dimensionalise bulk gas content in kg/m3

Source code in seaice3p/params/convert.py 
73
74
75
def convert_from_dimensional_bulk_gas(self, dimensional_bulk_gas):
    """Non dimensionalise bulk gas content in kg/m3"""
    return dimensional_bulk_gas / self.gas_density

convert_from_dimensional_bulk_salinity(dimensional_bulk_salinity)

Non dimensionalise bulk salinity in g/kg

Source code in seaice3p/params/convert.py 
63
64
65
66
67
def convert_from_dimensional_bulk_salinity(self, dimensional_bulk_salinity):
    """Non dimensionalise bulk salinity in g/kg"""
    return (
        dimensional_bulk_salinity - self.ocean_salinity
    ) / self.salinity_difference

convert_from_dimensional_dissolved_gas(dimensional_dissolved_gas)

convert from dissolved gas in kg(gas)/kg(liquid) to dimensionless

Source code in seaice3p/params/convert.py 
93
94
95
def convert_from_dimensional_dissolved_gas(self, dimensional_dissolved_gas):
    """convert from dissolved gas in kg(gas)/kg(liquid) to dimensionless"""
    return dimensional_dissolved_gas / self.saturation_concentration

convert_from_dimensional_grid(dimensional_grid)

Non dimensionalise domain depths in meters

Source code in seaice3p/params/convert.py 
47
48
49
def convert_from_dimensional_grid(self, dimensional_grid):
    """Non dimensionalise domain depths in meters"""
    return dimensional_grid / self.lengthscale

convert_from_dimensional_heat_flux(dimensional_heat_flux)

convert from heat flux in W/m2 to dimensionless units

Source code in seaice3p/params/convert.py 
110
111
112
113
114
115
116
def convert_from_dimensional_heat_flux(self, dimensional_heat_flux):
    """convert from heat flux in W/m2 to dimensionless units"""
    return (
        dimensional_heat_flux
        * self.lengthscale
        / (self.liquid_thermal_conductivity * self.temperature_difference)
    )

convert_from_dimensional_heating(dimensional_heating)

convert from heating rate in W/m3 to dimensionless units

Source code in seaice3p/params/convert.py 
102
103
104
105
106
107
108
def convert_from_dimensional_heating(self, dimensional_heating):
    """convert from heating rate in W/m3 to dimensionless units"""
    return (
        dimensional_heating
        * self.lengthscale**2
        / (self.liquid_thermal_conductivity * self.temperature_difference)
    )

convert_from_dimensional_temperature(dimensional_temperature)

Non dimensionalise temperature in deg C

Source code in seaice3p/params/convert.py 
35
36
37
38
39
def convert_from_dimensional_temperature(self, dimensional_temperature):
    """Non dimensionalise temperature in deg C"""
    return (
        dimensional_temperature - self.ocean_freezing_temperature
    ) / self.temperature_difference

convert_from_dimensional_time(dimensional_time)

Non dimensionalise time in days

Source code in seaice3p/params/convert.py 
55
56
57
def convert_from_dimensional_time(self, dimensional_time):
    """Non dimensionalise time in days"""
    return dimensional_time / self.time_scale

convert_to_dimensional_bulk_gas(bulk_gas)

Convert dimensionless bulk gas content to kg/m3

Source code in seaice3p/params/convert.py 
77
78
79
def convert_to_dimensional_bulk_gas(self, bulk_gas):
    """Convert dimensionless bulk gas content to kg/m3"""
    return self.gas_density * bulk_gas

convert_to_dimensional_bulk_salinity(bulk_salinity)

Convert non dimensional bulk salinity to g/kg

Source code in seaice3p/params/convert.py 
69
70
71
def convert_to_dimensional_bulk_salinity(self, bulk_salinity):
    """Convert non dimensional bulk salinity to g/kg"""
    return self.salinity_difference * bulk_salinity + self.ocean_salinity

convert_to_dimensional_dissolved_gas(dissolved_gas)

convert from non dimensional dissolved gas to dimensional dissolved gas in kg(gas)/kg(liquid)

Source code in seaice3p/params/convert.py 
 97
 98
 99
100
def convert_to_dimensional_dissolved_gas(self, dissolved_gas):
    """convert from non dimensional dissolved gas to dimensional dissolved gas in
    kg(gas)/kg(liquid)"""
    return self.saturation_concentration * dissolved_gas

convert_to_dimensional_grid(grid)

Get domain depths in meters from non dimensional values

Source code in seaice3p/params/convert.py 
51
52
53
def convert_to_dimensional_grid(self, grid):
    """Get domain depths in meters from non dimensional values"""
    return self.lengthscale * grid

convert_to_dimensional_heat_flux(heat_flux)

convert from dimensionless heat flux to heat flux in W/m2

Source code in seaice3p/params/convert.py 
118
119
120
121
122
123
124
def convert_to_dimensional_heat_flux(self, heat_flux):
    """convert from dimensionless heat flux to heat flux in W/m2"""
    return (
        heat_flux
        * (self.liquid_thermal_conductivity * self.temperature_difference)
        / self.lengthscale
    )

convert_to_dimensional_temperature(temperature)

get temperature in deg C from non dimensional temperature

Source code in seaice3p/params/convert.py 
41
42
43
44
45
def convert_to_dimensional_temperature(self, temperature):
    """get temperature in deg C from non dimensional temperature"""
    return (
        self.temperature_difference * temperature + self.ocean_freezing_temperature
    )

convert_to_dimensional_time(time)

Convert non dimensional time into time in days since start of simulation

Source code in seaice3p/params/convert.py 
59
60
61
def convert_to_dimensional_time(self, time):
    """Convert non dimensional time into time in days since start of simulation"""
    return self.time_scale * time

UniformInitialConditions dataclass

values for bottom (ocean) boundary

Source code in seaice3p/params/dimensional/initial_conditions.py 
6
7
8
9
@serde(type_check=coerce)
@dataclass(frozen=True)
class UniformInitialConditions:
    """values for bottom (ocean) boundary"""

YearlyForcing dataclass

Yearly sinusoidal temperature forcing

Source code in seaice3p/params/forcing.py 
45
46
47
48
49
50
51
52
@serde(type_check=coerce)
@dataclass(frozen=True)
class YearlyForcing:
    """Yearly sinusoidal temperature forcing"""

    offset: float = -1.0
    amplitude: float = 0.75
    period: float = 4.0

get_config(dimensional_params)

Return a Config object for the simulation.

physical parameters and Darcy law parameters are calculated from the dimensional input. You can modify the numerical parameters and boundary conditions and forcing provided for the simulation.

Source code in seaice3p/params/params.py 
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
def get_config(dimensional_params: DimensionalParams) -> Config:
    """Return a Config object for the simulation.

    physical parameters and Darcy law parameters are calculated from the dimensional
    input. You can modify the numerical parameters and boundary conditions and
    forcing provided for the simulation."""
    physical_params = get_dimensionless_physical_params(dimensional_params)
    initial_conditions_config = get_dimensionless_initial_conditions_config(
        dimensional_params
    )
    brine_convection_params = get_dimensionless_brine_convection_params(
        dimensional_params
    )
    bubble_params = get_dimensionless_bubble_params(dimensional_params)
    forcing_config = get_dimensionless_forcing_config(dimensional_params)
    ocean_forcing_config = get_dimensionless_ocean_forcing_config(dimensional_params)
    return Config(
        name=dimensional_params.name,
        physical_params=physical_params,
        initial_conditions_config=initial_conditions_config,
        brine_convection_params=brine_convection_params,
        bubble_params=bubble_params,
        forcing_config=forcing_config,
        ocean_forcing_config=ocean_forcing_config,
        numerical_params=dimensional_params.numerical_params,
        scales=dimensional_params.scales,
        total_time=dimensional_params.total_time,
        savefreq=dimensional_params.savefreq,
    )

enthalpy_method

enthalpy_method

Module containing enthalpy method to calculate state variables from bulk enthalpy, bulk salinity and bulk gas.

phase_boundaries

Module for calculating the phase boundaries needed for the enthalpy method.

calculates the phase boundaries neglecting the gas fraction so that

.. math:: \phi_s + \phi_l = 1

equations

RJW14

Module to calculate the sink terms for conservation equations when using the Rees Jones and Worster 2014 brine drainage parameterisation.

These terms represent loss through the brine channels and need to be added in the convecting region when using this parameterisation

brine_channel_sink_terms

_DISEQ_brine_convection_sink(state_BCs, cfg, grids)

TODO: check the sink terms for bulk_dissolved_gas and gas fraction

For now neglect the coupling of bubbles to the horizontal or vertical flow

Source code in seaice3p/equations/RJW14/brine_channel_sink_terms.py 
45
46
47
48
49
50
51
52
53
54
def _DISEQ_brine_convection_sink(state_BCs: DISEQStateBCs, cfg, grids) -> NDArray:
    """TODO: check the sink terms for bulk_dissolved_gas and gas fraction

    For now neglect the coupling of bubbles to the horizontal or vertical flow
    """
    heat_sink = _calculate_heat_sink(state_BCs, cfg, grids)
    salt_sink = _calculate_salt_sink(state_BCs, cfg, grids)
    bulk_dissolved_gas_sink = _calculate_bulk_dissolved_gas_sink(state_BCs, cfg, grids)
    gas_fraction_sink = np.zeros_like(heat_sink)
    return np.hstack((heat_sink, salt_sink, bulk_dissolved_gas_sink, gas_fraction_sink))
_EQM_brine_convection_sink(state_BCs, cfg, grids)

TODO: check the sink terms for bulk_dissolved_gas and gas fraction

For now neglect the coupling of bubbles to the horizontal or vertical flow

Source code in seaice3p/equations/RJW14/brine_channel_sink_terms.py 
34
35
36
37
38
39
40
41
42
def _EQM_brine_convection_sink(state_BCs: EQMStateBCs, cfg, grids) -> NDArray:
    """TODO: check the sink terms for bulk_dissolved_gas and gas fraction

    For now neglect the coupling of bubbles to the horizontal or vertical flow
    """
    heat_sink = _calculate_heat_sink(state_BCs, cfg, grids)
    salt_sink = _calculate_salt_sink(state_BCs, cfg, grids)
    gas_sink = _calculate_gas_sink(state_BCs, cfg, grids)
    return np.hstack((heat_sink, salt_sink, gas_sink))
_calculate_bulk_dissolved_gas_sink(state_BCs, cfg, grids)

This is for the DISEQ model

Source code in seaice3p/equations/RJW14/brine_channel_sink_terms.py 
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
def _calculate_bulk_dissolved_gas_sink(state_BCs, cfg: Config, grids):
    """This is for the DISEQ model"""
    liquid_fraction = state_BCs.liquid_fraction[1:-1]
    liquid_salinity = state_BCs.liquid_salinity[1:-1]
    dissolved_gas = state_BCs.dissolved_gas[1:-1]
    center_grid = grids.centers
    edge_grid = grids.edges

    if isinstance(cfg.brine_convection_params, NoBrineConvection):
        return np.zeros_like(liquid_fraction)

    sink = calculate_brine_channel_sink(
        liquid_fraction, liquid_salinity, center_grid, edge_grid, cfg
    )

    return sink * cfg.physical_params.expansion_coefficient * dissolved_gas
_calculate_gas_sink(state_BCs, cfg, grids)

This is for the EQM model

TODO: fix bug in bubble coupling to flow

Source code in seaice3p/equations/RJW14/brine_channel_sink_terms.py 
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
def _calculate_gas_sink(state_BCs, cfg: Config, grids):
    """This is for the EQM model

    TODO: fix bug in bubble coupling to flow
    """
    liquid_fraction = state_BCs.liquid_fraction[1:-1]
    liquid_salinity = state_BCs.liquid_salinity[1:-1]
    dissolved_gas = state_BCs.dissolved_gas[1:-1]
    gas_fraction = state_BCs.gas_fraction[1:-1]
    center_grid = grids.centers
    edge_grid = grids.edges

    if isinstance(cfg.brine_convection_params, NoBrineConvection):
        return np.zeros_like(liquid_fraction)

    sink = calculate_brine_channel_sink(
        liquid_fraction, liquid_salinity, center_grid, edge_grid, cfg
    )

    dissolved_gas_term = cfg.physical_params.expansion_coefficient * dissolved_gas

    if cfg.brine_convection_params.couple_bubble_to_horizontal_flow:
        if isinstance(cfg.bubble_params, MonoBubbleParams):
            lag_factor = calculate_mono_lag_factor(state_BCs.liquid_fraction, cfg)
        elif isinstance(cfg.bubble_params, PowerLawBubbleParams):
            lag_factor = calculate_power_law_lag_factor(state_BCs.liquid_fraction, cfg)
        else:
            raise NotImplementedError

        bubble_term = 2 * gas_fraction * geometric(lag_factor) / liquid_fraction
        is_frozen_solid = liquid_fraction == 0.0
        bubble_term[is_frozen_solid] = 0
    else:
        bubble_term = np.zeros_like(liquid_fraction)

    return sink * (dissolved_gas_term + bubble_term)

brine_drainage

Module to calculate the Rees Jones and Worster 2014 parameterisation for brine convection velocity and the strenght of the sink term.

Exports the functions:

calculate_brine_convection_liquid_velocity To be used in velocities module when using brine convection parameterisation.

calculate_brine_channel_sink To be used to add sink terms to conservation equations when using brine convection parameterisation.

calculate_Rayleigh(cell_centers, edge_grid, liquid_salinity, liquid_fraction, cfg)

Calculate the local Rayleigh number for brine convection as

.. math:: \text{Ra}(z) = \text{Ra}_S K(z) (z+h) \Theta_l

:param cell_centers: The vertical coordinates of cell centers. :type cell_centers: Numpy Array shape (I,) :param edge_grid: The vertical coordinate positions of the edge grid. :type edge_grid: Numpy Array (size I+1) :param liquid_salinity: liquid salinity on center grid :type liquid_salinity: Numpy Array shape (I,) :param liquid_fraction: liquid fraction on center grid :type liquid_fraction: Numpy Array (size I) :param cfg: Configuration object for the simulation. :type cfg: seaice3p.params.Config :return: Array of shape (I,) of Rayleigh number at cell centers

Source code in seaice3p/equations/RJW14/brine_drainage.py 
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
def calculate_Rayleigh(
    cell_centers, edge_grid, liquid_salinity, liquid_fraction, cfg: Config
):
    r"""Calculate the local Rayleigh number for brine convection as

    .. math:: \text{Ra}(z) = \text{Ra}_S K(z) (z+h) \Theta_l

    :param cell_centers: The vertical coordinates of cell centers.
    :type cell_centers: Numpy Array shape (I,)
    :param edge_grid: The vertical coordinate positions of the edge grid.
    :type edge_grid: Numpy Array (size I+1)
    :param liquid_salinity: liquid salinity on center grid
    :type liquid_salinity: Numpy Array shape (I,)
    :param liquid_fraction: liquid fraction on center grid
    :type liquid_fraction: Numpy Array (size I)
    :param cfg: Configuration object for the simulation.
    :type cfg: seaice3p.params.Config
    :return: Array of shape (I,) of Rayleigh number at cell centers
    """
    Rayleigh_salt = cfg.brine_convection_params.Rayleigh_salt
    ice_depth = calculate_ice_ocean_boundary_depth(liquid_fraction, edge_grid)
    averaged_permeabilities = np.array(
        [
            calculate_integrated_mean_permeability(
                z=z,
                liquid_fraction=liquid_fraction,
                ice_depth=ice_depth,
                cell_centers=cell_centers,
                cfg=cfg,
            )
            for z in cell_centers
        ]
    )
    return (
        Rayleigh_salt
        * (ice_depth + cell_centers)
        * averaged_permeabilities
        * liquid_salinity
    )
calculate_brine_channel_sink(liquid_fraction, liquid_salinity, center_grid, edge_grid, cfg)

Calculate the sink term due to brine channels.

.. math:: \text{sink} = \mathcal{A}

in the convecting region. Zero elsewhere.

NOTE: If no ice is present or if no convecting region exists returns zero

:param liquid_fraction: liquid fraction on center grid :type liquid_fraction: Numpy Array of shape (I,) :param liquid_salinity: liquid salinity on center grid :type liquid_salinity: Numpy Array of shape (I,) :param center_grid: vertical coordinate of center grid :type center_grid: Numpy Array of shape (I,) :param edge_grid: Vertical coordinates of cell edges :type edge_grid: Numpy Array of shape (I+1,) :param cfg: Configuration object for the simulation. :type cfg: seaice3p.params.Config :return: Strength of the sink term due to brine channels on the center grid.

Source code in seaice3p/equations/RJW14/brine_drainage.py 
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
def calculate_brine_channel_sink(
    liquid_fraction, liquid_salinity, center_grid, edge_grid, cfg: Config
):
    r"""Calculate the sink term due to brine channels.

    .. math:: \text{sink} = \mathcal{A}

    in the convecting region. Zero elsewhere.

    NOTE: If no ice is present or if no convecting region exists returns zero

    :param liquid_fraction: liquid fraction on center grid
    :type liquid_fraction: Numpy Array of shape (I,)
    :param liquid_salinity: liquid salinity on center grid
    :type liquid_salinity: Numpy Array of shape (I,)
    :param center_grid: vertical coordinate of center grid
    :type center_grid: Numpy Array of shape (I,)
    :param edge_grid: Vertical coordinates of cell edges
    :type edge_grid: Numpy Array of shape (I+1,)
    :param cfg: Configuration object for the simulation.
    :type cfg: seaice3p.params.Config
    :return: Strength of the sink term due to brine channels on the center grid.
    """
    ice_depth = calculate_ice_ocean_boundary_depth(liquid_fraction, edge_grid)
    Rayleigh_number = calculate_Rayleigh(
        center_grid, edge_grid, liquid_salinity, liquid_fraction, cfg
    )
    convecting_region_height = get_convecting_region_height(
        Rayleigh_number, edge_grid, cfg
    )
    brine_channel_strength = calculate_brine_channel_strength(
        Rayleigh_number, ice_depth, convecting_region_height, cfg
    )

    sink = np.zeros_like(center_grid)

    # No ice present
    if ice_depth == 0:
        return sink

    # ice present but no convection occuring
    if np.isnan(convecting_region_height):
        return sink

    # Make liquid vertical velocity continuous at bottom of the ice

    is_convecting_ice = (center_grid >= -ice_depth) & (
        center_grid <= convecting_region_height
    )
    is_liquid = center_grid < -ice_depth

    sink[is_convecting_ice] = brine_channel_strength
    sink[is_liquid] = 0

    return sink
calculate_brine_channel_strength(Rayleigh_number, ice_depth, convecting_region_height, cfg)

Calculate the brine channel strength in the convecting region as

.. math:: \mathcal{A} = \frac{\alpha \text{Ra}_e}{(h+z_c)^2}

the effective Rayleigh number multiplied by a tuning parameter (Rees Jones and Worster 2014) over the convecting region thickness squared.

:param Rayleigh_number: local Rayleigh number on center grid :type Rayleigh_number: Numpy Array of shape (I,) :param ice_depth: depth of ice (positive) :type ice_depth: float :param convecting_region_height: position of the convecting region boundary (negative) :type convecting_region_height: float :param cfg: Configuration object for the simulation. :type cfg: seaice3p.params.Config :return: Brine channel strength parameter

Source code in seaice3p/equations/RJW14/brine_drainage.py 
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
def calculate_brine_channel_strength(
    Rayleigh_number, ice_depth, convecting_region_height, cfg: Config
):
    r"""Calculate the brine channel strength in the convecting region as

    .. math:: \mathcal{A} = \frac{\alpha \text{Ra}_e}{(h+z_c)^2}

    the effective Rayleigh number multiplied by a tuning parameter (Rees Jones and
    Worster 2014) over the convecting region thickness squared.

    :param Rayleigh_number: local Rayleigh number on center grid
    :type Rayleigh_number: Numpy Array of shape (I,)
    :param ice_depth: depth of ice (positive)
    :type ice_depth: float
    :param convecting_region_height: position of the convecting region boundary (negative)
    :type convecting_region_height: float
    :param cfg: Configuration object for the simulation.
    :type cfg: seaice3p.params.Config
    :return: Brine channel strength parameter
    """
    convection_strength = cfg.brine_convection_params.convection_strength
    if ice_depth == 0:
        return 0

    if np.isnan(convecting_region_height):
        return 0

    convecting_layer_thickness = ice_depth + convecting_region_height
    effective_Rayleigh = get_effective_Rayleigh_number(Rayleigh_number, cfg)
    return convection_strength * effective_Rayleigh / convecting_layer_thickness**2
calculate_brine_convection_liquid_velocity(liquid_fraction, liquid_salinity, center_grid, edge_grid, cfg)

Calculate the vertical liquid Darcy velocity from Rees Jones and Worster 2014

.. math:: W_l = \mathcal{A} (z_c - z)

in the convecting region. The velocity is stagnant above the convecting region. The velocity is constant in the liquid region and continuous at the interface.

NOTE: If no ice is present or if no convecting region exists returns zero velocity

:param liquid_fraction: liquid fraction on center grid :type liquid_fraction: Numpy Array of shape (I,) :param liquid_salinity: liquid salinity on center grid :type liquid_salinity: Numpy Array of shape (I,) :param center_grid: vertical coordinate of center grid :type center_grid: Numpy Array of shape (I,) :param edge_grid: Vertical coordinates of cell edges :type edge_grid: Numpy Array of shape (I+1,) :param cfg: Configuration object for the simulation. :type cfg: seaice3p.params.Config :return: Liquid darcy velocity on the edge grid.

Source code in seaice3p/equations/RJW14/brine_drainage.py 
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
def calculate_brine_convection_liquid_velocity(
    liquid_fraction, liquid_salinity, center_grid, edge_grid, cfg: Config
):
    r"""Calculate the vertical liquid Darcy velocity from Rees Jones and Worster 2014

    .. math:: W_l = \mathcal{A} (z_c - z)

    in the convecting region. The velocity is stagnant above the convecting region.
    The velocity is constant in the liquid region and continuous at the interface.

    NOTE: If no ice is present or if no convecting region exists returns zero velocity

    :param liquid_fraction: liquid fraction on center grid
    :type liquid_fraction: Numpy Array of shape (I,)
    :param liquid_salinity: liquid salinity on center grid
    :type liquid_salinity: Numpy Array of shape (I,)
    :param center_grid: vertical coordinate of center grid
    :type center_grid: Numpy Array of shape (I,)
    :param edge_grid: Vertical coordinates of cell edges
    :type edge_grid: Numpy Array of shape (I+1,)
    :param cfg: Configuration object for the simulation.
    :type cfg: seaice3p.params.Config
    :return: Liquid darcy velocity on the edge grid.
    """
    ice_depth = calculate_ice_ocean_boundary_depth(liquid_fraction, edge_grid)
    Rayleigh_number = calculate_Rayleigh(
        center_grid, edge_grid, liquid_salinity, liquid_fraction, cfg
    )
    convecting_region_height = get_convecting_region_height(
        Rayleigh_number, edge_grid, cfg
    )
    brine_channel_strength = calculate_brine_channel_strength(
        Rayleigh_number, ice_depth, convecting_region_height, cfg
    )

    Wl = np.zeros_like(edge_grid)

    # No ice present
    if ice_depth == 0:
        return Wl

    # ice present but no convection occuring
    if np.isnan(convecting_region_height):
        return Wl

    # Make liquid vertical velocity continuous at bottom of the ice
    ocean_velocity = brine_channel_strength * (ice_depth + convecting_region_height)

    is_convecting_ice = (edge_grid >= -ice_depth) & (
        edge_grid <= convecting_region_height
    )
    is_liquid = edge_grid < -ice_depth

    Wl[is_convecting_ice] = brine_channel_strength * (
        convecting_region_height - edge_grid[is_convecting_ice]
    )
    Wl[is_liquid] = ocean_velocity

    return Wl
calculate_integrated_mean_permeability(z, liquid_fraction, ice_depth, cell_centers, cfg)

Calculate the harmonic mean permeability from the base of the ice up to the cell containing the specified z value using the expression of ReesJones2014.

.. math:: K(z) = (\frac{1}{h+z}\int_{-h}^{z} \frac{1}{\Pi(\phi_l(z'))}dz')^{-1}

:param z: height to integrate permeability up to :type z: float :param liquid_fraction: liquid fraction on the center grid :type liquid_fraction: Numpy Array shape (I,) :param ice_depth: positive depth position of ice ocean interface :type ice_depth: float :param cell_centers: cell center positions :type cell_centers: Numpy Array of shape (I,) :param cfg: Configuration object for the simulation. :type cfg: seaice3p.params.Config :return: permeability averaged from base of the ice up to given z value

Source code in seaice3p/equations/RJW14/brine_drainage.py 
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
def calculate_integrated_mean_permeability(
    z, liquid_fraction, ice_depth, cell_centers, cfg: Config
):
    r"""Calculate the harmonic mean permeability from the base of the ice up to the
    cell containing the specified z value using the expression of ReesJones2014.

    .. math:: K(z) = (\frac{1}{h+z}\int_{-h}^{z} \frac{1}{\Pi(\phi_l(z'))}dz')^{-1}

    :param z: height to integrate permeability up to
    :type z: float
    :param liquid_fraction: liquid fraction on the center grid
    :type liquid_fraction: Numpy Array shape (I,)
    :param ice_depth: positive depth position of ice ocean interface
    :type ice_depth: float
    :param cell_centers: cell center positions
    :type cell_centers: Numpy Array of shape (I,)
    :param cfg: Configuration object for the simulation.
    :type cfg: seaice3p.params.Config
    :return: permeability averaged from base of the ice up to given z value
    """
    if z < -ice_depth:
        return 0
    step = cfg.numerical_params.step
    ice_mask = (cell_centers > -ice_depth) & (cell_centers <= z)
    permeabilities = (
        calculate_permeability(liquid_fraction[ice_mask], cfg)
        / liquid_fraction[ice_mask].size
    )
    harmonic_mean = hmean(permeabilities)
    return (ice_depth + z + step / 2) * harmonic_mean / step
calculate_permeability(liquid_fraction, cfg)

Calculate the absolute permeability as a function of liquid fraction

.. math:: \Pi(\phi_l) = \phi_l^3

Alternatively if the porosity threshold flag is true

.. math:: \Pi(\phi_l) = \phi_l^2 (\phi_l - \phi_c)

:param liquid_fraction: liquid fraction :type liquid_fraction: Numpy Array :param cfg: Configuration object for the simulation. :type cfg: seaice3p.params.Config :return: permeability on the same grid as liquid fraction

Source code in seaice3p/equations/RJW14/brine_drainage.py 
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
def calculate_permeability(liquid_fraction, cfg: Config):
    r"""Calculate the absolute permeability as a function of liquid fraction

    .. math:: \Pi(\phi_l) = \phi_l^3

    Alternatively if the porosity threshold flag is true

    .. math:: \Pi(\phi_l) = \phi_l^2 (\phi_l - \phi_c)

    :param liquid_fraction: liquid fraction
    :type liquid_fraction: Numpy Array
    :param cfg: Configuration object for the simulation.
    :type cfg: seaice3p.params.Config
    :return: permeability on the same grid as liquid fraction
    """
    if cfg.bubble_params.porosity_threshold:
        cutoff = cfg.bubble_params.porosity_threshold_value
        step_function = np.heaviside(liquid_fraction - cutoff, 0)
        return liquid_fraction**2 * (liquid_fraction - cutoff) * step_function
    return liquid_fraction**3
get_convecting_region_height(Rayleigh_number, edge_grid, cfg)

Calculate the height of the convecting region as the top edge of the highest cell in the domain for which the quantity

.. math:: \text{Ra}(z) - \text{Ra}_c

is greater than or equal to zero.

NOTE: if no convecting region exists return np.nan

:param Rayleigh_number: local rayleigh number on center grid :type Rayleigh_number: Numpy Array of shape (I,) :param edge_grid: The vertical coordinate positions of the edge grid. :type edge_grid: Numpy Array (size I+1) :param cfg: Configuration object for the simulation. :type cfg: seaice3p.params.Config :return: Edge grid value at convecting boundary.

Source code in seaice3p/equations/RJW14/brine_drainage.py 
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
def get_convecting_region_height(Rayleigh_number, edge_grid, cfg: Config):
    r"""Calculate the height of the convecting region as the top edge of the highest
    cell in the domain for which the quantity

    .. math:: \text{Ra}(z) - \text{Ra}_c

    is greater than or equal to zero.

    NOTE: if no convecting region exists return np.nan

    :param Rayleigh_number: local rayleigh number on center grid
    :type Rayleigh_number: Numpy Array of shape (I,)
    :param edge_grid: The vertical coordinate positions of the edge grid.
    :type edge_grid: Numpy Array (size I+1)
    :param cfg: Configuration object for the simulation.
    :type cfg: seaice3p.params.Config
    :return: Edge grid value at convecting boundary.
    """
    Rayleigh_critical = cfg.brine_convection_params.Rayleigh_critical
    if np.all(Rayleigh_number - Rayleigh_critical < 0):
        return np.nan
    indices = np.where(Rayleigh_number >= Rayleigh_critical)
    return edge_grid[indices[0][-1] + 1]
get_effective_Rayleigh_number(Rayleigh_number, cfg)

Calculate the effective Rayleigh Number as the maximum of

.. math:: \text{Ra}(z) - \text{Ra}_c

in the convecting region.

NOTE: if no convecting region exists returns 0.

:param Rayleigh_number: local rayleigh number on center grid :type Rayleigh_number: Numpy Array of shape (I,) :param cfg: Configuration object for the simulation. :type cfg: seaice3p.params.Config :return: Effective Rayleigh number.

Source code in seaice3p/equations/RJW14/brine_drainage.py 
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
def get_effective_Rayleigh_number(Rayleigh_number, cfg: Config):
    r"""Calculate the effective Rayleigh Number as the maximum of

    .. math:: \text{Ra}(z) - \text{Ra}_c

    in the convecting region.

    NOTE: if no convecting region exists returns 0.

    :param Rayleigh_number: local rayleigh number on center grid
    :type Rayleigh_number: Numpy Array of shape (I,)
    :param cfg: Configuration object for the simulation.
    :type cfg: seaice3p.params.Config
    :return: Effective Rayleigh number.
    """
    Rayleigh_critical = cfg.brine_convection_params.Rayleigh_critical
    return np.max(
        np.where(
            Rayleigh_number >= Rayleigh_critical, Rayleigh_number - Rayleigh_critical, 0
        )
    )

equations

_prevent_gas_rise_into_saturated_cell(Vg, state_BCs, cfg)

Modify the gas interstitial velocity to prevent bubble rise into a cell which is already theoretically saturated with gas.

From the state with boundary conditions calculate the gas and solid fraction in the cells (except at lower ghost cell). If any of these are such that there is more gas fraction than pore space available then set gas insterstitial velocity to zero on the edge below. Make sure the very top boundary velocity is not changed as we want to always alow flux to the atmosphere regardless of the boundary conditions imposed.

:param Vg: gas insterstitial velocity on cell edges :type Vg: Numpy array (size I+1) :param state_BCs: state of system with boundary conditions :type state_BCs: seaice3p.state.StateBCs :return: filtered gas interstitial velocities on edges to prevent gas rise into a fully gas saturated cell

Source code in seaice3p/equations/equations.py 
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
def _prevent_gas_rise_into_saturated_cell(
    Vg, state_BCs: StateBCs, cfg: Config
) -> NDArray:
    """Modify the gas interstitial velocity to prevent bubble rise into a cell which
    is already theoretically saturated with gas.

    From the state with boundary conditions calculate the gas and solid fraction in the
    cells (except at lower ghost cell). If any of these are such that there is more gas
    fraction than pore space available then set gas insterstitial velocity to zero on
    the edge below. Make sure the very top boundary velocity is not changed as we want
    to always alow flux to the atmosphere regardless of the boundary conditions imposed.

    :param Vg: gas insterstitial velocity on cell edges
    :type Vg: Numpy array (size I+1)
    :param state_BCs: state of system with boundary conditions
    :type state_BCs: seaice3p.state.StateBCs
    :return: filtered gas interstitial velocities on edges to prevent gas rise into a
        fully gas saturated cell

    """
    # Prevent gas rising into already gas saturated cell
    gas_fraction_above = state_BCs.gas_fraction[1:]
    solid_fraction_above = 1 - state_BCs.liquid_fraction[1:]
    filtered_Vg = np.where(gas_fraction_above + solid_fraction_above >= 1, 0, Vg)

    if cfg.bubble_params.escape_ice_surface:
        # Allow gas to leave top boundary
        filtered_Vg[-1] = Vg[-1]
    else:
        # impermeable top boundary
        filtered_Vg[-1] = 0

    return filtered_Vg

flux

Module for calculating the fluxes using upwind scheme

bulk_dissolved_gas_flux

calculate the flux terms for the dissolved gas equation in DISEQ model

gas_fraction_flux

Calculate gas phase fluxes for disequilibrium model

heat_flux

calculate_conductive_heat_flux(state_BCs, D_g, cfg)

Calculate conductive heat flux as

.. math:: -[(\phi_l + \lambda \phi_s) \frac{\partial \theta}{\partial z}]

:param temperature: temperature including ghost cells :type temperature: Numpy Array of size I+2 :param D_g: difference matrix for ghost grid :type D_g: Numpy Array :param cfg: Simulation configuration :type cfg: seaice3p.params.Config :return: conductive heat flux

Source code in seaice3p/equations/flux/heat_flux.py 
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
def calculate_conductive_heat_flux(state_BCs, D_g, cfg):
    r"""Calculate conductive heat flux as

    .. math:: -[(\phi_l + \lambda \phi_s) \frac{\partial \theta}{\partial z}]

    :param temperature: temperature including ghost cells
    :type temperature: Numpy Array of size I+2
    :param D_g: difference matrix for ghost grid
    :type D_g: Numpy Array
    :param cfg: Simulation configuration
    :type cfg: seaice3p.params.Config
    :return: conductive heat flux

    """
    temperature = state_BCs.temperature
    edge_liquid_fraction = geometric(state_BCs.liquid_fraction)
    edge_solid_fraction = 1 - edge_liquid_fraction
    conductivity = calculate_conductivity(cfg, edge_solid_fraction)
    return -conductivity * np.matmul(D_g, temperature)
pure_liquid_switch(liquid_fraction)

Take the liquid fraction and return a smoothed switch that is equal to 1 in a pure liquid region and goes to zero rapidly outside of this

Source code in seaice3p/equations/flux/heat_flux.py 
 8
 9
10
11
12
def pure_liquid_switch(liquid_fraction: NDArray | float) -> NDArray | float:
    """Take the liquid fraction and return a smoothed switch that is equal to 1 in a
    pure liquid region and goes to zero rapidly outside of this"""
    SCALE = 1e-2
    return np.exp((liquid_fraction - 1) / SCALE)

salt_flux

calculate_diffusive_salt_flux(liquid_salinity, liquid_fraction, D_g, cfg)

Take liquid salinity and liquid fraction on ghost grid and interpolate liquid fraction geometrically

Source code in seaice3p/equations/flux/salt_flux.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
def calculate_diffusive_salt_flux(liquid_salinity, liquid_fraction, D_g, cfg: Config):
    """Take liquid salinity and liquid fraction on ghost grid and interpolate liquid
    fraction geometrically"""
    lewis_salt = cfg.physical_params.lewis_salt
    edge_liquid_fraction = geometric(liquid_fraction)
    # In pure liquid phase enhanced eddy diffusivity of dissolved salt
    salt_diffusivity = edge_liquid_fraction * (
        (1 / lewis_salt)
        + cfg.physical_params.eddy_diffusivity_ratio
        * pure_liquid_switch(edge_liquid_fraction)
    )
    return -salt_diffusivity * np.matmul(D_g, liquid_salinity)

nucleation

_DISEQ_nucleation(state_BCs, cfg)

implement nucleation term

Source code in seaice3p/equations/nucleation.py 
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
def _DISEQ_nucleation(state_BCs: DISEQStateBCs, cfg: Config) -> NDArray:
    """implement nucleation term"""
    zeros = np.zeros_like(state_BCs.enthalpy[1:-1])
    chi = cfg.physical_params.expansion_coefficient
    Da = cfg.physical_params.damkohler_number
    centers = np.s_[1:-1]
    bulk_dissolved_gas = state_BCs.bulk_dissolved_gas[centers]
    liquid_fraction = state_BCs.liquid_fraction[centers]
    saturation = chi * liquid_fraction
    gas_fraction = state_BCs.gas_fraction[centers]

    is_saturated = bulk_dissolved_gas > saturation
    nucleation = np.full_like(bulk_dissolved_gas, np.nan)
    nucleation[is_saturated] = Da * (
        bulk_dissolved_gas[is_saturated] - saturation[is_saturated]
    )
    nucleation[~is_saturated] = -Da * gas_fraction[~is_saturated]

    return np.hstack(
        (
            zeros,
            zeros,
            -nucleation,
            nucleation,
        )
    )

_EQM_nucleation(state_BCs, cfg)

implement nucleation term

Source code in seaice3p/equations/nucleation.py 
21
22
23
24
def _EQM_nucleation(state_BCs: EQMStateBCs, cfg: Config) -> NDArray:
    """implement nucleation term"""
    zeros = np.zeros_like(state_BCs.enthalpy[1:-1])
    return np.hstack((zeros, zeros, zeros))

radiative_heating

Calculate internal shortwave radiative heating due to oil droplets

_calculate_non_dimensional_shortwave_heating(state_bcs, cfg, grids)

Calculate internal shortwave heating due to oil droplets on center grid

Assumes a configuration with the RadForcing object as the forcing config is passed.

Source code in seaice3p/equations/radiative_heating.py 
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
def _calculate_non_dimensional_shortwave_heating(
    state_bcs: StateBCs, cfg: Config, grids: Grids
) -> NDArray:
    """Calculate internal shortwave heating due to oil droplets on center grid

    Assumes a configuration with the RadForcing object as the forcing config is
    passed."""
    # If we don't have radiative forcing then just return array of zeros for heating
    if not (
        isinstance(cfg.forcing_config, RadForcing)
        or isinstance(cfg.forcing_config, ERA5Forcing)
    ):
        return np.zeros_like(grids.centers)

    if isinstance(cfg.forcing_config.oil_heating, DimensionalNoHeating):
        return np.zeros_like(grids.centers)

    incident_SW_in_W_m2 = get_SW_forcing(state_bcs.time, cfg)
    # If incident shortwave is small then optimize by not running the two-stream model
    if incident_SW_in_W_m2 <= 0.5:
        return np.zeros_like(grids.centers)

    dimensionless_incident_SW = cfg.scales.convert_from_dimensional_heat_flux(
        incident_SW_in_W_m2
    )

    spectral_irradiances = run_two_stream_model(state_bcs, cfg, grids)
    integrated_irradiance = oi.integrate_over_SW(spectral_irradiances)

    dz_dF_net = grids.D_e @ integrated_irradiance.net_irradiance
    return dimensionless_incident_SW * dz_dF_net

get_radiative_heating(cfg, grids)

Calculate internal shortwave heating source for enthalpy equation.

if the RadForcing object is given as the forcing config then calculates internal heating based on the object given in the configuration for oil_heating.

If another forcing is chosen then just returns a function to create an array of zeros as no internal heating is calculated.

Source code in seaice3p/equations/radiative_heating.py 
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
def get_radiative_heating(cfg: Config, grids: Grids) -> Callable[[StateBCs], NDArray]:
    """Calculate internal shortwave heating source for enthalpy equation.

    if the RadForcing object is given as the forcing config then calculates internal
    heating based on the object given in the configuration for oil_heating.

    If another forcing is chosen then just returns a function to create an array of
    zeros as no internal heating is calculated.
    """
    fun_map = {
        EQMPhysicalParams: _EQM_radiative_heating,
        DISEQPhysicalParams: _DISEQ_radiative_heating,
    }

    def radiative_heating(state_BCs: StateBCs) -> NDArray:
        return fun_map[type(cfg.physical_params)](state_BCs, cfg, grids)

    return radiative_heating

velocities

Module to calculate Darcy velocities.

The liquid Darcy velocity must be parameterised.

The gas Darcy velocity is calculated as gas_fraction x interstitial bubble velocity

Interstitial bubble velocity is found by a steady state Stoke's flow calculation. We have implemented two cases mono: All bubbles nucleate and remain the same size power_law: A power law bubble size distribution with fixed max and min.

bubble_parameters

calculate_bubble_size_fraction(bubble_radius_scaled, liquid_fraction, cfg)

Takes bubble radius scaled and liquid fraction on edges and calculates the bubble size fraction as

.. math:: \lambda = \Lambda / (\phi_l^q + \text{reg})

Returns the bubble size fraction on the edge grid.

Source code in seaice3p/equations/velocities/bubble_parameters.py 
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
def calculate_bubble_size_fraction(bubble_radius_scaled, liquid_fraction, cfg: Config):
    r"""Takes bubble radius scaled and liquid fraction on edges and calculates the
    bubble size fraction as

    .. math:: \lambda = \Lambda / (\phi_l^q + \text{reg})

    Returns the bubble size fraction on the edge grid.
    """
    exponent = cfg.bubble_params.pore_throat_scaling
    reg = cfg.numerical_params.regularisation
    effective_tube_radius = liquid_fraction**exponent + reg
    return bubble_radius_scaled / effective_tube_radius

mono_distribution

calculate_lag_function(bubble_size_fraction)

Calculate lag function from bubble size fraction on edge grid as

.. math:: G(\lambda) = 1 - \lambda / 2

for 0<lambda<1. Edge cases are given by G(0)=1 and G(1) = 0.5 for values outside this range.

Source code in seaice3p/equations/velocities/mono_distribution.py 
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
def calculate_lag_function(bubble_size_fraction):
    r"""Calculate lag function from bubble size fraction on edge grid as

    .. math:: G(\lambda) = 1 - \lambda / 2

    for 0<lambda<1. Edge cases are given by G(0)=1 and G(1) = 0.5 for values outside
    this range.
    """
    lag = np.full_like(bubble_size_fraction, np.nan)
    intermediate = (bubble_size_fraction < 1) & (bubble_size_fraction >= 0)
    large = bubble_size_fraction >= 1
    lag[bubble_size_fraction < 0] = 1
    lag[intermediate] = 1 - 0.5 * bubble_size_fraction[intermediate]
    lag[large] = 0.5
    return lag
calculate_mono_lag_factor(liquid_fraction, cfg)

Take liquid fraction on the ghost grid and calculate the lag factor for a mono bubble size distribution as

.. math:: I_2 = G(\lambda)

returns lag factor on the edge grid

Source code in seaice3p/equations/velocities/mono_distribution.py 
71
72
73
74
75
76
77
78
79
80
81
82
83
def calculate_mono_lag_factor(liquid_fraction, cfg: Config):
    r"""Take liquid fraction on the ghost grid and calculate the lag factor
    for a mono bubble size distribution as

    .. math:: I_2 = G(\lambda)

    returns lag factor on the edge grid
    """
    bubble_radius_scaled = cfg.bubble_params.bubble_radius_scaled
    bubble_size_fraction = calculate_bubble_size_fraction(
        bubble_radius_scaled, geometric(liquid_fraction), cfg
    )
    return calculate_lag_function(bubble_size_fraction)
calculate_mono_wall_drag_factor(liquid_fraction, cfg)

Take liquid fraction on the ghost grid and calculate the wall drag factor for a mono bubble size distribution as

.. math:: I_1 = \frac{\lambda^2}{K(\lambda)}

returns wall drag factor on the edge grid

Source code in seaice3p/equations/velocities/mono_distribution.py 
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
def calculate_mono_wall_drag_factor(liquid_fraction, cfg: Config):
    r"""Take liquid fraction on the ghost grid and calculate the wall drag factor
    for a mono bubble size distribution as

    .. math:: I_1 = \frac{\lambda^2}{K(\lambda)}

    returns wall drag factor on the edge grid
    """
    bubble_radius_scaled = cfg.bubble_params.bubble_radius_scaled
    bubble_size_fraction = calculate_bubble_size_fraction(
        bubble_radius_scaled, geometric(liquid_fraction), cfg
    )
    drag_function = calculate_wall_drag_function(bubble_size_fraction, cfg)
    drag_factor = drag_function * bubble_size_fraction**2
    return drag_factor
calculate_wall_drag_function(bubble_size_fraction, cfg)

Calculate wall drag function from bubble size fraction on edge grid as

.. math:: \frac{1}{K(\lambda)} = (1 - \lambda)^r

in the power law case or in the Haberman case from the paper

.. math:: \frac{1}{K(\lambda)} = \frac{1 -1.5\lambda + 1.5\lambda^5 - \lambda^6}{1+1.5\lambda^5}

for 0<lambda<1. Edge cases are given by K(0)=1 and K(1) = 0 for values outside this range.

Source code in seaice3p/equations/velocities/mono_distribution.py 
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
def calculate_wall_drag_function(bubble_size_fraction, cfg: Config):
    r"""Calculate wall drag function from bubble size fraction on edge grid as

    .. math:: \frac{1}{K(\lambda)} = (1 - \lambda)^r

    in the power law case or in the Haberman case from the paper

    .. math:: \frac{1}{K(\lambda)} = \frac{1 -1.5\lambda + 1.5\lambda^5 - \lambda^6}{1+1.5\lambda^5}

    for 0<lambda<1. Edge cases are given by K(0)=1 and K(1) = 0 for values outside
    this range.
    """
    drag = np.full_like(bubble_size_fraction, np.nan)
    intermediate = (bubble_size_fraction < 1) & (bubble_size_fraction >= 0)
    large = bubble_size_fraction >= 1
    drag[bubble_size_fraction < 0] = 1
    drag[intermediate] = (
        1
        - 1.5 * bubble_size_fraction[intermediate]
        + 1.5 * bubble_size_fraction[intermediate] ** 5
        - bubble_size_fraction[intermediate] ** 6
    ) / (1 + 1.5 * bubble_size_fraction[intermediate] ** 5)
    drag[large] = 0
    return drag

power_law_distribution

calculate_lag_integrand(bubble_size_fraction, cfg)

Scalar function to calculate lag integrand for polydispersive case.

Bubble size fraction is given as a scalar input to calculate

.. math:: \lambda^{3-p} G(\lambda)

Source code in seaice3p/equations/velocities/power_law_distribution.py 
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
def calculate_lag_integrand(bubble_size_fraction: float, cfg: Config):
    r"""Scalar function to calculate lag integrand for polydispersive case.

    Bubble size fraction is given as a scalar input to calculate

    .. math:: \lambda^{3-p} G(\lambda)

    """
    power_law = cfg.bubble_params.bubble_distribution_power
    if bubble_size_fraction < 0:
        return 0
    elif (bubble_size_fraction >= 0) and (bubble_size_fraction < 1):
        return (1 - 0.5 * bubble_size_fraction) * (
            bubble_size_fraction ** (3 - power_law)
        )
    else:
        return 0.5
calculate_power_law_lag_factor(liquid_fraction, cfg)

Take liquid fraction on the ghost grid and calculate the lag factor for power law bubble size distribution.

Return on edge grid

Source code in seaice3p/equations/velocities/power_law_distribution.py 
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
def calculate_power_law_lag_factor(liquid_fraction, cfg: Config):
    r"""Take liquid fraction on the ghost grid and calculate the lag factor
    for power law bubble size distribution.

    Return on edge grid
    """
    minimum_size_fractions = calculate_bubble_size_fraction(
        cfg.bubble_params.minimum_bubble_radius_scaled,
        geometric(liquid_fraction),
        cfg,
    )
    maximum_size_fractions = calculate_bubble_size_fraction(
        cfg.bubble_params.maximum_bubble_radius_scaled,
        geometric(liquid_fraction),
        cfg,
    )
    lag_factor = np.full_like(minimum_size_fractions, np.nan)
    for i, (min, max) in enumerate(zip(minimum_size_fractions, maximum_size_fractions)):
        lag_factor[i] = calculate_lag_integral(min, max, cfg)
    return lag_factor
calculate_power_law_wall_drag_factor(liquid_fraction, cfg)

Take liquid fraction on the ghost grid and calculate the wall drag factor for power law bubble size distribution.

Return on edge grid

Source code in seaice3p/equations/velocities/power_law_distribution.py 
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
def calculate_power_law_wall_drag_factor(liquid_fraction, cfg: Config):
    r"""Take liquid fraction on the ghost grid and calculate the wall drag factor
    for power law bubble size distribution.

    Return on edge grid
    """
    minimum_size_fractions = calculate_bubble_size_fraction(
        cfg.bubble_params.minimum_bubble_radius_scaled,
        geometric(liquid_fraction),
        cfg,
    )
    maximum_size_fractions = calculate_bubble_size_fraction(
        cfg.bubble_params.maximum_bubble_radius_scaled,
        geometric(liquid_fraction),
        cfg,
    )
    drag_factor = np.full_like(minimum_size_fractions, np.nan)
    for i, (min, max) in enumerate(zip(minimum_size_fractions, maximum_size_fractions)):
        drag_factor[i] = calculate_wall_drag_integral(min, max, cfg)
    return drag_factor
calculate_volume_integrand(bubble_size_fraction, cfg)

Scalar function to calculate the integrand for volume under a power law bubble size distribution given as

.. math:: \lambda^{3-p}

in terms of the bubble size fraction.

Source code in seaice3p/equations/velocities/power_law_distribution.py 
58
59
60
61
62
63
64
65
66
67
def calculate_volume_integrand(bubble_size_fraction: float, cfg: Config):
    r"""Scalar function to calculate the integrand for volume under a power law
    bubble size distribution given as

    .. math:: \lambda^{3-p}

    in terms of the bubble size fraction.
    """
    p = cfg.bubble_params.bubble_distribution_power
    return bubble_size_fraction ** (3 - p)
calculate_wall_drag_integrand(bubble_size_fraction, cfg)

Scalar function to calculate wall drag integrand for polydispersive case.

Bubble size fraction is given as a scalar input to calculate

.. math:: \frac{\lambda^{5-p}}{K(\lambda)}

where the wall drag enhancement funciton K can be given by a power law fit or taken from the Haberman paper.

Source code in seaice3p/equations/velocities/power_law_distribution.py 
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
def calculate_wall_drag_integrand(bubble_size_fraction: float, cfg: Config):
    r"""Scalar function to calculate wall drag integrand for polydispersive case.

    Bubble size fraction is given as a scalar input to calculate

    .. math:: \frac{\lambda^{5-p}}{K(\lambda)}

    where the wall drag enhancement funciton K can be given by a power law fit
    or taken from the Haberman paper.
    """
    power_law = cfg.bubble_params.bubble_distribution_power
    if bubble_size_fraction < 0:
        return 0
    elif (bubble_size_fraction >= 0) and (bubble_size_fraction < 1):
        return (
            (
                1
                - 1.5 * bubble_size_fraction
                + 1.5 * bubble_size_fraction**5
                - bubble_size_fraction**6
            )
            / (1 + 1.5 * bubble_size_fraction**5)
        ) * (bubble_size_fraction ** (5 - power_law))
    else:
        return 0

velocities

calculate_gas_interstitial_velocity(liquid_fraction, liquid_darcy_velocity, wall_drag_factor, lag_factor, cfg)

Calculate Vg from liquid fraction on the ghost frid and liquid interstitial velocity

.. math:: V_g = \mathcal{B} (\phi_l^{2q} I_1) + U_0 I_2

Return Vg on edge grid

Source code in seaice3p/equations/velocities/velocities.py 
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
def calculate_gas_interstitial_velocity(
    liquid_fraction,
    liquid_darcy_velocity,
    wall_drag_factor,
    lag_factor,
    cfg: Config,
):
    r"""Calculate Vg from liquid fraction on the ghost frid and liquid interstitial velocity

    .. math:: V_g = \mathcal{B} (\phi_l^{2q} I_1) + U_0 I_2

    Return Vg on edge grid
    """
    B = cfg.bubble_params.B
    exponent = cfg.bubble_params.pore_throat_scaling

    REGULARISATION = 1e-10
    liquid_interstitial_velocity = (
        liquid_darcy_velocity * 2 / (geometric(liquid_fraction) + REGULARISATION)
    )

    viscosity_factor = (
        2
        * (1 + cfg.physical_params.gas_viscosity_ratio)
        / (2 + 3 * cfg.physical_params.gas_viscosity_ratio)
    )
    Vg = (
        viscosity_factor
        * B
        * wall_drag_factor
        * geometric(liquid_fraction) ** (2 * exponent)
        + liquid_interstitial_velocity * lag_factor
    )

    # apply a porosity cutoff to the gas interstitial velocity if necking occurs below
    # critical porosity.
    if cfg.bubble_params.porosity_threshold:
        return Vg * np.heaviside(
            geometric(liquid_fraction) - cfg.bubble_params.porosity_threshold_value,
            0,
        )

    return Vg
calculate_liquid_darcy_velocity(liquid_fraction, liquid_salinity, center_grid, edge_grid, cfg)

Calculate liquid Darcy velocity either using brine convection parameterisation or as stagnant

:param liquid_fraction: liquid fraction on ghost grid :type liquid_fraction: Numpy Array (size I+2) :param liquid_salinity: liquid salinity on ghost grid :type liquid_salinity: Numpy Array (size I+2) :param center_grid: vertical coordinates of cell centers :type center_grid: Numpy Array of shape (I,) :param edge_grid: Vertical coordinates of cell edges :type edge_grid: Numpy Array (size I+1) :param cfg: simulation configuration object :type cfg: seaice3p.params.Config :return: liquid darcy velocity on edge grid

Source code in seaice3p/equations/velocities/velocities.py 
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
def calculate_liquid_darcy_velocity(
    liquid_fraction, liquid_salinity, center_grid, edge_grid, cfg: Config
):
    r"""Calculate liquid Darcy velocity either using brine convection parameterisation
    or as stagnant


    :param liquid_fraction: liquid fraction on ghost grid
    :type liquid_fraction: Numpy Array (size I+2)
    :param liquid_salinity: liquid salinity on ghost grid
    :type liquid_salinity: Numpy Array (size I+2)
    :param center_grid: vertical coordinates of cell centers
    :type center_grid: Numpy Array of shape (I,)
    :param edge_grid: Vertical coordinates of cell edges
    :type edge_grid: Numpy Array (size I+1)
    :param cfg: simulation configuration object
    :type cfg: seaice3p.params.Config
    :return: liquid darcy velocity on edge grid
    """
    if isinstance(cfg.brine_convection_params, NoBrineConvection):
        return np.zeros_like(geometric(liquid_fraction))

    Wl = calculate_brine_convection_liquid_velocity(
        liquid_fraction[1:-1], liquid_salinity[1:-1], center_grid, edge_grid, cfg
    )
    return Wl
calculate_velocities(state_BCs, cfg)

Inputs on ghost grid, outputs on edge grid

needs the simulation config, liquid fraction, liquid salinity and grids

Source code in seaice3p/equations/velocities/velocities.py 
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
def calculate_velocities(state_BCs, cfg: Config):
    """Inputs on ghost grid, outputs on edge grid

    needs the simulation config, liquid fraction, liquid salinity and grids
    """
    liquid_fraction = state_BCs.liquid_fraction
    liquid_salinity = state_BCs.liquid_salinity
    center_grid, edge_grid = (
        Grids(cfg.numerical_params.I).centers,
        Grids(cfg.numerical_params.I).edges,
    )

    match cfg.bubble_params:
        case MonoBubbleParams():
            wall_drag_factor = calculate_mono_wall_drag_factor(liquid_fraction, cfg)
            lag_factor = calculate_mono_lag_factor(liquid_fraction, cfg)
        case PowerLawBubbleParams():
            wall_drag_factor = calculate_power_law_wall_drag_factor(
                liquid_fraction, cfg
            )
            lag_factor = calculate_power_law_lag_factor(liquid_fraction, cfg)
        case _:
            raise NotImplementedError

    # check if we want to couple the bubble to fluid motion in the vertical
    if not isinstance(cfg.brine_convection_params, NoBrineConvection):
        if not cfg.brine_convection_params.couple_bubble_to_vertical_flow:
            lag_factor = np.zeros_like(wall_drag_factor)

    Wl = calculate_liquid_darcy_velocity(
        liquid_fraction, liquid_salinity, center_grid, edge_grid, cfg
    )
    Vg = calculate_gas_interstitial_velocity(
        liquid_fraction, Wl, wall_drag_factor, lag_factor, cfg
    )
    V = calculate_frame_velocity(cfg)
    return Vg, Wl, V

example

Script to run a simulation starting with dimensional parameters and plot output

main(data_directory, frames_directory, simulation_dimensional_params)

Generate non dimensional simulation config and save along with dimensional config then run simulation and save data.

Source code in seaice3p/example.py 
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
def main(
    data_directory: Path,
    frames_directory: Path,
    simulation_dimensional_params: DimensionalParams,
):
    """Generate non dimensional simulation config and save along with dimensional
    config then run simulation and save data.
    """

    print(f"seaice3p version {__version__}")

    cfg = create_and_save_config(data_directory, simulation_dimensional_params)
    solve(cfg, data_directory, verbosity_level=1)

    # Analysis load simulation data
    # plot:
    # gas_fraction
    # salt
    # temperature
    # solid_fraction
    # save as frames in frames/gas_fraction etc...
    simulation_name = simulation_dimensional_params.name
    DIMENSIONAL_CONFIG_DATA_PATH = data_directory / f"{simulation_name}_dimensional.yml"
    results = load_simulation(
        DATA_DIRECTORY / "example_dimensional.yml",
        DATA_DIRECTORY / "example.npz",
        is_dimensional=True,
    )

    scales = results.cfg.scales
    dimensional_grid = scales.convert_to_dimensional_grid(results.grids.centers)
    dimensional_times = scales.convert_to_dimensional_time(results.times)

    GAS_FRACTION_DIR = frames_directory / "gas_fraction/"
    GAS_FRACTION_DIR.mkdir(exist_ok=True, parents=True)

    TEMPERATURE_DIR = frames_directory / "temperature/"
    TEMPERATURE_DIR.mkdir(exist_ok=True, parents=True)

    SOLID_FRAC_DIR = frames_directory / "solid_fraction/"
    SOLID_FRAC_DIR.mkdir(exist_ok=True, parents=True)

    BULK_AIR_DIR = frames_directory / "bulk_air/"
    BULK_AIR_DIR.mkdir(exist_ok=True, parents=True)

    BULK_SALT_DIR = frames_directory / "bulk_salt/"
    BULK_SALT_DIR.mkdir(exist_ok=True, parents=True)

    for n, state in enumerate(results.states):
        plt.figure(figsize=(5, 5))
        plt.plot(
            state.gas_fraction,
            dimensional_grid,
            "g*--",
        )
        plt.title(f"{dimensional_times[n]:.0f} days")
        plt.xlabel("gas fraction")
        plt.ylabel("depth [m]")
        plt.savefig(GAS_FRACTION_DIR / f"gas_fraction{n}.pdf")
        plt.close()

        plt.figure(figsize=(5, 5))
        plt.plot(
            scales.convert_to_dimensional_bulk_salinity(state.salt),
            dimensional_grid,
            "b*--",
        )
        plt.title(f"{dimensional_times[n]:.0f} days")
        plt.xlabel("bulk salinity [g/kg]")
        plt.ylabel("depth [m]")
        plt.savefig(BULK_SALT_DIR / f"bulk_salt{n}.pdf")
        plt.close()

        plt.figure(figsize=(5, 5))
        dimensional_temperature = scales.convert_to_dimensional_temperature(
            state.temperature
        )
        plt.plot(
            dimensional_temperature,
            dimensional_grid,
            "r*--",
        )
        plt.title(f"{dimensional_times[n]:.0f} days")
        plt.xlabel("temperature [deg C]")
        plt.ylabel("depth [m]")
        plt.savefig(TEMPERATURE_DIR / f"temperature{n}.pdf")
        plt.close()

        plt.figure(figsize=(5, 5))
        plt.plot(
            state.solid_fraction,
            dimensional_grid,
            "m*--",
        )
        plt.title(f"{dimensional_times[n]:.0f} days")
        plt.xlabel("solid fraction")
        plt.ylabel("depth [m]")
        plt.savefig(SOLID_FRAC_DIR / f"solid_fraction{n}.pdf")
        plt.close()

        plt.figure(figsize=(5, 5))
        dimensional_bulk_air = scales.convert_to_dimensional_bulk_gas(state.gas)
        argon_micromole_per_liter = (
            scales.convert_dimensional_bulk_air_to_argon_content(dimensional_bulk_air)
        )
        plt.plot(
            argon_micromole_per_liter,
            dimensional_grid,
            "m*--",
        )
        plt.title(f"{dimensional_times[n]:.0f} days")
        plt.xlabel("bulk argon [micromole/L]")
        plt.ylabel("depth [m]")
        plt.savefig(BULK_AIR_DIR / f"bulk_air{n}.pdf")
        plt.close()

    for converter, unit, attr in zip(
        [
            scales.convert_to_dimensional_temperature,
            scales.convert_to_dimensional_bulk_salinity,
            lambda x: x,
            lambda x: x,
            scales.convert_to_dimensional_bulk_gas,
        ],
        ["[deg C]", "[g/kg]", "", "", "[kg/m3]"],
        ["temperature", "salt", "gas_fraction", "solid_fraction", "bulk_gas"],
    ):
        plt.figure()
        plt.contourf(
            dimensional_times, dimensional_grid, converter(getattr(results, attr))
        )
        plt.colorbar()
        plt.title(f"{attr} {unit}")
        plt.xlabel("time [days]")
        plt.ylabel("depth [m]")
        plt.savefig(data_directory / f"contours_{attr}.pdf")

forcing

boundary_conditions

Module to provide functions to add boundary conditions to each quantity on the centered grid that needs to be on the ghost grid for the upwind scheme.

_dissolved_gas_BCs(dissolved_gas_centers, cfg)

Add ghost cells with BCs to center quantity

Source code in seaice3p/forcing/boundary_conditions.py 
90
91
92
93
94
def _dissolved_gas_BCs(dissolved_gas_centers, cfg: Config):
    """Add ghost cells with BCs to center quantity"""
    return add_ghost_cells(
        dissolved_gas_centers, bottom=cfg.ocean_forcing_config.ocean_gas_sat, top=1
    )

_enthalpy_BCs(enthalpy_centers, cfg, bottom_temperature)

Add ghost cells with BCs to center quantity

Source code in seaice3p/forcing/boundary_conditions.py 
134
135
136
137
138
def _enthalpy_BCs(enthalpy_centers, cfg: Config, bottom_temperature):
    """Add ghost cells with BCs to center quantity"""
    return add_ghost_cells(
        enthalpy_centers, bottom=bottom_temperature, top=enthalpy_centers[-1]
    )

_gas_BCs(gas_centers, cfg)

Add ghost cells with BCs to center quantity

Source code in seaice3p/forcing/boundary_conditions.py 
111
112
113
114
115
def _gas_BCs(gas_centers, cfg: Config):
    """Add ghost cells with BCs to center quantity"""
    chi = cfg.physical_params.expansion_coefficient
    far_gas_sat = cfg.ocean_forcing_config.ocean_gas_sat
    return add_ghost_cells(gas_centers, bottom=chi * far_gas_sat, top=chi)

_gas_fraction_BCs(gas_fraction_centers, cfg)

Add ghost cells with BCs to center quantity

Source code in seaice3p/forcing/boundary_conditions.py 
 97
 98
 99
100
101
102
103
104
105
106
107
108
def _gas_fraction_BCs(gas_fraction_centers, cfg: Config):
    """Add ghost cells with BCs to center quantity"""
    if isinstance(cfg.initial_conditions_config, OilInitialConditions):
        return add_ghost_cells(
            gas_fraction_centers,
            bottom=cfg.initial_conditions_config.initial_oil_volume_fraction,
            top=0,
        )
    else:
        return add_ghost_cells(
            gas_fraction_centers, bottom=gas_fraction_centers[0], top=0
        )

_liquid_fraction_BCs(liquid_fraction_centers)

Add ghost cells to liquid fraction such that top and bottom boundaries take the same value as the top and bottom cell center

Source code in seaice3p/forcing/boundary_conditions.py 
146
147
148
149
150
151
152
153
def _liquid_fraction_BCs(liquid_fraction_centers):
    """Add ghost cells to liquid fraction such that top and bottom boundaries take the
    same value as the top and bottom cell center"""
    return add_ghost_cells(
        liquid_fraction_centers,
        bottom=liquid_fraction_centers[0],
        top=liquid_fraction_centers[-1],
    )

_liquid_salinity_BCs(liquid_salinity_centers, cfg)

Add ghost cells with BCs to center quantity

Source code in seaice3p/forcing/boundary_conditions.py 
118
119
120
121
122
def _liquid_salinity_BCs(liquid_salinity_centers, cfg: Config):
    """Add ghost cells with BCs to center quantity"""
    return add_ghost_cells(
        liquid_salinity_centers, bottom=0, top=liquid_salinity_centers[-1]
    )

_salt_BCs(salt_centers, cfg)

Add ghost cells with BCs to center quantity

Source code in seaice3p/forcing/boundary_conditions.py 
141
142
143
def _salt_BCs(salt_centers, cfg: Config):
    """Add ghost cells with BCs to center quantity"""
    return add_ghost_cells(salt_centers, bottom=0, top=salt_centers[-1])

_temperature_BCs(state, cfg)

Add ghost cells with BCs to center quantity

Note this needs the current time as well as top temperature is forced.

Source code in seaice3p/forcing/boundary_conditions.py 
125
126
127
128
129
130
131
def _temperature_BCs(state: StateFull, cfg: Config):
    """Add ghost cells with BCs to center quantity

    Note this needs the current time as well as top temperature is forced."""
    far_temp = get_bottom_temperature_forcing(state, cfg)
    top_temp = get_temperature_forcing(state, cfg)
    return add_ghost_cells(state.temperature, bottom=far_temp, top=top_temp)

radiative_forcing

Module for providing surface radiative forcing to simulation.

Currently only total surface shortwave irradiance (integrated over entire shortwave part of the spectrum) is provided and this is used to calculate internal radiative heating.

Unlike temperature forcing this provides dimensional forcing

_constant_LW_forcing(time, cfg)

Returns constant surface longwave downwelling irradiance in W/m2 integrated over the entire longwave spectrum

Source code in seaice3p/forcing/radiative_forcing.py 
37
38
39
40
41
def _constant_LW_forcing(time, cfg: Config):
    """Returns constant surface longwave downwelling irradiance in W/m2 integrated
    over the entire longwave spectrum
    """
    return cfg.forcing_config.LW_forcing.LW_irradiance

_constant_SW_forcing(time, cfg)

Returns constant surface shortwave downwelling irradiance in W/m2 integrated over the entire shortwave spectrum

Source code in seaice3p/forcing/radiative_forcing.py 
21
22
23
24
25
def _constant_SW_forcing(time, cfg: Config):
    """Returns constant surface shortwave downwelling irradiance in W/m2 integrated
    over the entire shortwave spectrum
    """
    return cfg.forcing_config.SW_forcing.SW_irradiance

surface_energy_balance

surface_energy_balance

Module to compute the surface heat flux from geophysical energy balance

following [1]

Refs: [1] P. D. Taylor and D. L. Feltham, ‘A model of melt pond evolution on sea ice’, J. Geophys. Res., vol. 109, no. C12, p. 2004JC002361, Dec. 2004, doi: 10.1029/2004JC002361.

_calculate_total_heat_flux(cfg, time, top_cell_is_ice, surface_temp, temp_gradient, top_cell_conductivity)

Takes non-dimensional surface temperature and returns non-dimensional heat flux

Source code in seaice3p/forcing/surface_energy_balance/surface_energy_balance.py 
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
def _calculate_total_heat_flux(
    cfg: Config,
    time: float,
    top_cell_is_ice: bool,
    surface_temp: float,
    temp_gradient: float,
    top_cell_conductivity: float,
) -> float:
    """Takes non-dimensional surface temperature and returns non-dimensional heat flux"""
    if isinstance(cfg.forcing_config, ERA5Forcing):
        dimensional_temperature_gradient = (
            cfg.scales.temperature_difference * temp_gradient / cfg.scales.lengthscale
        )
        surface_temp_K = (
            _convert_non_dim_temperature_to_kelvin(cfg, surface_temp)
            + (top_cell_conductivity / cfg.physical_params.snow_conductivity_ratio)
            * cfg.forcing_config.get_snow_depth(time)
            * dimensional_temperature_gradient
        )
    else:
        surface_temp_K = _convert_non_dim_temperature_to_kelvin(cfg, surface_temp)
    emissivity = _calculate_emissivity(cfg, top_cell_is_ice)
    dimensional_heat_flux = (
        get_LW_forcing(time, cfg)
        - emissivity * STEFAN_BOLTZMANN * surface_temp_K**4
        + calculate_sensible_heat_flux(cfg, time, top_cell_is_ice, surface_temp_K)
        + calculate_latent_heat_flux(cfg, time, top_cell_is_ice, surface_temp_K)
    )
    return cfg.scales.convert_from_dimensional_heat_flux(dimensional_heat_flux)
find_ghost_cell_temperature(state, cfg)

Returns non dimensional ghost cell temperature such that surface heat flux is the sum of incoming LW, outgoing LW, sensible and latent heat fluxes. The SW heat flux is determined in the radiative heating term.

Source code in seaice3p/forcing/surface_energy_balance/surface_energy_balance.py 
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
def find_ghost_cell_temperature(state: StateFull, cfg: Config) -> float:
    """Returns non dimensional ghost cell temperature such that surface heat flux
    is the sum of incoming LW, outgoing LW, sensible and latent heat fluxes.
    The SW heat flux is determined in the radiative heating term."""
    if state.solid_fraction[-1] == 0:
        top_cell_is_ice = False
    else:
        top_cell_is_ice = True

    def residual(ghost_cell_temperature: float) -> float:
        surface_temperature = 0.5 * (ghost_cell_temperature + state.temperature[-1])
        temp_gradient = (1 / cfg.numerical_params.step) * (
            ghost_cell_temperature - state.temperature[-1]
        )
        conductivity = calculate_conductivity(cfg, state.solid_fraction[-1])
        return conductivity * temp_gradient - _calculate_total_heat_flux(
            cfg,
            state.time,
            top_cell_is_ice,
            surface_temperature,
            temp_gradient,
            conductivity,
        )

    initial_guess = state.temperature[-1]
    solution = fsolve(residual, initial_guess)[0]
    return solution

turbulent_heat_flux

Module to compute the turbulent atmospheric sensible and latent heat fluxes

All temperatures are in Kelvin in this module

Refs: [1] P. D. Taylor and D. L. Feltham, ‘A model of melt pond evolution on sea ice’, J. Geophys. Res., vol. 109, no. C12, p. 2004JC002361, Dec. 2004, doi: 10.1029/2004JC002361.

[2] E. E. Ebert and J. A. Curry, ‘An intermediate one-dimensional thermodynamic sea ice model for investigating ice-atmosphere interactions’, Journal of Geophysical Research: Oceans, vol. 98, no. C6, pp. 10085–10109, 1993, doi: 10.1029/93JC00656.

_calculate_bulk_transfer_coefficient(cfg, top_cell_is_ice, time, surface_temp)

Calculation of bulk transfer coeff from [2]

Source code in seaice3p/forcing/surface_energy_balance/turbulent_heat_flux.py 
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
def _calculate_bulk_transfer_coefficient(
    cfg: Config, top_cell_is_ice: bool, time: float, surface_temp: float
) -> float:
    """Calculation of bulk transfer coeff from [2]"""
    if top_cell_is_ice:
        CT0 = 1.3e-3
    else:
        CT0 = 1.0e-3
    BPRIME = 20
    C = 1961 * BPRIME * CT0
    ref_air_temp = _calculate_ref_air_temp(cfg, time)
    ref_windspeed = _calculate_ref_windspeed(cfg, time)
    ref_height = cfg.forcing_config.turbulent_flux.ref_height
    Richardson = (
        GRAVITY
        * (ref_air_temp - surface_temp)
        * ref_height
        / (ref_air_temp * ref_windspeed**2)
    )
    if Richardson < 0:
        frac = 2 * BPRIME * Richardson / (1 + C * np.sqrt(np.abs(Richardson)))
        return CT0 * (1 - frac)

    return CT0 * (1 + BPRIME * Richardson) ** (-2)
_calculate_ref_air_temp(cfg, time)

return air temperature at reference level above the ice in Kelvin

in the configuration the air temperature is given in deg C

Source code in seaice3p/forcing/surface_energy_balance/turbulent_heat_flux.py 
36
37
38
39
40
41
42
43
44
45
46
def _calculate_ref_air_temp(cfg: Config, time: float) -> float:
    """return air temperature at reference level above the ice in Kelvin

    in the configuration the air temperature is given in deg C
    """
    if isinstance(cfg.forcing_config, RadForcing):
        return cfg.forcing_config.turbulent_flux.air_temp + 273.15
    elif isinstance(cfg.forcing_config, ERA5Forcing):
        return cfg.forcing_config.get_2m_temp(time) + 273.15
    else:
        raise NotImplementedError
_calculate_ref_atmospheric_pressure(cfg, time)

return atmospheric pressure at reference level above the ice in KPa

Source code in seaice3p/forcing/surface_energy_balance/turbulent_heat_flux.py 
59
60
61
62
63
64
65
66
def _calculate_ref_atmospheric_pressure(cfg: Config, time: float) -> float:
    """return atmospheric pressure at reference level above the ice in KPa"""
    if isinstance(cfg.forcing_config, RadForcing):
        return cfg.forcing_config.turbulent_flux.atm_pressure
    elif isinstance(cfg.forcing_config, ERA5Forcing):
        return cfg.forcing_config.get_ATM(time)
    else:
        raise NotImplementedError
_calculate_ref_specific_humidity(cfg, time)

return specific humidity at reference level above the ice

Source code in seaice3p/forcing/surface_energy_balance/turbulent_heat_flux.py 
49
50
51
52
53
54
55
56
def _calculate_ref_specific_humidity(cfg: Config, time: float) -> float:
    """return specific humidity at reference level above the ice"""
    if isinstance(cfg.forcing_config, RadForcing):
        return cfg.forcing_config.turbulent_flux.specific_humidity
    elif isinstance(cfg.forcing_config, ERA5Forcing):
        return cfg.forcing_config.get_spec_hum(time)
    else:
        raise NotImplementedError
_calculate_ref_windspeed(cfg, time)

return windspeed at reference level above the ice

Source code in seaice3p/forcing/surface_energy_balance/turbulent_heat_flux.py 
26
27
28
29
30
31
32
33
def _calculate_ref_windspeed(cfg: Config, time: float) -> float:
    """return windspeed at reference level above the ice"""
    if isinstance(cfg.forcing_config, RadForcing):
        return cfg.forcing_config.turbulent_flux.windspeed
    elif isinstance(cfg.forcing_config, ERA5Forcing):
        return cfg.forcing_config.get_windspeed(time)
    else:
        raise NotImplementedError("No windspeed for this forcing configuration")
_calculate_surface_specific_humidity(cfg, time, surface_temp)

Following expression given in [1]

Source code in seaice3p/forcing/surface_energy_balance/turbulent_heat_flux.py 
 95
 96
 97
 98
 99
100
101
102
103
104
105
def _calculate_surface_specific_humidity(
    cfg: Config, time: float, surface_temp: float
) -> float:
    """Following expression given in [1]"""
    water_vapor_partial_pressure = 2.53e8 * np.exp(-(5420 / surface_temp))
    atm_pressure = _calculate_ref_atmospheric_pressure(cfg, time)
    return (
        0.622
        * water_vapor_partial_pressure
        / (atm_pressure - 0.378 * water_vapor_partial_pressure)
    )
calculate_latent_heat_flux(cfg, time, top_cell_is_ice, surface_temp)

Calculate latent heat flux from [2]

Source code in seaice3p/forcing/surface_energy_balance/turbulent_heat_flux.py 
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
def calculate_latent_heat_flux(
    cfg: Config, time: float, top_cell_is_ice: bool, surface_temp: float
) -> float:
    """Calculate latent heat flux from [2]"""
    air_density = cfg.forcing_config.turbulent_flux.air_density
    air_latent_heat_of_vaporisation = (
        cfg.forcing_config.turbulent_flux.air_latent_heat_of_vaporisation
    )
    windspeed = _calculate_ref_windspeed(cfg, time)
    ref_specific_humidity = _calculate_ref_specific_humidity(cfg, time)
    bulk_transfer_coeff = _calculate_bulk_transfer_coefficient(
        cfg, top_cell_is_ice, time, surface_temp
    )
    surface_specific_humidity = _calculate_surface_specific_humidity(
        cfg, time, surface_temp
    )
    return (
        air_density
        * air_latent_heat_of_vaporisation
        * bulk_transfer_coeff
        * windspeed
        * (ref_specific_humidity - surface_specific_humidity)
    )
calculate_sensible_heat_flux(cfg, time, top_cell_is_ice, surface_temp)

Calculate sensible heat flux from [2]

Source code in seaice3p/forcing/surface_energy_balance/turbulent_heat_flux.py 
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
def calculate_sensible_heat_flux(
    cfg: Config, time: float, top_cell_is_ice: bool, surface_temp: float
) -> float:
    """Calculate sensible heat flux from [2]"""
    air_density = cfg.forcing_config.turbulent_flux.air_density
    air_heat_capacity = cfg.forcing_config.turbulent_flux.air_heat_capacity
    ref_air_temp = _calculate_ref_air_temp(cfg, time)
    windspeed = _calculate_ref_windspeed(cfg, time)
    bulk_transfer_coeff = _calculate_bulk_transfer_coefficient(
        cfg, top_cell_is_ice, time, surface_temp
    )
    return (
        air_density
        * air_heat_capacity
        * bulk_transfer_coeff
        * windspeed
        * (ref_air_temp - surface_temp)
    )

temperature_forcing

Module for providing surface temperature forcing to simulation.

Note that the barrow temperature data is read in from a file if needed by the simulation configuration.

_Robin_forcing(state, cfg)

Returns non dimensional ghost cell temperature such that surface heat flux is given by Robin boundary condition

Source code in seaice3p/forcing/temperature_forcing.py 
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
def _Robin_forcing(state: StateFull, cfg: Config):
    """Returns non dimensional ghost cell temperature such that surface heat flux
    is given by Robin boundary condition"""

    def residual(ghost_cell_temperature: float) -> float:
        surface_temperature = 0.5 * (ghost_cell_temperature + state.temperature[-1])
        temp_gradient = (1 / cfg.numerical_params.step) * (
            ghost_cell_temperature - state.temperature[-1]
        )
        return calculate_conductivity(
            cfg, state.solid_fraction[-1]
        ) * temp_gradient - cfg.forcing_config.biot * (
            cfg.forcing_config.restoring_temperature - surface_temperature
        )

    initial_guess = state.temperature[-1]
    solution = fsolve(residual, initial_guess)[0]
    return solution

_barrow_ocean_temperature_forcing(state, cfg)

Take non dimensional time and return non dimensional ocean temperature at the Barrow site in 2009.

For this to work you must have created the configuration cfg from dimensional parameters as it must have the conversion scales object.

Source code in seaice3p/forcing/temperature_forcing.py 
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
def _barrow_ocean_temperature_forcing(state: StateFull, cfg: Config) -> float:
    """Take non dimensional time and return non dimensional ocean temperature at
    the Barrow site in 2009.

    For this to work you must have created the configuration cfg from dimensional
    parameters as it must have the conversion scales object.
    """
    time_in_days = cfg.scales.convert_to_dimensional_time(state.time)
    dimensional_temperature = _dimensional_barrow_ocean_temperature_forcing(
        time_in_days, cfg
    )
    temperature = cfg.scales.convert_from_dimensional_temperature(
        dimensional_temperature
    )
    return temperature

_barrow_temperature_forcing(state, cfg)

Take non dimensional time and return non dimensional air/snow/ice temperature at the Barrow site in 2009.

For this to work you must have created the configuration cfg from dimensional parameters as it must have the conversion scales object.

Source code in seaice3p/forcing/temperature_forcing.py 
71
72
73
74
75
76
77
78
79
80
81
82
83
def _barrow_temperature_forcing(state: StateFull, cfg: Config):
    """Take non dimensional time and return non dimensional air/snow/ice temperature at
    the Barrow site in 2009.

    For this to work you must have created the configuration cfg from dimensional
    parameters as it must have the conversion scales object.
    """
    time_in_days = cfg.scales.convert_to_dimensional_time(state.time)
    dimensional_temperature = _dimensional_barrow_temperature_forcing(time_in_days, cfg)
    temperature = cfg.scales.convert_from_dimensional_temperature(
        dimensional_temperature
    )
    return temperature

_dimensional_barrow_ocean_temperature_forcing(time_in_days, cfg)

Take time in days and linearly interp 2009 Barrow ocean temperature data to get temperature in degrees Celsius.

Source code in seaice3p/forcing/temperature_forcing.py 
110
111
112
113
114
115
116
117
118
def _dimensional_barrow_ocean_temperature_forcing(
    time_in_days: float, cfg: Config
) -> float:
    """Take time in days and linearly interp 2009 Barrow ocean temperature data to get
    temperature in degrees Celsius.
    """
    barrow_ocean_days = cfg.ocean_forcing_config.barrow_ocean_days
    barrow_bottom_temp = cfg.ocean_forcing_config.barrow_bottom_temp
    return np.interp(time_in_days, barrow_ocean_days, barrow_bottom_temp, right=np.nan)

_dimensional_barrow_temperature_forcing(time_in_days, cfg)

Take time in days and linearly interp 2009 Barrow air/snow/ice temperature data to get temperature in degrees Celsius.

Source code in seaice3p/forcing/temperature_forcing.py 
62
63
64
65
66
67
68
def _dimensional_barrow_temperature_forcing(time_in_days, cfg: Config):
    """Take time in days and linearly interp 2009 Barrow air/snow/ice temperature data to get
    temperature in degrees Celsius.
    """
    barrow_days = cfg.forcing_config.barrow_days
    barrow_top_temp = cfg.forcing_config.barrow_top_temp
    return np.interp(time_in_days, barrow_days, barrow_top_temp, right=np.nan)

grids

Module providing functions to initialise the different grids and interpolate quantities between them.

Grids dataclass

Class initialised from number of grid cells to contain:

grid cell width, center, edge and ghost grids and difference matrices

Source code in seaice3p/grids.py 
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
@dataclass(frozen=True)
class Grids:
    """Class initialised from number of grid cells to contain:

    grid cell width, center, edge and ghost grids and difference matrices
    """

    number_of_cells: int

    @cached_property
    def step(self) -> float:
        """Grid cell width"""
        return 1 / self.number_of_cells

    @cached_property
    def centers(self) -> NDArray:
        """Center grid"""
        return np.array(
            [-1 + (2 * i + 1) * self.step / 2 for i in range(self.number_of_cells)]
        )

    @cached_property
    def edges(self) -> NDArray:
        """Edge grid"""
        return np.array([-1 + i * self.step for i in range(self.number_of_cells + 1)])

    @cached_property
    def ghosts(self) -> NDArray:
        """Ghost grid"""
        return np.concatenate(
            (np.array([-1 - self.step / 2]), self.centers, np.array([self.step / 2]))
        )

    @cached_property
    def D_e(self) -> NDArray:
        """Difference matrix to differentiate edge grid quantities to the center grid"""
        return get_difference_matrix(self.number_of_cells, self.step)

    @cached_property
    def D_g(self) -> NDArray:
        """Difference matrix to differentiate ghost grid quantities to the edge grid"""
        return get_difference_matrix(self.number_of_cells + 1, self.step)

D_e cached property

Difference matrix to differentiate edge grid quantities to the center grid

D_g cached property

Difference matrix to differentiate ghost grid quantities to the edge grid

centers cached property

Center grid

edges cached property

Edge grid

ghosts cached property

Ghost grid

step cached property

Grid cell width

add_ghost_cells(centers, bottom, top)

Add specified bottom and top value to center grid

:param centers: numpy array on centered grid (size I). :type centers: Numpy array :param bottom: bottom value placed at index 0. :type bottom: float :param top: top value placed at index -1. :type top: float :return: numpy array on ghost grid (size I+2).

Source code in seaice3p/grids.py 
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
def add_ghost_cells(centers, bottom, top):
    """Add specified bottom and top value to center grid

    :param centers: numpy array on centered grid (size I).
    :type centers: Numpy array
    :param bottom: bottom value placed at index 0.
    :type bottom: float
    :param top: top value placed at index -1.
    :type top: float
    :return: numpy array on ghost grid (size I+2).
    """
    return np.concatenate((np.array([bottom]), centers, np.array([top])))

average(points)

Returns arithmetic mean of adjacent points in an array

takes ghosts -> edges -> centers

Source code in seaice3p/grids.py 
79
80
81
82
83
84
85
86
def average(points: NDArray) -> NDArray:
    """Returns arithmetic mean of adjacent points in an array

    takes ghosts -> edges -> centers
    """
    upper = points[1:]
    lower = points[:-1]
    return 0.5 * (upper + lower)

calculate_ice_ocean_boundary_depth(liquid_fraction, edge_grid)

Calculate the depth of the ice ocean boundary as the edge position of the first cell from the bottom to be not completely liquid. I.e the first time the liquid fraction goes below 1.

If the ice has made it to the bottom of the domain raise an error.

If the domain is completely liquid set h=0.

NOTE: depth is a positive quantity and our grid coordinate increases from -1 at the bottom of the domain to 0 at the top.

:param liquid_fraction: liquid fraction on center grid :type liquid_fraction: Numpy Array (size I) :param edge_grid: The vertical coordinate positions of the edge grid. :type edge_grid: Numpy Array (size I+1) :return: positive depth value of ice ocean interface

Source code in seaice3p/grids.py 
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
def calculate_ice_ocean_boundary_depth(liquid_fraction, edge_grid):
    r"""Calculate the depth of the ice ocean boundary as the edge position of the
    first cell from the bottom to be not completely liquid. I.e the first time the
    liquid fraction goes below 1.

    If the ice has made it to the bottom of the domain raise an error.

    If the domain is completely liquid set h=0.

    NOTE: depth is a positive quantity and our grid coordinate increases from -1 at the
    bottom of the domain to 0 at the top.

    :param liquid_fraction: liquid fraction on center grid
    :type liquid_fraction: Numpy Array (size I)
    :param edge_grid: The vertical coordinate positions of the edge grid.
    :type edge_grid: Numpy Array (size I+1)
    :return: positive depth value of ice ocean interface
    """
    # locate index on center grid where liquid fraction first drops below 1
    index = np.argmax(liquid_fraction < 1)

    # if domain is completely liquid set h=0
    if np.all(liquid_fraction == 1):
        index = edge_grid.size - 1

    # raise error if bottom of domain freezes
    if index == 0:
        raise ValueError("Ice ocean interface has reached bottom of domain")

    # ice interface is at bottom edge of first frozen cell
    depth = (-1) * edge_grid[index]
    return depth

geometric(ghosts)

Returns geometric mean of the first dimension of an array

Source code in seaice3p/grids.py 
72
73
74
75
76
def geometric(ghosts):
    """Returns geometric mean of the first dimension of an array"""
    upper_ghosts = ghosts[1:]
    lower_ghosts = ghosts[:-1]
    return np.sqrt(upper_ghosts * lower_ghosts)

initial_conditions

Module to provide initial state of bulk enthalpy, bulk salinity and bulk gas for the simulation.

_apply_value_in_ice_layer(depth_of_ice, ice_value, liquid_value, grid)

assume that top part of domain contains mushy ice of given depth and lower part of domain is liquid. This function returns output on the given grid where the ice part of the domain takes one value and the liquid a different.

This is useful for initialising the barrow simulation where we have an initial ice layer.

Source code in seaice3p/initial_conditions.py 
46
47
48
49
50
51
52
53
54
55
def _apply_value_in_ice_layer(depth_of_ice, ice_value, liquid_value, grid):
    """assume that top part of domain contains mushy ice of given depth and lower part
    of domain is liquid. This function returns output on the given grid where the ice
    part of the domain takes one value and the liquid a different.

    This is useful for initialising the barrow simulation where we have an initial ice
    layer.
    """
    output = np.where(grid > -depth_of_ice, ice_value, liquid_value)
    return output

_get_barrow_initial_conditions(cfg)

initialise domain with an initial ice layer of given temperature and bulk salinity. These values are hard coded in from Moreau paper to match barrow study. They also assume that the initial ice layer has 1/5 the saturation amount in pure liquid of dissolved gas to account for previous gas loss.

Initialise with bulk gas being (1/5) in ice and saturation in liquid. Bulk salinity is 5.92 g/kg in ice and ocean value in liquid. Enthalpy is calculated by inverting temperature relation in ice and ocean. Ice temperature is given as -8.15 degC and ocean is the far value from boundary config.

Source code in seaice3p/initial_conditions.py 
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
def _get_barrow_initial_conditions(cfg: Config):
    """initialise domain with an initial ice layer of given temperature and bulk
    salinity. These values are hard coded in from Moreau paper to match barrow study.
    They also assume that the initial ice layer has 1/5 the saturation amount in pure
    liquid of dissolved gas to account for previous gas loss.

    Initialise with bulk gas being (1/5) in ice and saturation in liquid.
    Bulk salinity is 5.92 g/kg in ice and ocean value in liquid.
    Enthalpy is calculated by inverting temperature relation in ice and ocean.
    Ice temperature is given as -8.15 degC and ocean is the far value from boundary
    config.
    """
    far_gas_sat = cfg.ocean_forcing_config.ocean_gas_sat
    ICE_DEPTH = cfg.scales.convert_from_dimensional_grid(0.7)

    # if we are going to have brine convection ice will desalinate on its own
    if not isinstance(cfg.brine_convection_params, NoBrineConvection):
        SALT_IN_ICE = 0
    else:
        SALT_IN_ICE = cfg.scales.convert_from_dimensional_bulk_salinity(5.92)

    BOTTOM_TEMP = cfg.scales.convert_from_dimensional_temperature(-1.8)
    BOTTOM_SALT = 0
    TEMP_IN_ICE = cfg.scales.convert_from_dimensional_temperature(-8.15)

    chi = cfg.physical_params.expansion_coefficient

    centers = Grids(cfg.numerical_params.I).centers
    salt = _apply_value_in_ice_layer(
        ICE_DEPTH, ice_value=SALT_IN_ICE, liquid_value=BOTTOM_SALT, grid=centers
    )
    gas = _apply_value_in_ice_layer(
        ICE_DEPTH,
        ice_value=cfg.initial_conditions_config.Barrow_initial_bulk_gas_in_ice * chi,
        liquid_value=chi * far_gas_sat,
        grid=centers,
    )

    temp = _apply_value_in_ice_layer(
        ICE_DEPTH, ice_value=TEMP_IN_ICE, liquid_value=BOTTOM_TEMP, grid=centers
    )
    solid_fraction_in_mush = (salt + temp) / (
        temp - cfg.physical_params.concentration_ratio
    )
    enthalpy = _apply_value_in_ice_layer(
        ICE_DEPTH,
        ice_value=temp - solid_fraction_in_mush * cfg.physical_params.stefan_number,
        liquid_value=temp,
        grid=centers,
    )

    return _pack_initial_state(cfg, enthalpy, salt, gas)

_get_oil_initial_conditions(cfg)

initialise domain with an initial ice layer of given temperature and bulk salinity given by values in the configuration.

This is an idealised initial condition to investigate the impact of shortwave radiative forcing on melting bare ice

Source code in seaice3p/initial_conditions.py 
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
def _get_oil_initial_conditions(cfg: Config):
    """initialise domain with an initial ice layer of given temperature and bulk
    salinity given by values in the configuration.

    This is an idealised initial condition to investigate the impact of shortwave
    radiative forcing on melting bare ice
    """
    ICE_DEPTH = cfg.initial_conditions_config.initial_ice_depth

    # Initialise with a constant bulk salinity in ice
    SALT_IN_ICE = cfg.initial_conditions_config.initial_ice_bulk_salinity

    BOTTOM_TEMP = cfg.initial_conditions_config.initial_ocean_temperature
    BOTTOM_SALT = 0
    TEMP_IN_ICE = cfg.initial_conditions_config.initial_ice_temperature

    INITIAL_OIL_VOLUME_FRACTION = (
        cfg.initial_conditions_config.initial_oil_volume_fraction
    )

    centers = Grids(cfg.numerical_params.I).centers
    salt = _apply_value_in_ice_layer(
        ICE_DEPTH, ice_value=SALT_IN_ICE, liquid_value=BOTTOM_SALT, grid=centers
    )
    gas = np.where(
        centers < -cfg.initial_conditions_config.initial_oil_free_depth,
        INITIAL_OIL_VOLUME_FRACTION,
        0,
    )

    temp = _apply_value_in_ice_layer(
        ICE_DEPTH, ice_value=TEMP_IN_ICE, liquid_value=BOTTOM_TEMP, grid=centers
    )
    solid_fraction_in_mush = (salt + temp) / (
        temp - cfg.physical_params.concentration_ratio
    )
    enthalpy = _apply_value_in_ice_layer(
        ICE_DEPTH,
        ice_value=temp - solid_fraction_in_mush * cfg.physical_params.stefan_number,
        liquid_value=temp,
        grid=centers,
    )

    return _pack_initial_state(cfg, enthalpy, salt, gas)

_get_previous_simulation_final_state(cfg)

Generate initial state from the final state of a saved simulation

:returns: initial solution arrays on ghost grid (enthalpy, salt, gas)

Source code in seaice3p/initial_conditions.py 
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
def _get_previous_simulation_final_state(cfg: Config):
    """Generate initial state from the final state of a saved simulation

    :returns: initial solution arrays on ghost grid (enthalpy, salt, gas)
    """
    with np.load(cfg.initial_conditions_config.data_path) as data:
        match cfg.physical_params:
            case EQMPhysicalParams():
                enthalpy = data["arr_1"]
                salt = data["arr_2"]
                bulk_gas = data["arr_3"]

                return EQMState(0, enthalpy[:, -1], salt[:, -1], bulk_gas[:, -1])

            case DISEQPhysicalParams():
                enthalpy = data["arr_1"]
                salt = data["arr_2"]
                bulk_dissolved_gas = data["arr_3"]
                gas_fraction = data["arr_4"]

                return DISEQState(
                    0,
                    enthalpy[:, -1],
                    salt[:, -1],
                    bulk_dissolved_gas[:, -1],
                    gas_fraction[:, -1],
                )

            case _:
                raise NotImplementedError

_get_uniform_initial_conditions(cfg)

Generate uniform initial solution on the ghost grid

:returns: initial solution arrays on ghost grid (enthalpy, salt, gas)

Source code in seaice3p/initial_conditions.py 
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
def _get_uniform_initial_conditions(cfg: Config):
    """Generate uniform initial solution on the ghost grid

    :returns: initial solution arrays on ghost grid (enthalpy, salt, gas)
    """
    chi = cfg.physical_params.expansion_coefficient

    bottom_temp = cfg.ocean_forcing_config.ocean_temp
    bottom_bulk_salinity = 0
    bottom_dissolved_gas = cfg.ocean_forcing_config.ocean_gas_sat
    bottom_bulk_gas = bottom_dissolved_gas * chi

    # Initialise uniform enthalpy assuming completely liquid initial domain
    enthalpy = np.full((cfg.numerical_params.I,), bottom_temp)
    salt = np.full_like(enthalpy, bottom_bulk_salinity)
    gas = np.full_like(enthalpy, bottom_bulk_gas)

    return _pack_initial_state(cfg, enthalpy, salt, gas)

load

DISEQResults dataclass

Bases: _BaseResults

Source code in seaice3p/load.py 
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
@dataclass
class DISEQResults(_BaseResults):
    bulk_dissolved_gas: NDArray
    gas_fraction: NDArray

    def _get_state(self, time: float) -> DISEQStateFull:
        index = self._get_index(time)
        state = DISEQState(
            self.times[index],
            self.enthalpy[:, index],
            self.salt[:, index],
            self.bulk_dissolved_gas[:, index],
            self.gas_fraction[:, index],
        )

        enthalpy_method = get_enthalpy_method(self.cfg)
        return enthalpy_method(state)

    @property
    def bulk_gas(self) -> NDArray:
        """Dimensionless bulk gas the same as the EQM model"""
        return self.bulk_dissolved_gas + self.gas_fraction

bulk_gas property

Dimensionless bulk gas the same as the EQM model

_BaseResults dataclass

Source code in seaice3p/load.py 
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
@dataclass
class _BaseResults:
    cfg: Config
    dcfg: None | DimensionalParams
    times: NDArray
    enthalpy: NDArray
    salt: NDArray

    def __post_init__(self):
        self.states = list(map(self._get_state, self.times))

        boundary_conditions = get_boundary_conditions(self.cfg)
        self.states_bcs = list(map(boundary_conditions, self.states))

        self.grids = Grids(self.cfg.numerical_params.I)

    def _get_state(self, time: float) -> StateFull:
        raise NotImplementedError

    def _get_index(self, time: float) -> int:
        return np.argmin(np.abs(self.times - time))

    def _is_ice(self, time: float) -> NDArray:
        """Boolean mask True where ice is present on center grid cells at given
        non-dimensional time"""
        return self.liquid_fraction[:, self._get_index(time)] < 1

    def _top_cell_is_ice(self, time: float) -> bool:
        """Return True if top cell is ice or False if liquid"""
        return self._is_ice(time)[-1]

    @property
    def dates(self) -> List[datetime]:
        if hasattr(self.cfg.forcing_config, "start_date"):
            days = self.times * self.cfg.scales.time_scale
            start_date = datetime.strptime(
                self.cfg.forcing_config.start_date, "%Y-%m-%d"
            )
            return [timedelta(days=day) + start_date for day in days]
        else:
            raise AttributeError("forcing has no start date")

    @property
    def solid_fraction(self) -> NDArray:
        return _get_array_data("solid_fraction", self.states)

    @property
    def liquid_fraction(self) -> NDArray:
        return _get_array_data("liquid_fraction", self.states)

    @property
    def temperature(self) -> NDArray:
        return _get_array_data("temperature", self.states)

    @property
    def liquid_salinity(self) -> NDArray:
        return _get_array_data("liquid_salinity", self.states)

    @property
    def dissolved_gas(self) -> NDArray:
        return _get_array_data("dissolved_gas", self.states)

    @property
    def dimensional_meltpond_onset_time(self) -> float:
        """Get meltpond onset time from start of simulation in days"""
        top_liquid_fraction = self.liquid_fraction[-1, :]
        times_with_meltpond = self.times[top_liquid_fraction == 1]
        if times_with_meltpond.size == 0:
            return np.nan
        return self.cfg.scales.convert_to_dimensional_time(times_with_meltpond[0])

    @property
    def oil_mass_ratio(self) -> NDArray:
        """in ng/g"""
        return convert_gas_fraction_to_oil_mass_ratio(
            self.gas_fraction, self.cfg.scales.gas_density, self.cfg.scales.ice_density
        )

    @property
    def bulk_argon(self) -> NDArray:
        """in mircomole Ar/L"""
        scales = self.cfg.scales
        return scales.convert_dimensional_bulk_air_to_argon_content(
            scales.convert_to_dimensional_bulk_gas(self.bulk_gas)
        )

    def get_spectral_irradiance(self, time: float) -> oi.SixBandSpectralIrradiance:
        if not (
            isinstance(self.cfg.forcing_config, RadForcing)
            or isinstance(self.cfg.forcing_config, ERA5Forcing)
        ):
            raise TypeError("Simulation was not run with radiative forcing")

        return run_two_stream_model(
            self.states_bcs[self._get_index(time)], self.cfg, self.grids
        )

    def total_albedo(self, time: float) -> float:
        """Total albedo including the effect of the surface scattering layer if present,
        if not present then the penetration fraction is 1 and so we regain just albedo
        calculated from the two stream radiative transfer model"""
        spec_irrad = self.get_spectral_irradiance(time)
        return oi.integrate_over_SW(spec_irrad).albedo

    def total_transmittance(self, time: float) -> float:
        """Total spectrally integrated transmittance"""
        spec_irrad = self.get_spectral_irradiance(time)
        return oi.integrate_over_SW(spec_irrad).transmittance

    def ice_ocean_boundary(self, time: float) -> float:
        index = self._get_index(time)
        liquid_fraction = self.liquid_fraction[:, index]

        # if no ice then no boundary
        if np.all(liquid_fraction == 1):
            return np.nan

        is_ice_centers = self._is_ice(time)
        is_ice_edges = np.hstack((is_ice_centers, is_ice_centers[-1]))
        return self.grids.edges[is_ice_edges][0]

    def ice_meltpond_boundary(self, time: float) -> float:
        index = self._get_index(time)
        liquid_fraction = self.liquid_fraction[:, index]

        # if no ice then no meltpond
        if np.all(liquid_fraction == 1):
            return np.nan

        is_ice_centers = self._is_ice(time)
        is_ice_edges = np.hstack((is_ice_centers[0], is_ice_centers))
        return self.grids.edges[is_ice_edges][-1]

    def ice_thickness(self, time: float) -> float:
        index = self._get_index(time)
        liquid_fraction = self.liquid_fraction[:, index]

        # if no ice no thickness
        if np.all(liquid_fraction == 1):
            return 0
        return self.ice_meltpond_boundary(time) - self.ice_ocean_boundary(time)

    def integrated_solid_fraction(self, time: float) -> float:
        return _integrate(
            self.solid_fraction[:, self._get_index(time)],
            self.cfg.numerical_params.step,
        )

    @property
    def corrected_solid_fraction(self) -> NDArray:
        """Adjusted so that corrected_solid_fraction + corrected_liquid_fraction + gas_fraction = 1"""
        corrected_solid_fraction = np.empty_like(self.solid_fraction)
        is_frozen = self.liquid_fraction == 0

        corrected_solid_fraction[is_frozen] = (
            self.solid_fraction[is_frozen] - self.gas_fraction[is_frozen]
        )
        corrected_solid_fraction[~is_frozen] = self.solid_fraction[~is_frozen]
        if np.any(corrected_solid_fraction < 0):
            raise ValueError("Corrected solid fraction is negative")
        return corrected_solid_fraction

    @property
    def corrected_liquid_fraction(self) -> NDArray:
        """Adjusted so that corrected_solid_fraction + corrected_liquid_fraction + gas_fraction = 1"""
        corrected_liquid_fraction = np.empty_like(self.liquid_fraction)
        is_frozen = self.liquid_fraction == 0

        corrected_liquid_fraction[is_frozen] = self.liquid_fraction[is_frozen]
        corrected_liquid_fraction[~is_frozen] = (
            self.liquid_fraction[~is_frozen] - self.gas_fraction[~is_frozen]
        )
        if np.any(corrected_liquid_fraction < 0):
            raise ValueError("Corrected liquid fraction is negative")
        return corrected_liquid_fraction

    @property
    def dimensional_salinity_dependent_liquid_density(self) -> NDArray:
        reference_liquid_density = self.cfg.scales.liquid_density
        return reference_liquid_density * (
            1
            + self.cfg.scales.haline_contraction_coefficient
            * self.cfg.scales.salinity_difference
            * self.liquid_salinity
        )

    @property
    def dimensional_bulk_density(self) -> NDArray:
        return (
            (self.corrected_solid_fraction * self.cfg.scales.ice_density)
            + (
                self.corrected_liquid_fraction
                * self.dimensional_salinity_dependent_liquid_density
            )
            + (self.gas_fraction * self.cfg.scales.gas_density)
        )

    def dimensional_ice_average_bulk_density(self, time: float) -> float:
        index = self._get_index(time)
        is_ice = self._is_ice(time)
        bulk_density = self.dimensional_bulk_density[is_ice, index]
        if bulk_density.size == 0:
            return np.nan

        return np.mean(bulk_density)

    def total_bulk_gas_content(self, time: float) -> float:
        """To get dimensional bulk gas in domain multiply by
        gas_density * lengthscale
        """
        index = self._get_index(time)
        return _integrate(self.bulk_gas[:, index], self.cfg.numerical_params.step)

    def ice_bulk_gas_content(self, time: float) -> float:
        """To get dimensional bulk gas in ice multiply by
        gas_density * lengthscale
        """
        index = self._get_index(time)
        is_ice = self._is_ice(time)
        return _integrate(self.bulk_gas[is_ice, index], self.cfg.numerical_params.step)

    def surface_temp_K(self, time: float) -> float:
        """Return surface temperature in K"""
        index = self._get_index(time)
        ghost_cell_temp = find_ghost_cell_temperature(self.states[index], self.cfg)
        return _convert_non_dim_temperature_to_kelvin(
            self.cfg, 0.5 * (ghost_cell_temp + self.temperature[-1, index])
        )

    def sensible_heat_flux(self, time: float) -> float:
        """W/m2"""
        return calculate_sensible_heat_flux(
            self.cfg, time, self._top_cell_is_ice(time), self.surface_temp_K(time)
        )

    def latent_heat_flux(self, time: float) -> float:
        """W/m2"""
        return calculate_latent_heat_flux(
            self.cfg, time, self._top_cell_is_ice(time), self.surface_temp_K(time)
        )

    def emitted_LW(self, time: float) -> float:
        """W/m2 radiated away from ice surface"""
        emissivity = _calculate_emissivity(self.cfg, self._top_cell_is_ice(time))
        return emissivity * STEFAN_BOLTZMANN * self.surface_temp_K(time) ** 4

    def net_LW(self, time: float) -> float:
        """W/m2 net into ice"""
        incident_LW = get_LW_forcing(time, self.cfg)
        return incident_LW - self.emitted_LW(time)

    def surface_heat_flux(self, time: float) -> float:
        """W/m2 net into ice"""
        return (
            self.net_LW(time)
            + self.sensible_heat_flux(time)
            + self.latent_heat_flux(time)
        )

bulk_argon property

in mircomole Ar/L

corrected_liquid_fraction property

Adjusted so that corrected_solid_fraction + corrected_liquid_fraction + gas_fraction = 1

corrected_solid_fraction property

Adjusted so that corrected_solid_fraction + corrected_liquid_fraction + gas_fraction = 1

dimensional_meltpond_onset_time property

Get meltpond onset time from start of simulation in days

oil_mass_ratio property

in ng/g

_is_ice(time)

Boolean mask True where ice is present on center grid cells at given non-dimensional time

Source code in seaice3p/load.py 
60
61
62
63
def _is_ice(self, time: float) -> NDArray:
    """Boolean mask True where ice is present on center grid cells at given
    non-dimensional time"""
    return self.liquid_fraction[:, self._get_index(time)] < 1

_top_cell_is_ice(time)

Return True if top cell is ice or False if liquid

Source code in seaice3p/load.py 
65
66
67
def _top_cell_is_ice(self, time: float) -> bool:
    """Return True if top cell is ice or False if liquid"""
    return self._is_ice(time)[-1]

emitted_LW(time)

W/m2 radiated away from ice surface

Source code in seaice3p/load.py 
279
280
281
282
def emitted_LW(self, time: float) -> float:
    """W/m2 radiated away from ice surface"""
    emissivity = _calculate_emissivity(self.cfg, self._top_cell_is_ice(time))
    return emissivity * STEFAN_BOLTZMANN * self.surface_temp_K(time) ** 4

ice_bulk_gas_content(time)

To get dimensional bulk gas in ice multiply by gas_density * lengthscale

Source code in seaice3p/load.py 
251
252
253
254
255
256
257
def ice_bulk_gas_content(self, time: float) -> float:
    """To get dimensional bulk gas in ice multiply by
    gas_density * lengthscale
    """
    index = self._get_index(time)
    is_ice = self._is_ice(time)
    return _integrate(self.bulk_gas[is_ice, index], self.cfg.numerical_params.step)

latent_heat_flux(time)

W/m2

Source code in seaice3p/load.py 
273
274
275
276
277
def latent_heat_flux(self, time: float) -> float:
    """W/m2"""
    return calculate_latent_heat_flux(
        self.cfg, time, self._top_cell_is_ice(time), self.surface_temp_K(time)
    )

net_LW(time)

W/m2 net into ice

Source code in seaice3p/load.py 
284
285
286
287
def net_LW(self, time: float) -> float:
    """W/m2 net into ice"""
    incident_LW = get_LW_forcing(time, self.cfg)
    return incident_LW - self.emitted_LW(time)

sensible_heat_flux(time)

W/m2

Source code in seaice3p/load.py 
267
268
269
270
271
def sensible_heat_flux(self, time: float) -> float:
    """W/m2"""
    return calculate_sensible_heat_flux(
        self.cfg, time, self._top_cell_is_ice(time), self.surface_temp_K(time)
    )

surface_heat_flux(time)

W/m2 net into ice

Source code in seaice3p/load.py 
289
290
291
292
293
294
295
def surface_heat_flux(self, time: float) -> float:
    """W/m2 net into ice"""
    return (
        self.net_LW(time)
        + self.sensible_heat_flux(time)
        + self.latent_heat_flux(time)
    )

surface_temp_K(time)

Return surface temperature in K

Source code in seaice3p/load.py 
259
260
261
262
263
264
265
def surface_temp_K(self, time: float) -> float:
    """Return surface temperature in K"""
    index = self._get_index(time)
    ghost_cell_temp = find_ghost_cell_temperature(self.states[index], self.cfg)
    return _convert_non_dim_temperature_to_kelvin(
        self.cfg, 0.5 * (ghost_cell_temp + self.temperature[-1, index])
    )

total_albedo(time)

Total albedo including the effect of the surface scattering layer if present, if not present then the penetration fraction is 1 and so we regain just albedo calculated from the two stream radiative transfer model

Source code in seaice3p/load.py 
135
136
137
138
139
140
def total_albedo(self, time: float) -> float:
    """Total albedo including the effect of the surface scattering layer if present,
    if not present then the penetration fraction is 1 and so we regain just albedo
    calculated from the two stream radiative transfer model"""
    spec_irrad = self.get_spectral_irradiance(time)
    return oi.integrate_over_SW(spec_irrad).albedo

total_bulk_gas_content(time)

To get dimensional bulk gas in domain multiply by gas_density * lengthscale

Source code in seaice3p/load.py 
244
245
246
247
248
249
def total_bulk_gas_content(self, time: float) -> float:
    """To get dimensional bulk gas in domain multiply by
    gas_density * lengthscale
    """
    index = self._get_index(time)
    return _integrate(self.bulk_gas[:, index], self.cfg.numerical_params.step)

total_transmittance(time)

Total spectrally integrated transmittance

Source code in seaice3p/load.py 
142
143
144
145
def total_transmittance(self, time: float) -> float:
    """Total spectrally integrated transmittance"""
    spec_irrad = self.get_spectral_irradiance(time)
    return oi.integrate_over_SW(spec_irrad).transmittance

oil_mass

convert_gas_fraction_to_oil_mass_ratio(gas_fraction, oil_density, ice_density)

Convert gas (oil) volume fraction to oil mass ratio in ng/g

Source code in seaice3p/oil_mass.py 
4
5
6
7
8
def convert_gas_fraction_to_oil_mass_ratio(
    gas_fraction: NDArray, oil_density: float, ice_density: float
) -> NDArray:
    """Convert gas (oil) volume fraction to oil mass ratio in ng/g"""
    return gas_fraction * 1e9 * oil_density / ice_density

convert_oil_mass_ratio_to_gas_fraction(oil_mass_ratio, oil_density, ice_density)

Convert oil mass ratio in ng/g to gas (oil) volume fraction

Source code in seaice3p/oil_mass.py 
11
12
13
14
15
def convert_oil_mass_ratio_to_gas_fraction(
    oil_mass_ratio: NDArray, oil_density: float, ice_density: float
) -> NDArray:
    """Convert oil mass ratio in ng/g to gas (oil) volume fraction"""
    return oil_mass_ratio * 1e-9 * ice_density / oil_density

params

BRW09InitialConditions dataclass

values for bottom (ocean) boundary

Source code in seaice3p/params/dimensional/initial_conditions.py 
12
13
14
15
16
17
@serde(type_check=coerce)
@dataclass(frozen=True)
class BRW09InitialConditions:
    """values for bottom (ocean) boundary"""

    Barrow_initial_bulk_gas_in_ice: float = 1 / 5

CubicLiquidus dataclass

Cubic fit to liquidus to give liquidus salinity in terms of temperature

S = a0 + a1 T + a2 T^2 + a3 T^3

defaults are taken from Notz PhD thesis for fit to Assur seawater data

Source code in seaice3p/params/dimensional/water.py 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
@serde(type_check=coerce)
@dataclass(frozen=True)
class CubicLiquidus:
    """Cubic fit to liquidus to give liquidus salinity in terms of temperature

    S = a0 + a1 T + a2 T^2 + a3 T^3

    defaults are taken from Notz PhD thesis for fit to Assur seawater data
    """

    eutectic_temperature: float = -21.1  # deg Celsius
    a0: float = -1.2
    a1: float = -21.8
    a2: float = -0.919
    a3: float = -0.0178

    def get_liquidus_salinity(self, temperature):
        return (
            self.a0
            + self.a1 * temperature
            + self.a2 * temperature**2
            + self.a3 * temperature**3
        )

    def get_liquidus_temperature(self, salinity):
        temperature = fsolve(
            lambda x: salinity - self.get_liquidus_salinity(x),
            np.full_like(salinity, -2),
        )
        if temperature.size == 1:
            return temperature[0]
        else:
            return temperature

DimensionalBRW09OceanForcing dataclass

Ocean temperature provided by Barrow 2009 data at 2.4m and specify ocean fixed gas saturation state

Source code in seaice3p/params/dimensional/ocean_forcing.py 
56
57
58
59
60
61
62
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalBRW09OceanForcing:
    """Ocean temperature provided by Barrow 2009 data at 2.4m and specify ocean
    fixed gas saturation state"""

    pass

DimensionalConstantTurbulentFlux dataclass

Parameters for calculating the turbulent surface sensible and latent heat fluxes

NOTE: If you are running a simulation with ERA5 reanalysis forcing you must set the ref_height=2m as this is the appropriate value for the atmospheric reanalysis quantities

The windspeed given here will only be used with ERA5 forcing if the windspeed key is set to None in the forcing_data_file_keys dictionary.

Source code in seaice3p/params/dimensional/forcing.py 
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalConstantTurbulentFlux:
    """Parameters for calculating the turbulent surface sensible and latent heat
    fluxes

    NOTE: If you are running a simulation with ERA5 reanalysis forcing you must set
    the ref_height=2m as this is the appropriate value for the atmospheric reanalysis
    quantities

    The windspeed given here will only be used with ERA5 forcing if the windspeed key
    is set to None in the forcing_data_file_keys dictionary.
    """

    ref_height: float = 10  # m
    windspeed: float = 5  # m/s
    air_temp: float = 0  # deg C
    specific_humidity: float = 3.6e-3  # kg water / kg air
    atm_pressure: float = 101.325  # KPa

    air_density: float = 1.275  # kg/m3
    air_heat_capacity: float = 1005  # J/kg K
    air_latent_heat_of_vaporisation: float = 2.501e6  # J/kg

DimensionalERA5Forcing

read ERA5 data from netCDF file located at data_path.

Simulation will take atmospheric forcings from the start date specified in the string format YYYY-MM-DD

forcing_data_file_keys is a mapping of the descriptive names of the forcing data to be provided to the simulationa and the values are the corresponding strings giving the name of that variable in the netCDF file. The default values are the ERA5 variable names and the SnowModel-LG snow depth name.

Note that if you pass "sd" for the snow depth the simulation will assume you have provided snow depth in m of water equivalent and you must provide a snow density for the conversion.

If you pass None for the snow depth the simulation will procede with no snow layer. If you pass None for the windspeed the simulation will use the constant windspeed defined in the turbulent flux forcing parameters.

Source code in seaice3p/params/dimensional/forcing.py 
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
@serde(type_check=coerce)
class DimensionalERA5Forcing:
    """read ERA5 data from netCDF file located at data_path.

    Simulation will take atmospheric forcings from the start date specified in the
    string format YYYY-MM-DD

    forcing_data_file_keys is a mapping of the descriptive names of the
    forcing data to be provided to the simulationa and the values are the corresponding
    strings giving the name of that variable in the netCDF file.
    The default values are the ERA5 variable names and the SnowModel-LG snow depth name.

    Note that if you pass "sd" for the snow depth the simulation will assume you have
    provided snow depth in m of water equivalent and you must provide a snow density for
    the conversion.

    If you pass None for the snow depth the simulation will procede with no snow layer.
    If you pass None for the windspeed the simulation will use the constant windspeed
    defined in the turbulent flux forcing parameters.
    """

    data_path: Path
    start_date: str  # YYYY-MM-DD
    forcing_data_file_keys: ERA5FileKeys = ERA5FileKeys()
    SW_forcing: DimensionalSWForcing = DimensionalConstantSWForcing()
    LW_forcing: DimensionalLWForcing = DimensionalConstantLWForcing()
    turbulent_flux: DimensionalTurbulentFlux = DimensionalConstantTurbulentFlux()
    oil_heating: DimensionalOilHeating = DimensionalBackgroundOilHeating()

DimensionalFixedHeatFluxOceanForcing dataclass

Provides constant ocean heat flux at the bottom of the domain

Parameters:

Name Type Description Default
ocean_heat_flux float

The constant heat flux at the bottom of the domain in W/m2

 1
Source code in seaice3p/params/dimensional/ocean_forcing.py 
14
15
16
17
18
19
20
21
22
23
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalFixedHeatFluxOceanForcing:
    """Provides constant ocean heat flux at the bottom of the domain

    Args:
        ocean_heat_flux: The constant heat flux at the bottom of the domain in W/m2
    """

    ocean_heat_flux: float = 1

DimensionalFixedTempOceanForcing dataclass

Fixed temperature and gas saturation ocean boundary condition

Source code in seaice3p/params/dimensional/ocean_forcing.py 
 6
 7
 8
 9
10
11
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalFixedTempOceanForcing:
    """Fixed temperature and gas saturation ocean boundary condition"""

    ocean_temp: float = -1

DimensionalMonoBubbleParams dataclass

Bases: DimensionalBaseBubbleParams

Source code in seaice3p/params/dimensional/bubble.py 
15
16
17
18
19
20
21
22
23
24
25
26
27
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalMonoBubbleParams(DimensionalBaseBubbleParams):
    bubble_radius: float = 1e-3  # bubble radius in m

    @property
    def bubble_radius_scaled(self):
        r"""calculate the bubble radius divided by the pore scale

        .. math:: \Lambda = R_B / R_0

        """
        return self.bubble_radius / self.pore_radius

bubble_radius_scaled property

calculate the bubble radius divided by the pore scale

.. math:: \Lambda = R_B / R_0

DimensionalMonthlyHeatFluxOceanForcing dataclass

Provides constant ocean heat flux at the bottom of the domain in each month

Proivde an average monthly ocean heat flux with the entries i=0, 1, 2, 3, ...., 11 in the tuple corresponding to the months January, February, March, April, ...., December

Parameters:

Name Type Description Default
monthly_ocean_heat_flux Tuple[float, ...]

Tuple of ocean heat flux values in each month in W/m2

 (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0)
Source code in seaice3p/params/dimensional/ocean_forcing.py 
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalMonthlyHeatFluxOceanForcing:
    """Provides constant ocean heat flux at the bottom of the domain in each month

    Proivde an average monthly ocean heat flux with the entries
    i=0, 1, 2, 3, ...., 11
    in the tuple corresponding to the months
    January, February, March, April, ...., December

    Args:
        monthly_ocean_heat_flux: Tuple of ocean heat flux values in each month in W/m2
    """

    monthly_ocean_heat_flux: Tuple[float, ...] = (
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
    )

DimensionalParams dataclass

Contains all dimensional parameters needed to calculate non dimensional numbers.

To see the units each input should have look at the comment next to the default value.

Source code in seaice3p/params/dimensional/dimensional.py 
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalParams:
    """Contains all dimensional parameters needed to calculate non dimensional numbers.

    To see the units each input should have look at the comment next to the default
    value.
    """

    name: str
    total_time_in_days: float
    savefreq_in_days: float
    lengthscale: float

    gas_params: DimensionalEQMGasParams | DimensionalDISEQGasParams
    bubble_params: DimensionalMonoBubbleParams | DimensionalPowerLawBubbleParams
    brine_convection_params: DimensionalRJW14Params | NoBrineConvection
    forcing_config: DimensionalRadForcing | DimensionalBRW09Forcing | DimensionalConstantForcing | DimensionalYearlyForcing | DimensionalRobinForcing | DimensionalERA5Forcing
    ocean_forcing_config: DimensionalBRW09OceanForcing | DimensionalFixedTempOceanForcing | DimensionalFixedHeatFluxOceanForcing | DimensionalMonthlyHeatFluxOceanForcing
    initial_conditions_config: DimensionalOilInitialConditions | UniformInitialConditions | BRW09InitialConditions | PreviousSimulation

    water_params: DimensionalWaterParams = DimensionalWaterParams()
    numerical_params: NumericalParams = NumericalParams()
    frame_velocity_dimensional: float = 0  # velocity of frame in m/day
    gravity: float = 9.81  # m/s2

    @property
    def damkohler_number(self):
        r"""Return damkohler number as ratio of thermal timescale to nucleation
        timescale
        """
        if isinstance(self.gas_params, DimensionalEQMGasParams):
            return None

        return (
            (self.lengthscale**2) / self.water_params.thermal_diffusivity
        ) / self.gas_params.nucleation_timescale

    @property
    def total_time(self):
        """calculate the total time in non dimensional units for the simulation"""
        return self.total_time_in_days / self.scales.time_scale

    @property
    def savefreq(self):
        """calculate the save frequency in non dimensional time"""
        return self.savefreq_in_days / self.scales.time_scale

    @property
    def frame_velocity(self):
        """calculate the frame velocity in non dimensional units"""
        return self.frame_velocity_dimensional / self.scales.velocity_scale

    @property
    def B(self):
        r"""calculate the non dimensional scale for buoyant rise of gas bubbles as

        .. math:: \mathcal{B} = \frac{\rho_l g R_0^2 h}{3 \mu \kappa}

        """
        stokes_velocity = (
            (self.water_params.liquid_density - self.gas_params.gas_density)
            * self.gravity
            * self.bubble_params.pore_radius**2
            / (3 * self.water_params.liquid_viscosity)
        )
        velocity_scale_in_m_per_second = (
            self.water_params.thermal_diffusivity / self.lengthscale
        )
        return stokes_velocity / velocity_scale_in_m_per_second

    @property
    def Rayleigh_salt(self):
        r"""Calculate the haline Rayleigh number as

        .. math:: \text{Ra}_S = \frac{\rho_l g \beta \Delta S H K_0}{\kappa \mu}

        """
        match self.brine_convection_params:
            case DimensionalRJW14Params():
                return (
                    self.water_params.liquid_density
                    * self.gravity
                    * self.water_params.haline_contraction_coefficient
                    * self.water_params.salinity_difference
                    * self.lengthscale
                    * self.brine_convection_params.reference_permeability
                    / (
                        self.water_params.thermal_diffusivity
                        * self.water_params.liquid_viscosity
                    )
                )
            case NoBrineConvection():
                return None

    @property
    def expansion_coefficient(self):
        r"""calculate

        .. math:: \chi = \rho_l \xi_{\text{sat}} / \rho_g

        """
        return (
            self.water_params.liquid_density
            * self.gas_params.saturation_concentration
            / self.gas_params.gas_density
        )

    @property
    def lewis_gas(self):
        r"""Calculate the lewis number for dissolved gas, return np.inf if there is no
        dissolved gas diffusion.

        .. math:: \text{Le}_\xi = \kappa / D_\xi

        """
        if self.gas_params.gas_diffusivity == 0:
            return np.inf

        return self.water_params.thermal_diffusivity / self.gas_params.gas_diffusivity

    @property
    def scales(self):
        """return a Scales object used for converting between dimensional and non
        dimensional variables."""
        return Scales(
            self.lengthscale,
            self.water_params.thermal_diffusivity,
            self.water_params.liquid_thermal_conductivity,
            self.water_params.ocean_salinity,
            self.water_params.salinity_difference,
            self.water_params.ocean_freezing_temperature,
            self.water_params.temperature_difference,
            self.gas_params.gas_density,
            self.water_params.liquid_density,
            self.water_params.ice_density,
            self.gas_params.saturation_concentration,
            self.bubble_params.pore_radius,
            self.water_params.haline_contraction_coefficient,
        )

    def save(self, directory: Path):
        """save this object to a yaml file in the specified directory.

        The name will be the name given with _dimensional appended to distinguish it
        from a saved non-dimensional configuration."""
        with open(directory / f"{self.name}_dimensional.yml", "w") as outfile:
            outfile.write(to_yaml(self))

    @classmethod
    def load(cls, path):
        """load this object from a yaml configuration file."""
        with open(path, "r") as infile:
            yaml = infile.read()
        return from_yaml(cls, yaml)

B property

calculate the non dimensional scale for buoyant rise of gas bubbles as

.. math:: \mathcal{B} = \frac{\rho_l g R_0^2 h}{3 \mu \kappa}

Rayleigh_salt property

Calculate the haline Rayleigh number as

.. math:: \text{Ra}_S = \frac{\rho_l g \beta \Delta S H K_0}{\kappa \mu}

damkohler_number property

Return damkohler number as ratio of thermal timescale to nucleation timescale

expansion_coefficient property

calculate

.. math:: \chi = \rho_l \xi_{\text{sat}} / \rho_g

frame_velocity property

calculate the frame velocity in non dimensional units

lewis_gas property

Calculate the lewis number for dissolved gas, return np.inf if there is no dissolved gas diffusion.

.. math:: \text{Le}\xi = \kappa / D\xi

savefreq property

calculate the save frequency in non dimensional time

scales property

return a Scales object used for converting between dimensional and non dimensional variables.

total_time property

calculate the total time in non dimensional units for the simulation

load(path) classmethod

load this object from a yaml configuration file.

Source code in seaice3p/params/dimensional/dimensional.py 
198
199
200
201
202
203
@classmethod
def load(cls, path):
    """load this object from a yaml configuration file."""
    with open(path, "r") as infile:
        yaml = infile.read()
    return from_yaml(cls, yaml)

save(directory)

save this object to a yaml file in the specified directory.

The name will be the name given with _dimensional appended to distinguish it from a saved non-dimensional configuration.

Source code in seaice3p/params/dimensional/dimensional.py 
190
191
192
193
194
195
196
def save(self, directory: Path):
    """save this object to a yaml file in the specified directory.

    The name will be the name given with _dimensional appended to distinguish it
    from a saved non-dimensional configuration."""
    with open(directory / f"{self.name}_dimensional.yml", "w") as outfile:
        outfile.write(to_yaml(self))

DimensionalPowerLawBubbleParams dataclass

Bases: DimensionalBaseBubbleParams

Source code in seaice3p/params/dimensional/bubble.py 
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalPowerLawBubbleParams(DimensionalBaseBubbleParams):
    bubble_distribution_power: float = 1.5
    minimum_bubble_radius: float = 1e-6
    maximum_bubble_radius: float = 1e-3

    @property
    def minimum_bubble_radius_scaled(self):
        r"""calculate the bubble radius divided by the pore scale

        .. math:: \Lambda = R_B / R_0

        """
        return self.minimum_bubble_radius / self.pore_radius

    @property
    def maximum_bubble_radius_scaled(self):
        r"""calculate the bubble radius divided by the pore scale

        .. math:: \Lambda = R_B / R_0

        """
        return self.maximum_bubble_radius / self.pore_radius

maximum_bubble_radius_scaled property

calculate the bubble radius divided by the pore scale

.. math:: \Lambda = R_B / R_0

minimum_bubble_radius_scaled property

calculate the bubble radius divided by the pore scale

.. math:: \Lambda = R_B / R_0

DimensionalRobinForcing dataclass

This forcing imposes a Robin boundary condition of the form surface_heat_flux=heat_transfer_coefficient * (restoring_temp - surface_temp)

Source code in seaice3p/params/dimensional/forcing.py 
164
165
166
167
168
169
170
171
172
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalRobinForcing:
    """This forcing imposes a Robin boundary condition of the form
    surface_heat_flux=heat_transfer_coefficient * (restoring_temp - surface_temp)
    """

    heat_transfer_coefficient: float = 6.3  # W/m2K
    restoring_temperature: float = -30  # deg C

DimensionalWaterParams dataclass

Source code in seaice3p/params/dimensional/water.py 
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalWaterParams:
    liquid_density: float = 1028  # kg/m3
    ice_density: float = 916  # kg/m3
    ocean_salinity: float = 34  # g/kg
    liquidus: LinearLiquidus | CubicLiquidus = LinearLiquidus()
    latent_heat: float = 334e3  # latent heat of fusion for ice in J/kg
    liquid_specific_heat_capacity: float = 4184  # J/kg degC
    solid_specific_heat_capacity: float = 2009  # J/kg degC
    liquid_thermal_conductivity: float = 0.54  # water thermal conductivity in W/m deg C
    solid_thermal_conductivity: float = 2.22  # ice thermal conductivity in W/m deg C
    snow_thermal_conductivity: float = 0.31  # snow thermal conductivity in W/m deg C
    snow_density: float = 150  # snow density kg/m3

    eddy_diffusivity: float = 0

    salt_diffusivity: float = 0  # molecular diffusivity of salt in water in m2/s
    # used to calculate Rayleigh number for convection and density contraction in liquid equation of state
    haline_contraction_coefficient: float = 7.5e-4  # 1/ppt

    # calculated from moreau et al 2014 value of kinematic viscosity for sewater 2.7e-6
    # dynamic liquid_viscosity = 2.7e-6 * liquid_density
    liquid_viscosity: float = 2.78e-3  # dynamic liquid viscosity in Pa.s

    @property
    def eutectic_salinity(self):
        if isinstance(self.liquidus, LinearLiquidus):
            return self.liquidus.eutectic_salinity
        if isinstance(self.liquidus, CubicLiquidus):
            return self.liquidus.get_liquidus_salinity(
                self.liquidus.eutectic_temperature
            )

        raise NotImplementedError

    @property
    def eutectic_temperature(self):
        if isinstance(self.liquidus, LinearLiquidus) or isinstance(
            self.liquidus, CubicLiquidus
        ):
            return self.liquidus.eutectic_temperature

        raise NotImplementedError

    @property
    def salinity_difference(self):
        r"""calculate difference between eutectic salinity and typical ocean salinity

        .. math:: \Delta S = S_E - S_i

        """
        return self.eutectic_salinity - self.ocean_salinity

    @property
    def ocean_freezing_temperature(self):
        """calculate salinity dependent freezing temperature using linear liquidus with
        ocean salinity

        .. math:: T_i = T_L(S_i) = T_E S_i / S_E

        or using a cubic fit for the liquidus curve

        """
        if isinstance(self.liquidus, LinearLiquidus):
            return (
                self.eutectic_temperature * self.ocean_salinity / self.eutectic_salinity
            )
        if isinstance(self.liquidus, CubicLiquidus):
            return self.liquidus.get_liquidus_temperature(self.ocean_salinity)

        raise NotImplementedError

    @property
    def temperature_difference(self):
        r"""calculate

        .. math:: \Delta T = T_i - T_E

        """
        return self.ocean_freezing_temperature - self.eutectic_temperature

    @property
    def concentration_ratio(self):
        r"""Calculate concentration ratio as

        .. math:: \mathcal{C} = S_i / \Delta S

        """
        return self.ocean_salinity / self.salinity_difference

    @property
    def stefan_number(self):
        r"""calculate Stefan number

        .. math:: \text{St} = L / c_p \Delta T

        """
        return self.latent_heat / (
            self.temperature_difference * self.liquid_specific_heat_capacity
        )

    @property
    def thermal_diffusivity(self):
        r"""Return thermal diffusivity in m2/s

        .. math:: \kappa = \frac{k}{\rho_l c_p}

        """
        return self.liquid_thermal_conductivity / (
            self.liquid_density * self.liquid_specific_heat_capacity
        )

    @property
    def conductivity_ratio(self):
        r"""Calculate the ratio of solid to liquid thermal conductivity

        .. math:: \lambda = \frac{k_s}{k_l}

        """
        return self.solid_thermal_conductivity / self.liquid_thermal_conductivity

    @property
    def specific_heat_ratio(self):
        r"""Calculate the ratio of solid to liquid specific heat capacities

        .. math:: \lambda = \frac{c_{p,s}}{c_{p,l}}

        """
        return self.solid_specific_heat_capacity / self.liquid_specific_heat_capacity

    @property
    def eddy_diffusivity_ratio(self):
        r"""Calculate the ratio of eddy diffusivity to thermal diffusivity in
        the liquid phase

        .. math:: \lambda = \frac{\kappa_\text{turbulent}}{\kappa_l}

        """
        return self.eddy_diffusivity / self.thermal_diffusivity

    @property
    def snow_conductivity_ratio(self):
        r"""Calculate the ratio of snow to liquid thermal conductivity

        .. math:: \lambda = \frac{k_{sn}}{k_l}

        """
        return self.snow_thermal_conductivity / self.liquid_thermal_conductivity

    @property
    def lewis_salt(self):
        r"""Calculate the lewis number for salt, return np.inf if there is no salt
        diffusion.

        .. math:: \text{Le}_S = \kappa / D_s

        """
        if self.salt_diffusivity == 0:
            return np.inf

        return self.thermal_diffusivity / self.salt_diffusivity

concentration_ratio property

Calculate concentration ratio as

.. math:: \mathcal{C} = S_i / \Delta S

conductivity_ratio property

Calculate the ratio of solid to liquid thermal conductivity

.. math:: \lambda = \frac{k_s}{k_l}

eddy_diffusivity_ratio property

Calculate the ratio of eddy diffusivity to thermal diffusivity in the liquid phase

.. math:: \lambda = \frac{\kappa_\text{turbulent}}{\kappa_l}

lewis_salt property

Calculate the lewis number for salt, return np.inf if there is no salt diffusion.

.. math:: \text{Le}_S = \kappa / D_s

ocean_freezing_temperature property

calculate salinity dependent freezing temperature using linear liquidus with ocean salinity

.. math:: T_i = T_L(S_i) = T_E S_i / S_E

or using a cubic fit for the liquidus curve

salinity_difference property

calculate difference between eutectic salinity and typical ocean salinity

.. math:: \Delta S = S_E - S_i

snow_conductivity_ratio property

Calculate the ratio of snow to liquid thermal conductivity

.. math:: \lambda = \frac{k_{sn}}{k_l}

specific_heat_ratio property

Calculate the ratio of solid to liquid specific heat capacities

.. math:: \lambda = \frac{c_{p,s}}{c_{p,l}}

stefan_number property

calculate Stefan number

.. math:: \text{St} = L / c_p \Delta T

temperature_difference property

calculate

.. math:: \Delta T = T_i - T_E

thermal_diffusivity property

Return thermal diffusivity in m2/s

.. math:: \kappa = \frac{k}{\rho_l c_p}

NoBrineConvection dataclass

No brine convection

Source code in seaice3p/params/dimensional/convection.py 
5
6
7
8
@serde(type_check=coerce)
@dataclass(frozen=True)
class NoBrineConvection:
    """No brine convection"""

NumericalParams dataclass

parameters needed for discretisation and choice of numerical method

Source code in seaice3p/params/dimensional/numerical.py 
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
@serde(type_check=coerce)
@dataclass(frozen=True)
class NumericalParams:
    """parameters needed for discretisation and choice of numerical method"""

    I: int = 50
    regularisation: float = 1e-6
    solver_choice: str = "RK23"  # scipy.integrate.solve_IVP solver choice

    @property
    def step(self):
        return 1 / self.I

UniformInitialConditions dataclass

values for bottom (ocean) boundary

Source code in seaice3p/params/dimensional/initial_conditions.py 
6
7
8
9
@serde(type_check=coerce)
@dataclass(frozen=True)
class UniformInitialConditions:
    """values for bottom (ocean) boundary"""

bubble

BaseBubbleParams dataclass

Not to be used directly but provides parameters for bubble model in sea ice common to other bubble parameter objects.

Source code in seaice3p/params/bubble.py 
10
11
12
13
14
15
16
17
18
19
20
21
@serde(type_check=coerce)
@dataclass(frozen=True)
class BaseBubbleParams:
    """Not to be used directly but provides parameters for bubble model in sea ice
    common to other bubble parameter objects.
    """

    B: float = 100
    pore_throat_scaling: float = 0.46
    porosity_threshold: bool = False
    porosity_threshold_value: float = 0.024
    escape_ice_surface: bool = True

MonoBubbleParams dataclass

Bases: BaseBubbleParams

Parameters for population of identical spherical bubbles.

Source code in seaice3p/params/bubble.py 
24
25
26
27
28
29
@serde(type_check=coerce)
@dataclass(frozen=True)
class MonoBubbleParams(BaseBubbleParams):
    """Parameters for population of identical spherical bubbles."""

    bubble_radius_scaled: float = 1.0

PowerLawBubbleParams dataclass

Bases: BaseBubbleParams

Parameters for population of bubbles following a power law size distribution between a minimum and maximum radius.

Source code in seaice3p/params/bubble.py 
32
33
34
35
36
37
38
39
40
41
@serde(type_check=coerce)
@dataclass(frozen=True)
class PowerLawBubbleParams(BaseBubbleParams):
    """Parameters for population of bubbles following a power law size distribution
    between a minimum and maximum radius.
    """

    bubble_distribution_power: float = 1.5
    minimum_bubble_radius_scaled: float = 1e-3
    maximum_bubble_radius_scaled: float = 1

convection

RJW14Params dataclass

Parameters for the RJW14 parameterisation of brine convection

Source code in seaice3p/params/convection.py 
 6
 7
 8
 9
10
11
12
13
14
15
16
@serde(type_check=coerce)
@dataclass(frozen=True)
class RJW14Params:
    """Parameters for the RJW14 parameterisation of brine convection"""

    Rayleigh_salt: float = 44105
    Rayleigh_critical: float = 2.9
    convection_strength: float = 0.13
    couple_bubble_to_horizontal_flow: bool = False
    couple_bubble_to_vertical_flow: bool = False
    advective_heat_flux_in_ocean: bool = True

convert

Scales dataclass

Source code in seaice3p/params/convert.py 
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
@serde(type_check=coerce)
@dataclass(frozen=True)
class Scales:
    lengthscale: float  # domain height in m
    thermal_diffusivity: float  # m2/s
    liquid_thermal_conductivity: float  # W/m deg C
    ocean_salinity: float  # g/kg
    salinity_difference: float  # g/kg
    ocean_freezing_temperature: float  # deg C
    temperature_difference: float  # deg C
    gas_density: float  # kg/m3
    liquid_density: float  # kg/m3
    ice_density: float  # kg/m3
    saturation_concentration: float  # kg(gas)/kg(liquid)
    pore_radius: float  # m
    haline_contraction_coefficient: float  # 1/ppt

    @property
    def time_scale(self):
        """in days"""
        return SECONDS_TO_DAYS * self.lengthscale**2 / self.thermal_diffusivity

    @property
    def velocity_scale(self):
        """in m /day"""
        return self.lengthscale / self.time_scale

    def convert_from_dimensional_temperature(self, dimensional_temperature):
        """Non dimensionalise temperature in deg C"""
        return (
            dimensional_temperature - self.ocean_freezing_temperature
        ) / self.temperature_difference

    def convert_to_dimensional_temperature(self, temperature):
        """get temperature in deg C from non dimensional temperature"""
        return (
            self.temperature_difference * temperature + self.ocean_freezing_temperature
        )

    def convert_from_dimensional_grid(self, dimensional_grid):
        """Non dimensionalise domain depths in meters"""
        return dimensional_grid / self.lengthscale

    def convert_to_dimensional_grid(self, grid):
        """Get domain depths in meters from non dimensional values"""
        return self.lengthscale * grid

    def convert_from_dimensional_time(self, dimensional_time):
        """Non dimensionalise time in days"""
        return dimensional_time / self.time_scale

    def convert_to_dimensional_time(self, time):
        """Convert non dimensional time into time in days since start of simulation"""
        return self.time_scale * time

    def convert_from_dimensional_bulk_salinity(self, dimensional_bulk_salinity):
        """Non dimensionalise bulk salinity in g/kg"""
        return (
            dimensional_bulk_salinity - self.ocean_salinity
        ) / self.salinity_difference

    def convert_to_dimensional_bulk_salinity(self, bulk_salinity):
        """Convert non dimensional bulk salinity to g/kg"""
        return self.salinity_difference * bulk_salinity + self.ocean_salinity

    def convert_from_dimensional_bulk_gas(self, dimensional_bulk_gas):
        """Non dimensionalise bulk gas content in kg/m3"""
        return dimensional_bulk_gas / self.gas_density

    def convert_to_dimensional_bulk_gas(self, bulk_gas):
        """Convert dimensionless bulk gas content to kg/m3"""
        return self.gas_density * bulk_gas

    def convert_dimensional_bulk_air_to_argon_content(self, dimensional_bulk_gas):
        """Convert kg/m3 of air to micromole of Argon per Liter of ice"""
        mass_ratio_of_argon_in_air = 0.01288
        micromoles_of_argon_in_a_kilogram_of_argon = 1 / (3.9948e-8)
        liters_in_a_meter_cubed = 1e3
        return (
            dimensional_bulk_gas
            * mass_ratio_of_argon_in_air
            * micromoles_of_argon_in_a_kilogram_of_argon
            / liters_in_a_meter_cubed
        )

    def convert_from_dimensional_dissolved_gas(self, dimensional_dissolved_gas):
        """convert from dissolved gas in kg(gas)/kg(liquid) to dimensionless"""
        return dimensional_dissolved_gas / self.saturation_concentration

    def convert_to_dimensional_dissolved_gas(self, dissolved_gas):
        """convert from non dimensional dissolved gas to dimensional dissolved gas in
        kg(gas)/kg(liquid)"""
        return self.saturation_concentration * dissolved_gas

    def convert_from_dimensional_heating(self, dimensional_heating):
        """convert from heating rate in W/m3 to dimensionless units"""
        return (
            dimensional_heating
            * self.lengthscale**2
            / (self.liquid_thermal_conductivity * self.temperature_difference)
        )

    def convert_from_dimensional_heat_flux(self, dimensional_heat_flux):
        """convert from heat flux in W/m2 to dimensionless units"""
        return (
            dimensional_heat_flux
            * self.lengthscale
            / (self.liquid_thermal_conductivity * self.temperature_difference)
        )

    def convert_to_dimensional_heat_flux(self, heat_flux):
        """convert from dimensionless heat flux to heat flux in W/m2"""
        return (
            heat_flux
            * (self.liquid_thermal_conductivity * self.temperature_difference)
            / self.lengthscale
        )
time_scale property

in days

velocity_scale property

in m /day

convert_dimensional_bulk_air_to_argon_content(dimensional_bulk_gas)

Convert kg/m3 of air to micromole of Argon per Liter of ice

Source code in seaice3p/params/convert.py 
81
82
83
84
85
86
87
88
89
90
91
def convert_dimensional_bulk_air_to_argon_content(self, dimensional_bulk_gas):
    """Convert kg/m3 of air to micromole of Argon per Liter of ice"""
    mass_ratio_of_argon_in_air = 0.01288
    micromoles_of_argon_in_a_kilogram_of_argon = 1 / (3.9948e-8)
    liters_in_a_meter_cubed = 1e3
    return (
        dimensional_bulk_gas
        * mass_ratio_of_argon_in_air
        * micromoles_of_argon_in_a_kilogram_of_argon
        / liters_in_a_meter_cubed
    )
convert_from_dimensional_bulk_gas(dimensional_bulk_gas)

Non dimensionalise bulk gas content in kg/m3

Source code in seaice3p/params/convert.py 
73
74
75
def convert_from_dimensional_bulk_gas(self, dimensional_bulk_gas):
    """Non dimensionalise bulk gas content in kg/m3"""
    return dimensional_bulk_gas / self.gas_density
convert_from_dimensional_bulk_salinity(dimensional_bulk_salinity)

Non dimensionalise bulk salinity in g/kg

Source code in seaice3p/params/convert.py 
63
64
65
66
67
def convert_from_dimensional_bulk_salinity(self, dimensional_bulk_salinity):
    """Non dimensionalise bulk salinity in g/kg"""
    return (
        dimensional_bulk_salinity - self.ocean_salinity
    ) / self.salinity_difference
convert_from_dimensional_dissolved_gas(dimensional_dissolved_gas)

convert from dissolved gas in kg(gas)/kg(liquid) to dimensionless

Source code in seaice3p/params/convert.py 
93
94
95
def convert_from_dimensional_dissolved_gas(self, dimensional_dissolved_gas):
    """convert from dissolved gas in kg(gas)/kg(liquid) to dimensionless"""
    return dimensional_dissolved_gas / self.saturation_concentration
convert_from_dimensional_grid(dimensional_grid)

Non dimensionalise domain depths in meters

Source code in seaice3p/params/convert.py 
47
48
49
def convert_from_dimensional_grid(self, dimensional_grid):
    """Non dimensionalise domain depths in meters"""
    return dimensional_grid / self.lengthscale
convert_from_dimensional_heat_flux(dimensional_heat_flux)

convert from heat flux in W/m2 to dimensionless units

Source code in seaice3p/params/convert.py 
110
111
112
113
114
115
116
def convert_from_dimensional_heat_flux(self, dimensional_heat_flux):
    """convert from heat flux in W/m2 to dimensionless units"""
    return (
        dimensional_heat_flux
        * self.lengthscale
        / (self.liquid_thermal_conductivity * self.temperature_difference)
    )
convert_from_dimensional_heating(dimensional_heating)

convert from heating rate in W/m3 to dimensionless units

Source code in seaice3p/params/convert.py 
102
103
104
105
106
107
108
def convert_from_dimensional_heating(self, dimensional_heating):
    """convert from heating rate in W/m3 to dimensionless units"""
    return (
        dimensional_heating
        * self.lengthscale**2
        / (self.liquid_thermal_conductivity * self.temperature_difference)
    )
convert_from_dimensional_temperature(dimensional_temperature)

Non dimensionalise temperature in deg C

Source code in seaice3p/params/convert.py 
35
36
37
38
39
def convert_from_dimensional_temperature(self, dimensional_temperature):
    """Non dimensionalise temperature in deg C"""
    return (
        dimensional_temperature - self.ocean_freezing_temperature
    ) / self.temperature_difference
convert_from_dimensional_time(dimensional_time)

Non dimensionalise time in days

Source code in seaice3p/params/convert.py 
55
56
57
def convert_from_dimensional_time(self, dimensional_time):
    """Non dimensionalise time in days"""
    return dimensional_time / self.time_scale
convert_to_dimensional_bulk_gas(bulk_gas)

Convert dimensionless bulk gas content to kg/m3

Source code in seaice3p/params/convert.py 
77
78
79
def convert_to_dimensional_bulk_gas(self, bulk_gas):
    """Convert dimensionless bulk gas content to kg/m3"""
    return self.gas_density * bulk_gas
convert_to_dimensional_bulk_salinity(bulk_salinity)

Convert non dimensional bulk salinity to g/kg

Source code in seaice3p/params/convert.py 
69
70
71
def convert_to_dimensional_bulk_salinity(self, bulk_salinity):
    """Convert non dimensional bulk salinity to g/kg"""
    return self.salinity_difference * bulk_salinity + self.ocean_salinity
convert_to_dimensional_dissolved_gas(dissolved_gas)

convert from non dimensional dissolved gas to dimensional dissolved gas in kg(gas)/kg(liquid)

Source code in seaice3p/params/convert.py 
 97
 98
 99
100
def convert_to_dimensional_dissolved_gas(self, dissolved_gas):
    """convert from non dimensional dissolved gas to dimensional dissolved gas in
    kg(gas)/kg(liquid)"""
    return self.saturation_concentration * dissolved_gas
convert_to_dimensional_grid(grid)

Get domain depths in meters from non dimensional values

Source code in seaice3p/params/convert.py 
51
52
53
def convert_to_dimensional_grid(self, grid):
    """Get domain depths in meters from non dimensional values"""
    return self.lengthscale * grid
convert_to_dimensional_heat_flux(heat_flux)

convert from dimensionless heat flux to heat flux in W/m2

Source code in seaice3p/params/convert.py 
118
119
120
121
122
123
124
def convert_to_dimensional_heat_flux(self, heat_flux):
    """convert from dimensionless heat flux to heat flux in W/m2"""
    return (
        heat_flux
        * (self.liquid_thermal_conductivity * self.temperature_difference)
        / self.lengthscale
    )
convert_to_dimensional_temperature(temperature)

get temperature in deg C from non dimensional temperature

Source code in seaice3p/params/convert.py 
41
42
43
44
45
def convert_to_dimensional_temperature(self, temperature):
    """get temperature in deg C from non dimensional temperature"""
    return (
        self.temperature_difference * temperature + self.ocean_freezing_temperature
    )
convert_to_dimensional_time(time)

Convert non dimensional time into time in days since start of simulation

Source code in seaice3p/params/convert.py 
59
60
61
def convert_to_dimensional_time(self, time):
    """Convert non dimensional time into time in days since start of simulation"""
    return self.time_scale * time

dimensional

bubble

DimensionalMonoBubbleParams dataclass

Bases: DimensionalBaseBubbleParams

Source code in seaice3p/params/dimensional/bubble.py 
15
16
17
18
19
20
21
22
23
24
25
26
27
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalMonoBubbleParams(DimensionalBaseBubbleParams):
    bubble_radius: float = 1e-3  # bubble radius in m

    @property
    def bubble_radius_scaled(self):
        r"""calculate the bubble radius divided by the pore scale

        .. math:: \Lambda = R_B / R_0

        """
        return self.bubble_radius / self.pore_radius
bubble_radius_scaled property

calculate the bubble radius divided by the pore scale

.. math:: \Lambda = R_B / R_0

DimensionalPowerLawBubbleParams dataclass

Bases: DimensionalBaseBubbleParams

Source code in seaice3p/params/dimensional/bubble.py 
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalPowerLawBubbleParams(DimensionalBaseBubbleParams):
    bubble_distribution_power: float = 1.5
    minimum_bubble_radius: float = 1e-6
    maximum_bubble_radius: float = 1e-3

    @property
    def minimum_bubble_radius_scaled(self):
        r"""calculate the bubble radius divided by the pore scale

        .. math:: \Lambda = R_B / R_0

        """
        return self.minimum_bubble_radius / self.pore_radius

    @property
    def maximum_bubble_radius_scaled(self):
        r"""calculate the bubble radius divided by the pore scale

        .. math:: \Lambda = R_B / R_0

        """
        return self.maximum_bubble_radius / self.pore_radius
maximum_bubble_radius_scaled property

calculate the bubble radius divided by the pore scale

.. math:: \Lambda = R_B / R_0

minimum_bubble_radius_scaled property

calculate the bubble radius divided by the pore scale

.. math:: \Lambda = R_B / R_0

convection

NoBrineConvection dataclass

No brine convection

Source code in seaice3p/params/dimensional/convection.py 
5
6
7
8
@serde(type_check=coerce)
@dataclass(frozen=True)
class NoBrineConvection:
    """No brine convection"""

dimensional

Dimensional parameters required to run a simulation and convert output to dimensional variables.

The DimensionalParams class contains all the dimensional parameters needed to produce a simulation configuration.

The Scales class contains all the dimensional parameters required to convert simulation output between physical and non-dimensional variables.

DimensionalParams dataclass

Contains all dimensional parameters needed to calculate non dimensional numbers.

To see the units each input should have look at the comment next to the default value.

Source code in seaice3p/params/dimensional/dimensional.py 
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalParams:
    """Contains all dimensional parameters needed to calculate non dimensional numbers.

    To see the units each input should have look at the comment next to the default
    value.
    """

    name: str
    total_time_in_days: float
    savefreq_in_days: float
    lengthscale: float

    gas_params: DimensionalEQMGasParams | DimensionalDISEQGasParams
    bubble_params: DimensionalMonoBubbleParams | DimensionalPowerLawBubbleParams
    brine_convection_params: DimensionalRJW14Params | NoBrineConvection
    forcing_config: DimensionalRadForcing | DimensionalBRW09Forcing | DimensionalConstantForcing | DimensionalYearlyForcing | DimensionalRobinForcing | DimensionalERA5Forcing
    ocean_forcing_config: DimensionalBRW09OceanForcing | DimensionalFixedTempOceanForcing | DimensionalFixedHeatFluxOceanForcing | DimensionalMonthlyHeatFluxOceanForcing
    initial_conditions_config: DimensionalOilInitialConditions | UniformInitialConditions | BRW09InitialConditions | PreviousSimulation

    water_params: DimensionalWaterParams = DimensionalWaterParams()
    numerical_params: NumericalParams = NumericalParams()
    frame_velocity_dimensional: float = 0  # velocity of frame in m/day
    gravity: float = 9.81  # m/s2

    @property
    def damkohler_number(self):
        r"""Return damkohler number as ratio of thermal timescale to nucleation
        timescale
        """
        if isinstance(self.gas_params, DimensionalEQMGasParams):
            return None

        return (
            (self.lengthscale**2) / self.water_params.thermal_diffusivity
        ) / self.gas_params.nucleation_timescale

    @property
    def total_time(self):
        """calculate the total time in non dimensional units for the simulation"""
        return self.total_time_in_days / self.scales.time_scale

    @property
    def savefreq(self):
        """calculate the save frequency in non dimensional time"""
        return self.savefreq_in_days / self.scales.time_scale

    @property
    def frame_velocity(self):
        """calculate the frame velocity in non dimensional units"""
        return self.frame_velocity_dimensional / self.scales.velocity_scale

    @property
    def B(self):
        r"""calculate the non dimensional scale for buoyant rise of gas bubbles as

        .. math:: \mathcal{B} = \frac{\rho_l g R_0^2 h}{3 \mu \kappa}

        """
        stokes_velocity = (
            (self.water_params.liquid_density - self.gas_params.gas_density)
            * self.gravity
            * self.bubble_params.pore_radius**2
            / (3 * self.water_params.liquid_viscosity)
        )
        velocity_scale_in_m_per_second = (
            self.water_params.thermal_diffusivity / self.lengthscale
        )
        return stokes_velocity / velocity_scale_in_m_per_second

    @property
    def Rayleigh_salt(self):
        r"""Calculate the haline Rayleigh number as

        .. math:: \text{Ra}_S = \frac{\rho_l g \beta \Delta S H K_0}{\kappa \mu}

        """
        match self.brine_convection_params:
            case DimensionalRJW14Params():
                return (
                    self.water_params.liquid_density
                    * self.gravity
                    * self.water_params.haline_contraction_coefficient
                    * self.water_params.salinity_difference
                    * self.lengthscale
                    * self.brine_convection_params.reference_permeability
                    / (
                        self.water_params.thermal_diffusivity
                        * self.water_params.liquid_viscosity
                    )
                )
            case NoBrineConvection():
                return None

    @property
    def expansion_coefficient(self):
        r"""calculate

        .. math:: \chi = \rho_l \xi_{\text{sat}} / \rho_g

        """
        return (
            self.water_params.liquid_density
            * self.gas_params.saturation_concentration
            / self.gas_params.gas_density
        )

    @property
    def lewis_gas(self):
        r"""Calculate the lewis number for dissolved gas, return np.inf if there is no
        dissolved gas diffusion.

        .. math:: \text{Le}_\xi = \kappa / D_\xi

        """
        if self.gas_params.gas_diffusivity == 0:
            return np.inf

        return self.water_params.thermal_diffusivity / self.gas_params.gas_diffusivity

    @property
    def scales(self):
        """return a Scales object used for converting between dimensional and non
        dimensional variables."""
        return Scales(
            self.lengthscale,
            self.water_params.thermal_diffusivity,
            self.water_params.liquid_thermal_conductivity,
            self.water_params.ocean_salinity,
            self.water_params.salinity_difference,
            self.water_params.ocean_freezing_temperature,
            self.water_params.temperature_difference,
            self.gas_params.gas_density,
            self.water_params.liquid_density,
            self.water_params.ice_density,
            self.gas_params.saturation_concentration,
            self.bubble_params.pore_radius,
            self.water_params.haline_contraction_coefficient,
        )

    def save(self, directory: Path):
        """save this object to a yaml file in the specified directory.

        The name will be the name given with _dimensional appended to distinguish it
        from a saved non-dimensional configuration."""
        with open(directory / f"{self.name}_dimensional.yml", "w") as outfile:
            outfile.write(to_yaml(self))

    @classmethod
    def load(cls, path):
        """load this object from a yaml configuration file."""
        with open(path, "r") as infile:
            yaml = infile.read()
        return from_yaml(cls, yaml)
B property

calculate the non dimensional scale for buoyant rise of gas bubbles as

.. math:: \mathcal{B} = \frac{\rho_l g R_0^2 h}{3 \mu \kappa}

Rayleigh_salt property

Calculate the haline Rayleigh number as

.. math:: \text{Ra}_S = \frac{\rho_l g \beta \Delta S H K_0}{\kappa \mu}

damkohler_number property

Return damkohler number as ratio of thermal timescale to nucleation timescale

expansion_coefficient property

calculate

.. math:: \chi = \rho_l \xi_{\text{sat}} / \rho_g

frame_velocity property

calculate the frame velocity in non dimensional units

lewis_gas property

Calculate the lewis number for dissolved gas, return np.inf if there is no dissolved gas diffusion.

.. math:: \text{Le}\xi = \kappa / D\xi

savefreq property

calculate the save frequency in non dimensional time

scales property

return a Scales object used for converting between dimensional and non dimensional variables.

total_time property

calculate the total time in non dimensional units for the simulation

load(path) classmethod

load this object from a yaml configuration file.

Source code in seaice3p/params/dimensional/dimensional.py 
198
199
200
201
202
203
@classmethod
def load(cls, path):
    """load this object from a yaml configuration file."""
    with open(path, "r") as infile:
        yaml = infile.read()
    return from_yaml(cls, yaml)
save(directory)

save this object to a yaml file in the specified directory.

The name will be the name given with _dimensional appended to distinguish it from a saved non-dimensional configuration.

Source code in seaice3p/params/dimensional/dimensional.py 
190
191
192
193
194
195
196
def save(self, directory: Path):
    """save this object to a yaml file in the specified directory.

    The name will be the name given with _dimensional appended to distinguish it
    from a saved non-dimensional configuration."""
    with open(directory / f"{self.name}_dimensional.yml", "w") as outfile:
        outfile.write(to_yaml(self))

forcing

DimensionalConstantTurbulentFlux dataclass

Parameters for calculating the turbulent surface sensible and latent heat fluxes

NOTE: If you are running a simulation with ERA5 reanalysis forcing you must set the ref_height=2m as this is the appropriate value for the atmospheric reanalysis quantities

The windspeed given here will only be used with ERA5 forcing if the windspeed key is set to None in the forcing_data_file_keys dictionary.

Source code in seaice3p/params/dimensional/forcing.py 
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalConstantTurbulentFlux:
    """Parameters for calculating the turbulent surface sensible and latent heat
    fluxes

    NOTE: If you are running a simulation with ERA5 reanalysis forcing you must set
    the ref_height=2m as this is the appropriate value for the atmospheric reanalysis
    quantities

    The windspeed given here will only be used with ERA5 forcing if the windspeed key
    is set to None in the forcing_data_file_keys dictionary.
    """

    ref_height: float = 10  # m
    windspeed: float = 5  # m/s
    air_temp: float = 0  # deg C
    specific_humidity: float = 3.6e-3  # kg water / kg air
    atm_pressure: float = 101.325  # KPa

    air_density: float = 1.275  # kg/m3
    air_heat_capacity: float = 1005  # J/kg K
    air_latent_heat_of_vaporisation: float = 2.501e6  # J/kg
DimensionalERA5Forcing

read ERA5 data from netCDF file located at data_path.

Simulation will take atmospheric forcings from the start date specified in the string format YYYY-MM-DD

forcing_data_file_keys is a mapping of the descriptive names of the forcing data to be provided to the simulationa and the values are the corresponding strings giving the name of that variable in the netCDF file. The default values are the ERA5 variable names and the SnowModel-LG snow depth name.

Note that if you pass "sd" for the snow depth the simulation will assume you have provided snow depth in m of water equivalent and you must provide a snow density for the conversion.

If you pass None for the snow depth the simulation will procede with no snow layer. If you pass None for the windspeed the simulation will use the constant windspeed defined in the turbulent flux forcing parameters.

Source code in seaice3p/params/dimensional/forcing.py 
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
@serde(type_check=coerce)
class DimensionalERA5Forcing:
    """read ERA5 data from netCDF file located at data_path.

    Simulation will take atmospheric forcings from the start date specified in the
    string format YYYY-MM-DD

    forcing_data_file_keys is a mapping of the descriptive names of the
    forcing data to be provided to the simulationa and the values are the corresponding
    strings giving the name of that variable in the netCDF file.
    The default values are the ERA5 variable names and the SnowModel-LG snow depth name.

    Note that if you pass "sd" for the snow depth the simulation will assume you have
    provided snow depth in m of water equivalent and you must provide a snow density for
    the conversion.

    If you pass None for the snow depth the simulation will procede with no snow layer.
    If you pass None for the windspeed the simulation will use the constant windspeed
    defined in the turbulent flux forcing parameters.
    """

    data_path: Path
    start_date: str  # YYYY-MM-DD
    forcing_data_file_keys: ERA5FileKeys = ERA5FileKeys()
    SW_forcing: DimensionalSWForcing = DimensionalConstantSWForcing()
    LW_forcing: DimensionalLWForcing = DimensionalConstantLWForcing()
    turbulent_flux: DimensionalTurbulentFlux = DimensionalConstantTurbulentFlux()
    oil_heating: DimensionalOilHeating = DimensionalBackgroundOilHeating()
DimensionalRobinForcing dataclass

This forcing imposes a Robin boundary condition of the form surface_heat_flux=heat_transfer_coefficient * (restoring_temp - surface_temp)

Source code in seaice3p/params/dimensional/forcing.py 
164
165
166
167
168
169
170
171
172
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalRobinForcing:
    """This forcing imposes a Robin boundary condition of the form
    surface_heat_flux=heat_transfer_coefficient * (restoring_temp - surface_temp)
    """

    heat_transfer_coefficient: float = 6.3  # W/m2K
    restoring_temperature: float = -30  # deg C

initial_conditions

BRW09InitialConditions dataclass

values for bottom (ocean) boundary

Source code in seaice3p/params/dimensional/initial_conditions.py 
12
13
14
15
16
17
@serde(type_check=coerce)
@dataclass(frozen=True)
class BRW09InitialConditions:
    """values for bottom (ocean) boundary"""

    Barrow_initial_bulk_gas_in_ice: float = 1 / 5
UniformInitialConditions dataclass

values for bottom (ocean) boundary

Source code in seaice3p/params/dimensional/initial_conditions.py 
6
7
8
9
@serde(type_check=coerce)
@dataclass(frozen=True)
class UniformInitialConditions:
    """values for bottom (ocean) boundary"""

numerical

NumericalParams dataclass

parameters needed for discretisation and choice of numerical method

Source code in seaice3p/params/dimensional/numerical.py 
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
@serde(type_check=coerce)
@dataclass(frozen=True)
class NumericalParams:
    """parameters needed for discretisation and choice of numerical method"""

    I: int = 50
    regularisation: float = 1e-6
    solver_choice: str = "RK23"  # scipy.integrate.solve_IVP solver choice

    @property
    def step(self):
        return 1 / self.I

ocean_forcing

DimensionalBRW09OceanForcing dataclass

Ocean temperature provided by Barrow 2009 data at 2.4m and specify ocean fixed gas saturation state

Source code in seaice3p/params/dimensional/ocean_forcing.py 
56
57
58
59
60
61
62
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalBRW09OceanForcing:
    """Ocean temperature provided by Barrow 2009 data at 2.4m and specify ocean
    fixed gas saturation state"""

    pass
DimensionalFixedHeatFluxOceanForcing dataclass

Provides constant ocean heat flux at the bottom of the domain

Parameters:

Name Type Description Default
ocean_heat_flux float

The constant heat flux at the bottom of the domain in W/m2

 1
Source code in seaice3p/params/dimensional/ocean_forcing.py 
14
15
16
17
18
19
20
21
22
23
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalFixedHeatFluxOceanForcing:
    """Provides constant ocean heat flux at the bottom of the domain

    Args:
        ocean_heat_flux: The constant heat flux at the bottom of the domain in W/m2
    """

    ocean_heat_flux: float = 1
DimensionalFixedTempOceanForcing dataclass

Fixed temperature and gas saturation ocean boundary condition

Source code in seaice3p/params/dimensional/ocean_forcing.py 
 6
 7
 8
 9
10
11
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalFixedTempOceanForcing:
    """Fixed temperature and gas saturation ocean boundary condition"""

    ocean_temp: float = -1
DimensionalMonthlyHeatFluxOceanForcing dataclass

Provides constant ocean heat flux at the bottom of the domain in each month

Proivde an average monthly ocean heat flux with the entries i=0, 1, 2, 3, ...., 11 in the tuple corresponding to the months January, February, March, April, ...., December

Parameters:

Name Type Description Default
monthly_ocean_heat_flux Tuple[float, ...]

Tuple of ocean heat flux values in each month in W/m2

 (1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0)
Source code in seaice3p/params/dimensional/ocean_forcing.py 
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalMonthlyHeatFluxOceanForcing:
    """Provides constant ocean heat flux at the bottom of the domain in each month

    Proivde an average monthly ocean heat flux with the entries
    i=0, 1, 2, 3, ...., 11
    in the tuple corresponding to the months
    January, February, March, April, ...., December

    Args:
        monthly_ocean_heat_flux: Tuple of ocean heat flux values in each month in W/m2
    """

    monthly_ocean_heat_flux: Tuple[float, ...] = (
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
    )

water

CubicLiquidus dataclass

Cubic fit to liquidus to give liquidus salinity in terms of temperature

S = a0 + a1 T + a2 T^2 + a3 T^3

defaults are taken from Notz PhD thesis for fit to Assur seawater data

Source code in seaice3p/params/dimensional/water.py 
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
@serde(type_check=coerce)
@dataclass(frozen=True)
class CubicLiquidus:
    """Cubic fit to liquidus to give liquidus salinity in terms of temperature

    S = a0 + a1 T + a2 T^2 + a3 T^3

    defaults are taken from Notz PhD thesis for fit to Assur seawater data
    """

    eutectic_temperature: float = -21.1  # deg Celsius
    a0: float = -1.2
    a1: float = -21.8
    a2: float = -0.919
    a3: float = -0.0178

    def get_liquidus_salinity(self, temperature):
        return (
            self.a0
            + self.a1 * temperature
            + self.a2 * temperature**2
            + self.a3 * temperature**3
        )

    def get_liquidus_temperature(self, salinity):
        temperature = fsolve(
            lambda x: salinity - self.get_liquidus_salinity(x),
            np.full_like(salinity, -2),
        )
        if temperature.size == 1:
            return temperature[0]
        else:
            return temperature
DimensionalWaterParams dataclass
Source code in seaice3p/params/dimensional/water.py 
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
@serde(type_check=coerce)
@dataclass(frozen=True)
class DimensionalWaterParams:
    liquid_density: float = 1028  # kg/m3
    ice_density: float = 916  # kg/m3
    ocean_salinity: float = 34  # g/kg
    liquidus: LinearLiquidus | CubicLiquidus = LinearLiquidus()
    latent_heat: float = 334e3  # latent heat of fusion for ice in J/kg
    liquid_specific_heat_capacity: float = 4184  # J/kg degC
    solid_specific_heat_capacity: float = 2009  # J/kg degC
    liquid_thermal_conductivity: float = 0.54  # water thermal conductivity in W/m deg C
    solid_thermal_conductivity: float = 2.22  # ice thermal conductivity in W/m deg C
    snow_thermal_conductivity: float = 0.31  # snow thermal conductivity in W/m deg C
    snow_density: float = 150  # snow density kg/m3

    eddy_diffusivity: float = 0

    salt_diffusivity: float = 0  # molecular diffusivity of salt in water in m2/s
    # used to calculate Rayleigh number for convection and density contraction in liquid equation of state
    haline_contraction_coefficient: float = 7.5e-4  # 1/ppt

    # calculated from moreau et al 2014 value of kinematic viscosity for sewater 2.7e-6
    # dynamic liquid_viscosity = 2.7e-6 * liquid_density
    liquid_viscosity: float = 2.78e-3  # dynamic liquid viscosity in Pa.s

    @property
    def eutectic_salinity(self):
        if isinstance(self.liquidus, LinearLiquidus):
            return self.liquidus.eutectic_salinity
        if isinstance(self.liquidus, CubicLiquidus):
            return self.liquidus.get_liquidus_salinity(
                self.liquidus.eutectic_temperature
            )

        raise NotImplementedError

    @property
    def eutectic_temperature(self):
        if isinstance(self.liquidus, LinearLiquidus) or isinstance(
            self.liquidus, CubicLiquidus
        ):
            return self.liquidus.eutectic_temperature

        raise NotImplementedError

    @property
    def salinity_difference(self):
        r"""calculate difference between eutectic salinity and typical ocean salinity

        .. math:: \Delta S = S_E - S_i

        """
        return self.eutectic_salinity - self.ocean_salinity

    @property
    def ocean_freezing_temperature(self):
        """calculate salinity dependent freezing temperature using linear liquidus with
        ocean salinity

        .. math:: T_i = T_L(S_i) = T_E S_i / S_E

        or using a cubic fit for the liquidus curve

        """
        if isinstance(self.liquidus, LinearLiquidus):
            return (
                self.eutectic_temperature * self.ocean_salinity / self.eutectic_salinity
            )
        if isinstance(self.liquidus, CubicLiquidus):
            return self.liquidus.get_liquidus_temperature(self.ocean_salinity)

        raise NotImplementedError

    @property
    def temperature_difference(self):
        r"""calculate

        .. math:: \Delta T = T_i - T_E

        """
        return self.ocean_freezing_temperature - self.eutectic_temperature

    @property
    def concentration_ratio(self):
        r"""Calculate concentration ratio as

        .. math:: \mathcal{C} = S_i / \Delta S

        """
        return self.ocean_salinity / self.salinity_difference

    @property
    def stefan_number(self):
        r"""calculate Stefan number

        .. math:: \text{St} = L / c_p \Delta T

        """
        return self.latent_heat / (
            self.temperature_difference * self.liquid_specific_heat_capacity
        )

    @property
    def thermal_diffusivity(self):
        r"""Return thermal diffusivity in m2/s

        .. math:: \kappa = \frac{k}{\rho_l c_p}

        """
        return self.liquid_thermal_conductivity / (
            self.liquid_density * self.liquid_specific_heat_capacity
        )

    @property
    def conductivity_ratio(self):
        r"""Calculate the ratio of solid to liquid thermal conductivity

        .. math:: \lambda = \frac{k_s}{k_l}

        """
        return self.solid_thermal_conductivity / self.liquid_thermal_conductivity

    @property
    def specific_heat_ratio(self):
        r"""Calculate the ratio of solid to liquid specific heat capacities

        .. math:: \lambda = \frac{c_{p,s}}{c_{p,l}}

        """
        return self.solid_specific_heat_capacity / self.liquid_specific_heat_capacity

    @property
    def eddy_diffusivity_ratio(self):
        r"""Calculate the ratio of eddy diffusivity to thermal diffusivity in
        the liquid phase

        .. math:: \lambda = \frac{\kappa_\text{turbulent}}{\kappa_l}

        """
        return self.eddy_diffusivity / self.thermal_diffusivity

    @property
    def snow_conductivity_ratio(self):
        r"""Calculate the ratio of snow to liquid thermal conductivity

        .. math:: \lambda = \frac{k_{sn}}{k_l}

        """
        return self.snow_thermal_conductivity / self.liquid_thermal_conductivity

    @property
    def lewis_salt(self):
        r"""Calculate the lewis number for salt, return np.inf if there is no salt
        diffusion.

        .. math:: \text{Le}_S = \kappa / D_s

        """
        if self.salt_diffusivity == 0:
            return np.inf

        return self.thermal_diffusivity / self.salt_diffusivity
concentration_ratio property

Calculate concentration ratio as

.. math:: \mathcal{C} = S_i / \Delta S

conductivity_ratio property

Calculate the ratio of solid to liquid thermal conductivity

.. math:: \lambda = \frac{k_s}{k_l}

eddy_diffusivity_ratio property

Calculate the ratio of eddy diffusivity to thermal diffusivity in the liquid phase

.. math:: \lambda = \frac{\kappa_\text{turbulent}}{\kappa_l}

lewis_salt property

Calculate the lewis number for salt, return np.inf if there is no salt diffusion.

.. math:: \text{Le}_S = \kappa / D_s

ocean_freezing_temperature property

calculate salinity dependent freezing temperature using linear liquidus with ocean salinity

.. math:: T_i = T_L(S_i) = T_E S_i / S_E

or using a cubic fit for the liquidus curve

salinity_difference property

calculate difference between eutectic salinity and typical ocean salinity

.. math:: \Delta S = S_E - S_i

snow_conductivity_ratio property

Calculate the ratio of snow to liquid thermal conductivity

.. math:: \lambda = \frac{k_{sn}}{k_l}

specific_heat_ratio property

Calculate the ratio of solid to liquid specific heat capacities

.. math:: \lambda = \frac{c_{p,s}}{c_{p,l}}

stefan_number property

calculate Stefan number

.. math:: \text{St} = L / c_p \Delta T

temperature_difference property

calculate

.. math:: \Delta T = T_i - T_E

thermal_diffusivity property

Return thermal diffusivity in m2/s

.. math:: \kappa = \frac{k}{\rho_l c_p}

forcing

BRW09Forcing

Surface and ocean temperature data loaded from thermistor temperature record during the Barrow 2009 field study.

Source code in seaice3p/params/forcing.py 
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
@serde(type_check=coerce)
class BRW09Forcing:
    """Surface and ocean temperature data loaded from thermistor temperature record
    during the Barrow 2009 field study.
    """

    Barrow_top_temperature_data_choice: str = "air"

    def __post_init__(self):
        """populate class attributes with barrow dimensional air temperature
        and time in days (with missing values filtered out).

        Note the metadata explaining how to use the barrow temperature data is also
        in seaice3p/forcing_data. The indices corresponding to days and air temp are
        hard coded in as class variables.
        """
        DATA_INDICES = {
            "time": 0,
            "air": 8,
            "bottom_snow": 18,
            "top_ice": 19,
        }
        data = np.genfromtxt(
            Path(__file__).parent.parent / "forcing_data/BRW09.txt", delimiter="\t"
        )
        top_temp_index = DATA_INDICES[self.Barrow_top_temperature_data_choice]
        time_index = DATA_INDICES["time"]

        barrow_top_temp = data[:, top_temp_index]
        barrow_days = data[:, time_index] - data[0, time_index]
        barrow_top_temp, barrow_days = _filter_missing_values(
            barrow_top_temp, barrow_days
        )

        self.barrow_top_temp = barrow_top_temp
        self.barrow_days = barrow_days
__post_init__()

populate class attributes with barrow dimensional air temperature and time in days (with missing values filtered out).

Note the metadata explaining how to use the barrow temperature data is also in seaice3p/forcing_data. The indices corresponding to days and air temp are hard coded in as class variables.

Source code in seaice3p/params/forcing.py 
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
def __post_init__(self):
    """populate class attributes with barrow dimensional air temperature
    and time in days (with missing values filtered out).

    Note the metadata explaining how to use the barrow temperature data is also
    in seaice3p/forcing_data. The indices corresponding to days and air temp are
    hard coded in as class variables.
    """
    DATA_INDICES = {
        "time": 0,
        "air": 8,
        "bottom_snow": 18,
        "top_ice": 19,
    }
    data = np.genfromtxt(
        Path(__file__).parent.parent / "forcing_data/BRW09.txt", delimiter="\t"
    )
    top_temp_index = DATA_INDICES[self.Barrow_top_temperature_data_choice]
    time_index = DATA_INDICES["time"]

    barrow_top_temp = data[:, top_temp_index]
    barrow_days = data[:, time_index] - data[0, time_index]
    barrow_top_temp, barrow_days = _filter_missing_values(
        barrow_top_temp, barrow_days
    )

    self.barrow_top_temp = barrow_top_temp
    self.barrow_days = barrow_days

ConstantForcing dataclass

Constant temperature forcing

Source code in seaice3p/params/forcing.py 
37
38
39
40
41
42
@serde(type_check=coerce)
@dataclass(frozen=True)
class ConstantForcing:
    """Constant temperature forcing"""

    constant_top_temperature: float = -1.5

ERA5Forcing

Forcing parameters for simulation forced with atmospheric variables from reanalysis data in netCDF file located at data_path.

Never create this object directly but instead initialise from a dimensional simulation configuration as we must pass it the simulation timescale to correctly read the atmospheric variables from the netCDF file.

Source code in seaice3p/params/forcing.py 
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
@serde(type_check=coerce)
class ERA5Forcing:
    """Forcing parameters for simulation forced with atmospheric variables
    from reanalysis data in netCDF file located at data_path.

    Never create this object directly but instead initialise from a dimensional
    simulation configuration as we must pass it the simulation timescale to correctly
    read the atmospheric variables from the netCDF file.
    """

    data_path: Path
    start_date: str
    timescale_in_days: float
    forcing_data_file_keys: ERA5FileKeys = ERA5FileKeys()
    snow_density: Optional[float] = None
    SW_forcing: DimensionalSWForcing = DimensionalConstantSWForcing()
    LW_forcing: DimensionalLWForcing = DimensionalConstantLWForcing()
    turbulent_flux: DimensionalTurbulentFlux = DimensionalConstantTurbulentFlux()
    oil_heating: DimensionalOilHeating = DimensionalBackgroundOilHeating()

    def __post_init__(self):
        data = xr.open_dataset(self.data_path)
        DATES = getattr(data, self.forcing_data_file_keys.time).to_numpy()
        DIMLESS_TIMES = (1 / self.timescale_in_days) * np.array(
            [
                (date - np.datetime64(self.start_date)) / np.timedelta64(1, "D")
                for date in DATES
            ]
        )

        # convert to deg C
        T2M = (
            getattr(data, self.forcing_data_file_keys.temperature_at_2m_in_K).to_numpy()
            - 273.15
        )
        D2M = (
            getattr(data, self.forcing_data_file_keys.dewpoint_at_2m_in_K).to_numpy()
            - 273.15
        )

        LW = getattr(
            data, self.forcing_data_file_keys.longwave_radiation_in_W_m2
        ).to_numpy()
        SW = getattr(
            data, self.forcing_data_file_keys.shortwave_radiation_in_W_m2
        ).to_numpy()

        # convert to KPa
        ATM = (
            getattr(data, self.forcing_data_file_keys.surface_pressure_in_Pa).to_numpy()
            / 1e3
        )

        wind_key = self.forcing_data_file_keys.windspeed_at_2m_in_m_s
        if wind_key is None:
            WIND = np.full_like(DIMLESS_TIMES, self.turbulent_flux.windspeed)
        else:
            WIND = getattr(data, wind_key).to_numpy()

        # Calculate specific humidity in kg/kg from dewpoint temperature
        SPEC_HUM = _calculate_specific_humidity(ATM, D2M)

        snow_key = self.forcing_data_file_keys.snow_depth_in_m
        # if ERA5 standard short name for snow depth in m of water equivalent use snow
        # density to convert to m of snow
        if snow_key == "sd":
            if self.snow_density is None:
                raise ValueError("No snow density provided")
            SNOW_DEPTH = getattr(data, "sd").to_numpy() * (1000 / self.snow_density)

        # If snow key is another name assume snow depth is just in m of snow
        elif snow_key is not None:
            SNOW_DEPTH = getattr(data, snow_key).to_numpy()

        # If snow key is None assume no snow
        else:
            SNOW_DEPTH = np.zeros_like(DIMLESS_TIMES)

        # Provide functions to interpolate forcing data at non-dimensional times
        # during simulation
        self.get_2m_temp = partial(
            np.interp, xp=DIMLESS_TIMES, fp=T2M, left=np.nan, right=np.nan
        )
        self.get_LW = partial(
            np.interp, xp=DIMLESS_TIMES, fp=LW, left=np.nan, right=np.nan
        )
        self.get_SW = partial(
            np.interp, xp=DIMLESS_TIMES, fp=SW, left=np.nan, right=np.nan
        )
        self.get_ATM = partial(
            np.interp, xp=DIMLESS_TIMES, fp=ATM, left=np.nan, right=np.nan
        )
        self.get_spec_hum = partial(
            np.interp, xp=DIMLESS_TIMES, fp=SPEC_HUM, left=np.nan, right=np.nan
        )
        self.get_snow_depth = partial(
            np.interp, xp=DIMLESS_TIMES, fp=SNOW_DEPTH, left=np.nan, right=np.nan
        )
        self.get_windspeed = partial(
            np.interp, xp=DIMLESS_TIMES, fp=WIND, left=np.nan, right=np.nan
        )

RadForcing dataclass

Forcing parameters for radiative transfer simulation with oil drops

we have not implemented the non-dimensionalisation for these parameters yet and so we just pass the dimensional values directly to the simulation

Source code in seaice3p/params/forcing.py 
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
@serde(type_check=coerce)
@dataclass(frozen=True)
class RadForcing:
    """Forcing parameters for radiative transfer simulation with oil drops

    we have not implemented the non-dimensionalisation for these parameters yet
    and so we just pass the dimensional values directly to the simulation"""

    SW_forcing: DimensionalSWForcing = DimensionalConstantSWForcing()
    LW_forcing: DimensionalLWForcing = DimensionalConstantLWForcing()
    turbulent_flux: DimensionalTurbulentFlux = DimensionalConstantTurbulentFlux()
    oil_heating: DimensionalOilHeating = DimensionalBackgroundOilHeating()

RobinForcing dataclass

Dimensionless forcing parameters for Robin boundary condition

Source code in seaice3p/params/forcing.py 
221
222
223
224
225
226
227
@serde(type_check=coerce)
@dataclass(frozen=True)
class RobinForcing:
    """Dimensionless forcing parameters for Robin boundary condition"""

    biot: float = 12
    restoring_temperature: float = -1.3

YearlyForcing dataclass

Yearly sinusoidal temperature forcing

Source code in seaice3p/params/forcing.py 
45
46
47
48
49
50
51
52
@serde(type_check=coerce)
@dataclass(frozen=True)
class YearlyForcing:
    """Yearly sinusoidal temperature forcing"""

    offset: float = -1.0
    amplitude: float = 0.75
    period: float = 4.0

_calculate_specific_humidity(pressure, dewpoint)

Take ERA5 data and return specific humidity at 2m in kg/kg

Source code in seaice3p/params/forcing.py 
210
211
212
213
214
215
216
217
218
def _calculate_specific_humidity(pressure: NDArray, dewpoint: NDArray) -> NDArray:
    """Take ERA5 data and return specific humidity at 2m in kg/kg"""
    return (
        specific_humidity_from_dewpoint(
            pressure * metpyunits.kPa, dewpoint * metpyunits.degC
        )
        .to("kg/kg")
        .magnitude
    )

_filter_missing_values(air_temp, days)

Filter out missing values are recorded as 9999

Source code in seaice3p/params/forcing.py 
31
32
33
34
def _filter_missing_values(air_temp, days):
    """Filter out missing values are recorded as 9999"""
    is_missing = np.abs(air_temp) > 100
    return air_temp[~is_missing], days[~is_missing]

initial_conditions

OilInitialConditions dataclass

values for bottom (ocean) boundary

Source code in seaice3p/params/initial_conditions.py 
12
13
14
15
16
17
18
19
20
21
22
23
@serde(type_check=coerce)
@dataclass(frozen=True)
class OilInitialConditions:
    """values for bottom (ocean) boundary"""

    # Non dimensional parameters for summer initial conditions
    initial_ice_depth: float = 0.5
    initial_ocean_temperature: float = -0.05
    initial_ice_temperature: float = -0.1
    initial_oil_volume_fraction: float = 1e-7
    initial_ice_bulk_salinity: float = -0.1
    initial_oil_free_depth: float = 0

ocean_forcing

BRW09OceanForcing

Ocean temperature provided by Barrow 2009 data at 2.4m and specify ocean fixed gas saturation state

Source code in seaice3p/params/ocean_forcing.py 
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
@serde(type_check=coerce)
class BRW09OceanForcing:
    """Ocean temperature provided by Barrow 2009 data at 2.4m and specify ocean
    fixed gas saturation state"""

    ocean_gas_sat: float = 1.0

    def __post_init__(self):
        """populate class attributes with barrow dimensional ocean temperature
        and time in days (with missing values filtered out).

        Note the metadata explaining how to use the barrow temperature data is also
        in seaice3p/forcing_data.
        """
        data = np.genfromtxt(
            Path(__file__).parent.parent / "forcing_data/BRW09.txt", delimiter="\t"
        )
        ocean_temp_index = 43
        time_index = 0

        barrow_bottom_temp = data[:, ocean_temp_index]
        barrow_ocean_days = data[:, time_index] - data[0, time_index]
        barrow_bottom_temp, barrow_ocean_days = _filter_missing_values(
            barrow_bottom_temp, barrow_ocean_days
        )

        self.barrow_bottom_temp = barrow_bottom_temp
        self.barrow_ocean_days = barrow_ocean_days
__post_init__()

populate class attributes with barrow dimensional ocean temperature and time in days (with missing values filtered out).

Note the metadata explaining how to use the barrow temperature data is also in seaice3p/forcing_data.

Source code in seaice3p/params/ocean_forcing.py 
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
def __post_init__(self):
    """populate class attributes with barrow dimensional ocean temperature
    and time in days (with missing values filtered out).

    Note the metadata explaining how to use the barrow temperature data is also
    in seaice3p/forcing_data.
    """
    data = np.genfromtxt(
        Path(__file__).parent.parent / "forcing_data/BRW09.txt", delimiter="\t"
    )
    ocean_temp_index = 43
    time_index = 0

    barrow_bottom_temp = data[:, ocean_temp_index]
    barrow_ocean_days = data[:, time_index] - data[0, time_index]
    barrow_bottom_temp, barrow_ocean_days = _filter_missing_values(
        barrow_bottom_temp, barrow_ocean_days
    )

    self.barrow_bottom_temp = barrow_bottom_temp
    self.barrow_ocean_days = barrow_ocean_days

FixedHeatFluxOceanForcing dataclass

Provides constant dimensionless ocean heat flux at the bottom of the domain and fixed gas saturation state.

Source code in seaice3p/params/ocean_forcing.py 
27
28
29
30
31
32
33
34
@serde(type_check=coerce)
@dataclass(frozen=True)
class FixedHeatFluxOceanForcing:
    """Provides constant dimensionless ocean heat flux at the bottom of the domain and fixed gas
    saturation state."""

    ocean_heat_flux: float = 1
    ocean_gas_sat: float = 1.0

FixedTempOceanForcing dataclass

Fixed temperature and gas saturation ocean boundary condition

Source code in seaice3p/params/ocean_forcing.py 
18
19
20
21
22
23
24
@serde(type_check=coerce)
@dataclass(frozen=True)
class FixedTempOceanForcing:
    """Fixed temperature and gas saturation ocean boundary condition"""

    ocean_temp: float = 0.1
    ocean_gas_sat: float = 1.0

MonthlyHeatFluxOceanForcing

Provides constant dimensionless ocean heat flux at the bottom of the domain in each month

and ocean gas saturation state.

Proivde an average monthly ocean heat flux with the entries i=0, 1, 2, 3, ...., 11 in the tuple corresponding to the months January, February, March, April, ...., December

Parameters:

Name Type Description Default
monthly_ocean_heat_flux

Tuple of dimensionless ocean heat flux values in

 required
Source code in seaice3p/params/ocean_forcing.py 
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
@serde(type_check=coerce)
class MonthlyHeatFluxOceanForcing:
    """Provides constant dimensionless ocean heat flux at the bottom of the domain in
    each month

    and ocean gas saturation state.

    Proivde an average monthly ocean heat flux with the entries
    i=0, 1, 2, 3, ...., 11
    in the tuple corresponding to the months
    January, February, March, April, ...., December

    Args:
        monthly_ocean_heat_flux: Tuple of dimensionless ocean heat flux values in
        each month
    """

    start_date: str
    timescale_in_days: float
    monthly_ocean_heat_flux: Tuple[float, ...] = (
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
        1.0,
    )
    ocean_gas_sat: float = 1.0

    def __post_init__(self):
        # Provide functions to interpolate day of year to monthly ocean heat flux
        self._interpolate_heat_flux = partial(
            np.interp,
            xp=np.array([15, 46, 74, 105, 135, 166, 196, 227, 258, 288, 319, 349]),
            fp=np.array(self.monthly_ocean_heat_flux),
            period=365,
        )

    def get_ocean_heat_flux(self, simulation_time: float) -> float:
        start_datetime = datetime.strptime(self.start_date, "%Y-%m-%d")
        current_datetime = start_datetime + timedelta(
            days=self.timescale_in_days * simulation_time
        )
        current_day = (
            current_datetime - datetime(start_datetime.year, 1, 1)
        ).total_seconds() / 86400
        return self._interpolate_heat_flux(current_day)

params

Classes containing parameters required to run a simulation

The config class contains all the parameters needed to run a simulation as well as methods to save and load this configuration to a yaml file.

Config dataclass

contains all information needed to run a simulation and save output

this config object can be saved and loaded to a yaml file.

Source code in seaice3p/params/params.py 
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
@serde(type_check=coerce)
@dataclass(frozen=True)
class Config:
    """contains all information needed to run a simulation and save output

    this config object can be saved and loaded to a yaml file."""

    name: str
    total_time: float
    savefreq: float

    physical_params: PhysicalParams
    bubble_params: BubbleParams
    brine_convection_params: BrineConvectionParams
    forcing_config: ForcingConfig
    ocean_forcing_config: OceanForcingConfig
    initial_conditions_config: InitialConditionsConfig
    numerical_params: NumericalParams = NumericalParams()
    scales: Scales | None = None

    def save(self, directory: Path):
        with open(directory / f"{self.name}.yml", "w") as outfile:
            outfile.write(to_yaml(self))

    @classmethod
    def load(cls, path):
        with open(path, "r") as infile:
            yaml = infile.read()
        return from_yaml(cls, yaml)

get_config(dimensional_params)

Return a Config object for the simulation.

physical parameters and Darcy law parameters are calculated from the dimensional input. You can modify the numerical parameters and boundary conditions and forcing provided for the simulation.

Source code in seaice3p/params/params.py 
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
def get_config(dimensional_params: DimensionalParams) -> Config:
    """Return a Config object for the simulation.

    physical parameters and Darcy law parameters are calculated from the dimensional
    input. You can modify the numerical parameters and boundary conditions and
    forcing provided for the simulation."""
    physical_params = get_dimensionless_physical_params(dimensional_params)
    initial_conditions_config = get_dimensionless_initial_conditions_config(
        dimensional_params
    )
    brine_convection_params = get_dimensionless_brine_convection_params(
        dimensional_params
    )
    bubble_params = get_dimensionless_bubble_params(dimensional_params)
    forcing_config = get_dimensionless_forcing_config(dimensional_params)
    ocean_forcing_config = get_dimensionless_ocean_forcing_config(dimensional_params)
    return Config(
        name=dimensional_params.name,
        physical_params=physical_params,
        initial_conditions_config=initial_conditions_config,
        brine_convection_params=brine_convection_params,
        bubble_params=bubble_params,
        forcing_config=forcing_config,
        ocean_forcing_config=ocean_forcing_config,
        numerical_params=dimensional_params.numerical_params,
        scales=dimensional_params.scales,
        total_time=dimensional_params.total_time,
        savefreq=dimensional_params.savefreq,
    )

physical

BasePhysicalParams dataclass

Not to be used directly but provides the common parameters for physical params objects

Source code in seaice3p/params/physical.py 
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
@serde(type_check=coerce)
@dataclass(frozen=True)
class BasePhysicalParams:
    """Not to be used directly but provides the common parameters for physical params
    objects
    """

    expansion_coefficient: float = 0.029
    concentration_ratio: float = 0.17
    stefan_number: float = 4.2
    lewis_salt: float = np.inf
    lewis_gas: float = np.inf
    frame_velocity: float = 0

    specific_heat_ratio: float = 0.5
    conductivity_ratio: float = 4.11
    eddy_diffusivity_ratio: float = 0
    snow_conductivity_ratio: float = 0.574

    # Option to change tolerable supersaturation
    tolerable_super_saturation_fraction: float = 1

    gas_viscosity_ratio: float = 0
    gas_bubble_eddy_diffusion: bool = False

    get_liquidus_temperature: Optional[Callable] = None
    get_liquidus_salinity: Optional[Callable] = None

DISEQPhysicalParams dataclass

Bases: BasePhysicalParams

non dimensional numbers for the mushy layer

Source code in seaice3p/params/physical.py 
50
51
52
53
54
55
56
@serde(type_check=coerce)
@dataclass(frozen=True)
class DISEQPhysicalParams(BasePhysicalParams):
    """non dimensional numbers for the mushy layer"""

    # only used in DISEQ model
    damkohler_number: float = 1

EQMPhysicalParams dataclass

Bases: BasePhysicalParams

non dimensional numbers for the mushy layer

Source code in seaice3p/params/physical.py 
44
45
46
47
@serde(type_check=coerce)
@dataclass(frozen=True)
class EQMPhysicalParams(BasePhysicalParams):
    """non dimensional numbers for the mushy layer"""

plot

script to visualise simulation data

usage: python -m seaice3p.plot "glob pattern to find npz files" Optional[True/False]

assumes the simulation configurations are to be found in the same directory as the data. If the simulation is a non-dimensional configuration file add the False option after the glob pattern.

run_simulation

Module to run the simulation on the given configuration with the appropriate solver.

Solve reduced model using scipy solve_ivp using RK23 solver.

Impose a maximum timestep constraint using courant number for thermal diffusion as this is an explicit method.

This solver uses adaptive timestepping which makes it a good choice for running simulations with large buoyancy driven gas bubble velocities and we save the output at intervals given by the savefreq parameter in configuration.

run_batch(list_of_cfg, directory, verbosity_level=0)

Run a batch of simulations from a list of configurations.

Each simulation name is logged, as well as if it successfully runs or crashes. Output from each simulation is saved in a .npz file.

:param list_of_cfg: list of configurations :type list_of_cfg: List[seaice3p.params.Config]

Source code in seaice3p/run_simulation.py 
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
def run_batch(list_of_cfg: List[Config], directory: Path, verbosity_level=0) -> None:
    """Run a batch of simulations from a list of configurations.

    Each simulation name is logged, as well as if it successfully runs or crashes.
    Output from each simulation is saved in a .npz file.

    :param list_of_cfg: list of configurations
    :type list_of_cfg: List[seaice3p.params.Config]

    """
    optprint = get_printer(verbosity_level, verbosity_threshold=1)
    for cfg in list_of_cfg:
        optprint(f"seaice3pv{__version__}: {cfg.name}")
        try:
            solve(cfg, directory, verbosity_level=verbosity_level)
        except Exception as e:
            optprint(f"{cfg.name} crashed")
            optprint(f"{e}")

state

disequilibrium_state

DISEQState dataclass

Contains the principal variables for solution with non-equilibrium gas phase. The total bulk gas is partitioned between dissolved gas and free phase gas with a finite nucleation rate (non dimensional damkohler number).

principal solution components: bulk enthalpy bulk salinity bulk dissolved gas gas fraction

all on the center grid.

Note: Define bulk dissolved gas for the system as

expansion_coefficient * liquid_fraction * dissolved_gas

so that this is different from the dissolved gas concentration and

bulk_gas = bulk_dissolved_gas + gas_fraction

in non-dimensional units.

Source code in seaice3p/state/disequilibrium_state.py 
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
@dataclass(frozen=True)
class DISEQState:
    """Contains the principal variables for solution with non-equilibrium gas phase.
    The total bulk gas is partitioned between dissolved gas and free phase gas with
    a finite nucleation rate (non dimensional damkohler number).

    principal solution components:
    bulk enthalpy
    bulk salinity
    bulk dissolved gas
    gas fraction

    all on the center grid.

    Note:
    Define bulk dissolved gas for the system as

    expansion_coefficient * liquid_fraction * dissolved_gas

    so that this is different from the dissolved gas concentration and

    bulk_gas = bulk_dissolved_gas + gas_fraction

    in non-dimensional units.
    """

    time: float
    enthalpy: NDArray
    salt: NDArray
    bulk_dissolved_gas: NDArray
    gas_fraction: NDArray

    @cached_property
    def gas(self) -> NDArray:
        """Calculate bulk gas content and use same attribute name as EQMState"""
        return self.bulk_dissolved_gas + self.gas_fraction
gas cached property

Calculate bulk gas content and use same attribute name as EQMState

DISEQStateBCs dataclass

Stores information needed for solution at one timestep with BCs on ghost cells as well

Initialiase the prime variables for the solver: enthalpy, bulk salinity and bulk air

Source code in seaice3p/state/disequilibrium_state.py 
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
@dataclass(frozen=True)
class DISEQStateBCs:
    """Stores information needed for solution at one timestep with BCs on ghost
    cells as well

    Initialiase the prime variables for the solver:
    enthalpy, bulk salinity and bulk air
    """

    time: float
    enthalpy: NDArray
    salt: NDArray

    temperature: NDArray
    liquid_salinity: NDArray
    dissolved_gas: NDArray
    liquid_fraction: NDArray
    bulk_dissolved_gas: NDArray
    gas_fraction: NDArray

DISEQStateFull dataclass

Contains all variables variables for solution with non-equilibrium gas phase after running the enthalpy method on DISEQSate. The total bulk gas is partitioned between dissolved gas and free phase gas with a finite nucleation rate (non dimensional damkohler number).

principal solution components: bulk enthalpy bulk salinity bulk dissolved gas gas fraction

enthalpy method variables: temperature liquid_fraction solid_fraction liquid_salinity dissolved_gas

all on the center grid.

Note: Define bulk dissolved gas for the system as

expansion_coefficient * liquid_fraction * dissolved_gas

so that this is different from the dissolved gas concentration and

bulk_gas = bulk_dissolved_gas + gas_fraction

in non-dimensional units.

Source code in seaice3p/state/disequilibrium_state.py 
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
@dataclass(frozen=True)
class DISEQStateFull:
    """Contains all variables variables for solution with non-equilibrium gas phase
    after running the enthalpy method on DISEQSate.
    The total bulk gas is partitioned between dissolved gas and free phase gas with
    a finite nucleation rate (non dimensional damkohler number).

    principal solution components:
    bulk enthalpy
    bulk salinity
    bulk dissolved gas
    gas fraction

    enthalpy method variables:
    temperature
    liquid_fraction
    solid_fraction
    liquid_salinity
    dissolved_gas

    all on the center grid.

    Note:
    Define bulk dissolved gas for the system as

    expansion_coefficient * liquid_fraction * dissolved_gas

    so that this is different from the dissolved gas concentration and

    bulk_gas = bulk_dissolved_gas + gas_fraction

    in non-dimensional units.
    """

    time: float
    enthalpy: NDArray
    salt: NDArray
    bulk_dissolved_gas: NDArray
    gas_fraction: NDArray

    temperature: NDArray
    liquid_fraction: NDArray
    solid_fraction: NDArray
    liquid_salinity: NDArray
    dissolved_gas: NDArray

    @cached_property
    def gas(self) -> NDArray:
        """Calculate bulk gas content and use same attribute name as EQMState"""
        return self.bulk_dissolved_gas + self.gas_fraction
gas cached property

Calculate bulk gas content and use same attribute name as EQMState

equilibrium_state

EQMState dataclass

Contains the principal variables for solution with equilibrium gas phase:

bulk enthalpy bulk salinity bulk gas

all on the center grid.

Source code in seaice3p/state/equilibrium_state.py 
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
@dataclass(frozen=True)
class EQMState:
    """Contains the principal variables for solution with equilibrium gas phase:

    bulk enthalpy
    bulk salinity
    bulk gas

    all on the center grid.
    """

    time: float
    enthalpy: NDArray
    salt: NDArray
    gas: NDArray

EQMStateBCs dataclass

Stores information needed for solution at one timestep with BCs on ghost cells as well

Initialiase the prime variables for the solver: enthalpy, bulk salinity and bulk air

Source code in seaice3p/state/equilibrium_state.py 
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
@dataclass(frozen=True)
class EQMStateBCs:
    """Stores information needed for solution at one timestep with BCs on ghost
    cells as well

    Initialiase the prime variables for the solver:
    enthalpy, bulk salinity and bulk air
    """

    time: float
    enthalpy: NDArray
    salt: NDArray
    gas: NDArray

    temperature: NDArray
    liquid_salinity: NDArray
    dissolved_gas: NDArray
    gas_fraction: NDArray
    liquid_fraction: NDArray

EQMStateFull dataclass

Contains all variables variables for solution with equilibrium gas phase after running the enthalpy method on EQMSate.

principal solution components: bulk enthalpy bulk salinity bulk gas

enthalpy method variables: temperature liquid_fraction solid_fraction liquid_salinity dissolved_gas gas_fraction

all on the center grid.

Source code in seaice3p/state/equilibrium_state.py 
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
@dataclass(frozen=True)
class EQMStateFull:
    """Contains all variables variables for solution with equilibrium gas phase
    after running the enthalpy method on EQMSate.

    principal solution components:
    bulk enthalpy
    bulk salinity
    bulk gas

    enthalpy method variables:
    temperature
    liquid_fraction
    solid_fraction
    liquid_salinity
    dissolved_gas
    gas_fraction

    all on the center grid.
    """

    time: float
    enthalpy: NDArray
    salt: NDArray
    gas: NDArray

    temperature: NDArray
    liquid_fraction: NDArray
    solid_fraction: NDArray
    liquid_salinity: NDArray
    dissolved_gas: NDArray
    gas_fraction: NDArray