PySDM.backends.impl_numba.methods.chemistry_methods
CPU implementation of backend methods for aqueous chemistry
1""" 2CPU implementation of backend methods for aqueous chemistry 3""" 4 5from collections import namedtuple 6 7import numba 8import numpy as np 9 10from PySDM.backends.impl_common.backend_methods import BackendMethods 11from PySDM.backends.impl_numba import conf 12from PySDM.backends.impl_numba.toms748 import toms748_solve 13from PySDM.dynamics.impl.chemistry_utils import ( 14 DIFFUSION_CONST, 15 DISSOCIATION_FACTORS, 16 GASEOUS_COMPOUNDS, 17 MASS_ACCOMMODATION_COEFFICIENTS, 18 EquilibriumConsts, 19 HenryConsts, 20 KineticConsts, 21 SpecificGravities, 22 k4, 23) 24from PySDM.physics.constants import K_H2O 25 26_MAX_ITER_QUITE_CLOSE = 8 27_MAX_ITER_DEFAULT = 32 28_REALY_CLOSE_THRESHOLD = 1e-6 29_QUITE_CLOSE_THRESHOLD = 1 30_QUITE_CLOSE_MULTIPLIER = 2 31 32_K = namedtuple("_K", ("NH3", "SO2", "HSO3", "HSO4", "HCO3", "CO2", "HNO3")) 33_conc = namedtuple("_conc", ("N_mIII", "N_V", "C_IV", "S_IV", "S_VI")) 34 35 36class ChemistryMethods(BackendMethods): 37 def __init__(self): 38 BackendMethods.__init__(self) 39 self.HENRY_CONST = HenryConsts(self.formulae) 40 self.KINETIC_CONST = KineticConsts(self.formulae) 41 self.EQUILIBRIUM_CONST = EquilibriumConsts(self.formulae) 42 self.specific_gravities = SpecificGravities(self.formulae.constants) 43 44 def dissolution( # pylint:disable=too-many-locals 45 self, 46 *, 47 n_cell, 48 n_threads, 49 cell_order, 50 cell_start_arg, 51 idx, 52 do_chemistry_flag, 53 mole_amounts, 54 env_mixing_ratio, 55 env_T, 56 env_p, 57 env_rho_d, 58 dissociation_factors, 59 timestep, 60 dv, 61 system_type, 62 droplet_volume, 63 multiplicity, 64 ): 65 assert n_cell == 1 66 assert n_threads == 1 67 for i in range(n_cell): 68 cell_id = cell_order[i] 69 70 cell_start = cell_start_arg[cell_id] 71 cell_end = cell_start_arg[cell_id + 1] 72 n_sd_in_cell = cell_end - cell_start 73 if n_sd_in_cell == 0: 74 continue 75 76 super_droplet_ids = numba.typed.List() 77 for sd_id in idx[cell_start:cell_end]: 78 if do_chemistry_flag.data[sd_id]: 79 super_droplet_ids.append(sd_id) 80 81 if len(super_droplet_ids) == 0: 82 return 83 84 for key, compound in GASEOUS_COMPOUNDS.items(): 85 ChemistryMethods.dissolution_body( 86 super_droplet_ids=super_droplet_ids, 87 mole_amounts=mole_amounts[key].data, 88 env_mixing_ratio=env_mixing_ratio[compound][cell_id : cell_id + 1], 89 henrysConstant=self.HENRY_CONST.HENRY_CONST[compound].at( 90 env_T[cell_id] 91 ), 92 env_p=env_p[cell_id], 93 env_T=env_T[cell_id], 94 env_rho_d=env_rho_d[cell_id], 95 timestep=timestep, 96 dv=dv, 97 droplet_volume=droplet_volume.data, 98 multiplicity=multiplicity.data, 99 system_type=system_type, 100 specific_gravity=self.specific_gravities[compound], 101 alpha=MASS_ACCOMMODATION_COEFFICIENTS[compound], 102 diffusion_const=DIFFUSION_CONST[compound], 103 dissociation_factor=dissociation_factors[compound].data, 104 radius=self.formulae.trivia.radius, 105 const=self.formulae.constants, 106 ) 107 108 @staticmethod 109 @numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False}}) 110 def dissolution_body( # pylint: disable=too-many-locals 111 *, 112 super_droplet_ids, 113 mole_amounts, 114 env_mixing_ratio, 115 henrysConstant, 116 env_p, 117 env_T, 118 env_rho_d, 119 timestep, 120 dv, 121 droplet_volume, 122 multiplicity, 123 system_type, 124 specific_gravity, 125 alpha, 126 diffusion_const, 127 dissociation_factor, 128 radius, 129 const, 130 ): 131 mole_amount_taken = 0 132 for i in super_droplet_ids: 133 Mc = specific_gravity * const.Md 134 Rc = const.R_str / Mc 135 cinf = env_p / env_T / (const.Rd / env_mixing_ratio[0] + Rc) / Mc 136 r_w = radius(volume=droplet_volume[i]) 137 v_avg = np.sqrt(8 * const.R_str * env_T / (np.pi * Mc)) 138 dt_over_scale = timestep / ( 139 4 * r_w / (3 * v_avg * alpha) + r_w**2 / (3 * diffusion_const) 140 ) 141 A_old = mole_amounts[i] / droplet_volume[i] 142 H_eff = henrysConstant * dissociation_factor[i] 143 A_new = (A_old + dt_over_scale * cinf) / ( 144 1 + dt_over_scale / H_eff / const.R_str / env_T 145 ) 146 new_mole_amount_per_real_droplet = A_new * droplet_volume[i] 147 assert new_mole_amount_per_real_droplet >= 0 148 149 mole_amount_taken += multiplicity[i] * ( 150 new_mole_amount_per_real_droplet - mole_amounts[i] 151 ) 152 mole_amounts[i] = new_mole_amount_per_real_droplet 153 delta_mr = mole_amount_taken * specific_gravity * const.Md / (dv * env_rho_d) 154 assert delta_mr <= env_mixing_ratio 155 if system_type == "closed": 156 env_mixing_ratio -= delta_mr 157 158 def oxidation( # pylint: disable=too-many-locals 159 self, 160 *, 161 n_sd, 162 cell_ids, 163 do_chemistry_flag, 164 k0, 165 k1, 166 k2, 167 k3, 168 K_SO2, 169 K_HSO3, 170 timestep, 171 droplet_volume, 172 pH, 173 dissociation_factor_SO2, 174 # output 175 moles_O3, 176 moles_H2O2, 177 moles_S_IV, 178 moles_S_VI, 179 ): 180 ChemistryMethods.oxidation_body( 181 n_sd=n_sd, 182 cell_ids=cell_ids.data, 183 do_chemistry_flag=do_chemistry_flag.data, 184 explicit_euler=self.formulae.trivia.explicit_euler, 185 pH2H=self.formulae.trivia.pH2H, 186 k0=k0.data, 187 k1=k1.data, 188 k2=k2.data, 189 k3=k3.data, 190 K_SO2=K_SO2.data, 191 K_HSO3=K_HSO3.data, 192 timestep=timestep, 193 droplet_volume=droplet_volume.data, 194 pH=pH.data, 195 dissociation_factor_SO2=dissociation_factor_SO2.data, 196 # output 197 moles_O3=moles_O3.data, 198 moles_H2O2=moles_H2O2.data, 199 moles_S_IV=moles_S_IV.data, 200 moles_S_VI=moles_S_VI.data, 201 ) 202 203 @staticmethod 204 @numba.njit(**conf.JIT_FLAGS) 205 def oxidation_body( # pylint: disable=too-many-locals 206 *, 207 n_sd, 208 cell_ids, 209 do_chemistry_flag, 210 explicit_euler, 211 pH2H, 212 k0, 213 k1, 214 k2, 215 k3, 216 K_SO2, 217 K_HSO3, 218 timestep, 219 droplet_volume, 220 pH, 221 dissociation_factor_SO2, 222 # output 223 moles_O3, 224 moles_H2O2, 225 moles_S_IV, 226 moles_S_VI, 227 ): 228 for i in numba.prange(n_sd): # pylint: disable=not-an-iterable 229 if not do_chemistry_flag[i]: 230 continue 231 232 cid = cell_ids[i] 233 H = pH2H(pH[i]) 234 SO2aq = moles_S_IV[i] / droplet_volume[i] / dissociation_factor_SO2[i] 235 236 # NB: This might not be entirely correct 237 # https://doi.org/10.1029/JD092iD04p04171 238 # https://doi.org/10.5194/acp-16-1693-2016 239 240 ozone = ( 241 ( 242 k0[cid] 243 + (k1[cid] * K_SO2[cid] / H) 244 + (k2[cid] * K_SO2[cid] * K_HSO3[cid] / H**2) 245 ) 246 * (moles_O3[i] / droplet_volume[i]) 247 * SO2aq 248 ) 249 peroxide = ( 250 k3[cid] 251 * K_SO2[cid] 252 / (1 + k4 * H) 253 * (moles_H2O2[i] / droplet_volume[i]) 254 * SO2aq 255 ) 256 dt_times_volume = timestep * droplet_volume[i] 257 258 dconc_dt_O3 = -ozone 259 dconc_dt_S_IV = -(ozone + peroxide) 260 dconc_dt_H2O2 = -peroxide 261 dconc_dt_S_VI = ozone + peroxide 262 263 if ( 264 moles_O3[i] + dconc_dt_O3 * dt_times_volume < 0 265 or moles_S_IV[i] + dconc_dt_S_IV * dt_times_volume < 0 266 or moles_S_VI[i] + dconc_dt_S_VI * dt_times_volume < 0 267 or moles_H2O2[i] + dconc_dt_H2O2 * dt_times_volume < 0 268 ): 269 continue 270 271 moles_O3[i] = explicit_euler(moles_O3[i], dt_times_volume, dconc_dt_O3) 272 moles_S_IV[i] = explicit_euler( 273 moles_S_IV[i], dt_times_volume, dconc_dt_S_IV 274 ) 275 moles_S_VI[i] = explicit_euler( 276 moles_S_VI[i], dt_times_volume, dconc_dt_S_VI 277 ) 278 moles_H2O2[i] = explicit_euler( 279 moles_H2O2[i], dt_times_volume, dconc_dt_H2O2 280 ) 281 282 def chem_recalculate_drop_data( 283 self, dissociation_factors, equilibrium_consts, cell_id, pH 284 ): 285 for i in range(len(pH)): 286 H = self.formulae.trivia.pH2H(pH.data[i]) 287 for key in DIFFUSION_CONST: 288 dissociation_factors[key].data[i] = DISSOCIATION_FACTORS[key]( 289 H, equilibrium_consts, cell_id.data[i] 290 ) 291 292 def chem_recalculate_cell_data( 293 self, equilibrium_consts, kinetic_consts, temperature 294 ): 295 for i in range(len(temperature)): 296 for key in equilibrium_consts: 297 equilibrium_consts[key].data[i] = ( 298 self.EQUILIBRIUM_CONST.EQUILIBRIUM_CONST[key].at( 299 temperature.data[i] 300 ) 301 ) 302 for key in kinetic_consts: 303 kinetic_consts[key].data[i] = self.KINETIC_CONST.KINETIC_CONST[key].at( 304 temperature.data[i] 305 ) 306 307 def equilibrate_H( 308 self, 309 *, 310 equilibrium_consts, 311 cell_id, 312 conc, 313 do_chemistry_flag, 314 pH, 315 H_min, 316 H_max, 317 ionic_strength_threshold, 318 rtol, 319 ): 320 ChemistryMethods.equilibrate_H_body( 321 within_tolerance=self.formulae.trivia.within_tolerance, 322 pH2H=self.formulae.trivia.pH2H, 323 H2pH=self.formulae.trivia.H2pH, 324 cell_id=cell_id.data, 325 conc=_conc( 326 N_mIII=conc.N_mIII.data, 327 N_V=conc.N_V.data, 328 C_IV=conc.C_IV.data, 329 S_IV=conc.S_IV.data, 330 S_VI=conc.S_VI.data, 331 ), 332 K=_K( 333 NH3=equilibrium_consts["K_NH3"].data, 334 SO2=equilibrium_consts["K_SO2"].data, 335 HSO3=equilibrium_consts["K_HSO3"].data, 336 HSO4=equilibrium_consts["K_HSO4"].data, 337 HCO3=equilibrium_consts["K_HCO3"].data, 338 CO2=equilibrium_consts["K_CO2"].data, 339 HNO3=equilibrium_consts["K_HNO3"].data, 340 ), 341 # output 342 do_chemistry_flag=do_chemistry_flag.data, 343 pH=pH.data, 344 # params 345 H_min=H_min, 346 H_max=H_max, 347 ionic_strength_threshold=ionic_strength_threshold, 348 rtol=rtol, 349 ) 350 351 @staticmethod 352 @numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False, "cache": False}}) 353 def equilibrate_H_body( # pylint: disable=too-many-arguments,too-many-locals 354 within_tolerance, 355 pH2H, 356 H2pH, 357 cell_id, 358 conc, 359 K, 360 do_chemistry_flag, 361 pH, 362 # params 363 H_min, 364 H_max, 365 ionic_strength_threshold, 366 rtol, 367 ): 368 for i, pH_i in enumerate(pH): 369 cid = cell_id[i] 370 args = ( 371 _conc( 372 N_mIII=conc.N_mIII[i], 373 N_V=conc.N_V[i], 374 C_IV=conc.C_IV[i], 375 S_IV=conc.S_IV[i], 376 S_VI=conc.S_VI[i], 377 ), 378 _K( 379 NH3=K.NH3[cid], 380 SO2=K.SO2[cid], 381 HSO3=K.HSO3[cid], 382 HSO4=K.HSO4[cid], 383 HCO3=K.HCO3[cid], 384 CO2=K.CO2[cid], 385 HNO3=K.HNO3[cid], 386 ), 387 ) 388 a = pH2H(pH_i) 389 fa = acidity_minfun(a, *args) 390 if abs(fa) < _REALY_CLOSE_THRESHOLD: 391 continue 392 b = np.nan 393 fb = np.nan 394 use_default_range = False 395 if abs(fa) < _QUITE_CLOSE_THRESHOLD: 396 b = a * _QUITE_CLOSE_MULTIPLIER 397 fb = acidity_minfun(b, *args) 398 if fa * fb > 0: 399 b = a 400 fb = fa 401 a = b / _QUITE_CLOSE_MULTIPLIER / _QUITE_CLOSE_MULTIPLIER 402 fa = acidity_minfun(a, *args) 403 if fa * fb > 0: 404 use_default_range = True 405 else: 406 use_default_range = True 407 if use_default_range: 408 a = H_min 409 b = H_max 410 fa = acidity_minfun(a, *args) 411 fb = acidity_minfun(b, *args) 412 max_iter = _MAX_ITER_DEFAULT 413 else: 414 max_iter = _MAX_ITER_QUITE_CLOSE 415 H, _iters_taken = toms748_solve( 416 acidity_minfun, 417 args, 418 a, 419 b, 420 fa, 421 fb, 422 rtol=rtol, 423 max_iter=max_iter, 424 within_tolerance=within_tolerance, 425 ) 426 assert _iters_taken != max_iter 427 pH[i] = H2pH(H) 428 ionic_strength = calc_ionic_strength(H, *args) 429 do_chemistry_flag[i] = ionic_strength <= ionic_strength_threshold 430 431 432@numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False}}) 433def calc_ionic_strength(H, conc, K): 434 # Directly adapted 435 # https://github.com/igfuw/libcloudphxx/blob/0b4e2455fba4f95c7387623fc21481a85e7b151f/src/impl/particles_impl_chem_strength.ipp#L50 436 # https://en.wikipedia.org/wiki/Ionic_strength 437 438 # H+ and OH- 439 water = H + K_H2O / H 440 441 # HSO4- and SO4 2- 442 cz_S_VI = H * conc.S_VI / (H + K.HSO4) + 4 * K.HSO4 * conc.S_VI / (H + K.HSO4) 443 444 # HCO3- and CO3 2- 445 cz_CO2 = K.CO2 * H * conc.C_IV / ( 446 H * H + K.CO2 * H + K.CO2 * K.HCO3 447 ) + 4 * K.CO2 * K.HCO3 * conc.C_IV / (H * H + K.CO2 * H + K.CO2 * K.HCO3) 448 449 # HSO3- and HSO4 2- 450 cz_SO2 = K.SO2 * H * conc.S_IV / ( 451 H * H + K.SO2 * H + K.SO2 * K.HSO3 452 ) + 4 * K.SO2 * K.HSO3 * conc.S_IV / (H * H + K.SO2 * H + K.SO2 * K.HSO3) 453 454 # NO3- 455 cz_HNO3 = K.HNO3 * conc.N_V / (H + K.HNO3) 456 457 # NH4+ 458 cz_NH3 = K.NH3 * H * conc.N_mIII / (K_H2O + K.NH3 * H) 459 460 return 0.5 * (water + cz_S_VI + cz_CO2 + cz_SO2 + cz_HNO3 + cz_NH3) 461 462 463@numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False}}) 464def acidity_minfun(H, conc, K): 465 ammonia = (conc.N_mIII * H * K.NH3) / (K_H2O + K.NH3 * H) 466 nitric = conc.N_V * K.HNO3 / (H + K.HNO3) 467 sulfous = ( 468 conc.S_IV * K.SO2 * (H + 2 * K.HSO3) / (H * H + H * K.SO2 + K.SO2 * K.HSO3) 469 ) 470 water = K_H2O / H 471 sulfuric = conc.S_VI * (H + 2 * K.HSO4) / (H + K.HSO4) 472 carbonic = ( 473 conc.C_IV * K.CO2 * (H + 2 * K.HCO3) / (H * H + H * K.CO2 + K.CO2 * K.HCO3) 474 ) 475 zero = H + ammonia - (nitric + sulfous + water + sulfuric + carbonic) 476 return zero
37class ChemistryMethods(BackendMethods): 38 def __init__(self): 39 BackendMethods.__init__(self) 40 self.HENRY_CONST = HenryConsts(self.formulae) 41 self.KINETIC_CONST = KineticConsts(self.formulae) 42 self.EQUILIBRIUM_CONST = EquilibriumConsts(self.formulae) 43 self.specific_gravities = SpecificGravities(self.formulae.constants) 44 45 def dissolution( # pylint:disable=too-many-locals 46 self, 47 *, 48 n_cell, 49 n_threads, 50 cell_order, 51 cell_start_arg, 52 idx, 53 do_chemistry_flag, 54 mole_amounts, 55 env_mixing_ratio, 56 env_T, 57 env_p, 58 env_rho_d, 59 dissociation_factors, 60 timestep, 61 dv, 62 system_type, 63 droplet_volume, 64 multiplicity, 65 ): 66 assert n_cell == 1 67 assert n_threads == 1 68 for i in range(n_cell): 69 cell_id = cell_order[i] 70 71 cell_start = cell_start_arg[cell_id] 72 cell_end = cell_start_arg[cell_id + 1] 73 n_sd_in_cell = cell_end - cell_start 74 if n_sd_in_cell == 0: 75 continue 76 77 super_droplet_ids = numba.typed.List() 78 for sd_id in idx[cell_start:cell_end]: 79 if do_chemistry_flag.data[sd_id]: 80 super_droplet_ids.append(sd_id) 81 82 if len(super_droplet_ids) == 0: 83 return 84 85 for key, compound in GASEOUS_COMPOUNDS.items(): 86 ChemistryMethods.dissolution_body( 87 super_droplet_ids=super_droplet_ids, 88 mole_amounts=mole_amounts[key].data, 89 env_mixing_ratio=env_mixing_ratio[compound][cell_id : cell_id + 1], 90 henrysConstant=self.HENRY_CONST.HENRY_CONST[compound].at( 91 env_T[cell_id] 92 ), 93 env_p=env_p[cell_id], 94 env_T=env_T[cell_id], 95 env_rho_d=env_rho_d[cell_id], 96 timestep=timestep, 97 dv=dv, 98 droplet_volume=droplet_volume.data, 99 multiplicity=multiplicity.data, 100 system_type=system_type, 101 specific_gravity=self.specific_gravities[compound], 102 alpha=MASS_ACCOMMODATION_COEFFICIENTS[compound], 103 diffusion_const=DIFFUSION_CONST[compound], 104 dissociation_factor=dissociation_factors[compound].data, 105 radius=self.formulae.trivia.radius, 106 const=self.formulae.constants, 107 ) 108 109 @staticmethod 110 @numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False}}) 111 def dissolution_body( # pylint: disable=too-many-locals 112 *, 113 super_droplet_ids, 114 mole_amounts, 115 env_mixing_ratio, 116 henrysConstant, 117 env_p, 118 env_T, 119 env_rho_d, 120 timestep, 121 dv, 122 droplet_volume, 123 multiplicity, 124 system_type, 125 specific_gravity, 126 alpha, 127 diffusion_const, 128 dissociation_factor, 129 radius, 130 const, 131 ): 132 mole_amount_taken = 0 133 for i in super_droplet_ids: 134 Mc = specific_gravity * const.Md 135 Rc = const.R_str / Mc 136 cinf = env_p / env_T / (const.Rd / env_mixing_ratio[0] + Rc) / Mc 137 r_w = radius(volume=droplet_volume[i]) 138 v_avg = np.sqrt(8 * const.R_str * env_T / (np.pi * Mc)) 139 dt_over_scale = timestep / ( 140 4 * r_w / (3 * v_avg * alpha) + r_w**2 / (3 * diffusion_const) 141 ) 142 A_old = mole_amounts[i] / droplet_volume[i] 143 H_eff = henrysConstant * dissociation_factor[i] 144 A_new = (A_old + dt_over_scale * cinf) / ( 145 1 + dt_over_scale / H_eff / const.R_str / env_T 146 ) 147 new_mole_amount_per_real_droplet = A_new * droplet_volume[i] 148 assert new_mole_amount_per_real_droplet >= 0 149 150 mole_amount_taken += multiplicity[i] * ( 151 new_mole_amount_per_real_droplet - mole_amounts[i] 152 ) 153 mole_amounts[i] = new_mole_amount_per_real_droplet 154 delta_mr = mole_amount_taken * specific_gravity * const.Md / (dv * env_rho_d) 155 assert delta_mr <= env_mixing_ratio 156 if system_type == "closed": 157 env_mixing_ratio -= delta_mr 158 159 def oxidation( # pylint: disable=too-many-locals 160 self, 161 *, 162 n_sd, 163 cell_ids, 164 do_chemistry_flag, 165 k0, 166 k1, 167 k2, 168 k3, 169 K_SO2, 170 K_HSO3, 171 timestep, 172 droplet_volume, 173 pH, 174 dissociation_factor_SO2, 175 # output 176 moles_O3, 177 moles_H2O2, 178 moles_S_IV, 179 moles_S_VI, 180 ): 181 ChemistryMethods.oxidation_body( 182 n_sd=n_sd, 183 cell_ids=cell_ids.data, 184 do_chemistry_flag=do_chemistry_flag.data, 185 explicit_euler=self.formulae.trivia.explicit_euler, 186 pH2H=self.formulae.trivia.pH2H, 187 k0=k0.data, 188 k1=k1.data, 189 k2=k2.data, 190 k3=k3.data, 191 K_SO2=K_SO2.data, 192 K_HSO3=K_HSO3.data, 193 timestep=timestep, 194 droplet_volume=droplet_volume.data, 195 pH=pH.data, 196 dissociation_factor_SO2=dissociation_factor_SO2.data, 197 # output 198 moles_O3=moles_O3.data, 199 moles_H2O2=moles_H2O2.data, 200 moles_S_IV=moles_S_IV.data, 201 moles_S_VI=moles_S_VI.data, 202 ) 203 204 @staticmethod 205 @numba.njit(**conf.JIT_FLAGS) 206 def oxidation_body( # pylint: disable=too-many-locals 207 *, 208 n_sd, 209 cell_ids, 210 do_chemistry_flag, 211 explicit_euler, 212 pH2H, 213 k0, 214 k1, 215 k2, 216 k3, 217 K_SO2, 218 K_HSO3, 219 timestep, 220 droplet_volume, 221 pH, 222 dissociation_factor_SO2, 223 # output 224 moles_O3, 225 moles_H2O2, 226 moles_S_IV, 227 moles_S_VI, 228 ): 229 for i in numba.prange(n_sd): # pylint: disable=not-an-iterable 230 if not do_chemistry_flag[i]: 231 continue 232 233 cid = cell_ids[i] 234 H = pH2H(pH[i]) 235 SO2aq = moles_S_IV[i] / droplet_volume[i] / dissociation_factor_SO2[i] 236 237 # NB: This might not be entirely correct 238 # https://doi.org/10.1029/JD092iD04p04171 239 # https://doi.org/10.5194/acp-16-1693-2016 240 241 ozone = ( 242 ( 243 k0[cid] 244 + (k1[cid] * K_SO2[cid] / H) 245 + (k2[cid] * K_SO2[cid] * K_HSO3[cid] / H**2) 246 ) 247 * (moles_O3[i] / droplet_volume[i]) 248 * SO2aq 249 ) 250 peroxide = ( 251 k3[cid] 252 * K_SO2[cid] 253 / (1 + k4 * H) 254 * (moles_H2O2[i] / droplet_volume[i]) 255 * SO2aq 256 ) 257 dt_times_volume = timestep * droplet_volume[i] 258 259 dconc_dt_O3 = -ozone 260 dconc_dt_S_IV = -(ozone + peroxide) 261 dconc_dt_H2O2 = -peroxide 262 dconc_dt_S_VI = ozone + peroxide 263 264 if ( 265 moles_O3[i] + dconc_dt_O3 * dt_times_volume < 0 266 or moles_S_IV[i] + dconc_dt_S_IV * dt_times_volume < 0 267 or moles_S_VI[i] + dconc_dt_S_VI * dt_times_volume < 0 268 or moles_H2O2[i] + dconc_dt_H2O2 * dt_times_volume < 0 269 ): 270 continue 271 272 moles_O3[i] = explicit_euler(moles_O3[i], dt_times_volume, dconc_dt_O3) 273 moles_S_IV[i] = explicit_euler( 274 moles_S_IV[i], dt_times_volume, dconc_dt_S_IV 275 ) 276 moles_S_VI[i] = explicit_euler( 277 moles_S_VI[i], dt_times_volume, dconc_dt_S_VI 278 ) 279 moles_H2O2[i] = explicit_euler( 280 moles_H2O2[i], dt_times_volume, dconc_dt_H2O2 281 ) 282 283 def chem_recalculate_drop_data( 284 self, dissociation_factors, equilibrium_consts, cell_id, pH 285 ): 286 for i in range(len(pH)): 287 H = self.formulae.trivia.pH2H(pH.data[i]) 288 for key in DIFFUSION_CONST: 289 dissociation_factors[key].data[i] = DISSOCIATION_FACTORS[key]( 290 H, equilibrium_consts, cell_id.data[i] 291 ) 292 293 def chem_recalculate_cell_data( 294 self, equilibrium_consts, kinetic_consts, temperature 295 ): 296 for i in range(len(temperature)): 297 for key in equilibrium_consts: 298 equilibrium_consts[key].data[i] = ( 299 self.EQUILIBRIUM_CONST.EQUILIBRIUM_CONST[key].at( 300 temperature.data[i] 301 ) 302 ) 303 for key in kinetic_consts: 304 kinetic_consts[key].data[i] = self.KINETIC_CONST.KINETIC_CONST[key].at( 305 temperature.data[i] 306 ) 307 308 def equilibrate_H( 309 self, 310 *, 311 equilibrium_consts, 312 cell_id, 313 conc, 314 do_chemistry_flag, 315 pH, 316 H_min, 317 H_max, 318 ionic_strength_threshold, 319 rtol, 320 ): 321 ChemistryMethods.equilibrate_H_body( 322 within_tolerance=self.formulae.trivia.within_tolerance, 323 pH2H=self.formulae.trivia.pH2H, 324 H2pH=self.formulae.trivia.H2pH, 325 cell_id=cell_id.data, 326 conc=_conc( 327 N_mIII=conc.N_mIII.data, 328 N_V=conc.N_V.data, 329 C_IV=conc.C_IV.data, 330 S_IV=conc.S_IV.data, 331 S_VI=conc.S_VI.data, 332 ), 333 K=_K( 334 NH3=equilibrium_consts["K_NH3"].data, 335 SO2=equilibrium_consts["K_SO2"].data, 336 HSO3=equilibrium_consts["K_HSO3"].data, 337 HSO4=equilibrium_consts["K_HSO4"].data, 338 HCO3=equilibrium_consts["K_HCO3"].data, 339 CO2=equilibrium_consts["K_CO2"].data, 340 HNO3=equilibrium_consts["K_HNO3"].data, 341 ), 342 # output 343 do_chemistry_flag=do_chemistry_flag.data, 344 pH=pH.data, 345 # params 346 H_min=H_min, 347 H_max=H_max, 348 ionic_strength_threshold=ionic_strength_threshold, 349 rtol=rtol, 350 ) 351 352 @staticmethod 353 @numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False, "cache": False}}) 354 def equilibrate_H_body( # pylint: disable=too-many-arguments,too-many-locals 355 within_tolerance, 356 pH2H, 357 H2pH, 358 cell_id, 359 conc, 360 K, 361 do_chemistry_flag, 362 pH, 363 # params 364 H_min, 365 H_max, 366 ionic_strength_threshold, 367 rtol, 368 ): 369 for i, pH_i in enumerate(pH): 370 cid = cell_id[i] 371 args = ( 372 _conc( 373 N_mIII=conc.N_mIII[i], 374 N_V=conc.N_V[i], 375 C_IV=conc.C_IV[i], 376 S_IV=conc.S_IV[i], 377 S_VI=conc.S_VI[i], 378 ), 379 _K( 380 NH3=K.NH3[cid], 381 SO2=K.SO2[cid], 382 HSO3=K.HSO3[cid], 383 HSO4=K.HSO4[cid], 384 HCO3=K.HCO3[cid], 385 CO2=K.CO2[cid], 386 HNO3=K.HNO3[cid], 387 ), 388 ) 389 a = pH2H(pH_i) 390 fa = acidity_minfun(a, *args) 391 if abs(fa) < _REALY_CLOSE_THRESHOLD: 392 continue 393 b = np.nan 394 fb = np.nan 395 use_default_range = False 396 if abs(fa) < _QUITE_CLOSE_THRESHOLD: 397 b = a * _QUITE_CLOSE_MULTIPLIER 398 fb = acidity_minfun(b, *args) 399 if fa * fb > 0: 400 b = a 401 fb = fa 402 a = b / _QUITE_CLOSE_MULTIPLIER / _QUITE_CLOSE_MULTIPLIER 403 fa = acidity_minfun(a, *args) 404 if fa * fb > 0: 405 use_default_range = True 406 else: 407 use_default_range = True 408 if use_default_range: 409 a = H_min 410 b = H_max 411 fa = acidity_minfun(a, *args) 412 fb = acidity_minfun(b, *args) 413 max_iter = _MAX_ITER_DEFAULT 414 else: 415 max_iter = _MAX_ITER_QUITE_CLOSE 416 H, _iters_taken = toms748_solve( 417 acidity_minfun, 418 args, 419 a, 420 b, 421 fa, 422 fb, 423 rtol=rtol, 424 max_iter=max_iter, 425 within_tolerance=within_tolerance, 426 ) 427 assert _iters_taken != max_iter 428 pH[i] = H2pH(H) 429 ionic_strength = calc_ionic_strength(H, *args) 430 do_chemistry_flag[i] = ionic_strength <= ionic_strength_threshold
def
dissolution( self, *, n_cell, n_threads, cell_order, cell_start_arg, idx, do_chemistry_flag, mole_amounts, env_mixing_ratio, env_T, env_p, env_rho_d, dissociation_factors, timestep, dv, system_type, droplet_volume, multiplicity):
45 def dissolution( # pylint:disable=too-many-locals 46 self, 47 *, 48 n_cell, 49 n_threads, 50 cell_order, 51 cell_start_arg, 52 idx, 53 do_chemistry_flag, 54 mole_amounts, 55 env_mixing_ratio, 56 env_T, 57 env_p, 58 env_rho_d, 59 dissociation_factors, 60 timestep, 61 dv, 62 system_type, 63 droplet_volume, 64 multiplicity, 65 ): 66 assert n_cell == 1 67 assert n_threads == 1 68 for i in range(n_cell): 69 cell_id = cell_order[i] 70 71 cell_start = cell_start_arg[cell_id] 72 cell_end = cell_start_arg[cell_id + 1] 73 n_sd_in_cell = cell_end - cell_start 74 if n_sd_in_cell == 0: 75 continue 76 77 super_droplet_ids = numba.typed.List() 78 for sd_id in idx[cell_start:cell_end]: 79 if do_chemistry_flag.data[sd_id]: 80 super_droplet_ids.append(sd_id) 81 82 if len(super_droplet_ids) == 0: 83 return 84 85 for key, compound in GASEOUS_COMPOUNDS.items(): 86 ChemistryMethods.dissolution_body( 87 super_droplet_ids=super_droplet_ids, 88 mole_amounts=mole_amounts[key].data, 89 env_mixing_ratio=env_mixing_ratio[compound][cell_id : cell_id + 1], 90 henrysConstant=self.HENRY_CONST.HENRY_CONST[compound].at( 91 env_T[cell_id] 92 ), 93 env_p=env_p[cell_id], 94 env_T=env_T[cell_id], 95 env_rho_d=env_rho_d[cell_id], 96 timestep=timestep, 97 dv=dv, 98 droplet_volume=droplet_volume.data, 99 multiplicity=multiplicity.data, 100 system_type=system_type, 101 specific_gravity=self.specific_gravities[compound], 102 alpha=MASS_ACCOMMODATION_COEFFICIENTS[compound], 103 diffusion_const=DIFFUSION_CONST[compound], 104 dissociation_factor=dissociation_factors[compound].data, 105 radius=self.formulae.trivia.radius, 106 const=self.formulae.constants, 107 )
@staticmethod
@numba.njit(**{**conf.JIT_FLAGS, **{'parallel': False}})
def
dissolution_body( *, super_droplet_ids, mole_amounts, env_mixing_ratio, henrysConstant, env_p, env_T, env_rho_d, timestep, dv, droplet_volume, multiplicity, system_type, specific_gravity, alpha, diffusion_const, dissociation_factor, radius, const):
109 @staticmethod 110 @numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False}}) 111 def dissolution_body( # pylint: disable=too-many-locals 112 *, 113 super_droplet_ids, 114 mole_amounts, 115 env_mixing_ratio, 116 henrysConstant, 117 env_p, 118 env_T, 119 env_rho_d, 120 timestep, 121 dv, 122 droplet_volume, 123 multiplicity, 124 system_type, 125 specific_gravity, 126 alpha, 127 diffusion_const, 128 dissociation_factor, 129 radius, 130 const, 131 ): 132 mole_amount_taken = 0 133 for i in super_droplet_ids: 134 Mc = specific_gravity * const.Md 135 Rc = const.R_str / Mc 136 cinf = env_p / env_T / (const.Rd / env_mixing_ratio[0] + Rc) / Mc 137 r_w = radius(volume=droplet_volume[i]) 138 v_avg = np.sqrt(8 * const.R_str * env_T / (np.pi * Mc)) 139 dt_over_scale = timestep / ( 140 4 * r_w / (3 * v_avg * alpha) + r_w**2 / (3 * diffusion_const) 141 ) 142 A_old = mole_amounts[i] / droplet_volume[i] 143 H_eff = henrysConstant * dissociation_factor[i] 144 A_new = (A_old + dt_over_scale * cinf) / ( 145 1 + dt_over_scale / H_eff / const.R_str / env_T 146 ) 147 new_mole_amount_per_real_droplet = A_new * droplet_volume[i] 148 assert new_mole_amount_per_real_droplet >= 0 149 150 mole_amount_taken += multiplicity[i] * ( 151 new_mole_amount_per_real_droplet - mole_amounts[i] 152 ) 153 mole_amounts[i] = new_mole_amount_per_real_droplet 154 delta_mr = mole_amount_taken * specific_gravity * const.Md / (dv * env_rho_d) 155 assert delta_mr <= env_mixing_ratio 156 if system_type == "closed": 157 env_mixing_ratio -= delta_mr
def
oxidation( self, *, n_sd, cell_ids, do_chemistry_flag, k0, k1, k2, k3, K_SO2, K_HSO3, timestep, droplet_volume, pH, dissociation_factor_SO2, moles_O3, moles_H2O2, moles_S_IV, moles_S_VI):
159 def oxidation( # pylint: disable=too-many-locals 160 self, 161 *, 162 n_sd, 163 cell_ids, 164 do_chemistry_flag, 165 k0, 166 k1, 167 k2, 168 k3, 169 K_SO2, 170 K_HSO3, 171 timestep, 172 droplet_volume, 173 pH, 174 dissociation_factor_SO2, 175 # output 176 moles_O3, 177 moles_H2O2, 178 moles_S_IV, 179 moles_S_VI, 180 ): 181 ChemistryMethods.oxidation_body( 182 n_sd=n_sd, 183 cell_ids=cell_ids.data, 184 do_chemistry_flag=do_chemistry_flag.data, 185 explicit_euler=self.formulae.trivia.explicit_euler, 186 pH2H=self.formulae.trivia.pH2H, 187 k0=k0.data, 188 k1=k1.data, 189 k2=k2.data, 190 k3=k3.data, 191 K_SO2=K_SO2.data, 192 K_HSO3=K_HSO3.data, 193 timestep=timestep, 194 droplet_volume=droplet_volume.data, 195 pH=pH.data, 196 dissociation_factor_SO2=dissociation_factor_SO2.data, 197 # output 198 moles_O3=moles_O3.data, 199 moles_H2O2=moles_H2O2.data, 200 moles_S_IV=moles_S_IV.data, 201 moles_S_VI=moles_S_VI.data, 202 )
@staticmethod
@numba.njit(**conf.JIT_FLAGS)
def
oxidation_body( *, n_sd, cell_ids, do_chemistry_flag, explicit_euler, pH2H, k0, k1, k2, k3, K_SO2, K_HSO3, timestep, droplet_volume, pH, dissociation_factor_SO2, moles_O3, moles_H2O2, moles_S_IV, moles_S_VI):
204 @staticmethod 205 @numba.njit(**conf.JIT_FLAGS) 206 def oxidation_body( # pylint: disable=too-many-locals 207 *, 208 n_sd, 209 cell_ids, 210 do_chemistry_flag, 211 explicit_euler, 212 pH2H, 213 k0, 214 k1, 215 k2, 216 k3, 217 K_SO2, 218 K_HSO3, 219 timestep, 220 droplet_volume, 221 pH, 222 dissociation_factor_SO2, 223 # output 224 moles_O3, 225 moles_H2O2, 226 moles_S_IV, 227 moles_S_VI, 228 ): 229 for i in numba.prange(n_sd): # pylint: disable=not-an-iterable 230 if not do_chemistry_flag[i]: 231 continue 232 233 cid = cell_ids[i] 234 H = pH2H(pH[i]) 235 SO2aq = moles_S_IV[i] / droplet_volume[i] / dissociation_factor_SO2[i] 236 237 # NB: This might not be entirely correct 238 # https://doi.org/10.1029/JD092iD04p04171 239 # https://doi.org/10.5194/acp-16-1693-2016 240 241 ozone = ( 242 ( 243 k0[cid] 244 + (k1[cid] * K_SO2[cid] / H) 245 + (k2[cid] * K_SO2[cid] * K_HSO3[cid] / H**2) 246 ) 247 * (moles_O3[i] / droplet_volume[i]) 248 * SO2aq 249 ) 250 peroxide = ( 251 k3[cid] 252 * K_SO2[cid] 253 / (1 + k4 * H) 254 * (moles_H2O2[i] / droplet_volume[i]) 255 * SO2aq 256 ) 257 dt_times_volume = timestep * droplet_volume[i] 258 259 dconc_dt_O3 = -ozone 260 dconc_dt_S_IV = -(ozone + peroxide) 261 dconc_dt_H2O2 = -peroxide 262 dconc_dt_S_VI = ozone + peroxide 263 264 if ( 265 moles_O3[i] + dconc_dt_O3 * dt_times_volume < 0 266 or moles_S_IV[i] + dconc_dt_S_IV * dt_times_volume < 0 267 or moles_S_VI[i] + dconc_dt_S_VI * dt_times_volume < 0 268 or moles_H2O2[i] + dconc_dt_H2O2 * dt_times_volume < 0 269 ): 270 continue 271 272 moles_O3[i] = explicit_euler(moles_O3[i], dt_times_volume, dconc_dt_O3) 273 moles_S_IV[i] = explicit_euler( 274 moles_S_IV[i], dt_times_volume, dconc_dt_S_IV 275 ) 276 moles_S_VI[i] = explicit_euler( 277 moles_S_VI[i], dt_times_volume, dconc_dt_S_VI 278 ) 279 moles_H2O2[i] = explicit_euler( 280 moles_H2O2[i], dt_times_volume, dconc_dt_H2O2 281 )
def
chem_recalculate_drop_data(self, dissociation_factors, equilibrium_consts, cell_id, pH):
283 def chem_recalculate_drop_data( 284 self, dissociation_factors, equilibrium_consts, cell_id, pH 285 ): 286 for i in range(len(pH)): 287 H = self.formulae.trivia.pH2H(pH.data[i]) 288 for key in DIFFUSION_CONST: 289 dissociation_factors[key].data[i] = DISSOCIATION_FACTORS[key]( 290 H, equilibrium_consts, cell_id.data[i] 291 )
def
chem_recalculate_cell_data(self, equilibrium_consts, kinetic_consts, temperature):
293 def chem_recalculate_cell_data( 294 self, equilibrium_consts, kinetic_consts, temperature 295 ): 296 for i in range(len(temperature)): 297 for key in equilibrium_consts: 298 equilibrium_consts[key].data[i] = ( 299 self.EQUILIBRIUM_CONST.EQUILIBRIUM_CONST[key].at( 300 temperature.data[i] 301 ) 302 ) 303 for key in kinetic_consts: 304 kinetic_consts[key].data[i] = self.KINETIC_CONST.KINETIC_CONST[key].at( 305 temperature.data[i] 306 )
def
equilibrate_H( self, *, equilibrium_consts, cell_id, conc, do_chemistry_flag, pH, H_min, H_max, ionic_strength_threshold, rtol):
308 def equilibrate_H( 309 self, 310 *, 311 equilibrium_consts, 312 cell_id, 313 conc, 314 do_chemistry_flag, 315 pH, 316 H_min, 317 H_max, 318 ionic_strength_threshold, 319 rtol, 320 ): 321 ChemistryMethods.equilibrate_H_body( 322 within_tolerance=self.formulae.trivia.within_tolerance, 323 pH2H=self.formulae.trivia.pH2H, 324 H2pH=self.formulae.trivia.H2pH, 325 cell_id=cell_id.data, 326 conc=_conc( 327 N_mIII=conc.N_mIII.data, 328 N_V=conc.N_V.data, 329 C_IV=conc.C_IV.data, 330 S_IV=conc.S_IV.data, 331 S_VI=conc.S_VI.data, 332 ), 333 K=_K( 334 NH3=equilibrium_consts["K_NH3"].data, 335 SO2=equilibrium_consts["K_SO2"].data, 336 HSO3=equilibrium_consts["K_HSO3"].data, 337 HSO4=equilibrium_consts["K_HSO4"].data, 338 HCO3=equilibrium_consts["K_HCO3"].data, 339 CO2=equilibrium_consts["K_CO2"].data, 340 HNO3=equilibrium_consts["K_HNO3"].data, 341 ), 342 # output 343 do_chemistry_flag=do_chemistry_flag.data, 344 pH=pH.data, 345 # params 346 H_min=H_min, 347 H_max=H_max, 348 ionic_strength_threshold=ionic_strength_threshold, 349 rtol=rtol, 350 )
@staticmethod
@numba.njit(**{**conf.JIT_FLAGS, **{'parallel': False, 'cache': False}})
def
equilibrate_H_body( within_tolerance, pH2H, H2pH, cell_id, conc, K, do_chemistry_flag, pH, H_min, H_max, ionic_strength_threshold, rtol):
352 @staticmethod 353 @numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False, "cache": False}}) 354 def equilibrate_H_body( # pylint: disable=too-many-arguments,too-many-locals 355 within_tolerance, 356 pH2H, 357 H2pH, 358 cell_id, 359 conc, 360 K, 361 do_chemistry_flag, 362 pH, 363 # params 364 H_min, 365 H_max, 366 ionic_strength_threshold, 367 rtol, 368 ): 369 for i, pH_i in enumerate(pH): 370 cid = cell_id[i] 371 args = ( 372 _conc( 373 N_mIII=conc.N_mIII[i], 374 N_V=conc.N_V[i], 375 C_IV=conc.C_IV[i], 376 S_IV=conc.S_IV[i], 377 S_VI=conc.S_VI[i], 378 ), 379 _K( 380 NH3=K.NH3[cid], 381 SO2=K.SO2[cid], 382 HSO3=K.HSO3[cid], 383 HSO4=K.HSO4[cid], 384 HCO3=K.HCO3[cid], 385 CO2=K.CO2[cid], 386 HNO3=K.HNO3[cid], 387 ), 388 ) 389 a = pH2H(pH_i) 390 fa = acidity_minfun(a, *args) 391 if abs(fa) < _REALY_CLOSE_THRESHOLD: 392 continue 393 b = np.nan 394 fb = np.nan 395 use_default_range = False 396 if abs(fa) < _QUITE_CLOSE_THRESHOLD: 397 b = a * _QUITE_CLOSE_MULTIPLIER 398 fb = acidity_minfun(b, *args) 399 if fa * fb > 0: 400 b = a 401 fb = fa 402 a = b / _QUITE_CLOSE_MULTIPLIER / _QUITE_CLOSE_MULTIPLIER 403 fa = acidity_minfun(a, *args) 404 if fa * fb > 0: 405 use_default_range = True 406 else: 407 use_default_range = True 408 if use_default_range: 409 a = H_min 410 b = H_max 411 fa = acidity_minfun(a, *args) 412 fb = acidity_minfun(b, *args) 413 max_iter = _MAX_ITER_DEFAULT 414 else: 415 max_iter = _MAX_ITER_QUITE_CLOSE 416 H, _iters_taken = toms748_solve( 417 acidity_minfun, 418 args, 419 a, 420 b, 421 fa, 422 fb, 423 rtol=rtol, 424 max_iter=max_iter, 425 within_tolerance=within_tolerance, 426 ) 427 assert _iters_taken != max_iter 428 pH[i] = H2pH(H) 429 ionic_strength = calc_ionic_strength(H, *args) 430 do_chemistry_flag[i] = ionic_strength <= ionic_strength_threshold
@numba.njit(**{**conf.JIT_FLAGS, **{'parallel': False}})
def
calc_ionic_strength(H, conc, K):
433@numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False}}) 434def calc_ionic_strength(H, conc, K): 435 # Directly adapted 436 # https://github.com/igfuw/libcloudphxx/blob/0b4e2455fba4f95c7387623fc21481a85e7b151f/src/impl/particles_impl_chem_strength.ipp#L50 437 # https://en.wikipedia.org/wiki/Ionic_strength 438 439 # H+ and OH- 440 water = H + K_H2O / H 441 442 # HSO4- and SO4 2- 443 cz_S_VI = H * conc.S_VI / (H + K.HSO4) + 4 * K.HSO4 * conc.S_VI / (H + K.HSO4) 444 445 # HCO3- and CO3 2- 446 cz_CO2 = K.CO2 * H * conc.C_IV / ( 447 H * H + K.CO2 * H + K.CO2 * K.HCO3 448 ) + 4 * K.CO2 * K.HCO3 * conc.C_IV / (H * H + K.CO2 * H + K.CO2 * K.HCO3) 449 450 # HSO3- and HSO4 2- 451 cz_SO2 = K.SO2 * H * conc.S_IV / ( 452 H * H + K.SO2 * H + K.SO2 * K.HSO3 453 ) + 4 * K.SO2 * K.HSO3 * conc.S_IV / (H * H + K.SO2 * H + K.SO2 * K.HSO3) 454 455 # NO3- 456 cz_HNO3 = K.HNO3 * conc.N_V / (H + K.HNO3) 457 458 # NH4+ 459 cz_NH3 = K.NH3 * H * conc.N_mIII / (K_H2O + K.NH3 * H) 460 461 return 0.5 * (water + cz_S_VI + cz_CO2 + cz_SO2 + cz_HNO3 + cz_NH3)
@numba.njit(**{**conf.JIT_FLAGS, **{'parallel': False}})
def
acidity_minfun(H, conc, K):
464@numba.njit(**{**conf.JIT_FLAGS, **{"parallel": False}}) 465def acidity_minfun(H, conc, K): 466 ammonia = (conc.N_mIII * H * K.NH3) / (K_H2O + K.NH3 * H) 467 nitric = conc.N_V * K.HNO3 / (H + K.HNO3) 468 sulfous = ( 469 conc.S_IV * K.SO2 * (H + 2 * K.HSO3) / (H * H + H * K.SO2 + K.SO2 * K.HSO3) 470 ) 471 water = K_H2O / H 472 sulfuric = conc.S_VI * (H + 2 * K.HSO4) / (H + K.HSO4) 473 carbonic = ( 474 conc.C_IV * K.CO2 * (H + 2 * K.HCO3) / (H * H + H * K.CO2 + K.CO2 * K.HCO3) 475 ) 476 zero = H + ammonia - (nitric + sulfous + water + sulfuric + carbonic) 477 return zero