Skip to content

API Reference

Simulate steadily solidifying three-phase mushy layers in a Hele-Shaw cell.

boundary_conditions

get_boundary_conditions(non_dimensional_params, bottom_variables, top_variables)

Function to return the boundary conditions for scipy solve_BVP.

The returned array is zero when the boundary conditions are satisfied.

Parameters:

Name Type Description Default
non_dimensional_params NonDimensionalParams

Non-dimensional parameters

 required
bottom_variables NDArray

Array of the solution variables evaluated at the bottom boundary.

 required
top_variables NDArray

Array of the solution variables evaluated at the top boundary.

 required

Returns:

Name Type Description
NDArray NDArray

residual of the boundary conditions
Source code in mush3p/boundary_conditions.py 
 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
def get_boundary_conditions(
    non_dimensional_params: NonDimensionalParams,
    bottom_variables: NDArray,
    top_variables: NDArray,
) -> NDArray:
    """Function to return the boundary conditions for scipy solve_BVP.

    The returned array is zero when the boundary conditions are satisfied.

    Args:
        non_dimensional_params (NonDimensionalParams): Non-dimensional parameters
        bottom_variables (NDArray): Array of the solution variables evaluated at the bottom boundary.
        top_variables (NDArray): Array of the solution variables evaluated at the top boundary.

    Returns:
        NDArray: residual of the boundary conditions
    """

    OPTIONS = {
        "full": BoundaryConditionsFull,
        "incompressible": BoundaryConditionsIncompressible,
        "reduced": BoundaryConditionsReduced,
        "instant": BoundaryConditionsInstant,
    }
    return OPTIONS[non_dimensional_params.model_choice](
        non_dimensional_params, bottom_variables, top_variables
    ).boundary_conditions

model

full

FullModel

Implement the full model equations and provide the ode_fun method to be used by scipy solve_BVP.

The gas fraction is calculated by solving a non-linear equation and is fully coupled to the liquid darcy flux. Solid, liquid and gas fractions sum to one. Gas density is calculated from the ideal gas law.

Parameters:

Name Type Description Default
params NonDimensionalParams

Non-dimensional parameters

 required
height NDArray

Height values

 required
temperature NDArray

Temperature values

 required
temperature_derivative NDArray

Temperature derivative values

 required
dissolved_gas_concentration NDArray

Dissolved gas concentration values

 required
hydrostatic_pressure NDArray

Hydrostatic pressure values

 required
frozen_gas_fraction NDArray

Frozen gas fraction

 required
mushy_layer_depth NDArray

Mushy layer depth

 required
Source code in mush3p/model/full.py 
 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
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
class FullModel:
    """Implement the full model equations and provide the ode_fun method to be used by
    scipy solve_BVP.

    The gas fraction is calculated by solving a non-linear equation and is fully coupled
    to the liquid darcy flux.
    Solid, liquid and gas fractions sum to one.
    Gas density is calculated from the ideal gas law.

    Args:
        params (NonDimensionalParams): Non-dimensional parameters
        height (NDArray): Height values
        temperature (NDArray): Temperature values
        temperature_derivative (NDArray): Temperature derivative values
        dissolved_gas_concentration (NDArray): Dissolved gas concentration values
        hydrostatic_pressure (NDArray): Hydrostatic pressure values
        frozen_gas_fraction (NDArray): Frozen gas fraction
        mushy_layer_depth (NDArray): Mushy layer depth"""

    def __init__(
        self,
        params: NonDimensionalParams,
        height: NDArray,
        temperature: NDArray,
        temperature_derivative: NDArray,
        dissolved_gas_concentration: NDArray,
        hydrostatic_pressure: NDArray,
        frozen_gas_fraction: NDArray,
        mushy_layer_depth: NDArray,
    ) -> None:
        self.params = params
        self.height = height
        self.temperature = temperature
        self.temperature_derivative = temperature_derivative
        self.dissolved_gas_concentration = dissolved_gas_concentration
        self.hydrostatic_pressure = hydrostatic_pressure
        self.frozen_gas_fraction = frozen_gas_fraction
        self.mushy_layer_depth = mushy_layer_depth

    @property
    def solid_salinity(self) -> NDArray:
        return np.full_like(self.temperature, -self.params.concentration_ratio)

    @property
    def liquid_salinity(self) -> NDArray:
        return -self.temperature

    @property
    def solid_fraction(self) -> NDArray:
        concentration_ratio = self.params.concentration_ratio
        return (
            -(1 - self.frozen_gas_fraction)
            * self.temperature
            / (concentration_ratio - self.temperature)
        )

    @property
    def liquid_fraction(self) -> NDArray:
        return calculate_liquid_fraction(self.gas_fraction, self.solid_fraction)

    @property
    def gas_darcy_velocity(
        self,
    ) -> NDArray:
        return calculate_gas_darcy_velocity(
            self.solid_fraction,
            self.gas_fraction,
            self.params,
        )

    @property
    def gas_density(
        self,
    ):
        return calculate_gas_density(
            self.height,
            self.mushy_layer_depth,
            self.temperature,
            self.hydrostatic_pressure,
            self.params,
        )

    @property
    def gas_fraction(
        self,
    ) -> Any:
        return calculate_gas_fraction(
            self.frozen_gas_fraction,
            self.solid_fraction,
            self.temperature,
            self.dissolved_gas_concentration,
            self.gas_density,
            self.params,
        )

    @property
    def liquid_darcy_velocity(self):
        return calculate_liquid_darcy_velocity(
            self.gas_fraction, self.frozen_gas_fraction
        )

    @property
    def permeability(self) -> NDArray:
        liquid_permeability_reciprocal = (
            1 - self.liquid_fraction
        ) ** 2 / self.liquid_fraction**3
        reference = self.params.hele_shaw_permeability_scaled
        return ((1 / reference) + liquid_permeability_reciprocal) ** (-1)

    @property
    def saturation_concentration(self) -> NDArray:
        return np.full_like(self.temperature, 1)

    @property
    def nucleation_rate(self) -> NDArray:
        indicator = np.where(
            self.dissolved_gas_concentration >= self.saturation_concentration, 1, 0
        )

        return (
            self.liquid_fraction
            * indicator
            * (self.dissolved_gas_concentration - self.saturation_concentration)
        )

    @property
    def solid_fraction_derivative(
        self,
    ) -> NDArray:
        concentration_ratio = self.params.concentration_ratio
        return (
            -concentration_ratio
            * (1 - self.frozen_gas_fraction)
            * self.temperature_derivative
            / (concentration_ratio - self.temperature) ** 2
        )

    @property
    def gas_fraction_derivative(self) -> NDArray:
        """Numerically approximate the derivative with finite difference."""
        return np.gradient(self.gas_fraction, self.height)

    @property
    def hydrostatic_pressure_derivative(
        self,
    ) -> NDArray:
        return -self.mushy_layer_depth * self.liquid_darcy_velocity / self.permeability

    @property
    def effective_heat_capacity(self):
        solid_specific_heat_capacity_ratio = (
            self.params.solid_specific_heat_capacity_ratio
        )
        gas_specific_heat_capacity_ratio = self.params.gas_specific_heat_capacity_ratio
        density_ratio = self.params.gas_density_ratio
        return (
            1
            - (1 - solid_specific_heat_capacity_ratio) * self.solid_fraction
            - (1 - gas_specific_heat_capacity_ratio * density_ratio * self.gas_density)
            * self.gas_fraction
        )

    @property
    def effective_thermal_conductivity(self):
        gas_conductivity_ratio = self.params.gas_conductivity_ratio
        solid_conductivity_ratio = self.params.solid_conductivity_ratio
        return (
            1
            - (1 - solid_conductivity_ratio) * self.solid_fraction
            - (1 - gas_conductivity_ratio) * self.gas_fraction
        )

    @property
    def temperature_second_derivative(
        self,
    ) -> NDArray:
        stefan_number = self.params.stefan_number
        gas_specific_heat_capacity_ratio = self.params.gas_specific_heat_capacity_ratio
        density_ratio = self.params.gas_density_ratio
        gas_conductivity_ratio = self.params.gas_conductivity_ratio
        solid_conductivity_ratio = self.params.solid_conductivity_ratio

        heat_capacity_term = (
            self.mushy_layer_depth
            * self.effective_heat_capacity
            * self.temperature_derivative
        )
        liquid_advection_term = (
            self.mushy_layer_depth
            * self.liquid_darcy_velocity
            * self.temperature_derivative
        )
        gas_advection_term = (
            self.mushy_layer_depth
            * density_ratio
            * self.gas_density
            * gas_specific_heat_capacity_ratio
            * self.gas_darcy_velocity
            * self.temperature_derivative
        )
        latent_heat_term = (
            -self.mushy_layer_depth * stefan_number * self.solid_fraction_derivative
        )
        conductivity_change_term = (
            (1 - solid_conductivity_ratio) * self.solid_fraction_derivative
            + (1 - gas_conductivity_ratio) * self.gas_fraction_derivative
        ) * self.temperature_derivative

        return (1 / self.effective_thermal_conductivity) * (
            heat_capacity_term
            + liquid_advection_term
            + gas_advection_term
            + latent_heat_term
            + conductivity_change_term
        )

    @property
    def dissolved_gas_concentration_derivative(
        self,
    ) -> NDArray:

        damkholer_number = self.params.damkholer_number
        freezing = self.dissolved_gas_concentration * self.solid_fraction_derivative
        dissolution = -damkholer_number * self.mushy_layer_depth * self.nucleation_rate

        return (freezing + dissolution) / (
            1 - self.frozen_gas_fraction - self.solid_fraction
        )

    @property
    def check_volume_fractions_sum_to_one(self):
        if (
            np.max(
                np.abs(
                    self.solid_fraction + self.liquid_fraction + self.gas_fraction - 1
                )
            )
            > VOLUME_SUM_TOLERANCE
        ):
            return False
        return True

    @property
    def ode_fun(self):

        if not self.check_volume_fractions_sum_to_one:
            raise ValueError("Volume fractions do not sum to 1")

        return np.vstack(
            (
                self.temperature_derivative,
                self.temperature_second_derivative,
                self.dissolved_gas_concentration_derivative,
                self.hydrostatic_pressure_derivative,
                np.zeros_like(self.temperature),
                np.zeros_like(self.temperature),
            )
        )
gas_fraction_derivative: NDArray property

Numerically approximate the derivative with finite difference.

incompressible

IncompressibleModel

Bases: FullModel

implement equations for the incompressible model.

These are identical to the full model but with the dimensionless gas density set to 1.0.

Source code in mush3p/model/incompressible.py 
 5
 6
 7
 8
 9
10
11
12
13
14
15
class IncompressibleModel(FullModel):
    """implement equations for the incompressible model.

    These are identical to the full model but with the dimensionless gas density set to 1.0.
    """

    @property
    def gas_density(
        self,
    ):
        return np.ones_like(self.temperature)

instant

The equations for solving the instant nucleation model

All quantities are calculated from the smaller set of variables: temperature temperature_derivative hydrostatic_pressure frozen_gas_fraction mushy_layer_depth

height (vertical coordinate)

InstantNucleationModel

Bases: ReducedModel

Implements the equations to solve the instant model.

This is an extension of the reduced model where the nucleation rate is infinite and and so any supersaturations is immediately converted to the gas phase. This means we do not need to solve the ODE for dissolved gas concentration in this case.

Parameters:

Name Type Description Default
params NonDimensionalParams

Non-dimensional parameters

 required
height NDArray

Height values

 required
temperature NDArray

Temperature values

 required
temperature_derivative NDArray

Temperature derivative values

 required
hydrostatic_pressure NDArray

Hydrostatic pressure values

 required
frozen_gas_fraction NDArray

Frozen gas fraction

 required
mushy_layer_depth NDArray

Mushy layer depth

 required
Source code in mush3p/model/instant.py 
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
class InstantNucleationModel(ReducedModel):
    """Implements the equations to solve the instant model.

    This is an extension of the reduced model where the nucleation rate is infinite and
    and so any supersaturations is immediately converted to the gas phase.
    This means we do not need to solve the ODE for dissolved gas concentration in this
    case.

    Args:
        params (NonDimensionalParams): Non-dimensional parameters
        height (NDArray): Height values
        temperature (NDArray): Temperature values
        temperature_derivative (NDArray): Temperature derivative values
        hydrostatic_pressure (NDArray): Hydrostatic pressure values
        frozen_gas_fraction (NDArray): Frozen gas fraction
        mushy_layer_depth (NDArray): Mushy layer depth"""

    def __init__(
        self,
        params: NonDimensionalParams,
        height: NDArray,
        temperature: NDArray,
        temperature_derivative: NDArray,
        hydrostatic_pressure: NDArray,
        frozen_gas_fraction: NDArray,
        mushy_layer_depth: NDArray,
    ) -> None:
        self.params = params
        self.height = height
        self.temperature = temperature
        self.temperature_derivative = temperature_derivative
        self.hydrostatic_pressure = hydrostatic_pressure
        self.frozen_gas_fraction = frozen_gas_fraction
        self.mushy_layer_depth = mushy_layer_depth

    @property
    def dissolved_gas_concentration(self):
        return np.minimum(
            np.ones_like(self.liquid_fraction),
            self.params.far_dissolved_concentration_scaled / self.liquid_fraction,
        )

    @property
    def ode_fun(self) -> Any:

        if not self.check_volume_fractions_sum_to_one:
            raise ValueError("Volume fractions do not sum to 1")

        return np.vstack(
            (
                self.temperature_derivative,
                self.temperature_second_derivative,
                self.hydrostatic_pressure_derivative,
                np.zeros_like(self.temperature),
                np.zeros_like(self.temperature),
            )
        )

reduced

ReducedModel

Bases: IncompressibleModel

Implement the equations for the reduced model

The reduced model is an approximation to the full model where the exsolved gas volume fraction is small and incompressible. Under this approximation scheme the solid and liquid volume fractions sum to 1 and the exsolution of gas does not drive a liquid flow.

Source code in mush3p/model/reduced.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
 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
class ReducedModel(IncompressibleModel):
    """Implement the equations for the reduced model

    The reduced model is an approximation to the full model where the exsolved gas
    volume fraction is small and incompressible.
    Under this approximation scheme the solid and liquid volume fractions sum to 1
    and the exsolution of gas does not drive a liquid flow.
    """

    @property
    def liquid_darcy_velocity(self) -> NDArray:
        return np.zeros_like(self.temperature)

    @property
    def solid_fraction(self) -> NDArray:
        concentration_ratio = self.params.concentration_ratio
        return self.temperature / (self.temperature - concentration_ratio)

    @property
    def liquid_fraction(self) -> NDArray:
        return 1 - self.solid_fraction

    @property
    def gas_fraction(
        self,
    ) -> Any:
        return calculate_gas_fraction(
            0,
            self.solid_fraction,
            self.temperature,
            self.dissolved_gas_concentration,
            1,
            self.params,
        )

    @property
    def solid_fraction_derivative(
        self,
    ) -> NDArray:
        concentration_ratio = self.params.concentration_ratio
        return (
            -concentration_ratio
            * self.temperature_derivative
            / ((self.temperature - concentration_ratio) ** 2)
        )

    @property
    def hydrostatic_pressure_derivative(
        self,
    ) -> NDArray:
        return np.zeros_like(self.temperature)

    @property
    def effective_heat_capacity(self):
        solid_specific_heat_capacity_ratio = (
            self.params.solid_specific_heat_capacity_ratio
        )
        return 1 - (1 - solid_specific_heat_capacity_ratio) * self.solid_fraction

    @property
    def effective_thermal_conductivity(self):
        solid_conductivity_ratio = self.params.solid_conductivity_ratio
        return 1 - (1 - solid_conductivity_ratio) * self.solid_fraction

    @property
    def temperature_second_derivative(
        self,
    ) -> NDArray:
        stefan_number = self.params.stefan_number
        solid_conductivity_ratio = self.params.solid_conductivity_ratio

        heat_capacity_term = (
            self.mushy_layer_depth
            * self.effective_heat_capacity
            * self.temperature_derivative
        )
        latent_heat_term = (
            -self.mushy_layer_depth * stefan_number * self.solid_fraction_derivative
        )
        conductivity_change_term = (
            (1 - solid_conductivity_ratio) * self.solid_fraction_derivative
        ) * self.temperature_derivative

        return (1 / self.effective_thermal_conductivity) * (
            heat_capacity_term + latent_heat_term + conductivity_change_term
        )

    @property
    def dissolved_gas_concentration_derivative(
        self,
    ) -> NDArray:

        damkholer_number = self.params.damkholer_number
        dissolution = -damkholer_number * self.mushy_layer_depth * self.nucleation_rate

        return (1 / self.liquid_fraction) * (
            self.dissolved_gas_concentration * self.solid_fraction_derivative
            + dissolution
        )

    @property
    def check_volume_fractions_sum_to_one(self):
        if (
            np.max(np.abs(self.solid_fraction + self.liquid_fraction - 1))
            > VOLUME_SUM_TOLERANCE
        ):
            return False
        return True

output

NonDimensionalResults

Class to store non-dimensional solution profiles and provide methods to linearly interpolate them with depth.

The class also provides a save method to serialize the non-dimensional parameters and solution arrays to a JSON file, and a load method to deserialize them back into a NonDimensionalResults object.

The class calcualtes all mushy layer arrays using the model specified in the non-dimensional parameters.

Parameters:

Name Type Description Default
non_dimensional_parameters NonDimensionalParams

Non-dimensional parameters

 required
temperature_array ndarray

Temperature profile

 required
temperature_derivative_array ndarray

Temperature derivative profile

 required
concentration_array ndarray

Concentration profile

 required
hydrostatic_pressure_array ndarray

Hydrostatic pressure profile

 required
frozen_gas_fraction float

Frozen gas fraction

 required
mushy_layer_depth float

Mushy layer depth

 required
height_array ndarray

Vertical grid points.

 required

Attributes:

Name Type Description
name str

Name of the simulation
params NonDimensionalParams

Non-dimensional parameters
temperature_array ndarray

Temperature profile
temperature_derivative_array ndarray

Temperature derivative profile
concentration_array ndarray

Concentration profile
hydrostatic_pressure_array ndarray

Hydrostatic pressure profile
frozen_gas_fraction float

Frozen gas fraction
mushy_layer_depth float

Mushy layer depth
height_array ndarray

Vertical grid points
solid_salinity_array ndarray

Solid salinity profile
liquid_salinity_array ndarray

Liquid salinity profile
solid_fraction_array ndarray

Solid fraction profile
liquid_fraction_array ndarray

Liquid fraction profile
gas_fraction_array ndarray

Gas fraction profile
gas_density_array ndarray

Gas density profile
liquid_darcy_velocity_array ndarray

Liquid Darcy velocity profile
gas_darcy_velocity_array ndarray

Gas Darcy velocity profile
Source code in mush3p/output.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
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
class NonDimensionalResults:
    """Class to store non-dimensional solution profiles and provide methods to linearly
    interpolate them with depth.

    The class also provides a save method to serialize the non-dimensional parameters
    and solution arrays to a JSON file, and a load method to deserialize them back into
    a NonDimensionalResults object.

    The class calcualtes all mushy layer arrays using the model specified in the
    non-dimensional parameters.

    Args:
        non_dimensional_parameters (NonDimensionalParams): Non-dimensional parameters
        temperature_array (np.ndarray): Temperature profile
        temperature_derivative_array (np.ndarray): Temperature derivative profile
        concentration_array (np.ndarray): Concentration profile
        hydrostatic_pressure_array (np.ndarray): Hydrostatic pressure profile
        frozen_gas_fraction (float): Frozen gas fraction
        mushy_layer_depth (float): Mushy layer depth
        height_array (np.ndarray): Vertical grid points.

    Attributes:
        name (str): Name of the simulation
        params (NonDimensionalParams): Non-dimensional parameters
        temperature_array (np.ndarray): Temperature profile
        temperature_derivative_array (np.ndarray): Temperature derivative profile
        concentration_array (np.ndarray): Concentration profile
        hydrostatic_pressure_array (np.ndarray): Hydrostatic pressure profile
        frozen_gas_fraction (float): Frozen gas fraction
        mushy_layer_depth (float): Mushy layer depth
        height_array (np.ndarray): Vertical grid points
        solid_salinity_array (np.ndarray): Solid salinity profile
        liquid_salinity_array (np.ndarray): Liquid salinity profile
        solid_fraction_array (np.ndarray): Solid fraction profile
        liquid_fraction_array (np.ndarray): Liquid fraction profile
        gas_fraction_array (np.ndarray): Gas fraction profile
        gas_density_array (np.ndarray): Gas density profile
        liquid_darcy_velocity_array (np.ndarray): Liquid Darcy velocity profile
        gas_darcy_velocity_array (np.ndarray): Gas Darcy velocity profile
    """

    def __init__(
        self,
        non_dimensional_parameters,
        temperature_array,
        temperature_derivative_array,
        concentration_array,
        hydrostatic_pressure_array,
        frozen_gas_fraction,
        mushy_layer_depth,
        height_array,
    ):
        self.name = non_dimensional_parameters.name
        self.params = non_dimensional_parameters

        self.temperature_array = np.array(temperature_array)
        self.temperature_derivative_array = np.array(temperature_derivative_array)
        self.concentration_array = np.array(concentration_array)
        self.hydrostatic_pressure_array = np.array(hydrostatic_pressure_array)
        self.frozen_gas_fraction = frozen_gas_fraction
        self.mushy_layer_depth = mushy_layer_depth
        self.height_array = np.array(height_array)

        # Calculate all mushy layer arrays
        if self.params.model_choice == "instant":
            model = MODEL_OPTIONS[self.params.model_choice](
                self.params,
                self.height_array,
                self.temperature_array,
                self.temperature_derivative_array,
                self.hydrostatic_pressure_array,
                self.frozen_gas_fraction,
                self.mushy_layer_depth,
            )
        else:
            model = MODEL_OPTIONS[self.params.model_choice](
                self.params,
                self.height_array,
                self.temperature_array,
                self.temperature_derivative_array,
                self.concentration_array,
                self.hydrostatic_pressure_array,
                self.frozen_gas_fraction,
                self.mushy_layer_depth,
            )
        self.solid_salinity_array = model.solid_salinity
        self.liquid_salinity_array = model.liquid_salinity
        self.solid_fraction_array = model.solid_fraction
        self.liquid_fraction_array = model.liquid_fraction
        self.gas_fraction_array = model.gas_fraction
        self.gas_density_array = model.gas_density
        self.liquid_darcy_velocity_array = model.liquid_darcy_velocity
        self.gas_darcy_velocity_array = model.gas_darcy_velocity

    def save(self, filename: str) -> None:
        data = {
            "non_dimensional_parameters": asdict(self.params),
            "temperature_array": self.temperature_array.tolist(),
            "temperature_derivative_array": self.temperature_derivative_array.tolist(),
            "concentration_array": self.concentration_array.tolist(),
            "hydrostatic_pressure_array": self.hydrostatic_pressure_array.tolist(),
            "frozen_gas_fraction": self.frozen_gas_fraction,
            "mushy_layer_depth": self.mushy_layer_depth,
            "height_array": self.height_array.tolist(),
        }
        with open(f"{filename}.json", "w") as fp:
            json.dump(data, fp, indent=4)

    @classmethod
    def load(cls, filename: str):
        with open(f"{filename}.json", "r") as fp:
            data = json.load(fp)
        params = data["non_dimensional_parameters"]
        data["non_dimensional_parameters"] = NonDimensionalParams(**params)
        return cls(**data)

    def liquid_salinity(self, height):
        return np.interp(
            height, self.height_array, self.liquid_salinity_array, right=np.nan
        )

    def temperature(self, height):
        liquid_darcy_velocity_at_bottom = self.liquid_darcy_velocity_array[0]
        liquid_range = np.linspace(-10, -1.1, 100)
        liquid_values = self.params.far_temperature_scaled * (
            1
            - np.exp(
                (1 + liquid_darcy_velocity_at_bottom)
                * self.mushy_layer_depth
                * (liquid_range + 1)
            )
        )
        return np.interp(
            height,
            np.hstack((liquid_range, self.height_array)),
            np.hstack((liquid_values, self.temperature_array)),
            left=self.params.far_temperature_scaled,
        )

    def concentration(self, height):
        return np.interp(
            height, self.height_array, self.concentration_array, right=np.nan
        )

    def hydrostatic_pressure(self, height):
        return np.interp(
            height, self.height_array, self.hydrostatic_pressure_array, right=np.nan
        )

    def solid_fraction(self, height):
        return np.interp(
            height,
            self.height_array,
            self.solid_fraction_array,
            right=1 - self.frozen_gas_fraction,
        )

    def liquid_fraction(self, height):
        return np.interp(height, self.height_array, self.liquid_fraction_array, right=0)

    def gas_fraction(self, height):
        return np.interp(
            height,
            self.height_array,
            self.gas_fraction_array,
            left=0,
            right=self.frozen_gas_fraction,
        )

    def liquid_darcy_velocity(self, height):
        return np.interp(
            height, self.height_array, self.liquid_darcy_velocity_array, right=np.nan
        )

    def gas_darcy_velocity(self, height):
        return np.interp(
            height,
            self.height_array,
            self.gas_darcy_velocity_array,
            left=np.nan,
            right=0,
        )

    def gas_density(self, height):
        gas_density_filtered = np.where(
            self.gas_fraction_array <= 0, np.nan, self.gas_density_array
        )
        return np.interp(height, self.height_array, gas_density_filtered, left=np.nan)

params

NonDimensionalParams dataclass

Non-dimensional parameters for the three phase mushy layer system.

Note: these can be initialised directly or by using the non_dimensionalise method of PhysicalParams.

This class also provides save and load methods to serialize and deserialize to JSON.

Parameters:

Name Type Description Default
name str

Name of the simulation

 required
model_choice str

Choice of model to use, either "full", "incompressible", "reduced" or "instant"

 required
concentration_ratio float

Ratio of the far field salinity to the salinity difference

 required
stefan_number float

Ratio of the latent heat to the sensible heat

 required
hele_shaw_permeability_scaled float

Ratio of the permeability of the Hele-Shaw cell to the reference permeability

 required
far_temperature_scaled float

Non-dimensionalised far field temperature

 required
damkholer_number float

Ratio of the thermal time scale to the nucleation time scale

 required
expansion_coefficient float

Relative volume increase on exsolving a saturation amount of dissolved gas

 required
stokes_rise_velocity_scaled float

Ratio of the bubble Stokes' rise velocity to the reference velocity

 required
bubble_radius_scaled float

Ratio of the bubble radius to the pore throat radius

 required
pore_throat_exponent float

Power law exponent for pore throat radius as a function of porosity

 required
far_dissolved_concentration_scaled float

Ratio of the far field dissolved gas concentration to the saturation concentration

 required
Source code in mush3p/params.py 
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
296
297
298
299
300
301
302
303
304
305
306
307
308
@dataclass
class NonDimensionalParams:
    """Non-dimensional parameters for the three phase mushy layer system.

    Note: these can be initialised directly or by using the non_dimensionalise method of PhysicalParams.

    This class also provides save and load methods to serialize and deserialize to JSON.

    Args:
        name (str): Name of the simulation
        model_choice (str): Choice of model to use, either "full", "incompressible", "reduced" or "instant"
        concentration_ratio (float): Ratio of the far field salinity to the salinity difference
        stefan_number (float): Ratio of the latent heat to the sensible heat
        hele_shaw_permeability_scaled (float): Ratio of the permeability of the Hele-Shaw cell to the reference permeability
        far_temperature_scaled (float): Non-dimensionalised far field temperature
        damkholer_number (float): Ratio of the thermal time scale to the nucleation time scale
        expansion_coefficient (float): Relative volume increase on exsolving a saturation amount of dissolved gas
        stokes_rise_velocity_scaled (float): Ratio of the bubble Stokes' rise velocity to the reference velocity
        bubble_radius_scaled (float): Ratio of the bubble radius to the pore throat radius
        pore_throat_exponent (float): Power law exponent for pore throat radius as a function of porosity
        far_dissolved_concentration_scaled (float): Ratio of the far field dissolved gas concentration to the saturation concentration
    """

    name: str
    model_choice: str

    # mushy layer params
    concentration_ratio: float
    stefan_number: float
    hele_shaw_permeability_scaled: float
    far_temperature_scaled: float
    solid_conductivity_ratio: float
    solid_specific_heat_capacity_ratio: float
    gas_specific_heat_capacity_ratio: float

    # gas params
    damkholer_number: float
    expansion_coefficient: float
    stokes_rise_velocity_scaled: float
    bubble_radius_scaled: float
    pore_throat_exponent: float
    far_dissolved_concentration_scaled: float
    gas_conductivity_ratio: float
    gas_density_ratio: float

    # compressible gas params
    hydrostatic_pressure_scale: float
    laplace_pressure_scale: float
    kelvin_conversion_temperature: float
    atmospheric_pressure_scaled: float

    @classmethod
    def load(cls, filename: str) -> NonDimensionalParams:
        params = json.load(open(f"{filename}.json"))
        return cls(**params)

    def save(self, filename: str) -> None:
        json.dump(asdict(self), open(f"{filename}.json", "w"), indent=4)

PhysicalParams dataclass

Dimensional parameters for the three phase mushy layer system.

This class implements the calculation of non-dimensional parameters for the system and provides the non_dimensionalise method to provide a NonDimensionalParams object.

This class also provides save and load methods to serialize and deserialize to JSON.

Parameters:

Name Type Description Default
name str

Name of the simulation

 required
model_choice str

Choice of model to use, either "full", "incompressible", "reduced" or "instant"

 'full'
bubble_radius float

Radius of the bubbles [m]

 0.001
nucleation_time_scale float

Time scale for bubble nucleation [s]

 1000
far_salinity float

Salinity of the far field [g/kg]

 35
far_temperature float

Temperature of the far field [degC]

 0.1
far_dissolved_gas_concentration float

Dissolved gas concentration in the far field [kg/kg]

 3.71e-05
reference_velocity float

Pulling speed [m/s], the default value is 1 micron/s

 1e-06
liquid_density float

Density of the liquid [kg/m3]

 1028
gravitational_acceleration float

Acceleration due to gravity [m/s2]

 9.81
liquid_dynamic_viscosity float

Dynamic viscosity of the liquid [kg/m s], the default value gives the same kinematic viscosity used in Moreau et al 2014.

 0.00278
surface_tension float

Surface tension of the liquid [N/m]

 0.077
liquidus_slope float

Slope of the liquidus curve [deg C / g/kg]

 0.07
eutectic_temperature float

Eutectic temperature [deg C]

 -21.2
latent_heat float

Latent heat of fusion [J/kg]

 333000.0
liquid_specific_heat_capacity float

Specific heat capacity of the liquid [J/kg degC]

 4209
solid_specific_heat_capacity float

Specific heat capacity of the solid [J/kg degC]

 2108
gas_specific_heat_capacity float

Specific heat capacity of the gas [J/kg degC]

 1004
liquid_thermal_conductivity float

Thermal conductivity of the liquid [W/m degC]

 0.523
solid_thermal_conductivity float

Thermal conductivity of the solid [W/m degC]

 2.22
gas_thermal_conductivity float

Thermal conductivity of the gas [W/m degC]

 0.02
hele_shaw_gap_width float

Width of the Hele-Shaw cell [m], the default value is the same value used in the experiments of Peppin et al 2007

 0.005
reference_permeability float

Reference permeability [m2], the defulat value is the value used in Rees Jones and Worster 2014

 1e-08
reference_pore_scale float

Fitted pore throat radius at zero solid fraction [m], default value from Maus et al 2021

 0.000195
pore_throat_exponent float

Exponent for pore throat radius power law as a function of porosity, default value from Maus et al 2021

 0.46
reference_saturation_concentration float

Saturation concentration for dissolved gas [kg/kg]

 3.71e-05
specific_gas_constant float

Specific gas constant [J/kg degK]

 286
atmospheric_pressure float

Atmospheric pressure [Pa]

 101000.0
Source code in mush3p/params.py 
  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
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
@dataclass
class PhysicalParams:
    """Dimensional parameters for the three phase mushy layer system.

    This class implements the calculation of non-dimensional parameters for the system
    and provides the non_dimensionalise method to provide a NonDimensionalParams object.

    This class also provides save and load methods to serialize and deserialize to JSON.

    Args:
        name (str): Name of the simulation
        model_choice (str): Choice of model to use, either "full", "incompressible", "reduced" or "instant"
        bubble_radius (float): Radius of the bubbles [m]
        nucleation_time_scale (float): Time scale for bubble nucleation [s]
        far_salinity (float): Salinity of the far field [g/kg]
        far_temperature (float): Temperature of the far field [degC]
        far_dissolved_gas_concentration (float): Dissolved gas concentration in the far field [kg/kg]
        reference_velocity (float): Pulling speed [m/s], the default value is 1 micron/s
        liquid_density (float): Density of the liquid [kg/m3]
        gravitational_acceleration (float): Acceleration due to gravity [m/s2]
        liquid_dynamic_viscosity (float): Dynamic viscosity of the liquid [kg/m s],
            the default value gives the same kinematic viscosity used in Moreau et al 2014.
        surface_tension (float): Surface tension of the liquid [N/m]
        liquidus_slope (float): Slope of the liquidus curve [deg C / g/kg]
        eutectic_temperature (float): Eutectic temperature [deg C]
        latent_heat (float): Latent heat of fusion [J/kg]
        liquid_specific_heat_capacity (float): Specific heat capacity of the liquid [J/kg degC]
        solid_specific_heat_capacity (float): Specific heat capacity of the solid [J/kg degC]
        gas_specific_heat_capacity (float): Specific heat capacity of the gas [J/kg degC]
        liquid_thermal_conductivity (float): Thermal conductivity of the liquid [W/m degC]
        solid_thermal_conductivity (float): Thermal conductivity of the solid [W/m degC]
        gas_thermal_conductivity (float): Thermal conductivity of the gas [W/m degC]
        hele_shaw_gap_width (float): Width of the Hele-Shaw cell [m], the default value
            is the same value used in the experiments of Peppin et al 2007
        reference_permeability (float): Reference permeability [m2], the defulat value
            is the value used in Rees Jones and Worster 2014
        reference_pore_scale (float): Fitted pore throat radius at zero solid fraction [m],
            default value from Maus et al 2021
        pore_throat_exponent (float): Exponent for pore throat radius power law as a function of porosity,
            default value from Maus et al 2021
        reference_saturation_concentration (float): Saturation concentration for dissolved gas [kg/kg]
        specific_gas_constant (float): Specific gas constant [J/kg degK]
        atmospheric_pressure (float): Atmospheric pressure [Pa]
    """

    name: str
    model_choice: str = "full"
    bubble_radius: float = 1e-3
    nucleation_time_scale: float = 1000
    far_salinity: float = 35
    far_temperature: float = 0.1
    far_dissolved_gas_concentration: float = 3.71e-5
    reference_velocity: float = 1e-6
    liquid_density: float = 1028
    gravitational_acceleration: float = 9.81
    liquid_dynamic_viscosity: float = 2.78e-3
    surface_tension: float = 77e-3
    liquidus_slope: float = 0.07
    eutectic_temperature: float = -21.2
    latent_heat: float = 333e3
    liquid_specific_heat_capacity: float = 4209
    solid_specific_heat_capacity: float = 2108
    gas_specific_heat_capacity: float = 1004
    liquid_thermal_conductivity: float = 0.523
    solid_thermal_conductivity: float = 2.22
    gas_thermal_conductivity: float = 2e-2
    hele_shaw_gap_width: float = 5e-3
    reference_permeability: float = 1e-8
    reference_pore_scale: float = 1.95e-4
    pore_throat_exponent: float = 0.46
    reference_saturation_concentration: float = 3.71e-5
    specific_gas_constant: float = 286
    atmospheric_pressure: float = 1.01e5

    @property
    def initial_temperature(self) -> float:
        """Liquidus freezing temperature of the liquid at salinty far_salinity"""
        return -self.liquidus_slope * self.far_salinity

    @property
    def eutectic_salinity(self) -> float:
        """Calculated eutectic salinity from linear liquidus relation"""
        return -self.eutectic_temperature / self.liquidus_slope

    @property
    def liquid_thermal_diffusivity(self) -> float:
        return self.liquid_thermal_conductivity / (
            self.liquid_density * self.liquid_specific_heat_capacity
        )

    @property
    def length_scale(self) -> float:
        return self.liquid_thermal_diffusivity / self.reference_velocity

    @property
    def time_scale(self) -> float:
        return self.liquid_thermal_diffusivity / self.reference_velocity**2

    @property
    def reference_gas_density(self) -> float:
        return self.atmospheric_pressure / (
            self.specific_gas_constant * (self.initial_temperature + CELSIUS_TO_KELVIN)
        )

    @property
    def gas_density_ratio(self) -> float:
        return self.reference_gas_density / self.liquid_density

    @property
    def pressure_scale(self) -> float:
        return (
            self.liquid_thermal_diffusivity
            * self.liquid_dynamic_viscosity
            / self.reference_permeability
        )

    @property
    def concentration_ratio(self) -> float:
        salinity_diff = self.eutectic_salinity - self.far_salinity
        return self.far_salinity / salinity_diff

    @property
    def stefan_number(self) -> float:
        temperature_diff = self.initial_temperature - self.eutectic_temperature
        return self.latent_heat / (
            temperature_diff * self.liquid_specific_heat_capacity
        )

    @property
    def hele_shaw_permeability_scaled(self) -> float:
        return self.hele_shaw_gap_width**2 / (12 * self.reference_permeability)

    @property
    def far_temperature_scaled(self) -> float:
        return (self.far_temperature - self.initial_temperature) / (
            self.initial_temperature - self.eutectic_temperature
        )

    @property
    def damkholer_number(self) -> float:
        return self.time_scale / self.nucleation_time_scale

    @property
    def expansion_coefficient(self) -> float:
        return (
            self.liquid_density * self.reference_saturation_concentration
        ) / self.reference_gas_density

    @property
    def stokes_rise_velocity_scaled(self) -> float:
        return (
            self.liquid_density
            * self.gravitational_acceleration
            * self.reference_pore_scale**2
        ) / (3 * self.liquid_dynamic_viscosity * self.reference_velocity)

    @property
    def bubble_radius_scaled(self) -> float:
        return self.bubble_radius / self.reference_pore_scale

    @property
    def far_dissolved_concentration_scaled(self) -> float:
        return (
            self.far_dissolved_gas_concentration
            / self.reference_saturation_concentration
        )

    @property
    def gas_conductivity_ratio(self) -> float:
        return self.gas_thermal_conductivity / self.liquid_thermal_conductivity

    @property
    def solid_conductivity_ratio(self) -> float:
        return self.solid_thermal_conductivity / self.liquid_thermal_conductivity

    @property
    def solid_specific_heat_capacity_ratio(self) -> float:
        return self.solid_specific_heat_capacity / self.liquid_specific_heat_capacity

    @property
    def gas_specific_heat_capacity_ratio(self) -> float:
        return self.gas_specific_heat_capacity / self.liquid_specific_heat_capacity

    @property
    def hydrostatic_pressure_scale(self) -> float:
        return (
            self.liquid_density
            * self.gravitational_acceleration
            * self.liquid_thermal_diffusivity
        ) / (self.atmospheric_pressure * self.reference_velocity)

    @property
    def laplace_pressure_scale(self) -> float:
        return (
            2 * self.surface_tension / (self.bubble_radius * self.atmospheric_pressure)
        )

    @property
    def kelvin_conversion_temperature(self) -> float:
        return (self.initial_temperature + CELSIUS_TO_KELVIN) / (
            self.initial_temperature - self.eutectic_temperature
        )

    @property
    def atmospheric_pressure_scaled(self) -> float:
        return self.atmospheric_pressure / self.pressure_scale

    def non_dimensionalise(self) -> NonDimensionalParams:
        non_dimensional_params: Dict[str, Any] = {
            "name": self.name,
            "model_choice": self.model_choice,
            "concentration_ratio": self.concentration_ratio,
            "stefan_number": self.stefan_number,
            "hele_shaw_permeability_scaled": self.hele_shaw_permeability_scaled,
            "far_temperature_scaled": self.far_temperature_scaled,
            "damkholer_number": self.damkholer_number,
            "expansion_coefficient": self.expansion_coefficient,
            "stokes_rise_velocity_scaled": self.stokes_rise_velocity_scaled,
            "bubble_radius_scaled": self.bubble_radius_scaled,
            "pore_throat_exponent": self.pore_throat_exponent,
            "far_dissolved_concentration_scaled": self.far_dissolved_concentration_scaled,
            "gas_conductivity_ratio": self.gas_conductivity_ratio,
            "solid_conductivity_ratio": self.solid_conductivity_ratio,
            "solid_specific_heat_capacity_ratio": self.solid_specific_heat_capacity_ratio,
            "gas_specific_heat_capacity_ratio": self.gas_specific_heat_capacity_ratio,
            "hydrostatic_pressure_scale": self.hydrostatic_pressure_scale,
            "laplace_pressure_scale": self.laplace_pressure_scale,
            "kelvin_conversion_temperature": self.kelvin_conversion_temperature,
            "atmospheric_pressure_scaled": self.atmospheric_pressure_scaled,
            "gas_density_ratio": self.gas_density_ratio,
        }
        return NonDimensionalParams(**non_dimensional_params)

    @classmethod
    def load(cls, filename: str) -> PhysicalParams:
        params = json.load(open(f"{filename}.json"))
        return cls(**params)

    def save(self, filename: str) -> None:
        json.dump(asdict(self), open(f"{filename}.json", "w"), indent=4)

eutectic_salinity: float property

Calculated eutectic salinity from linear liquidus relation

initial_temperature: float property

Liquidus freezing temperature of the liquid at salinty far_salinity

static_settings

get_initial_solution(model_choice)

Get intial solution for scipy solve_BVP which satisfies the boundary conditions.

Parameters:

Name Type Description Default
model_choice str

one of "full", "incompressible" "reduced" or "instant"

 required

Returns: NDArray: initial solution for scipy solve_BVP

Source code in mush3p/static_settings.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
def get_initial_solution(model_choice: str) -> NDArray:
    """Get intial solution for scipy solve_BVP which satisfies the boundary conditions.

    Args:
        model_choice (str): one of "full", "incompressible" "reduced" or "instant"
    Returns:
        NDArray: initial solution for scipy solve_BVP
    """
    if model_choice == "instant":
        return np.vstack(
            (
                INITIAL_TEMPERATURE,
                INITIAL_TEMPERATURE_DERIVATIVE,
                INITIAL_HYDROSTATIC_PRESSURE,
                INITIAL_FROZEN_GAS_FRACTION,
                INITIAL_MUSHY_LAYER_DEPTH,
            )
        )
    return np.vstack(
        (
            INITIAL_TEMPERATURE,
            INITIAL_TEMPERATURE_DERIVATIVE,
            INITIAL_DISSOLVED_GAS_CONCENTRATION,
            INITIAL_HYDROSTATIC_PRESSURE,
            INITIAL_FROZEN_GAS_FRACTION,
            INITIAL_MUSHY_LAYER_DEPTH,
        )
    )