Construction of hydrostatically stratified base state
Here we describe how ERF initializes base state values of density and potential temperature such that the density, pressure and potential temperature satisfy both the hydrostatic balance and the equation of state. The base-state density stored in ERF is the dry-air density \(\rho_{d,0}\). The only moisture variable carried in the base state is the water vapor mixing ratio \(q_{v,0}\); when total base-state density is needed, ERF forms \(\rho_{d,0}(1 + q_{v,0})\).
Users have the option to define a dry or moist background state.
Computation of the dry density
We express the total pressure \(p\) as:
By definition, we have:
This gives:
which, on rearranging, yields:
To obtain \(\theta_m\), consider the density of dry air:
Substituting for \(\rho_d\) from the above equation, we get:
which implies:
We also have:
where \(p_v\) is the partial pressure of water vapor, \(r = R_d/R_v \approx 0.622\), and \(q_v\) is the vapor mass mixing ratio—the ratio of the mass of vapor to the mass of dry air. Rearranging and using \(q_v \ll r\), we get:
which, on substitution in the equation for \(\frac{T_v}{T}\), gives:
As \(q_v \ll 1\), a binomial expansion, ignoring higher-order terms, gives:
Hence, the density of dry air is given by:
Initialization with a second-order integration of the hydrostatic equation
We have the hydrostatic equation given by
where \(\rho = \rho_d(1 + q_v)\), \(\rho_d\) is the dry density, and \(q_v\) is the water vapor mixing ratio carried in the base state. Using an average value of \(\rho\) for the integration from \(z = z(k-1)\) to \(z(k)\), we get
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 \(\rho(k) = f(p(k))\). Hence
Now, we define
This is a non-linear equation in \(p(k)\). Consider a Newton-Raphson iteration (where \(n\) denotes the iteration number) procedure
which implies
with the gradient being evaluated as
and the iteration update is given by
For the first cell (\(k=0\)), which is at a height of \(z = \frac{\Delta z}{2}\) from the base, we have
where \(p_0 = 1e5 \, \text{N/m}^2\) is the pressure at the base. Hence, we define
and the Newton-Raphson procedure is the same.
The base state on a refined mesh
With mesh refinement, the base state is constructed independently on each level, against that level’s own cell-centered heights \(z_{cc}\), by the same procedure described above. Running the same construction on every level, rather than interpolating the coarse answer, is what leaves each level in discrete hydrostatic balance on its own mesh.
Before that construction runs, the base state on a refined level is filled by conservative interpolation from its parent, and the physical boundary conditions for the base state are applied. That is what gives values to the part of the refined level that lies inside the domain but outside the refined grids, and to the ghost cells; the per-level construction then overwrites only the refined grids.
Two consequences are worth stating explicitly.
A refined level need not span the depth of the domain. The vertical integration is
causal upward: the value in a cell depends only on the cells at or below it. A refined
region that stops partway up the domain therefore gets exactly the base state it would
have had if it had been refined all the way to the model top – the values in the cells it
does contain are unchanged. What the integration does require is a valid starting value
in the lowest cell of each column, which is why the grids must not be decomposed in the
vertical (see Avoiding Decomposition in the Vertical Direction); amr.refine_grid_layout_z
defaults to 0 in ERF, so this holds unless it is overridden.
The base state does not depend on the coupling type. Under two-way coupling ERF averages \(\det J\) down from fine to coarse, so that the volume weighting used when averaging the state down telescopes and remains conservative. The cell-centered heights \(z_{cc}\) are deliberately not averaged down. Were they, the coarse heights would change after the coarse base state had already been built against them, leaving that base state out of hydrostatic balance with the heights the dynamics subsequently use for vertical gradients, Rayleigh damping and the sponge zones. Not averaging them also keeps \(z_{cc}\) consistent with the nodal heights \(z_{nd}\) it is derived from, which are not averaged down either.
Note that the base states of two levels will not be identical to one another wherever their terrain differs – for instance when a finer level resolves topography that its parent does not, either because it re-reads a terrain file at its own resolution or because it reads a nested wrfinput file. Each level is individually hydrostatic; the difference between them is set by the difference in surface elevation and is carried up the column as an essentially constant offset in pressure.
The base state in the lateral ghost cells
At a lateral boundary that is not periodic, the base state has to be given values in the ghost cells outside the domain. Where the boundary reflects – a symmetry or slip wall – the nodal mesh is extended by mirroring, so a ghost cell sits at exactly the height of the cell it reflects, and reflecting the base state with it is right.
Where the boundary is outflow or open, the nodal mesh is instead extrapolated past the domain, and a ghost cell does not sit at the height of the cell just inside it whenever the terrain has a slope where it meets the boundary. Copying the base state outward at constant index, as a zero-gradient condition would, therefore leaves values that belong to a different height than the one the mesh assigns to the cell: the base state there is neither the hydrostatic profile at that height nor in hydrostatic balance along the ghost column, while \(\det J\) and the metric terms in the same column are height-consistent.
ERF instead builds the base state in those ghost cells as the same reference atmosphere sampled at the height the mesh puts them at. Writing \(z_g\) for the height of the ghost cell and \(z_r\) for the height of the cell inside the domain it would otherwise have copied – the nearest cell that is inside the domain laterally, which for a corner means the nearest corner cell – and \(\delta z = z_g - z_r\):
\(p_0\) is carried across the height offset by the hydrostatic relation, using the scale height of the cell it came from,
\[p_0(z_g) = p_0(z_r)\, \exp\!\left(-\frac{g\, \delta z}{R_d T_v}\right), \qquad R_d T_v = \frac{p_0(z_r)}{\rho_0(z_r)}.\]\(\theta_0\) and \(q_{v,0}\) carry over unchanged. Reconstructing them at \(z_g\) would need the profile of the column above and below the cell, and the cells that can be read there are the ones the box in hand happens to hold, so the stencil – and with it the answer – would depend on where the grids are split in the vertical. Holding them fixed leaves an error of \((\partial\theta_0/\partial z)\,\delta z\), which for the offsets a terrain-fitted mesh produces is a few thousandths of a kelvin, against the several pascals in \(p_0\) that the transfer above removes.
\(\rho_0\) and \(\pi_0\) then follow from the equation of state, so that it holds exactly in the ghost cell.
The construction is local to each ghost cell, so it applies to any number of ghost layers and to the corners, and it reduces to the plain copy – bit for bit – wherever the two heights agree. That covers every constant-\(\Delta z\) mesh and any terrain that is flat where it meets the boundary, so runs that were not affected by the height offset in the first place are unchanged.
When terrain is built from a WRF input or metgrid file, the heights the file supplies live at w-faces, while ERF stores terrain at nodes. The nodal heights are obtained by inverting the four-node averaging operator exactly, so that the average of the four nodes around each w-face reproduces the file’s w-face height to round-off; the cell-centered heights, each being the average of the two w-faces that bracket it, then agree with the source as well. That operator is a tensor product of two one-dimensional two-point averages and is surjective, so an exact solution always exists; the remaining freedom is a sawtooth, which is fixed by minimizing the first-difference energy along each line. Should the exact heights fall outside the range of the source data, which very rough terrain can force, ERF warns and falls back to a regularized least-squares fit that stays bounded but no longer reproduces the source exactly.