1 Introduction

Planetary atmospheric boundary layer research has long been a popular field in atmospheric sciences (Mahrt, 1982, 1998; Tan et al., 2006; Baklanov et al., 2011). Studies have shown that anabatic and katabatic flows occurring in the boundary layer affect glacial mass variation, sea level rise, and climate change; thus, these processes are under the scientific spotlight (Grisogono and Oerlemans, 2001a; Barthélemy et al., 2012; Grisogono et al., 2015; Salvador et al., 2016). Katabatic and anabatic winds occur over a sloped surface. When the temperature of an air parcel near a sloping surface falls, a difference between it and the ambient temperature at the same elevation is established. The buoyancy force related to this temperature contrast projects in the along-slope direction and induces the katabatic wind. Its counterpart, i.e., anabatic wind, corresponds to heating over a sloping surface (Shapiro et al., 2012). Research on the modeling and theory of a stably stratified atmospheric boundary layer is an active field (Grisogono and Oerlemans, 2001b; Grisogono, 2003; Baklanov et al., 2011; Ingel, 2011; Grisogono et al., 2015). Flows over sloping terrain have been simulated with numerical weather prediction and climate models (Smith and Skyllingstad, 2005; Jeričević et al., 2010; De Ridder et al., 2015). However, the handling of the relevant lower boundary conditions and near-surface turbulence parameterization remains unresolved (Epifanio, 2007; Buzzi et al., 2011). Therefore, researchers continue to focus their attention on anabatic and katabatic flows (Shapiro et al., 2012).

Among the studies on anabatic and katabatic flows, simulation based on analytical modeling is the most fundamental. Most analytical models of katabatic flows are one dimensional; furthermore, the most common analytical model of thermally driven simple slope flows is the Prandtl model (Egger, 1990; Axelsen and van Dop, 2009). The strengths and weaknesses of the classical Prandtl model have been explored by Grisogono and Axelsen (2012), mainly in large eddy simulations based on limited observational data (van den Broeke, 1997). This line of work was improved by Mo (2013); however, serious limitations remain, not least insofar as most analytical models of katabatic flow are linear and the eddy thermal conductivity is assumed to be constant. The main disadvantage of a linear Prandtl model is that the thermodynamic equation does not consider the fact that the flow-related potential temperature gradient near the surface feeds back to the environmental potential temperature gradient (Grisogono et al., 2015).

The abovementioned defect of the classical Prandtl model was improved by Grisogono et al. (2015), who introduced a small parameter ε into the thermodynamic equilibrium equation that controls the transport from the flow-induced potential temperature gradient to the environment. The addition of this parameter makes the model more realistic, as it transforms from a linear model into a weakly nonlinear model. The scope of the small parame-ter $ \varepsilon $ is ${\rm{0}} \leqslant \varepsilon \leqslant {\rm{0}}{\rm{.01}}$ , according to Grisogono et al. (2015).

A second disadvantage of the classical Prandtl model was amended by assuming gradually varying eddy thermal conductivity K; that is, K is a function of height z (Grisogono and Oerlemans, 2001a,b, 2002), which is closer to reality. Shapiro et al. (2016) found that the time-varying eddy thermal conductivity K was potentially more significant than the vertical change when they studied the Great Plains nocturnal low-level jet. Because our research is based on the weakly nonlinear Prandtl model in the stable boundary layer, the vertical change in K is our main concern. However, the varying K(z) is generally obtained by an empirical parameterization scheme, which is an approximation of reality. Thus, it is worth studying how to precisely retrieve the varying eddy thermal conductivity.

As is well known, the problems in parameter retrieval are difficult to resolve because they are generally nonlinear and sometimes ill-posed. Indeed, there are still many important theoretical issues being investigated (Huang et al., 2006). Parameter retrieval via inversion is an approach that has been successfully applied in many disciplines, such as the geosciences, biological sciences, and atmospheric and oceanic sciences. Some previous studies on inverse problems can be found in Huang et al. (2005), Zhao and Huang (2012), and Yan and Huang (2014). Typically, parameter inversion is mainly tackled by using the variational method and Kalman filter me-thod, and there have been many studies conducted using boundary layer parameter inversion. For example, Hu et al. (2010) explored the treatment of model errors and uncertainties through simultaneous state and parameter estimation with an ensemble Kalman filter (EnKF). Nielsen-Gammon et al. (2010) evaluated the sensitivities of a planetary boundary layer scheme for parameter estimation with an EnKF. Huang et al. (2005) analyzed the retrieval of atmospheric and oceanic parameters with variational data assimilation and the relevant numerical calculation. Zhao and Huang (2012) retrieved atmos-pheric refractivity profiles from sea surface backscattered radar clutter with a variational method.

As for the retrieval of gradually varying K(z), this too is a typical parameter retrieval problem. In both the atmospheric sciences and the hydromechanics research, K(z) is a critical parameter and is usually derived from an empirical parameterization scheme. It is actually very difficult to obtain K because it is closely related to the dynamic thermal structure of the fluid. It would therefore be an important step to develop a method that can achieve a reasonable retrieval of K(z) in the weakly nonlinear Prandtl model, as this would be highly useful in improving theoretical analyses and numerical models.

This paper describes the retrieval of eddy thermal conductivity in the weakly nonlinear Prandtl model through the variational assimilation method. Because observational data are difficult to obtain, we demonstrate the effectiveness and reliability of the method through numerical experiments.

The rest of this paper is arranged as follows. Section 2 introduces the weakly nonlinear Prandtl model; Section 3 describes the variational assimilation scheme for the model; Section 4 reports the numerical experiment results; and finally, Section 5 provides a summary and discussion. Details of the Wentzel–Kramers–Brillouin (WKB) calculation steps are also included, in the Appendix.

2 Weakly nonlinear Prandtl model 2.1 Model

The classical Prandtl model has been validated and used by numerous scientists and engineers over the past few decades (Mahrt, 1982; Egger, 1990; Axelsen and van Dop, 2009). This linear model is imperfect, given that the actual physical process is generally nonlinear. In the past decade, the model has been modified to more accurately represent the real physical process. Among these improvements, introducing the small parameter $\varepsilon $ in the classical Prandtl model is the most recent and the best advancement (Grisogono et al., 2015). The introduction of the small parameter $\varepsilon $ can change the model into a weakly nonlinear one.

The weakly nonlinear Prandtl model is essentially a simplified, steady-state, one-dimensional model with governing equations in the tilted coordinate system and with a constant slope α, $\left| \alpha \right| \ll 1$ ; therefore, the quasi-hydrostatic approximation is satisfied (Mahrt, 1982; Grisogono and Axelsen, 2012). The momentum and thermodynamic equations are as follows (Denby, 1999):

$\begin{split}& g\frac{\theta }{{{\Theta _{\rm{0}}}}}\sin (\alpha) + K\Pr \frac{{{\partial ^2}u}}{{\partial {z^2}}} = 0, \\ & - (\Gamma + \varepsilon \frac{{\partial \theta }}{{\partial z}})u\sin (\alpha) + K\frac{{{\partial ^2}\theta }}{{\partial {z^2}}} = 0.\end{split}$

(1)

The term $\varepsilon \cdot \partial \theta /\partial z$ is the new term compared with the classical Prandtl model. In this paper, we use the same value, $\varepsilon \,{\rm{ = 0}}{\rm{.005}}$ , as in Grisogono et al. (2015). The classical Prandtl model becomes a weakly nonlinear form, where g is the gravitational acceleration, $\theta $ represents the potential temperature deviation from the background, K(z) is eddy thermal conductivity, and Pr is the turbulent Prandtl number (assumed to be constant in this paper). According to the conclusion in Grisogono et al. (2015), $K\!(z) = {K_0}(z/h)\exp (- {z^2}/2{h^2}) + {K_{\min }}$ , where K_min is typically tenfold of that for the molecular diffusivity of air [ ${K_{\min }} = 0.0002$ ; see Eq. (25) in Sun et al. (2013)]; h is the level where K(z) reaches the maximum, i.e., $\max(K\!(z)) = {K_0} \cdot \exp (- 1/2)$ ; u is the downslope wind; z is the coordinate that is orthogonal to the slope surface; α < 0 and α > 0 represent the anabatic and katabatic flows, respectively; and ${\Theta _0}$ is the reference temperature. In addition, $\Gamma = {\rm d}\Theta /{\rm d}{z^*}$ , where $z^* $ represents the true vertical direction. The slope coordinate system is shown in Fig. 1. We assume the system is steady state and one dimensional, so the local variation term in the equation is zero. The Boussinesq and rotation-free approximations and quasi-hydrostatic balance are included in this system.

2.2 WKB solutions

The forward problem of the weakly nonlinear Prandtl model is a weakly nonlinear two-point boundary value problem. Because of $z \in [0, + \infty)$ , the result is inaccurate when using the conventional numerical approach or spectral method, i.e., the two-point boundary value problem is mathematically ill-posed. Thus, we adopt the WKB method, which is the same as the method used in Grisogono et al. (2015), to solve the forward problem approximately. The basic idea of the method is described as follows. Assume

$\begin{array}{l}u = {u_0} + \varepsilon {u_1} + \varepsilon u_2^2 + \cdot \cdot \cdot, \\\theta = {\theta _0} + \varepsilon {\theta _1} + \varepsilon \theta _2^2 + \cdot \cdot \cdot, \end{array}$

(2)

where the subscripts represent the order of expansion, with zero-order flow variables being the solution for the classic Prandtl model, and the addition to first-order corrections being the solution for the weakly nonlinear mo-del. We maintain this expansion only for the first-order corrections. By substituting Eq. (2) into Eq. (1), the solution for the zero-order expansion is obtained, i.e., the solution for the classical Prandtl model (Mahrt, 1998; Grisogono and Oerlemans, 2002; Stiperski et al., 2007) is

$\frac{{{{\rm d}^4}({u_0}, {\theta _0})}}{{{\rm d}{z^4}}} + {N^2}\frac{{{{\sin }^2}(\alpha)}}{{\Pr K\!{{(z)}^2}}}({u_0}, {\theta _0}) = 0, $

(3)

where ${N^2} = \Gamma g/{\Theta _0}$ , and the boundary conditions are: ${\theta _{\rm{0}}}(0) = C$ , ${u_{\rm{0}}}(0) = 0$ ; ${\theta _{\rm{0}}}(z \to \infty) = 0$ , ${u_{\rm{0}}}(z \to \infty) = 0$ . The form of the solution is obtained as follows,

$\begin{array}{l}{u_0} = - C\mu \exp [ - I]\sin [I], \\{\theta _0} = C\exp [ - I]\cos [I], \end{array}$

(4)

where C is the surface potential temperature deficit in the case of katabatic flows; $\mu = {(g/{\Theta _0}\Gamma \Pr)^{\frac{1}{2}}}$ , which is a dimensional parameter in the unit of m °C ^–1 s^–1; $I(z) = {({\sigma _0}/2)^{\frac{1}{2}}}\int_0^z {K\!{{(z)}^{ - \frac{1}{2}}}{\rm d}z} $ ; and ${\sigma _0} = N\sin (\alpha){\Pr ^{ - \frac{1}{2}}}$ .

The solution for the first-order corrections $({u_1}, {\theta _1})$ is included in the following equation:

$\begin{split}& g\frac{{{\theta _1}}}{{{\Theta _0}}}\sin (\alpha) + K\Pr \frac{{{{\rm d}^2}{u_1}}}{{{\rm d}{z^2}}} = 0, \\& - (\Gamma {u_1} + \frac{{{\rm d}{\theta _0}}}{{{\rm d}z}}{u_0})\sin (\alpha) + K\frac{{{{\rm d}^2}{\theta _1}}}{{{\rm d}{z^2}}} = 0.\end{split}$

(5)

Similarly, as before, Eq. (5) can be transformed into:

$\begin{array}{c} \frac{{{{\rm d}^4}{u_1}}}{{{\rm d}{z^4}}} + {N^2}\frac{{{{\sin }^2}(\alpha)}}{{\Pr K\!{{(z)}^2}}}{u_1} = - \frac{{g{{\sin }^2}(\alpha)}}{{{\Theta _0}K\!{{(z)}^2}\Pr }}{u_0}\frac{{{\rm d}{\theta _0}}}{{{\rm d}z}}, \\\!\!\!\!\! \frac{{{{\rm d}^4}{\theta _1}}}{{{\rm d}{z^4}}} + {N^2}\frac{{{{\sin }^2}(\alpha)}}{{\Pr K\!{{(z)}^2}}}{\theta _1} = \frac{{\sin (\alpha)}}{{K\!(z)}}\frac{{{{\rm d}^2}}}{{{\rm d}{z^2}}}({u_0}\frac{{{\rm d}{\theta _0}}}{{{\rm d}z}}).\end{array}$

(6)

The associated boundary conditions are: ${\theta _1}(0) = 0$ , ${u_1}(0) = 0$ ; ${\theta _1}(z \to \infty) = 0$ , ${u_1}(z \to \infty) = 0$ .

For details of the calculation steps, please refer to the Appendix. The form of the solution is obtained as follows:

$\begin{aligned}{u_1} & = {u_A}\exp (- I)[ - \frac{1}{3}\sin (I) + \frac{2}{{15}}\cos (I)] \\ & + {u_A}\exp (- 2I)[\frac{1}{{30}}\sin (2I) - \frac{1}{{30}}\cos (2I) - \frac{1}{{10}}], \\{\theta _1} & = {\theta _A}\exp (- I)[ - \frac{1}{{15}}\sin (I) - \frac{1}{6}\cos (I)] \\ &+ {\theta _A}\exp (- 2I)[\frac{1}{{15}}\sin (2I) + \frac{1}{{15}}\cos (2I) + \frac{1}{{10}}], \end{aligned}$

(7)

where ${u_A} = {({\sigma _0}/2)^{\frac{1}{2}}}{C^2}(\mu /\Gamma)K\!{(z)^{ - \frac{1}{2}}}$ , ${\theta _A} = {(2/{\sigma _0})^{\frac{1}{2}}}{C^2}\mu \sin (\alpha)K\!{(z)^{ - \frac{1}{2}}}$ .

Thus, the total WKB approximate solutions, up to the first order, are as follows:

$\begin{array}{l}u = {u_0} + \varepsilon {u_1}, \\\theta = {\theta _0} + \varepsilon {\theta _1}.\end{array}$

(8)

3 Variational assimilation scheme

The idea of variational data assimilation was proposed by the atmospheric science community in the 1980s (Lorenc, 1988; Gao et al., 1995; Navon, 1998). The primary idea of variational assimilation technology is to combine the model and observation data to obtain the optimal model parameter or model initial value. The essence of this method is optimal control theory (Tosaka et al., 1999). After decades of development, variational assimilation has become a very mature method and has been implemented in operational models. Thus, we use the variational assimilation method to retrieve the eddy thermal conductivity in the weakly nonlinear Prandtl model for katabatic flows.

3.1 Defining the cost function

Let

$J[{K_0}, h] \!=\! \frac{1}{2}\int_0^H \!\! {[{{(u \!-\! {u^{{\rm {obs}}}})}^2} \!+\! \gamma {{(\theta \!-\! {\theta ^{\rm {obs}}})}^2}]{\rm d}z} {\rm{ = }}\min, $

(9)

where $({K_0}, h)$ is the parameter to be retrieved; $({u^{\rm {obs}}}, {\theta ^{\rm {obs}}})$ are the observation data of wind speed and potential temperature, provided in advance; H is the maximum height of measurement; and $\gamma $ is the dimension adjustment value (taken as $\gamma = 1$ in the dimensional analy-sis in this paper). Our aim is to identify the parameter for which the cost function is minimal; that is, the value calculated by the weakly nonlinear Prandtl model that is the closest to the observation data. Regardless of the method we choose to minimize the cost function (e.g., Newton’s iteration method, steepest descent method, conjugate gradient methods), the computation of ${\nabla _{{K_0}}}J$ and ${\nabla _h}J$ is necessary.

3.2 Variation of the function

According to the form of variation:

$\begin{aligned}\delta J & = \int_0^H {[(u - {u^{{\rm {obs}}}})\delta u + \gamma (\theta - {\theta ^{{\rm {obs}}}})\delta \theta ]} {{\rm d}z}\\ & = \int_0^H {[(u - {u^{{\rm {obs}}}})(\frac{{\partial u}}{{\partial I}}\delta I + \frac{{\partial u}}{{\partial K}}\delta K) + \gamma (\theta - {\theta ^{{\rm {obs}}}})(\frac{{\partial \theta }}{{\partial I}}\delta I + \frac{{\partial \theta }}{{\partial K}}\delta K)]} {{\rm d}z}\\ & = \int_0^H {(u - {u^{{\rm {obs}}}})[\frac{{\partial u}}{{\partial I}}(\frac{{\partial I}}{{\partial K}}, \delta K) + \frac{{\partial u}}{{\partial K}}\delta K]{{\rm d}z}} \\ &{\rm{ + }}\int_0^H {\gamma (\theta - {\theta ^{{\rm {obs}}}})[\frac{{\partial \theta }}{{\partial I}}(\frac{{\partial I}}{{\partial K}}, \delta K) + \frac{{\partial \theta }}{{\partial K}}\delta K]{{\rm d}z}} \\ & = \int_0^H {(u - {u^{{\rm {obs}}}})[\frac{{\partial u}}{{\partial I}}(\frac{{\partial I}}{{\partial K}}, \frac{{\partial K}}{{\partial {K_0}}}\delta {K_0} + \frac{{\partial K}}{{\partial h}}\delta h) + \frac{{\partial u}}{{\partial K}}(\frac{{\partial K}}{{\partial {K_0}}}\delta {K_0} + \frac{{\partial K}}{{\partial h}}\delta h)]{{\rm d}z}} \\ & + \int_0^H {\gamma (\theta - {\theta ^{{\rm {obs}}}})[\frac{{\partial \theta }}{{\partial I}}(\frac{{\partial I}}{{\partial K}}, \frac{{\partial K}}{{\partial {K_0}}}\delta {K_0} + \frac{{\partial K}}{{\partial h}}\delta h) + \frac{{\partial \theta }}{{\partial K}}(\frac{{\partial K}}{{\partial {K_0}}}\delta {K_0} + \frac{{\partial K}}{{\partial h}}\delta h)]} {{\rm d}z}.\end{aligned}$

(10)

Here, $u = u(I(K), K)$ , so $\delta u = \displaystyle\frac{{\partial u}}{{\partial I}}\delta I + \displaystyle\frac{{\partial u}}{{\partial K}}\delta K$ . Also, $I = {({\sigma _0}/2)^{\frac{1}{2}}}\int_0^z {K\!{{(z)}^{ - \frac{1}{2}}}{\rm d}z} $ , which is a functional expression, so $\delta I = (\displaystyle\frac{{\partial I}}{{\partial K}}, \delta K \! )$ = $\displaystyle \int_{0}^{H}\displaystyle{\frac{\partial I}{\partial K}}\cdot \delta K\text{d}z $ , according to the definition of variation $\delta J= {\nabla _{{K_0}}}J \cdot \delta {K_0} +{\nabla _h}J \cdot \delta h$ . Therefore,

$\begin{array}{l}{\nabla _{{K_0}}}J = \int_0^H [{(u - {u^{{\rm {obs}}}})(\frac{{\partial u}}{{\partial I}}(\frac{{\partial I}}{{\partial K}}, \frac{{\partial K}}{{\partial {K_0}}}) + \frac{{\partial u}}{{\partial K}}\frac{{\partial K}}{{\partial {K_0}}}) + \gamma (\theta - {\theta ^{{\rm {obs}}}})(\frac{{\partial \theta }}{{\partial I}}(\frac{{\partial I}}{{\partial K}}, \frac{{\partial K}}{{\partial {K_0}}}) + \frac{{\partial \theta }}{{\partial K}}\frac{{\partial K}}{{\partial {K_0}}})]{{\rm {d}}z}}, \\{\nabla _h}J = \int_0^H [{(u - {u^{{\rm {obs}}}})(\frac{{\partial u}}{{\partial I}}(\frac{{\partial I}}{{\partial K}}, \frac{{\partial K}}{{\partial h}}) + \frac{{\partial u}}{{\partial K}}\frac{{\partial K}}{{\partial h}}) + \gamma (\theta - {\theta ^{{\rm {obs}}}})(\frac{{\partial \theta }}{{\partial I}}(\frac{{\partial I}}{{\partial K}}, \frac{{\partial K}}{{\partial h}}) + \frac{{\partial \theta }}{{\partial K}}\frac{{\partial K}}{{\partial h}})]{{\rm {d}}z}} .\end{array}$

(11)

3.3 Numerical calculation procedure

In this paper, the steepest descent method is selected for numerical calculation. The procedure is as follows:

Step 1: Prescribe an initial value for $(K_0^0, {h^0})$ ;

Step 2: Calculate K(z) according to the formula $K\!(z) = {K_0}(z/h)\exp (- {z^2}/2{h^2}) + {K\!_{\min }}$ (K_min is known beforehand);

Step 3: obtain the solution for the weakly nonlinear Prandtl model through the WKB method and the known parameter;

Step 4: Compute the gradient $({\nabla _{{K_0}}}J, {\nabla _h}J)$ through Eq. (11);

Step 5: Set a step length, find the value of $(K_0^{n + 1}, {h^{n + 1}})$ , and calculate the cost function $J[{K_{\rm{0}}}, h]$ of the n+1 iterations. If it is smaller than the previous time, then go to Step 6. If not, adjust the step length until the value of the cost function is smaller than the previous time.

Step 6: Print the result if the cost function is smaller than a very small value; otherwise, replace $(K_0^0, {h^0})$ with $(K_0^{n + 1}, {h^{n + 1}})$ , and return to Step 2.

A flow chart illustrating the numerical calculation is provided in Fig. 2.

The initial value is ideally supposed to be logical. It should not be too far away from the truth. We take $K_0^0 = 1$ and ${h^0} = H/2$ , where H is the maximum height of simulated measurement. We assume H = 80 m. The initial iterative step length is one. If it is large, then the iteration step length is reduced by multiplying by 0.618. If δ is small enough, we assume δ = 10^–4.

4 Numerical experiment

To validate the robustness and reliability of the proposed method for parameter retrieval in the weakly nonlinear Prandtl model of katabatic flows, numerical testing is carried out through two experiments, each with two versions: observation with errors and without errors. Near-surface observation data with high accuracy are difficult to obtain; therefore, in this paper, we only use the numerical experiment with simulated observation data to test the proposed method.

However, when validating the proposed method with simulated observations, a problem emerges insofar as the observations are highly dependent on the model (Arnold and Dey, 1986). Considering the independence of the data to the model, we decide to take $\tilde K\!(z) =$ $ {K_0}(z/h)(1 \!+\! \eta \cdot z/h)\exp (- {z^2}/2{h^2} \cdot (1 \!+\! \eta \cdot z/h)) + {K\!_{\min }}$ [first-order Taylor series expansion of original K(z)] as the model of the simulated observation data. The introduction of nonlinear effects alleviates the dependency. The determination of the small parameter η should sati-sfy the following two requirements: the distance between the two curves is not too large; the nonlinear effects should be significant. After comparing the distribution of K(z) and $\tilde K(z)$ with height, we choose η = 0.05. To avoid the “inverse crime” property, grid sizes for the forward and inverse problem should not be the same (Kaipio and Somersalo, 2005). A gradual grid is chosen in the differential process. In the simulated observation data, the grid density decreases with height. This corresponds to the fact that the real observation density in the Pasterze experiment (PASTEX-94), Austria, decreases with height. The grid distance of the simulated observation in the numerical experiment is 0.3 m when the height is below 9 m; 0.6 m when the height is between 9 and 21 m; 1.2 m when the height is between 21 and 45 m; and 2.5 m when the height is between 45 and 80 m. Because the grid distance of the inverse problem is 0.25 m, the grid of the simulated observation is different from the grid of the inverse problem. The model dependency involves a further decrease. In future, adding observation noise into the experiment might reduce the model dependency.

4.1 Experiment Ⅰ 4.1.1 Case Ⅰ: $({u^{\rm {obs}}}, {\theta ^{\rm {obs}}})$

Assuming $({K^0}, h) = (0.4746, 30)$ , the values of the remaining parameters are as follows: $g = 9.8\,\, {\rm{m}}\,\, {{\rm{s}}^{ - 2}}$ ; ${\Theta _0} = 273.15\,\, {\rm{K}}$ ; $\Gamma {\rm{ = 0}}{\rm{.003}}\,\, {\rm{K}}\,\, {{\rm{m}}^{{\rm{ - 1}}}}$ ; Pr = 2; α = –5°;ε = 0.005; and C = –6°C (these are the typical values for the Prandtl model). The boundary conditions are: $\theta (0) = C$ , $u(0) = 0$ ; $\theta (z \to \infty) = 0$ , $u(z \to \infty) = 0$ . The simulation data $({u^{\rm {obs}}}, {\theta ^{\rm {obs}}})$ are obtained by the WKB method. Next, the parameter K(z) is retrieved in accordance with the numerical calculation procedure. We chose the root-mean-square error (RMSE) and relative root-mean-square (R-RMSE) as the error criterion. The result is shown in Fig. 3.

The error mainly arises in the region of maximum K. The maximum absolute error is 0.005, which is small compared with the maximum value of K. The maximum relative error is about 4%, and the mean relative error is about 1.7%. The R-RMSE is 1.9%, and the RMSE is 0.0027. The result demonstrates that the proposed me-thod is effective for observations without noise.

4.1.2 Case Ⅱ: $({u^{\rm {obs}}}, {\theta ^{\rm {obs}}})$ with random noise

The proposed variational assimilation scheme is robust; it can restrain the influence of the observation error. We illustrate this characteristic in the numerical experiment of Case II. The values of all the parameters are the same as in Case I. Assuming there is random noise with the observation, and the error disturbance has a normal distribution with mean of 0 and standard variance of 0.3. The random noises imposed on wind speed and potential temperature are shown in Fig. 4a, b. The random observation errors depicted by Fig. 4 are expected to represent the real observation errors to the most reasonable degree.

We use the observations of wind speed and potential temperature with errors to retrieve the parameter K(z). The retrieved profiles of K(z) are shown in Fig. 5. In Fig. 5, the mean relative error is 6.1%; the maximum relative error is 12%; R-RMSE is 7%; and RMSE is 0.0096. These numbers indicate that the variational assimilation scheme we used has decently restrained the influence of the observation error. Therefore, we conclude that the proposed method is capable of accommodating the observation noises and suitable to retrieve the parameter in the weakly nonlinear Prandtl model with observation errors.

4.2 Experiment Ⅱ

It is expected that an increase in the number of samples should make the variational assimilation scheme more persuasive. In this experiment, we assume that $({K^0}, h) = (0.{\rm{8}}, {\rm{15}})$ , and the remaining parameters and boundary conditions are the same as in Experiment I. The results are shown in Fig. 6. Figure 6a displays the inversion results without observation errors while Fig. 6b displays the inversion results with observation errors. The error disturbance is in accordance with the normal distribution (average value is 0; standard variance is 0.3)—the same as in Experiment I. Figure 6 shows that the accuracy of the results is high, especially for the situation with no observation errors. The accuracy will decrease when the observation has random noise.

The error distribution with height is shown in Fig. 7. The maximum relative error of the retrieval without observation noise is 25.7% and the mean relative error is 10.5% (Fig. 7a). Meanwhile, the relative error below 40 m is less than 10%, and the R-RMSE and RMSE is 13.7% and 0.0175, respectively (Fig. 7a). The maximum relative error of the retrieval with observation noise is 39.5%; the mean relative error is 18.5%; the relative error below 40 m is less than 20%; the R-RMSE is 19.8%; the RMSE is 0.022; and the maximum absolute error is 0.04 (Fig. 7b). The maximum relative errors appear in the upper part of both Fig. 7a and Fig. 7b, where exact eddy thermal conductivities are close to zero. The large relative errors above 40 m seem to be associated with small K(z) values at zero and its vicinity.

5 Conclusions and discussion

Research on the slope flows in the atmospheric boundary layer has been a popular field in atmospheric sciences. Research on katabatic flows provides guidance to the development of weather and climate prediction models. The weakly nonlinear Prandtl model is more efficient than the classical Prandtl model as it introduces the small parameter $\varepsilon $ and gradually varying eddy thermal conductivity K(z). The gradually varying K(z) has inestimable value, in both the fluid dynamics and the atmospheric sciences. This parameter is difficult to determine, so an empirical parameterization scheme is usually adopted to solve the eddy thermal conductivity. We utilize an idealized observation of u and θ, together with the variational assimilation technique, to solve the parameterization of K(z) in this paper. The coefficients in the parameterization scheme can be obtained precisely by using this method. The efficiency of the proposed method is validated via results from numerical experiments. The effects of the observation error can also be effectively reduced. All the above results demonstrate that the method is feasible and reliable.

Research on inversion problems tells us that obtaining a retrieval model parameter is an ill-posed problem. It is more difficult to retrieve the K(z) without a fixed form. If we want to obtain the vertical distribution of K(z), we must obtain the statistical background information of K(z) first, and then use the variational regularization method to solve K(z). The WKB method is selected in this paper, and the solution is only an approximate one. Further study of direct numerical calculation is required.

Acknowledgments. We thank Professor Branko Grisogono for the PASTEX-94 dataset. We also thank the reviewers and Editor Xiaoming Hu for their constructive criticism that contributes to significant improvements to the original manuscript and greater reliability in interpretation of the numerical experiment results.

Appendix

This appendix contains the details of the WKB calculation steps. The weakly nonlinear Prandtl model is expressed as follows (meaning of each expression is the same as detailed in Section 2):

$\begin{split}& g\frac{\theta }{{{\Theta _{\rm{0}}}}}\sin (\alpha) + K\Pr \frac{{{\partial ^2}u}}{{\partial {z^2}}} = 0, \\& - (\Gamma + \varepsilon \frac{{\partial \theta }}{{\partial z}})u\sin (\alpha) + K\frac{{{\partial ^2}\theta }}{{\partial {z^2}}} = 0.\end{split}$

(A1)

The expression of the global expansion of $(u, \theta)$ is:

$\begin{array}{l}u = {u_0} + \varepsilon {u_1} + \varepsilon u_2^2 + \cdot \cdot \cdot, \\\theta = {\theta _0} + \varepsilon {\theta _1} + \varepsilon \theta _2^2 + \cdot \cdot \cdot .\end{array}$

(A2)

Now, we substitute Eq. (A1) with Eq. (A2) and keep the expansion up to the first order:

$\begin{align} & g\frac{{{\theta }_{0}}}{{{\Theta }_{\text{0}}}}\sin (\alpha )+K\Pr \frac{{{{\rm d}}^{2}}{{u}_{0}}}{{\rm d}{{z}^{2}}} \\ \quad & +\varepsilon (g\frac{{{\theta }_{1}}}{{{\Theta }_{\text{0}}}}\sin (\alpha )+K\Pr \frac{{{{\rm d}}^{2}}{{u}_{1}}}{{\rm d}{{z}^{2}}})=0, \\ \quad & -\Gamma {{u}_{0}}\sin (\alpha )+K\frac{{{{\rm d}}^{2}}{{\theta }_{0}}}{{\rm d}{{z}^{2}}} \\ \quad & +\varepsilon (-(\Gamma {{u}_{1}}+\frac{{\rm d}{{\theta }_{0}}}{{\rm d}z}{{u}_{0}})\sin (\alpha )+K\frac{{{{\rm d}}^{2}}{{\theta }_{1}}}{{\rm d}{{z}^{2}}})=0. \\ \end{align}$

(A3)

A1. Zero-order expansion

The zero-order expansion reads

$\begin{split}& g\frac{{{\theta _0}}}{{{\Theta _{\rm{0}}}}}\sin (\alpha) + K\Pr \frac{{{{\rm d}^2}{u_0}}}{{{\rm d}{z^2}}} = 0, \\ & - \Gamma {u_0}\sin (\alpha) + K\frac{{{{\rm d}^2}{\theta _0}}}{{{\rm d}{z^2}}} = 0.\end{split}$

(A4)

System (A4) can be rewritten to a single governing equation of the fourth order by the elimination method. The terms dK/dz and ${{\rm d}^2}K \!\!/{\rm d}{z^2}$ are ignorable according to Grisogono (2003). Since K(z) is gradually varying with z, i.e., $K\!(Z) = K\!(\varepsilon z)$ , so ${\rm d}K\!(Z)/{\rm d}z = \left({{\rm d}K\!/{\rm d}Z} \right)\left({{\rm d}Z\!/{\rm d}z} \right) = \varepsilon \left({{\rm d}K\!\!/{\rm d}Z} \right)$ , ${{\rm d}^2}K\!(Z)/{\rm d}{z^2} = {\varepsilon ^2}\left({{{\rm d}^2}K\!/{\rm d}{Z^2}} \right)$ .

As $ \varepsilon $ is a small parameter, it makes sense to proceed neglecting ${\rm d}K\!(Z)/{\rm d}z$ and ${{\rm d}^2}K\!(Z)/{\rm d}{z^2}$ :

$\frac{{{{\rm d}^4}({u_0}, {\theta _0})}}{{{\rm d}{z^4}}} + {N^2}\frac{{{{\sin }^2}(\alpha)}}{{\Pr K\!{{(z)}^2}}}({u_0}, {\theta _0}) = 0.$

(A5)

If $K\!(z) = {\rm constant}$ in System (A5), the solution is already known from numerous previous studies (Egger, 1990; Grisogono and Oerlemans, 2001a). It is a homogeneous linear differential equation whose solution can be obtained by the method of undetermined coefficients with boundary conditions. It reads:

$\begin{array}{l}{u_0} = - C\mu \exp (- \frac{z}{{{h_{\rm p}}}})\sin (\frac{z}{{{h_{\rm p}}}}), \\{\theta _0} = C\exp (- \frac{z}{{{h_{\rm p}}}})\cos (\frac{z}{{{h_{\rm p}}}}), \end{array}$

(A6)

where C is the surface potential temperature deficit in the case of katabatic flows. The opposite, C < 0, is the temperature surplus, for anabatic flows. The parameter $\mu = {(g/{\Theta _0}\Gamma \Pr)^{\frac{1}{2}}}$ is a dimensional parameter (units: m °C s ^–1); ${h_{\rm p}} = \sqrt 2 /\sigma $ is the characteristic depth of the Prandtl layer; and σ is the characteristic “wave number” expressed as $\sigma = {[N\sin (\alpha)/(K{\Pr ^{\frac{1}{2}}})]^{\frac{1}{2}}}$ , where N is the background or free-flow buoyancy frequency (N² = $ \Gamma g/{\Theta _0}$ ). The boundary conditions are the same as detailed in Section 2.

If K(z) is the gradually varying term in Eq. (A5), then the argument in Eq. (A6), $z/{h_{\rm p}}$ , becomes an integral, i.e., $z/{h_{\rm{p}}} \to \int {\rm{d}} \left( {z/{h_{\rm{p}}}} \right) = I(z)$ . Thus,

$I(z) = {({\sigma _0}/2)^{\frac{1}{2}}}\int_0^z {K\!{{(z)}^{ - \frac{1}{2}}}{\rm d}z}, $

(A7)

where σ₀ is a factor in I (z) and ${\sigma _0} = N\sin (\alpha){\Pr ^{ - \frac{1}{2}}}$ . In this way, Eq. (A6) generalizes to:

$\begin{array}{l}{u_0} = - C\mu \exp (- I)\sin (I), \\{\theta _0} = C\exp (- I)\cos (I).\end{array}$

(A8)

Note that the term ${\rm d}K\!(Z)/{\rm d}z$ is neglected.

A2. First-order correction

Next, the first-order correction is analyzed. The corresponding system at the order of ε¹ is

$\begin{split}& g\frac{{{\theta _1}}}{{{\Theta _0}}}\sin (\alpha) + K\Pr \frac{{{{\rm d}^2}{u_1}}}{{{\rm d}{z^2}}} = 0, \\& - (\Gamma {u_1} + \frac{{{\rm d}{\theta _0}}}{{{\rm d}z}}{u_0})\sin (\alpha) + K\frac{{{{\rm d}^2}{\theta _1}}}{{{\rm d}{z^2}}} = 0.\end{split}$

(A9)

The unknowns are $({u_1}, {\theta _1})$ . We assume that ${\rm d}K\!(Z)/{\rm d}z{\rm{ = 0}}$ and ${{\rm d}^2}K\!/{\rm d}{z^2}{\rm{ = 0}}$ . Similar to Eq. (A4), Eq. (A9) also becomes the damped oscillator equation:

$\begin{array}{l}\frac{{{{\rm d}^4}{u_1}}}{{{\rm d}{z^4}}} + {N^2}\frac{{{{\sin }^2}(\alpha)}}{{\Pr K\!{{(z)}^2}}}{u_1} = - \frac{{g{{\sin }^2}(\alpha)}}{{{\Theta _0}K\!{{(z)}^2}\Pr }}{u_0}\frac{{{\rm d}{\theta _0}}}{{{\rm d}z}}, \\\frac{{{{\rm d}^4}{\theta _1}}}{{{\rm d}{z^4}}} + {N^2}\frac{{{{\sin }^2}(\alpha)}}{{\Pr K\!{{(z)}^2}}}{\theta _1} = \frac{{\sin (\alpha)}}{{K\!(z)}}\frac{{{{\rm d}^2}}}{{{\rm d}{z^2}}}({u_0}\frac{{{\rm d}{\theta _0}}}{{{\rm d}z}}).\end{array}$

(A10)

Equation (A10) comprises inhomogeneous ordinary differential equations whose solutions consist of a general solution and a special solution; the knowns are $({u_0}, {\theta _0})$ . The differences between Eqs. (A10) and (A5) are the boundary conditions and forcing term on the right-hand side. Thus, the form of the general solution of Eq. (A10) is similar to Eq. (A5). The forcing terms satisfy the form $f(z) = {e^{\lambda z}}[{P_l}(z)\cos \omega z + {P_n}(z)\sin \omega z]$ , where ${P_l}(z)$ and ${P_n}(z)$ are the l-order polynomial and n-order polynomial, respectively. Thus, the form of the special solution is ${u^*} = {z^k}{e^{\lambda z}}[R_m^{(1)}(z)\cos \omega z + R_m^{(2)}(z)\sin \omega z]$ , where u^* is similar to θ^*; and $R_m^{(1)}(z)$ and $R_m^{(2)}(z)$ are the m-order polynomials ( $m = \max \{ l, n\} $ ). The method of undetermined coefficients is applied to obtain the solution of Eq. (A10), which reads:

$\begin{aligned}{u_1} & = {u_A}\exp (- I)[ - \frac{1}{3}\sin (I) + \frac{2}{{15}}\cos (I)]\\ & + {u_A}\exp (- 2I)[\frac{1}{{30}}\sin (2I) - \frac{1}{{30}}\cos (2I) - \frac{1}{{10}}], \\ {\theta _1} & = {\theta _A}\exp (- I)[ - \frac{1}{{15}}\sin (I) - \frac{1}{6}\cos (I)] \quad \quad \quad (\rm {A}11)\\ & + {\theta _A}\exp (- 2I)[\frac{1}{{15}}\sin (2I) + \frac{1}{{15}}\cos (2I) + \frac{1}{{10}}].\end{aligned}$

Here, ${u_A} = {({\sigma _0}/2)^{\frac{1}{2}}}{C^2}(\mu /\Gamma)K{(z)^{ - \frac{1}{2}}}$ , ${\theta _A} = {(2/{\sigma _0})^{\frac{1}{2}}}{C^2}\mu \sin (\alpha)K{(z)^{ - \frac{1}{2}}}$ .

Thus, the total WKB approximate solutions, up to the first order, are as follows:

$\begin{array}{l}u = {u_0} + \varepsilon {u_1}, \\\theta = {\theta _0} + \varepsilon {\theta _1}.\end{array}$