diff --git a/Docs/sphinx_doc/index.rst b/Docs/sphinx_doc/index.rst index 62e87c7ce..5e57dfa3d 100644 --- a/Docs/sphinx_doc/index.rst +++ b/Docs/sphinx_doc/index.rst @@ -49,6 +49,7 @@ In addition to this documentation, there is API documentation for ERF generated theory/DryEquations.rst theory/WetEquations.rst + theory/Initialization.rst theory/Buoyancy.rst theory/Microphysics.rst theory/DNSvsLES.rst diff --git a/Docs/sphinx_doc/theory/Buoyancy.rst b/Docs/sphinx_doc/theory/Buoyancy.rst index 6213c379f..dbf289ac2 100644 --- a/Docs/sphinx_doc/theory/Buoyancy.rst +++ b/Docs/sphinx_doc/theory/Buoyancy.rst @@ -109,3 +109,60 @@ which removes the need to compute horizontal averages of these quantities. Thi We note that this version of the buoyancy force matches that given in Marat F. Khairoutdinov and David A. Randall's paper (J. Atm Sciences, 607, 1983) if we neglect :math:`\frac{p^\prime}{\bar{p_0}}`. + +Type 4 +------ +This expression for buoyancy is from `khairoutdinov2003cloud`_ and `bryan2002benchmark`_. + +.. _`khairoutdinov2003cloud`: https://journals.ametsoc.org/view/journals/atsc/60/4/1520-0469_2003_060_0607_crmota_2.0.co_2.xml +.. _`bryan2002benchmark`: https://journals.ametsoc.org/view/journals/mwre/130/12/1520-0493_2002_130_2917_absfmn_2.0.co_2.xml + +.. math:: + + \begin{equation} + \mathbf{B} = \rho'\mathbf{g} \approx -\rho\Bigg(\frac{T'}{T} + 0.61 q_v' - q_c - q_p - \frac{p'}{p}\Bigg), + \end{equation} + +The derivation follows. The total density is given by :math:`\rho = \rho_d(1 + q_v + q_c + q_p)`, which can be written as + +.. math:: + + \rho = \frac{p (1 + q_v + q_c + q_p)}{R_dT\Bigg(1 + \cfrac{R_v}{R_d}q_v\Bigg)} + +This can be written using binomial expansion as + +.. math:: + + \begin{align*} + \rho &= \frac{p}{R_dT} (1 + q_v + q_c + q_p)\Bigg(1 + \frac{R_v}{R_d}q_v\Bigg)^{-1} \\ + &= \frac{p}{R_dT} (1 + q_v + q_c + q_p)\Bigg(1 - \frac{R_v}{R_d}q_v + O(q_v^2)\Bigg) \\ + &= \frac{p}{R_dT}\Bigg(1 + q_v + q_c + q_p - \frac{R_v}{R_d}q_v + \text{H.O.T. such as } O(q_v^2) + O(q_vq_c)\Bigg) \\ + &\approx \frac{p}{R_dT}\Bigg(1 + q_v + q_c + q_p - \frac{R_v}{R_d}q_v\Bigg) + \end{align*} + +Taking log on both sides, we get + +.. math:: + + \log{\rho} = \log{p} - \log{R_d} - \log{T} + \log(1 - 0.61 q_v + q_c + q_p) + +Taking derivative gives + +.. math:: + + \frac{\rho'}{\rho} = \frac{p'}{p} - \frac{T'}{T} + \frac{(-0.61 q_v' + q_c' + q_p')}{(1 - 0.61 q_v + q_c + q_p)} + +Using :math:`- 0.61 q_v + q_c + q_p \ll 1`, we have + +.. math:: + + \frac{\rho'}{\rho} = \frac{p'}{p} - \frac{T'}{T} + (-0.61 q_v' + q_c' + q_p') + +Since the background values of cloud water and precipitate mass mixing ratios -- :math:`q_c` and :math:`q_p` are zero, we have :math:`q_c' = q_c` and :math:`q_p' = q_p`. Hence, we have + +.. math:: + + \begin{equation} + \rho'\approx -\rho\Bigg(\frac{T'}{T} + 0.61 q_v' - q_c - q_p - \frac{p'}{p}\Bigg), + \end{equation} + diff --git a/Docs/sphinx_doc/theory/Initialization.rst b/Docs/sphinx_doc/theory/Initialization.rst new file mode 100644 index 000000000..10a7e1a50 --- /dev/null +++ b/Docs/sphinx_doc/theory/Initialization.rst @@ -0,0 +1,149 @@ + + .. role:: cpp(code) + :language: c++ + + .. role:: f(code) + :language: fortran + +.. _Initialization: + +Initialization +================== + +To initialize a background (base) state for a simulation with a hydrostatic atmosphere, the hydrostatic equation balancing the pressure gradient +and gravity must be satisfied. This section details the procedure for the initialization of the background state. The procedure is similar for the +cases with and without moisture, the only difference being that the density for the cases with moisture has to be the total density +:math:`\rho = \rho_d(1 + q_t)`, where :math:`\rho_d` is the dry density, and :math:`q_t` is the total mass mixing ratio -- water vapor and liquid water, instead +of the dry density :math:`\rho_d` for cases without moisture. + +Computation of the dry density +------------------------------- +We express the total pressure :math:`p` as: + +.. math:: + p = \rho_d R_d T_v. + +By definition, we have: + +.. math:: + T_v = \theta_m\left(\frac{p}{p_0}\right)^{R_d/C_p}, + +.. math:: + T = \theta_d\left(\frac{p}{p_0}\right)^{R_d/C_p}. + +This gives: + +.. math:: + p = \rho_d R_d \theta_m\left(\frac{p}{p_0}\right)^{R_d/C_p}, + +which, on rearranging, yields: + +.. math:: + p = p_0\left(\frac{\rho_d R_d \theta_m}{p_0}\right)^\gamma. + +To obtain :math:`\theta_m`, consider the density of dry air: + +.. math:: + \rho_d = \frac{p - p_v}{R_d T}. + +Substituting for :math:`\rho_d` from the above equation, we get: + +.. math:: + \frac{p}{T_v} = \frac{p - p_v}{T}, + +which implies: + +.. math:: + \frac{T_v}{T} = \frac{p}{p - p_v} = \frac{1}{1-\left(\cfrac{p_v}{p}\right)}. + +We also have: + +.. math:: + q_v = \frac{\rho_v}{\rho_d} = \frac{p_v}{R_v T}\frac{R_d T}{p-p_v} = \frac{r p_v}{p - p_v}, + +where :math:`p_v` is the partial pressure of water vapor, :math:`r = R_d/R_v \approx 0.622`, and :math:`q_v` is the vapor mass mixing ratio—the ratio of the +mass of vapor to the mass of dry air. Rearranging and using :math:`q_v \ll r`, we get: + +.. math:: + \frac{p_v}{p} = \frac{1}{1 + \left(\cfrac{r}{q_v}\right)} \approx \frac{q_v}{r}, + +which, on substitution in the equation for :math:`\frac{T_v}{T}`, gives: + +.. math:: + \frac{T_v}{T} = \frac{1}{1 - \left(\cfrac{q_v}{r}\right)}. + +As :math:`q_v \ll 1`, a binomial expansion, ignoring higher-order terms, gives: + +.. math:: + T_v = T\left(1 + \frac{R_v}{R_d}q_v\right). + +Hence, the density of dry air is given by: + +.. math:: + \rho_d = \frac{p}{R_d T_v} = \frac{p}{R_d T\left(1 + \cfrac{R_v}{R_d}q_v\right)}. + + +Initialization with a second-order integration of the hydrostatic equation +---------------------------------------------------------------------------- + +We have the hydrostatic equation given by + +.. math:: + + \frac{\partial p}{\partial z} = -\rho g, + +where :math:`\rho = \rho_d(1 + q_t)`, :math:`\rho_d` is the dry density, and :math:`q_t` is the total mass mixing ratio -- water vapor and liquid water. Using an average value of :math:`\rho` for the integration from :math:`z = z(k-1)` to :math:`z(k)`, we get + +.. math:: + + p(k) = p(k-1) - \frac{(\rho(k-1) + \rho(k))}{2} g\Delta z. + +The density at a point is a function of the pressure, potential temperature, and relative humidity. The latter two quantities are computed using user-specified profiles, and hence, for simplicity, we write :math:`\rho(k) = f(p(k))`. Hence + +.. math:: + + p(k) = p(k-1) - \frac{\rho(k-1)}{2}g\Delta z - \frac{f(p(k))}{2}g\Delta z. + +Now, we define + +.. math:: + + F(p(k)) \equiv p(k) - p(k-1) + \frac{\rho(k-1)}{2}g\Delta z + \frac{f(p(k))}{2}g\Delta z = 0. + +This is a non-linear equation in :math:`p(k)`. Consider a Newton-Raphson iteration (where :math:`n` denotes the iteration number) procedure + +.. math:: + + F(p+\delta p) \approx F(p) + \delta p \frac{\partial F}{\partial p} = 0, + +which implies + +.. math:: + + \delta p = -\frac{F}{F'}, + +with the gradient being evaluated as + +.. math:: + + F' = \frac{F(p+\epsilon) - F(p)}{\epsilon}, + +and the iteration update is given by + +.. math:: + + p^{n+1} = p^n + \delta p. + +For the first cell (:math:`k=0`), which is at a height of :math:`z = \frac{\Delta z}{2}` from the base, we have + +.. math:: + + p(0) = p_0 - \rho(0)g\frac{\Delta z}{2}, + +where :math:`p_0 = 1e5 \, \text{N/m}^2` is the pressure at the base. Hence, we define + +.. math:: + + F(p(0)) \equiv p(0) - p_0 + \rho(0)g\frac{\Delta z}{2}, + +and the Newton-Raphson procedure is the same. diff --git a/Docs/sphinx_doc/theory/Microphysics.rst b/Docs/sphinx_doc/theory/Microphysics.rst index a644d86d4..05607a0b4 100644 --- a/Docs/sphinx_doc/theory/Microphysics.rst +++ b/Docs/sphinx_doc/theory/Microphysics.rst @@ -7,8 +7,12 @@ .. _Microphysics: +Microphysics model +==================== + Kessler Microphysics model -=========================== +--------------------------- +The Kessler microphysics model is a simple version of cloud microphysics which has precipitation only in the form of rain. Hence :math:`q_p = q_r`. Governing equations for the microphysical quantities for Kessler microphysics from `gabervsek2012dry`_ are .. math:: @@ -24,14 +28,78 @@ Governing equations for the microphysical quantities for Kessler microphysics fr where :math:`C_c` is the rate of condensation of water vapor to cloud water, :math:`E_c` is the rate of evaporation of cloud water to water vapor, :math:`A_c` is the autoconversion of cloud water to rain, :math:`K_c` is the accretion of cloud water to rain drops, :math:`E_r` is the evaporation of rain to water vapor and :math:`F_r` is the sedimentation of rain. The parametrization used is given in `klemp1978simulation`_, and is given -below. Note that in all the equations, :math:`p` is specified in millibars and :math:`\overline{\rho}` is specified in g cm :math:`^{-3}`. +below. Note that in all the equations, :math:`p` is specified in millibars and :math:`\overline{\rho}` is specified in g cm :math:`^{-3}`. The parametrization +of the source terms are given below. .. _`gabervsek2012dry`: https://journals.ametsoc.org/view/journals/mwre/140/10/mwr-d-11-00144.1.xml + +.. raw:: latex + + \newgeometry{left=2cm,right=2cm,top=2cm,bottom=2cm} + +Rate of condensation of water vapor/evaporation of cloud water +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +From `klemp1978simulation`_, we have the following expressions. + .. _`klemp1978simulation`: https://journals.ametsoc.org/view/journals/atsc/35/6/1520-0469_1978_035_1070_tsotdc_2_0_co_2.xml +If the air is not saturated, i.e. :math:`q_v > q_{vs}` + +.. math:: + C_c = \frac{q_v - q_{vs}}{1 + \cfrac{q_{vs}^*4093L}{C_p(T-36)^2}} + +If the air is not saturated, i.e. :math:`q_v < q_{vs}`, then cloud water evaporates to water vapor at the rate + +.. math:: + E_c = \frac{q_{vs} - q_v}{1 + \cfrac{q_{vs}^*4093L}{C_p(T-36)^2}} + +Rate of autoconversion of cloud water into rain +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The rate of autoconversion of cloud water into rain is given by + +.. math:: + A_c = k_1(q_c - a) + +where :math:`k_1 = 0.001` s\ :sup:`-1`, :math:`a = 0.001` kg kg\ :sup:`-1`. + +Rate of accretion of cloud water into rain water drops +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The rate of accretion of cloud water into rain water drops is given by + +.. math:: + K_c = k_2q_cq_r^{0.875} + +where :math:`k_2= 2.2` s\ :sup:`-1`. + +The rate of evaporation of rain into water vapor +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The rate of evaporation of rain into water vapor is given by + +.. math:: + E_r = \cfrac{1}{\overline{\rho}}\cfrac{(1- q_v/q_s)C(\overline{\rho}q_r)^{0.525}}{5.4\times10^5 + 2.55\times10^6/(\overline{p}q_s)}, + +where the ventilation factor :math:`C = 1.6 + 124.9(\overline{\rho}q_r)^{0.2046}`. + +Terminal fall velocity of rain +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The terminal fall velocity of rain is given by + +.. math:: + V = 3634(\overline{\rho}q_r)^{0.1346}\Bigg(\cfrac{\overline{\rho}}{\rho_0}\Bigg)^{-\frac{1}{2}}~\text{[cm/s]} + +.. raw:: latex + + \restoregeometry + + Single Moment Microphysics Model -=================================== +---------------------------------- The conversion rates among the moist hydrometeors are parameterized assuming that .. math:: @@ -54,7 +122,7 @@ Marat F. Khairoutdinov and David A. Randall's (J. Atm Sciences, 607, 1983). The implementation of microphysics model in ERF is similar to the that in the SAM code (http://rossby.msrc.sunysb.edu/~marat/SAM.html) Accretion ------------------- +~~~~~~~~~~~~ There are several different type of accretional growth mechanisms that need to be included; these describe the interaction of water vapor and cloud water with rain water. @@ -69,21 +137,21 @@ The accretion of raindrops by accretion of cloud water is Q_{racw} = \frac{\pi E_{RW}n_{0R}\alpha q_{c}\Gamma(3+b)}{4\lambda_{R}^{3+b}}(\frac{\rho_{0}}{\rho})^{1/2} The bergeron Process ------------------------- +~~~~~~~~~~~~~~~~~~~~~~ The cloud water transform to snow by deposition and rimming can be written as .. math:: Q_{sfw} = N_{150}\left(\alpha_{1}m_{150}^{\alpha_{2}}+\pi E_{iw}\rho q_{c}R_{150}^{2}U_{150}\right) Autoconversion ------------------------- +~~~~~~~~~~~~~~~~~~~~~~ The collision and coalescence of cloud water to from raindrops is parameterized as following .. math:: Q_{raut} = \rho\left(q_{c}-q_{c0}\right)^{2}\left[1.2 \times 10^{-4}+{1.569 \times 10^{-12}N_{1}/[D_{0}(q_{c}-q_{c0})]}\right]^{-1} Evaporation ------------------------- +~~~~~~~~~~~~~~~~~~~~~~ The evaporation rate of rain is .. math:: diff --git a/Docs/sphinx_doc/theory/WetEquations.rst b/Docs/sphinx_doc/theory/WetEquations.rst index 38288e484..27e75e28f 100644 --- a/Docs/sphinx_doc/theory/WetEquations.rst +++ b/Docs/sphinx_doc/theory/WetEquations.rst @@ -138,4 +138,3 @@ In this set of equations, the subgrid turbulent parameterization effects are inc :math:`\mathbf{F}` stands for the external force, and :math:`Q` and :math:`F_Q` represent the mass and energy transformation of water vapor to/from water through condensation/evaporation, which is determined by the microphysics parameterization processes. :math:`\mathbf{B}` is the buoyancy force, which is defined in :ref:`Buoyancy `. - diff --git a/Exec/SquallLine_2D/inputs_moisture b/Exec/SquallLine_2D/inputs_moisture index 2e2757fe5..f82dd6612 100644 --- a/Exec/SquallLine_2D/inputs_moisture +++ b/Exec/SquallLine_2D/inputs_moisture @@ -70,7 +70,7 @@ erf.molec_diff_type = "ConstantAlpha" erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities erf.dynamicViscosity = 200.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s erf.alpha_T = 00.0 # [m^2/s] -erf.alpha_C = 25.0 +erf.alpha_C = 100.0 erf.moisture_model = "Kessler" erf.use_moist_background = true diff --git a/Exec/SquallLine_2D/inputs_moisture_Gabersek b/Exec/SquallLine_2D/inputs_moisture_Gabersek new file mode 100644 index 000000000..ad5b642b4 --- /dev/null +++ b/Exec/SquallLine_2D/inputs_moisture_Gabersek @@ -0,0 +1,91 @@ +# ------------------ INPUTS TO MAIN PROGRAM ------------------- +max_step = 18000 +stop_time = 90000.0 + +amrex.fpe_trap_invalid = 1 + +fabarray.mfiter_tile_size = 2048 1024 2048 + +# PROBLEM SIZE & GEOMETRY +geometry.prob_lo = -60000. 0. 0. +geometry.prob_hi = 60000. 400. 20000. +amr.n_cell = 384 4 128 # dx=dy=dz=100 m + +# periodic in x to match WRF setup +# - as an alternative, could use symmetry at x=0 and outflow at x=25600 +geometry.is_periodic = 1 1 0 +#xlo.type = "Outflow" +#xhi.type = "Outflow" +zlo.type = "SlipWall" +zhi.type = "Outflow" + +erf.sponge_strength = 2.0 +#erf.use_zhi_sponge_damping = true +erf.zhi_sponge_start = 12000.0 + +erf.sponge_density = 1.2 +erf.sponge_x_velocity = 0.0 +erf.sponge_y_velocity = 0.0 +erf.sponge_z_velocity = 0.0 + +# TIME STEP CONTROL +erf.use_native_mri = 1 +erf.fixed_dt = 0.5 # fixed time step [s] -- Straka et al 1993 +erf.fixed_fast_dt = 0.25 # fixed time step [s] -- Straka et al 1993 +#erf.no_substepping = 1 + +# DIAGNOSTICS & VERBOSITY +erf.sum_interval = 1 # timesteps between computing mass +erf.v = 1 # verbosity in ERF.cpp +amr.v = 1 # verbosity in Amr.cpp + +# REFINEMENT / REGRIDDING +amr.max_level = 0 # maximum level number allowed + +# CHECKPOINT FILES +amr.check_file = chk # root name of checkpoint file +amr.check_int = 1000 # number of timesteps between checkpoints +#amr.restart = chk09000 + +# PLOTFILES +erf.plot_file_1 = plt # root name of plotfile +erf.plot_int_1 = 100 # number of timesteps between plotfiles +erf.plot_vars_1 = density rhotheta rhoQt rhoQp x_velocity y_velocity z_velocity pressure theta temp qt qp qv qc qi scalar pert_dens + +# SOLVER CHOICE +erf.use_gravity = true +erf.buoyancy_type = 1 +erf.use_coriolis = false +erf.use_rayleigh_damping = false + +#erf.les_type = "Smagorinsky" +erf.Cs = 0.25 +erf.les_type = "None" + +# +# diffusion coefficient from Straka, K = 75 m^2/s +# +erf.molec_diff_type = "ConstantAlpha" +#erf.molec_diff_type = "Constant" +erf.rho0_trans = 1.0 # [kg/m^3], used to convert input diffusivities +erf.dynamicViscosity = 100.0 # [kg/(m-s)] ==> nu = 75.0 m^2/s +erf.alpha_T = 00.0 # [m^2/s] +erf.alpha_C = 50.0 + +erf.moisture_model = "Kessler" +erf.use_moist_background = true + +erf.moistscal_horiz_adv_string = "Centered_2nd" +erf.moistscal_vert_adv_string = "Centered_2nd" + +# PROBLEM PARAMETERS (optional) +prob.z_tr = 12000.0 +prob.height = 1200.0 +prob.theta_0 = 300.0 +prob.theta_tr = 343.0 +prob.T_tr = 213.0 +prob.x_c = 0.0 +prob.z_c = 2000.0 +prob.x_r = 10000.0 +prob.z_r = 1500.0 +prob.theta_c = 3.0 diff --git a/Source/Microphysics/Kessler/Kessler.cpp b/Source/Microphysics/Kessler/Kessler.cpp index 529bb03e9..6c32a0ff2 100644 --- a/Source/Microphysics/Kessler/Kessler.cpp +++ b/Source/Microphysics/Kessler/Kessler.cpp @@ -197,8 +197,12 @@ void Kessler::AdvanceKessler() { } + if(std::fabs(fz_array(i,j,k+1)) < 1e-14) fz_array(i,j,k+1) = 0.0; + if(std::fabs(fz_array(i,j,k)) < 1e-14) fz_array(i,j,k) = 0.0; Real dq_sed = 1.0/rho_array(i,j,k)*(fz_array(i,j,k+1) - fz_array(i,j,k))/dz*dtn; - + if(std::fabs(dq_sed) < 1e-14)dq_sed = 0.0; + //dq_sed = 0.0; + qt_array(i,j,k) = qt_array(i,j,k) + dq_rain_to_vapor - dq_clwater_to_rain; qp_array(i,j,k) = qp_array(i,j,k) + dq_sed + dq_clwater_to_rain - dq_rain_to_vapor; qn_array(i,j,k) = qn_array(i,j,k) + dq_vapor_to_clwater - dq_clwater_to_vapor - dq_clwater_to_rain;