Evaporation scheme described in the paper #1023
Conversation
|
| # it instead allocated and returned the field. | ||
| # This needs to be updated everywhere, and then these functions | ||
| # can become e.g. surface_specific_humidity!(p.soil.q_sfc, ...) | ||
| ClimaLand.surface_specific_humidity(model, Y, p, T_sfc, ρ_sfc) |
There was a problem hiding this comment.
This also updates p.soil.ice_frac in place, and that seems confusing. I will move this into a different function and make it obvious (but am going to let the current long run keep going)
| @@ -960,18 +955,8 @@ function soil_turbulent_fluxes_at_a_point( | |||
| q_air::FT = | |||
There was a problem hiding this comment.
quite a large diff - may be better to look at file directly?
| porous media, Geophys. Res. Lett., 35, L19407, doi:10.1029/ | ||
| 2008GL035230. | ||
| """ | ||
| function soil_tortuosity(θ_l::FT, θ_i::FT, ν::FT) where {FT} |
There was a problem hiding this comment.
We had previously implement a lot of these equations, commented them out, and this adds them back to soil_hydrology_parameterizations.
|
Pre-turbulent fluxes ClimaLand.jl/src/standalone/Soil/energy_hydrology.jl Lines 843 to 846 in fcf1385 As noted, the call to q_soil also updates ice_frac (I will change this to be explicit): This uses Equation 74, q_soiil, Equation 73: ClimaLand.jl/src/standalone/Soil/energy_hydrology.jl Lines 1116 to 1130 in fcf1385 Within turbulent fluxes: Equation 65: Equation 66 Equation 67 Both have q_0 = q_soil, the latter has r_0 = r_soil. Equations 75/76: ClimaLand.jl/src/standalone/Soil/energy_hydrology.jl Lines 994 to 995 in fcf1385 Soil resistance Equation 79: Branch of 79: ClimaLand.jl/src/standalone/Soil/energy_hydrology.jl Lines 990 to 992 in fcf1385 Equation 82: Equation 80, with a worse name: Equation 81 + top branch of 79: |
| vg_n = FT(8.91)#FT(6.34) optimized values by root finding | ||
| vg_α = FT(3) #FT(7.339) optimized values by root finding |
There was a problem hiding this comment.
delete commented out values
| # it instead allocated and returned the field. | ||
| # This needs to be updated everywhere, and then these functions | ||
| # can become e.g. surface_specific_humidity!(p.soil.q_sfc, ...) | ||
| ClimaLand.surface_specific_humidity(model, Y, p, T_sfc, ρ_sfc) |
There was a problem hiding this comment.
| ClimaLand.surface_specific_humidity(model, Y, p, T_sfc, ρ_sfc) | |
| ClimaLand.surface_specific_humidity!(model, Y, p, T_sfc, ρ_sfc) |
This function should have a ! since it's modifying p in-place
| # can become e.g. surface_specific_humidity!(p.soil.q_sfc, ...) | ||
| ClimaLand.surface_specific_humidity(model, Y, p, T_sfc, ρ_sfc) | ||
| q_sfc = p.soil.q_sfc | ||
| ClimaLand.surface_resistance(model, Y, p, t) |
There was a problem hiding this comment.
| ClimaLand.surface_resistance(model, Y, p, t) | |
| ClimaLand.surface_resistance!(model, Y, p, t) |
here too (modifies p.soil.r_sfc)
| @@ -822,7 +822,6 @@ function ClimaLand.get_drivers(model::EnergyHydrology) | |||
| end | |||
|
|
|||
|
|
|||
|
|
|||
| function turbulent_fluxes!( | |||
There was a problem hiding this comment.
could you add a docstring for this function please?
| return S_w < S_c ? d_ds * (S_c - S_w) / S_c : FT(0) | ||
| ice_frac = p.soil.ice_frac | ||
| @. ice_frac = ice_fraction(θ_l_sfc, θ_i_sfc, ν_sfc, θ_r_sfc, _ρ_l, _ρ_i) | ||
| @. p.soil.q_sfc = |
There was a problem hiding this comment.
this function should be surface_specific_humidity because it's modifying p
| ) where {FT <: AbstractFloat, C, EP} | ||
|
|
||
| Computes turbulent surface fluxes for soil at a point on a surface given | ||
| (1) Surface state conditions (`T_sfc`, `θ_l_sfc`, `θ_i_sfc`) | ||
| (1) Surface state conditions (`T_sfc`, `q_sfc`, `ρ_sfc`, `ice_frac`) | ||
| (2) Surface properties, such as the topographical height of the surface (`h_sfc`), |
| ice_frac * Thermodynamics.q_vap_saturation_generic( | ||
| thermo_params, | ||
| T_sfc, | ||
| ρ_sfc, | ||
| Thermodynamics.Ice(), |
There was a problem hiding this comment.
Minor thing, but you could collect Thermodynamics.q_vap_saturation_generic in this expression to compute it only once
| Computes the resistance of the top of the soil column to | ||
| water vapor diffusion, as a function of the surface | ||
| volumetric liquid water fraction `θ_l`, the augmented | ||
| liquid water fraction `ϑ_l`, the volumetric ice water |
There was a problem hiding this comment.
It would be helpful to make sure all the function parameters are clearly identified in the docstring. For example, I don't understand how ϑ_l enters here, or what d_ds is.
| dsl::FT = dry_soil_layer_thickness(S_w, S_c, d_ds) | ||
| r_pore::FT = 2 * _σ * α / _ρ_liq / _grav | ||
| θ_safe = max(eps(FT), (θ_i + θ_l - θ_r)) | ||
| r_shell::FT = r_pore / _D_vapor / (4 * θ_safe) * (π - 2 * (θ_safe)^(1 / 2)) |
There was a problem hiding this comment.
| r_shell::FT = r_pore / _D_vapor / (4 * θ_safe) * (π - 2 * (θ_safe)^(1 / 2)) | |
| r_shell::FT = r_pore / _D_vapor / (4 * θ_safe) * (π - 2 * sqrt(θ_safe)) |
Floating point exponents can be much more computationally expensive
| θ_safe = max(eps(FT), (θ_i + θ_l - θ_r)) | ||
| r_shell::FT = r_pore / _D_vapor / (4 * θ_safe) * (π - 2 * (θ_safe)^(1 / 2)) | ||
| # This factor just damps r_shell to zero more quickly as x -> 0. numerical. | ||
| x = θ_safe / FT(0.001) |
There was a problem hiding this comment.
| x = θ_safe / FT(0.001) | |
| x = θ_safe / 1_000 |
potentially avoids a conversion and better preserves accuracy
| """ | ||
| dry_soil_layer_thickness(S_w::FT, S_c::FT, d_ds::FT)::FT where {FT} | ||
|
|
||
| Returns the maximum dry soil layer thickness that can develop under vapor flux; |
There was a problem hiding this comment.
It would be helpful to clearly identify the physical meaning of the input parameters
| """ | ||
| ice_fraction(θ_l::FT, θ_i::FT, ν::FT, θ_r::FT, ρ_l::FT, ρ_i::FT)::FT | ||
| where {FT} | ||
|
|
||
| Computes and returns the | ||
| fraction of the humidity in the pore space | ||
| that is due to sublimation. | ||
|
|
||
| f = θ_iρ_i/(θ_iρ_i+(θ_l-θ_r)ρ_l) | ||
| """ |
There was a problem hiding this comment.
"""
ice_fraction(θ_l::FT, θ_i::FT, ν::FT, θ_r::FT, ρ_l::FT, ρ_i::FT)::FT where {FT}
Return the fraction of the humidity in the pore space that is due to sublimation.
"""
There was a problem hiding this comment.
This doesn't render well on github, but you can add a math block so that it renders as equation int he docs
There was a problem hiding this comment.
```math
f = \frac{θ_i ρ_i}{θ_i ρ_i + (θ_l - θ_r) ρ_l}
```
| """ | ||
| ice_fraction(θ_l::FT, θ_i::FT, ν::FT, θ_r::FT, ρ_l::FT, ρ_i::FT)::FT | ||
| where {FT} | ||
|
|
||
| Computes and returns the | ||
| fraction of the humidity in the pore space | ||
| that is due to sublimation. | ||
|
|
||
| f = θ_iρ_i/(θ_iρ_i+(θ_l-θ_r)ρ_l) | ||
| """ |
There was a problem hiding this comment.
This doesn't render well on github, but you can add a math block so that it renders as equation int he docs
| """ | ||
| ice_fraction(θ_l::FT, θ_i::FT, ν::FT, θ_r::FT, ρ_l::FT, ρ_i::FT)::FT | ||
| where {FT} | ||
|
|
||
| Computes and returns the | ||
| fraction of the humidity in the pore space | ||
| that is due to sublimation. | ||
|
|
||
| f = θ_iρ_i/(θ_iρ_i+(θ_l-θ_r)ρ_l) | ||
| """ |
There was a problem hiding this comment.
```math
f = \frac{θ_i ρ_i}{θ_i ρ_i + (θ_l - θ_r) ρ_l}
```
| dsl::FT = dry_soil_layer_thickness(S_w, S_c, d_ds) | ||
| r_pore::FT = 2 * _σ * α / _ρ_liq / _grav | ||
| θ_safe = max(eps(FT), (θ_i + θ_l - θ_r)) | ||
| r_shell::FT = r_pore / _D_vapor / (4 * θ_safe) * (π - 2 * (θ_safe)^(1 / 2)) |
There was a problem hiding this comment.
| r_shell::FT = r_pore / _D_vapor / (4 * θ_safe) * (π - 2 * (θ_safe)^(1 / 2)) | |
| r_shell::FT = r_pore / (4 * _D_vapor * θ_safe) * (π - 2 * (θ_safe)^(1 / 2)) |
| S_w = effective_saturation(ν, θ_l + θ_i, θ_r) | ||
| τ_a = soil_tortuosity(θ_l, θ_i, ν) | ||
| dsl::FT = dry_soil_layer_thickness(S_w, S_c, d_ds) | ||
| r_pore::FT = 2 * _σ * α / _ρ_liq / _grav |
There was a problem hiding this comment.
| r_pore::FT = 2 * _σ * α / _ρ_liq / _grav | |
| r_pore::FT = 2 * _σ * α / (_ρ_liq * _grav) |
| _D_vapor = LP.D_vapor(earth_param_set) | ||
| _ρ_liq = LP.ρ_cloud_liq(earth_param_set) | ||
| _grav = LP.grav(earth_param_set) | ||
| _σ = FT(7.2e-2) # need to add to CP, surface tension of water N/m |
There was a problem hiding this comment.
Maybe call this _σ if there's already a Stafan-Boltzmann constant?
| """ | ||
| function soil_tortuosity(θ_l::FT, θ_i::FT, ν::FT) where {FT} | ||
| safe_θ_a = max(ν - θ_l - θ_i, eps(FT)) | ||
| return safe_θ_a^(FT(2.5)) / ν |
There was a problem hiding this comment.
| return safe_θ_a^(FT(2.5)) / ν | |
| return sqrt(safe_θ_a^5) / ν |
Purpose
Adds the evaporation scheme from the paper.
Equations
in turbulent fluxes, surface specific humidity, and ice fraction functions:
soil resistance:
Comparison to data
Plot produced by docs/tutorials/standalone/Soil/evaporation.jl
To-do
Look at long runs