Canoe Example#4：Juno MWR辐射传输模拟

本节解读一下juno_mwr_eq示例的具体实现和实验使用到的理论和假设。

情景描述

Juno Mission
朱诺（Juno）是一项由美国宇航局（NASA）进行的木星探测任务, 于2011年8月5日从佛罗里达州的卡纳维拉尔角空军基地发射,于2016年7月4日成功进入了木星轨道, 主要目标是深入研究木星的内部结构、大气层和磁场，以帮助科学家们了解这个巨大气体行星的形成和演化过程。
朱诺搭载了多种科学仪器，包括微波辐射计(MWR)、微波雷达、飞行器磁力计、高能粒子探测器等，用于测量木星大气的温度、气体成分、磁场和辐射环境等参数。
极化轨道：

朱诺以一个相对较小的轨道倾角环绕木星，被称为“极化轨道”，以减少朱诺受到木星辐射带的辐射影响，并允许更全面地研究木星的极区域。
Juno的轨道离心率较大，可以使其再近木点钻进木星辐射带进行大气辐射的观测，减少木星磁场的干扰。
NASA：Juno’s Orbit

NASA：Juno’s Orbit

已知木星大气特性：伽利略号探测器对木星的穿透性探测得到了一系列的probe in-situ测了结果，其中包括，大气气体成分主要是水和氨气；1 bar 高度的气温约为169 K。
模拟需要实现的步骤：

输入已知的探测结果和热力学参数；
基于绝热假设构建大气温湿度廓线和气体含量廓线；
辐射传输模拟，计算天顶微波亮温。

运行个例

构建项目和编译

bash

$ cd
(env) jihenghu@canopy:~ [11:31]$ mkdir -p Juno_MWR_Proj/build
(env) jihenghu@canopy:~ [11:31]$ cd Juno_MWR_Proj/build/
(env) jihenghu@canopy:build [11:32]$ cmake ~/canoe/ -DTASK=juno . 
(env) jihenghu@canopy:build [11:32]$ make -j8 ## 8核编译

准备所需要的文件：

bash

$ cd ..
(env) jihenghu@canopy:Juno_MWR_Proj [11:35]$ ln -s build/bin/juno_mwr.release .
(env) jihenghu@canopy:Juno_MWR_Proj [11:35]$ cp build/bin/juno_mwr.inp .
(env) jihenghu@canopy:Juno_MWR_Proj [11:36]$ cp build/bin/mwr_channels.yaml .

输入和配置文件juno_mwr.inp | mwr_channels.yaml我们会在后面介绍其具体意义。

运行

bash

(env) jihenghu@canopy:Juno_MWR_Proj [11:40]$ ./juno_mwr.release -i juno_mwr.inp 
######################################################
##              CANOE MODEL STARTS                  ##
######################################################
Log, "2023-11-15 11:40:57",        canoe, 2., "Installing monitor canoe"
Log, "2023-11-15 11:40:57",        canoe, 2.1., "Initialize IndexMap"
......
Log, "2023-11-15 11:40:59",         snap, 47.1., "Destroy Thermodynamics"
Log, "2023-11-15 11:40:59",        canoe, 48.1., "Destroy IndexMap"
Log, "2023-11-15 11:40:59",      outputs, 49.1., "Destroy MetadataTable"

Terminating on success
cpu time used  = 1.6321420000000000e+00 (s)

…点击查看完整日志…

检查输出：

bash

$ tree
.
├── build/
├── juno_mwr.inp
├── juno_mwr.out2.00000.nc
├── juno_mwr.out3.00000.nc
├── juno_mwr.out4.00000.nc
├── juno_mwr.out5.00000.nc
├── juno_mwr.release -> build/bin/juno_mwr.release
└── mwr_channels.yaml

$ ln -s build/bin/combine.py 
$ python3 combine.py -o test. ## 将输出合并成一个nc文件
$ tree
├── build/
├── combine.py -> build/bin/combine.py
├── combine_rules
├── juno_mwr.inp
├── juno_mwr.release -> build/bin/juno_mwr.release
├── juno_mwr-test-main.nc
└── mwr_channels.yam

结果分析

大气廓线输出：包括温度、位温、两种气体的浓度廓线。
juno-output-profiles.png
辐射传输模拟：包括6个通道的天顶亮温，临边昏暗和权重廓线
juno-output-radiances _24_.png

radiance _15_.png

Note:这里的权重指的是：

$\begin{array}{r} (1) & w (z) = - e^{- τ} d τ / d z \end{array}$

源码分析

比较特殊的头文件如下：

cpp

#include <snap/thermodynamics/vapors/sodium_vapors.hpp>
#include <astro/Jupiter/jup_fletcher16_cirs.hpp>
#include <astro/planets.hpp>
//radiative transfer
#include <harp/radiation.hpp>
#include <harp/radiation_band.hpp>
// tracer
#include <tracer/tracer.hpp>
// opacity
#include <opacity/Giants/microwave/mwr_absorbers.hpp>
// inversion
#include <inversion/profile_inversion.hpp>
// special includes
#include "juno_mwr_specs.hpp"

全局变量
Real grav, P0, T0, Tmin, clat;
Real xHe, xCH4, xH2S, xNa, xKCl, metallicity; // 痕量气体的含量

设定输出的诊断量，常规的温度和能量，加上几种condensible气体的相对湿度。这里只有H2O和NH3两种（NHYDRO=2）。

cpp

void MeshBlock::InitUserMeshBlockData(ParameterInput *pin) {
  AllocateUserOutputVariables(4 + NVAPOR);
  SetUserOutputVariableName(0, "temp");
  SetUserOutputVariableName(1, "theta");
  SetUserOutputVariableName(2, "thetav");
  SetUserOutputVariableName(3, "mse");
  for (int n = 1; n <= NVAPOR; ++n) {
    std::string name = "rh" + std::to_string(n);
    SetUserOutputVariableName(3 + n, name.c_str());
  }
}

NVAPOR,NHYDRO的定义在defs.hpp.in

输出前的计算

cpp

void MeshBlock::UserWorkBeforeOutput(ParameterInput *pin) {
  auto pthermo = Thermodynamics::GetInstance();

  for (int k = ks; k <= ke; ++k)
    for (int j = js; j <= je; ++j)
      for (int i = is; i <= ie; ++i) {
        user_out_var(0, k, j, i) = pthermo->GetTemp(this, k, j, i);
        user_out_var(1, k, j, i) = pthermo->PotentialTemp(this, P0, k, j, i);
        // theta_v
        user_out_var(2, k, j, i) =
            user_out_var(1, k, j, i) * pthermo->RovRd(this, k, j, i);
        // mse
        user_out_var(3, k, j, i) =
            pthermo->MoistStaticEnergy(this, grav * pcoord->x1v(i), k, j, i);
        for (int n = 1; n <= NVAPOR; ++n)
          user_out_var(3 + n, k, j, i) =
              pthermo->RelativeHumidity(this, n, k, j, i);
      }
}

初始化，读取适用于所有网格的设置：

cpp

void Mesh::InitUserMeshData(ParameterInput *pin) {
  grav = -pin->GetReal("hydro", "grav_acc1");
  P0 = pin->GetReal("mesh", "ReferencePressure");
  T0 = pin->GetReal("problem", "T1bar");

  Tmin = pin->GetReal("problem", "Tmin");
  clat = pin->GetOrAddReal("problem", "clat", 0.);

  xH2S = pin->GetReal("problem", "xH2S");

  metallicity = pin->GetOrAddReal("problem", "metallicity", 0.);

  xNa = pin->GetReal("problem", "xNa");
  xNa *= pow(10., metallicity);

  xKCl = pin->GetReal("problem", "xKCl");
  xKCl *= pow(10., metallicity);

  xHe = pin->GetReal("problem", "xHe");

  xCH4 = pin->GetReal("problem", "xCH4");
}

ProblemGenerator

这个个例比较复杂，我们逐块解释：

数据读取

cpp

Application::Logger app("main");  // 日志注册
app->Log("ProblemGenerator: juno"); 
auto pthermo = Thermodynamics::GetInstance(); // 获取热力学指针对象

// mesh limits
Real x1min = pmy_mesh->mesh_size.x1min;
Real x1max = pmy_mesh->mesh_size.x1max;
Real H0 = pcoord->GetPressureScaleHeight();  // 大气标高H0

这里读取了模拟的垂直坐标的范围和大气标高H0：
Input文件的范围是 x1min = -270.E3 到 x1max = 200.E3 的高度范围。
这里说明一下模拟使用的坐标系统：
$\begin{array}{r} (2) & z^{*} = H_{0} \cdot l n \frac{P}{P_{0}}; P_{0} = 1 b a r \end{array}$

对于等温大气： $z = z^{*} = H_{0} \cdot l n \frac{P}{P_{0}}$ ；
对于等熵大气： $z! = z^{*}$

cpp

// request temperature and pressure
app->Log("request T", T0);
app->Log("request P", P0);

// thermodynamic constants
Real gamma = pin->GetReal("hydro", "gamma");
Real Rd = pthermo->GetRd();
Real cp = gamma / (gamma - 1.) * Rd;

// index
auto pindex = IndexMap::GetInstance();
int iH2O = pindex->GetVaporId("H2O");
int iNH3 = pindex->GetVaporId("NH3");

app->Log("index of H2O", iH2O);
app->Log("index of NH3", iNH3);

迭代求算 $T_{s}$ , $P_{s}$ :

①计算地表气压和地表气温（干绝热模型）： $P_{s} = P_{0} \cdot e^{- x 1 m i n / H 0}$ ； $T_{s} = T_{0} \cdot {\frac{P_{s}}{P_{0}}}^{R_{d} / C_{p}}$ ;
②从地表开始构建一条绝热廓线到1bar的位置（湿绝热），获得是 $T_{0}^{'}$ ;
计算1bar处的温度偏差 $Δ T = T_{0}^{'} - T_{0}$ ，修正Ts： $T_{s}^{'} = T_{s} - Δ T$ ，
③自下而上再次构建湿绝热廓线，获得新的T_0”；
重复3和4，直到 $Δ T < 0.01 K$ 。

① ②两个过程并非简单的逆过程，否则2过程得到的 $T_{0}^{'} = T_{0}$ ，即迭代没有意义。实际上之所以这么做，是因为过程2在干绝热过程计算时，每一层混合气体的总 $C_{p}$ 随着气体组分变化。并且，单一组分气体的 $C_{p}$ 也在随温度变化的，而过程1使用了上下恒定的 $C_{p}$ 。
下图展示了最初几步① ② ③过程示意图：
求算Ts,Ps

求算Ts,Ps

稳定性条件

代码如下：

cpp

// estimate surface temperature and pressure
Real Ps = P0 * exp(-x1min / H0);
Real Ts = T0 * pow(Ps / P0, Rd / cp);

while (iter++ < max_iter) {
  // read in vapors
  air.w[iH2O] = xH2O;
  air.w[iNH3] = xNH3;
  air.w[IPR] = Ps;
  air.w[IDN] = Ts;

  // stop at just above P0
  for (int i = is; i <= ie; ++i) {
    pthermo->Extrapolate(&air, -dlnp / 2.,
                         Thermodynamics::Method::DryAdiabat);
    if (air.w[IPR] < P0) break;
  }

  // extrapolate down to where air is
  pthermo->Extrapolate(&air, log(P0 / air.w[IPR]),
                       Thermodynamics::Method::DryAdiabat);

  // make up for the difference
  Ts += T0 - air.w[IDN];
  if (std::abs(T0 - air.w[IDN]) < 0.01) break;

  app->Log("Iteration #", iter);
  app->Log("T", air.w[IDN]);
}

if (iter > max_iter) {
  throw RuntimeError("ProblemGenerator", "maximum iteration reached");
}

构建大气廓线：

cpp

// construct atmosphere from bottom up
air.ToMoleFraction();
for (int k = ks; k <= ke; ++k)
  for (int j = js; j <= je; ++j) {
    air.SetZero();
    air.w[iH2O] = xH2O;
    air.w[iNH3] = xNH3;
    air.w[IPR] = Ps;
    air.w[IDN] = Ts;

    int i = is;
    for (; i <= ie; ++i) {
      // check relative humidity
      Real rh = pthermo->RelativeHumidity(air, iNH3);
      air.w[iNH3] *= std::min(rh_max_nh3 / rh, 1.);

      AirParcelHelper::distribute_to_primitive(this, k, j, i, air);

      pthermo->Extrapolate(&air, -dlnp, Thermodynamics::Method::DryAdiabat);

      if (air.w[IDN] < Tmin) break;
    }

    // Replace adiabatic atmosphere with isothermal atmosphere if temperature
    // is too low
    pthermo->Extrapolate(&air, dlnp, Thermodynamics::Method::DryAdiabat);
    for (; i <= ie; ++i) {
      pthermo->Extrapolate(&air, -dlnp, Thermodynamics::Method::Isothermal);
      AirParcelHelper::distribute_to_primitive(this, k, j, i, air);
    }
  }

手动修正大气地不稳定度，使用adlnTdlnP，adlnNH3dlnP这个输入量来完成，会作用到Extrapolate函数上。

cpp

// if requested
// modify atmospheric stability
Real adlnTdlnP = pin->GetOrAddReal("problem", "adlnTdlnP", 0.);
Real adlnNH3dlnP = pin->GetOrAddReal("problem", "adlnNH3dlnP", 0.);

if (adlnTdlnP != 0.) {
  Real pmin = pin->GetOrAddReal("problem", "adlnTdlnP.pmin", 1.);
  Real pmax = pin->GetOrAddReal("problem", "adlnTdlnP.pmax", 20.);

  for (int k = ks; k <= ke; ++k)
    for (int j = js; j <= je; ++j) {
      int ibegin = find_pressure_level_lesser(pmax, phydro->w, k, j, is, ie);
      int iend = find_pressure_level_lesser(pmin, phydro->w, k, j, is, ie);

      auto &&air = AirParcelHelper::gather_from_primitive(this, k, j, ibegin);
      air.ToMoleFraction();

      for (int i = ibegin; i < iend; ++i) {
        pthermo->Extrapolate(&air, -dlnp, Thermodynamics::Method::DryAdiabat,
                             0., adlnTdlnP);
        AirParcelHelper::distribute_to_primitive(this, k, j, i + 1, air);
      }
    }
}

if (adlnNH3dlnP != 0.) {
  Real pmin = pin->GetOrAddReal("problem", "adlnNH3dlnP.pmin", 1.);
  Real pmax = pin->GetOrAddReal("problem", "adlnNH3dlnP.pmax", 20.);

  for (int k = ks; k <= ke; ++k)
    for (int j = js; j <= je; ++j) {
      int ibegin = find_pressure_level_lesser(pmax, phydro->w, k, j, is, ie);
      int iend = find_pressure_level_lesser(pmin, phydro->w, k, j, is, ie);

      auto &&air = AirParcelHelper::gather_from_primitive(this, k, j, ibegin);
      air.ToMoleFraction();

      for (int i = ibegin; i < iend; ++i) {
        pthermo->Extrapolate(&air, -dlnp, Thermodynamics::Method::DryAdiabat);
        air.w[iNH3] += adlnNH3dlnP * air.w[iNH3] * dlnp;
        auto rates = pthermo->TryEquilibriumTP_VaporCloud(air, iNH3);
        air.w[iNH3] += rates[0];
        AirParcelHelper::distribute_to_primitive(this, k, j, i + 1, air);
      }
    }
}

设置痕量气体和电子，碱金属。计算蒸汽压和浓度：

cpp

// set tracers, electron and Na
int ielec = pindex->GetTracerId("e-");
int iNa = pindex->GetTracerId("Na");
auto ptracer = pimpl->ptracer;

for (int k = ks; k <= ke; ++k)
  for (int j = js; j <= je; ++j)
    for (int i = is; i <= ie; ++i) {
      Real temp = pthermo->GetTemp(this, k, j, i);
      Real pH2S = xH2S * phydro->w(IPR, k, j, i);
      Real pNa = xNa * phydro->w(IPR, k, j, i);
      Real svp = sat_vapor_p_Na_H2S_Visscher(temp, pH2S);
      pNa = std::min(svp, pNa);

      ptracer->u(iNa, k, j, i) = pNa / (Constants::kBoltz * temp);
      ptracer->u(ielec, k, j, i) = saha_ionization_electron_density(
          temp, ptracer->u(iNa, k, j, i), 5.14);
    }

// primitive to conserved conversion (hydro)
peos->PrimitiveToConserved(phydro->w, pfield->bcc, phydro->u, pcoord, is, ie,
                           js, je, ks, ke);

// conserved to primitive conversion (tracer)
peos->PassiveScalarConservedToPrimitive(pscalars->s, phydro->u, pscalars->r,
                                        pscalars->r, pcoord, is, ie, js, je,
                                        ks, ke);

最后是辐射传输计算。

cpp

// Microwave radiative transfer needs temperatures at cell interfaces, which
// are interpolated from cell centered hydrodynamic variables. Normally, the
// boundary conditions are taken care of internally. But, since we call
// radiative tranfer directly in pgen, we would need to update the boundary
// conditions manually. The following lines of code updates the boundary
// conditions.
phydro->hbvar.SwapHydroQuantity(phydro->w, HydroBoundaryQuantity::prim);
pbval->ApplyPhysicalBoundaries(0., 0., pbval->bvars_main_int);

// calculate radiative transfer
auto prad = pimpl->prad;

for (int k = ks; k <= ke; ++k) {
  // run RT models
  app->Log("run microwave radiative transfer");
  for (int j = js; j <= je; ++j) prad->CalRadiance(this, k, j);
}

本节涉及到了使用湿绝热垂直递减率来构建大气温度廓线的方法，考虑了多种水凝物潜热释放的情形，具体请参考文献【Moist Adiabats with Multiple Condensing Species: A New Theory with Application to Giant-Planet Atmospheres】