Parameter estimation for land-surface models using Neural Physics

We propose a novel inverse-modelling approach that estimates the parameters of a simple land-surface model (LSM) by assimilating data into a differentiable, physics-based forward model formulated using convolutional operations. The governing equations are expressed within the Neural Physics framework, allowing direct gradient-based optimisation of time-dependent parameters without the need to derive and maintain adjoint formulations. The model parameters are estimated by minimising the mismatch between model predictions and synthetic or observational data. Although differentiability is enabled through machine-learning libraries, the forward model itself remains entirely physics-based and neither the forward model nor the parameter estimation procedure involve training.

To evaluate the approach, we first generate synthetic observations of soil temperature by running the forward model with known parameter values and subsequently treat these parameters as unknown in an inverse problem. We show that observations of soil temperature at a single depth are insufficient to reliably constrain the model parameters. Using observations at two depths, however, does yield reliable parameter estimates, although the individual contributions of latent and sensible heat fluxes cannot be distinguished.

We also apply the approach to urban flux tower data from Phoenix, United States, and show that the thermal conductivity, volumetric heat capacity and the combined sensible-latent heat transfer coefficient can be reliably estimated whilst using an observed value for the effective surface albedo. The resulting model accurately predicts the outgoing longwave radiation, conductive soil fluxes and the combined sensible-latent heat fluxes, demonstrating that the Neural Physics framework can be used to accurately determine the parameters of the particular LSM used here. This model is intentionally simple and does not include, for example, a subsurface moisture model. This simplicity facilitates an exploration of parameter identifiability, confounding and equifinality.

1 Introduction

The exchange of momentum and thermal energy between the land surface and the atmosphere is of great importance to atmospheric processes and thermal comfort (Oke, 2017). At present, such processes are typically modelled with land-surface models (LSMs) using meteorological forcing data and various parameters associated with the site. Many LSMs have been developed over the years, each differing in how they handle the urban surface, and in their modelling of various morphological, vegetation, hydrological and anthropogenic processes. Modern LSMs are capable of closely replicating the heterogeneous surface structure of cities (Lee and Lee, 2020; Krayenhoff and Voogt, 2007), and capturing the effects of various urban processes, such as urban heat storage and vegetation dynamics, using various soil, hydrological or vegetation models (Lipson et al., 2017; Arsenault et al., 2018).

However, larger complexity of LSMs does not always translate to better predictions (Grimmond et al., 2010, 2011; Lipson et al., 2024). The Urban-PLUMBER project (Lipson et al., 2024) evaluated the performance of 30 LSMs in estimating surface fluxes at a suburban site in Melbourne, Australia. The study found that models with highly complex urban schemes often do not perform as well as models with simpler urban schemes, as their performances are penalised by having relatively simple representation of hydrological/vegetation models. Most complex models were also unable to make efficient use of site-specific information.

Besides capturing the complexity of land-surface interactions, another key challenge in modelling urban surface fluxes lies in estimating (or calibrating) the relevant input parameters (Lipson et al., 2017; Chaney et al., 2016). LSMs typically have different sensitivities to changes in different parameters, making it difficult to perform data assimilation or improve model accuracy by manually updating parameter values, as a slight alteration can lead to large and unexpected changes in the final output (Massoud et al., 2019). Grimmond et al. (2011) found that poor estimation of parameter values can result in a significant decline in performance of LSMs. It is hence critical to find a reliable method of parameter estimation.

Raoult et al. (2024) presented an overview of the popular techniques for estimating the parameters of LSMs, which include statistical methods and inverse modelling. Statistical methods involve estimating the probability distributions of model parameters based on prior knowledge and observed data. Some examples include the Markov Chain Monte Carlo (MCMC) method (Hastings, 1970), which draws samples from posterior distributions to determine the likely values of parameters. However, statistical methods can be computationally intensive, and therefore have been mainly used in computationally inexpensive LSMs such as Simplified PnET (SIPNET) ecosystem model (Fer et al., 2018), or to estimate isolated processes (e.g. plant respiration; Jones et al., 2024). Inverse modelling employs optimisation algorithms, such as the gradient-descent method or genetic algorithms, to determine the optimal set of parameters that minimises the error between the model output and observational data. Such minimisation algorithms have been shown to have better overall optimisation efficiency than statistical methods (Raoult et al., 2024). However, inverse modelling algorithms, especially those based on the gradient-descent method, are prone to become stuck in local minima, leading to non-unique solutions that are dependent on the initial conditions (Kuppel et al., 2014; Raoult et al., 2016).

A number of researchers have begun exploring the coupling of machine-learning models with differentiable physics-based models to construct differentiable hybrid models (Shen et al., 2023). In these approaches, neural networks are embedded within a physics-based framework, which is often trained end-to-end, to infer unknown parameters or represent unknown physical relationships. Several studies focus on learning physical processes or closure relationships. For example, Fang and Gentine (2024) proposed a hybrid model to learn the impact of soil water stress on photosynthesis and evapotranspiration, outperforming traditional terrestrial biosphere models that use empirical soil water stress parametrisations. In a similar vein, ElGhawi et al. (2023) represent leaf surface and aerodynamic resistances using two neural networks within an evapotranspiration model, while Norouzi et al. (2026) use neural networks to separate adsorption and capillary effects on the soil water retention curve. Other work embeds neural networks within land-surface models, learning a snow process representation (Deck et al., 2026) and constitutive relationships for frozen soil (Ren et al., 2026). Hybrid models have also been applied to inverse problems and parameter estimation. For example, Ouyang et al. (2025) re-implemented the Xin'anjiang (XAJ) hydrological model in PyTorch representing the discretised equations for run-off and evapotranspiration by neural networks. An inverse problem is solved to learn the weights of an LSTM model which predicts the 15 parameters of the XAJ model. With limited data, their model gave better streamflow predictions than manually calibrating XAJ. A further example is JAX-CanVeg (Jiang et al., 2025), in which a multi-layer land-surface canopy model is augmented with a neural-network–based closure for leaf relative humidity, improving predictions of latent heat and net ecosystem exchange across four flux tower sites. These studies demonstrate the potential of differentiable hybrid models to improve the performance and calibration of LSMs. However, challenges still remain for these differentiable hybrid methods, including sensitivity to training data; tendency of predictions to become physically inconsistent; and computational cost, arising when physics-based models require many forward simulations for parameter optimisation, particularly in end-to-end differentiable settings (Fang and Gentine, 2024; Dukes et al., 2026).

Figure 1 Definition sketch of the domain with the governing equations and boundary conditions. Shown on the right are the nodes of the discretised problem.

A related class of differentiable but fully physics-based models has also been developed (Zhu et al., 2021; Bezgin et al., 2023; Chen et al., 2026). One such training-free approach, called Neural Physics, capitalises on the equivalence between discretisation stencils and convolutional operators, allowing partial differential equations and their solvers to be expressed as the convolutional layers of neural networks with analytically determined weights. While differentiability can also be achieved using general-purpose automatic differentiation frameworks (e.g. Julia, or JAX and NumPy-based implementations), Neural Physics provides an alternative to these by formulating numerical discretisations of physical laws using machine-learning primitives such as convolutional operations. The motivation for using Neural Physics in this work extends beyond differentiability. Among its principal advantages are its established and versatile framework, differentiability, reduced implementation complexity and support for a broad range of numerical methods and hardware architectures. The forward model has been applied to a wide range of problems including single-phase flow benchmarks (Chen et al., 2026), multiphase flow (Chen et al., 2024), shallow water dynamics (Chen et al., 2025) and applications in nuclear engineering (Phillips et al., 2023). It is also platform agnostic, being able to execute on CPUs, GPUs and AI accelerators (Cerebras, 2026; Graphcore, 2026) without modification. The end-to-end differentiability of Neural Physics enables the computation of sensitivities of objective functions with respect to model parameters, allowing standard optimisation tools from machine-learning libraries to be used for parameter estimation. Furthermore, the use of convolutional operations and other functions from machine-learning libraries abstracts away much of the implementation complexity associated with accelerator portability and adjoint development. Neural Physics and the associated NN4PDEs code also supports a broader class of numerical methods, including higher-order discretisations (Phillips et al., 2024); finite volume (Xue et al., 2026), finite element and finite difference discretisations (Chen et al., 2026) and multigrid solvers (Phillips et al., 2023). The latter can be exactly represented as a U-Net architectures (Phillips et al., 2023), which require no training and reduce the number of iterations needed for convergence compared with classical Jacobi-based schemes. Although not used in this work, the multigrid solver would reduce the computational burden in more complex LSMs. We hope to exploit these advantages for a range of inverse problems in future work.

The NN4PDEs code has recently been used to determine the conductivity of the subsurface given a set of steady-state observations (Li et al., 2026). Here, we extend Li et al.'s methodology to a transient application to determine the parameters of a simple land-surface model, and to study whether a parameter can be uniquely determined from the observations (parameter identifiability), whether the effects of different parameters on the model response can be distinguished from one another (confounding) and whether many different parameter combinations result in similar parameter values (equifinality). Although developed independently, our proposed AI-enabled physics-based approach is similar to the NoahPy model which has been developed to improve the representation of permafrost in LSMs (Tian et al., 2026). In that work, numerical discretisations are formulated as a recurrent neural network, enabling parameter estimation from observations at multiple soil depths. The model estimates four parameters using the Adam optimiser and was shown to converge more efficiently than traditional gradient-free optimisation methods.

The key contribution of this paper is the formulation of a land-surface model within the Neural Physics framework, an end-to-end, fully differentiable, physics-based system that enables efficient parameter estimation. We intentionally chose a simple LSM that still includes a surface energy balance and time-dependence, in order to isolate parameter identifiability, confounding and equifinality. The surface energy balance equation is coupled to a soil layer (Sect. 2.1), but does not incorporate a subsurface moisture model or explicitly represent buildings. By expressing the governing equations as convolutions within a machine-learning framework using Neural Physics (Sect. 2.2), we eliminate the need to derive and maintain adjoint models, which is a significant practical barrier in traditional LSM calibration workflows. Instead, we demonstrate that our approach can be applied directly using gradient-based optimisation to estimate parameters in time-dependent problems (Sect. 2.3), a process that can otherwise be tedious (Lipson et al., 2024). Using a synthetic dataset for which the correct input parameters are known, the sensitivity to the initial estimate and dependencies between parameters are investigated (Sect. 3.1). We show that simultaneous measurements at two depths are needed to obtain reliable estimates for the soil properties independent of the initial estimates (Sect. 3.2). Finally, we use the model to determine the site parameters for a 150 h period of the flux data recorded in the city of Phoenix, United States and evaluate the model on a subsequent 50 h period (Sect. 3.3). This approach provides insight into the degree of underdetermination in the system (e.g., identifiable parameter combinations and the role of observational configuration), which are often difficult to assess using conventional calibration methods. We comment on this in the results and conclusions (Sects. 3 and 4).

2 Methods

2.1 Land-surface model

The land-surface model considered comprises a soil layer coupled to a surface energy balance at the soil-atmosphere interface (Fig. 1). The surface energy balance equation is given by

$$\begin{array}{}\text{(1)}& \begin{array}{rl}G\left(t\right)=& (\mathrm{1}-\mathit{\alpha })K^{\downarrow }\left(t\right)+L^{\downarrow }\left(t\right)-\mathit{ϵ}\mathit{\sigma }{T}_s^{\mathrm{4}}\left(t\right)\\ & -h(\mathrm{1}+{\mathit{\beta }}^{-\mathrm{1}})\left(T_{\mathrm{s}}\right(t)-T_{\mathrm{a}}(t\left)\right),\end{array}\end{array}$$

where G is the ground heat flux, K^↓ and L^↓ are the downwelling shortwave and longwave radiation, respectively, α is the surface albedo, ϵ is the emissivity, σ = 5.67 × 10⁻⁸ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{4}}$ is the Stefan–Boltzmann constant, h is the heat transfer coefficient, T_a is the atmospheric temperature and T_s is the absolute surface temperature. The quantity β is the Bowen ratio, defined as

$$\begin{array}{}\text{(2)}& \mathit{\beta }={\displaystyle \frac{Q_{\mathrm{H}}}{Q_{\mathrm{E}}}},\end{array}$$

where Q_H is the sensible heat flux and Q_E is the latent heat flux. The enthalpy equation for the soil is given by

$$\begin{array}{}\text{(3)}& C{\displaystyle \frac{\partial T}{\partial t}}=\mathit{\lambda }{\displaystyle \frac{{\partial }^{\mathrm{2}}T}{\partial x^{\mathrm{2}}}},\end{array}$$

where x is a downward coordinate starting at the surface, T(x,t) is the soil temperature, C is the volumetric heat capacity and λ is the thermal conductivity. The boundary conditions are given by

$$\begin{array}{}\text{(4)}& -\mathit{\lambda }{\displaystyle \frac{\partial T}{\partial x}}(x=\mathrm{0},t)=G\left(t\right),T(x=d,t)=T_{\mathrm{b}},\end{array}$$

where d is the vertical extent of the soil layer and T_b is the fixed temperature at the bottom boundary of the soil layer. The initial condition is given by

$$\begin{array}{}\text{(5)}& T(x,t=\mathrm{0})=T_{\mathrm{b}}.\end{array}$$

2.2 Forward model

2.2.1 Discretisation

Time and space are discretised uniformly as tⁿ=nΔt and x_i=iΔx, respectively, where Δt is the time increment, $\mathrm{\Delta }x=d/N$ is the grid increment and N is the number of intervals. Using a central finite difference discretisation in space and a backward Euler discretisation in time, Eq. (3) becomes

$$\begin{array}{}\text{(6)}& \begin{array}{rl}& -r{T}_{i+\mathrm{1}}^{n+\mathrm{1}}+(\mathrm{1}+\mathrm{2}r)T_i^{n+\mathrm{1}}-r{T}_{i-\mathrm{1}}^{n+\mathrm{1}}=T_i^n\\ & \forall i\in \mathit{\{}\mathrm{1},\mathrm{\dots },N-\mathrm{1}\mathit{\}},\end{array}\end{array}$$

where $T_i^n=T(x_i,t^n)$ and $r=\mathit{\lambda }\mathrm{\Delta }t/\left(C\mathrm{\Delta }{x}^{\mathrm{2}}\right)$. The backward Euler scheme is chosen here because it is absolutely stable and allows for large time-steps, which is useful for the computational efficiency of the inverse model.

The surface boundary condition in Eq. (1) is imposed by choosing a second order central discretisation for the surface flux and combining this with the PDE at the boundary (Bradle, 2007), which after discretisation with an Euler backward time-integration method becomes:

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{T_{\mathrm{0}}^{n+\mathrm{1}}-T_{\mathrm{0}}^n}{\mathrm{\Delta }t}}={\displaystyle \frac{\mathrm{2}\mathit{\lambda }}{C\mathrm{\Delta }x}}\left({\displaystyle \frac{T_{\mathrm{1}}^{n+\mathrm{1}}-T_{\mathrm{0}}^{n+\mathrm{1}}}{\mathrm{\Delta }x}}+{\displaystyle \frac{G\left(t^{n+\mathrm{1}}\right)}{\mathit{\lambda }}}\right),\end{array}$$

where the ground heat flux G(tⁿ⁺¹) uses a mixed discretisation for the outgoing longwave radiation term:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}G\left(t^{n+\mathrm{1}}\right)=& (\mathrm{1}-\mathit{\alpha })K^{\downarrow }\left(t^{n+\mathrm{1}}\right)+L^{\downarrow }\left(t^{n+\mathrm{1}}\right)-\mathit{ϵ}\mathit{\sigma }{\left(T_{\mathrm{0}}^{n+\mathrm{1}}\right)}^{\mathrm{4}}\\ & -h(\mathrm{1}+{\mathit{\beta }}^{-\mathrm{1}})\left(T_{\mathrm{0}}^{n+\mathrm{1}}-T_{\mathrm{a}}\left(t^{n+\mathrm{1}}\right)\right).\end{array}\end{array}$$

The discretised bottom boundary condition is given by

$$\begin{array}{}\text{(9)}& T_N^{n+\mathrm{1}}=T_{\mathrm{b}}.\end{array}$$

In order to avoid inverting a matrix, the system Eqs. (6)–(9) is solved using Jacobi iterations (Greenbaum, 1997). For the PDE discretisation in Eq. (6), this implies an iteration scheme of the form

$$\begin{array}{}\text{(10)}& T_i^{n+\mathrm{1},\left(k\right)}={\displaystyle \frac{T_i^n+r\left(T_{i+\mathrm{1}}^{n+\mathrm{1},(k-\mathrm{1})}+T_{i-\mathrm{1}}^{n+\mathrm{1},(k-\mathrm{1})}\right)}{\mathrm{1}+\mathrm{2}r}},\end{array}$$

where k denotes the iteration number. The surface boundary condition (Eq. 8) is solved by linearising the outgoing longwave radiation term (i.e. $\mathit{ϵ}\mathit{\sigma }(T_{\mathrm{0}}^{n+\mathrm{1},(k-\mathrm{1})}{)}^{\mathrm{3}}{T}_{\mathrm{0}}^{n+\mathrm{1},\left(k\right)}$).

Figure 2 A schematic diagram of the inverse problem. Initial guesses for a set of parameters are fed into the forward model (here, a physics-based model expressed as convolutions using Neural Physics) resulting in a solution for temperature (blue arrows). The mismatch between temperature and the observations is minimised with gradient descent and the parameters are updated (orange arrows). This process continues until convergence.

2.2.2 Implementation of Forward Model with Neural Physics

The NN4PDEs code of Neural Physics uses a convolutional layer with pre-determined weights to implement the finite difference scheme described in Sect. 2.2.1. Equation (10) can be re-written as

$$\begin{array}{}\text{(11)}& T_i^{n+\mathrm{1},\left(k\right)}={\displaystyle \frac{T_i^n}{\mathrm{1}+\mathrm{2}r}}+\mathit{w}\ast {\mathit{T}}^{n+\mathrm{1},(k-\mathrm{1})}{\mathrm{|}}_i\end{array}$$

in which ∗ represents the discrete convolution, defined as

$$\begin{array}{}\text{(12)}& \mathit{w}\ast \mathit{T}{\mathrm{|}}_i\equiv \sum _{j=-\mathrm{1}}^{\mathrm{1}}{w}_j{T}_{i+j},\end{array}$$

where ${\mathit{T}}^{n+\mathrm{1},\left(k\right)}$ is a vector containing the nodal values of temperature at iteration k, and the vector w has values $[r/(\mathrm{1}+\mathrm{2}r),\mathrm{0},r/(\mathrm{1}+\mathrm{2}r\left)\right]$. Equations (11) and (12) apply to all interior nodes $i\in \mathit{\{}\mathrm{1},\mathrm{\dots },N-\mathrm{1}\mathit{\}}$. The discretisation in Eq. (11) is analogous to a convolutional layer in a neural network where the filter or kernel has weights given by w.

The PyTorch machine-learning library is used to build the forward model. To initialise the neural network, the convolution function and solver for the implicit time-stepping scheme are created. Within the solver, Jacobi iteration for interior nodes is performed by repeated application of the convolution function (Eq. 11). In each iteration, the boundary temperatures $T_{\mathrm{0}}^{n+\mathrm{1}}$ and $T_N^{n+\mathrm{1}}$ are updated and appended to the vector for interior nodes to give the full temperature vector. Once the error, which is taken as the maximum absolute difference between consecutive approximations, is smaller than the defined tolerance, the latest temperature vector is saved and used as the new guess for temperature at the next time level. The iteration procedure is repeated for subsequent time steps until the final time is reached. Convergence checks have been performed to ensure that the numerical scheme is stable and consistent (Huang, 2024).

2.3 Inverse model

The inverse problem aims to find a set of parameters p that could have given rise to a set of observations. First, an initial guess is made for the values of the parameters. Second, the forward model is solved through time. Third, a data mismatch functional is calculated, based on the difference between the model output and the observations:

$$\begin{array}{}\text{(13)}& J\left(p\right)={\displaystyle \frac{\mathrm{1}}{\left|\mathcal{I}\right|\left|\mathcal{N}\right|}}\sum _{i\in \mathcal{I}}\sum _{n\in \mathcal{N}}{\left(T\left(x_i,t^n;p\right)-\stackrel{\mathrm{̃}}{T}\left(x_i,t^n\right)\right)}^{\mathrm{2}},\end{array}$$

where $T(x_i,t^n;p)$ denotes the model prediction for parameters p, $\stackrel{\mathrm{̃}}{T}(x_i,t^n)$ denotes the observations, ℐ is a set of indices which indicate where the observations are taken and 𝒩 is a set of indices indicating when there are observations. The cardinality or size of these sets is represented by $\left|\mathcal{I}\right|$ and $\left|\mathcal{N}\right|$. The final stage is to minimise this mismatch with respect to the parameters:

$$\begin{array}{}\text{(14)}& {\mathit{p}}^{\ast }={\text{argmin}}_{p}J\left(p\right),\end{array}$$

where p^∗ is the optimised parameter vector. Classical optimisation methods such as stochastic gradient descent can be used to minimise the data mismatch and calculate an iterative update of the parameters based on the gradient (Fletcher, 1987; Nocedal and Wright, 2006). We use a modified version of stochastic gradient descent known as Adam, which scales the step length with estimates of the mean and variance of the gradients during the iterative process, in order to improve the speed and stability of the algorithm (Kingma and Ba, 2015). Each parameter has an associated step length, the initial value of which set according to the order of magnitude of the parameter. C has the largest initial step length of 10⁴, T_b has an initial step length of 1, and the remaining parameters have initial step lengths of 0.1. We note that in a machine-learning context, step length is more often referred to as learning rate. However, as we are an using optimisation method to minimise a functional without training a neural network, we use the term step length, in line with classical texts on optimisation (Fletcher, 1987; Nocedal and Wright, 2006).

The optimisation procedure involves solving the forward model, minimising the mismatch between the model output and the observations using the Adam method to give an improved estimation of parameter values that are consistent with the observations. These steps are repeated until the maximum number of iterations (defined as 150) is reached, or the solution has converged. The solution is considered converged at iteration n if both the relative change in the data mismatch functional J is less than 0.01% and the relative changes in all parameters p are less than 0.1%. This procedure is illustrated in Fig. 2. The optimisation process is expedited using Python's multi-processing library and the coding platform Google Colab. Each active session in Google Colab uses 2 CPUs, allowing a maximum of 10 trials to be run simultaneously with 5 active sessions.

Table 1 Input parameter ranges and values for the land-surface model. The “Value” column contains the parameter values used to create the synthetic dataset. Data sourced from Brutsaert (2005), Oke (2002), Laloui and Rotta Loria (2020), and Oke (2017).

3 Results

3.1 Observations at a single depth

To investigate the performance of the inverse model, we produce synthetic data using the forward model, so that the correct parameter values are known. We initially pursue a “naive” approach, which is to simply run the inverse model to determine the six main input parameters, namely p={α, β, h, C, λ, T_b}, without consideration of the physics of the underlying model.

Table 1 shows all the inputs that are needed to obtain model output, including the parameter values used to create the synthetic dataset. There is relatively little variation in the emissivity ϵ of longwave radiation, and we will therefore simply assume a fixed value of 0.95. The soil thickness must be chosen such that it is larger than the penetration thickness $d_{\mathit{\omega }}=\sqrt{\mathrm{8}{\mathit{\pi }}^{\mathrm{2}}\mathit{\lambda }/\left(C\mathit{\omega }\right)}$ (Carslaw and Jaeger, 1959), where ω is the angular frequency of a diurnal cycle. For the synthetic dataset, d_ω = 0.628 m, which means that d = 1 m is sufficiently deep to assume that the temperature remains unchanged at that depth. The soil layer is discretised with N=100 intervals. A large number of intervals is used to ensure that spatial discretisation errors will not affect the results.

The shortwave radiation varies according to

$$\begin{array}{}\text{(15)}& K^{\downarrow }\left(t\right)=S_{\mathrm{b}}\cos \left(\mathit{\theta }\right(t\left)\right).\end{array}$$

Here, θ is the solar zenith angle, which is the angle between the direction of the sun's rays and the axis perpendicular to the ground, and S_b is the direct solar irradiance when θ = 0. The solar zenith angle θ is paramaterised as

$$\mathit{\theta }\left(\widehat{t}\right)=\left\{\begin{array}{ll}-\frac{\mathit{\pi }}{\mathrm{2}}& \text{if }\widehat{t}⩽\mathrm{7};\\ -\frac{\mathit{\pi }}{\mathrm{2}}+\frac{\mathit{\pi }}{\mathrm{14}}(\widehat{t}-\mathrm{7})& \text{if }\mathrm{7}<\widehat{t}<\mathrm{21};\\ \frac{\mathit{\pi }}{\mathrm{2}}& \text{if }\widehat{t}⩾\mathrm{21}\end{array}\right.$$

where $\widehat{t}$ represents the hour in the day and is given by $t/\mathrm{3600}\text{mod}\mathrm{24}$. It is assumed that the sun rises at 07:00 LT, and sets at 21:00 LT Before 07:00 a.m. and after 21:00 LT, $\mathit{\theta }=-\frac{\mathit{\pi }}{\mathrm{2}}$ which gives zero K^↓. Between 07:00 LT and 21:00 LT, K^↓ varies sinusoidally and peaks at 14:00 LT. The longwave radiation L^↓ and atmospheric temperature T_a are assumed constant. The values assumed for S_b, L^↓ and T_a are stated in Table 1. The simulation is run for 100 h and the temperatures recorded over time at a depth of 5 cm are used as the observations to determine all six parameters with the inverse model.

In order to investigate whether the parameters to be determined are sensitive to the initial parameter values, 50 trials are performed, each starting from different initial parameter values. The parameters for trial i are initialised using $p_i=p_i^{min}+R(p_i^{max}-p_i^{min})$, where $p_i^{min}$ and $p_i^{max}$ are the minimum and maximum parameter values defined such that the range of the initial values is sufficiently large, but still lies within the typical ranges of the parameter (see Table 1), and R is a random number taken from a uniform probability distribution.

Figure 3 Scatter plots of initial vs final values for the synthetic data set of (a) surface albedo α, (b) heat transfer coefficient h, (c) Bowen ratio β, (d) $h(\mathrm{1}+{\mathit{\beta }}^{-\mathrm{1}})$, (e) Initial/interior temperature T_b, (f) volumetric heat capacity C, (g) thermal conductivity λ and (h) μ compared with their true values (black dotted line). The blue circles assimilate data from depth of 0.05 m and the orange circles assimilate data from two heights, 0.05 and 0 m (i.e., the surface).

The results from the 50 trials are shown in Fig. 3 in the blue circles. Scatter plots are shown for each parameter, with its initial value on the x-axis and its final value on the y-axis. Clearly, the ideal outcome is that the same final value is reached regardless of the initial parameter value, but this is not always the case. The results indicate that the algorithm is able to determine T_b and α. However, the value of h is strongly dependent on its initial value value. The parameters β, λ and C have a substantial spread in their optimal values, in a seemingly uncorrelated manner with their initial value. Hence, the conclusion from the “naive” parameter estimation is that it is not possible to obtain a reliable estimation of the parameter values by assuming the parameters are independent and using a single observation of the soil temperature profile.

The reason why this happened is because an optimiser converges to a global minimum only when the optimisation problem is convex. This condition will be violated when two parameters are dependent on each other. In particular, the heat transfer coefficient h and the Bowen ratio β occur only through the product $h(\mathrm{1}+{\mathit{\beta }}^{-\mathrm{1}})$ in Eq. (1). As a result, the optimisation problem is non-unique. Figure 4a shows the final values for h and β plotted against each other, showing a clear relation between the two quantities. Figure 3d shows that the product $h(\mathrm{1}+{\mathit{\beta }}^{-\mathrm{1}})$ has much less spread than h and β individually. Thus, the optimisation should use the parameter $h(\mathrm{1}+{\mathit{\beta }}^{-\mathrm{1}})$ instead of h and β separately.

Figure 4 Covariance plots of final parameter values for the synthetic dataset of (a) β against h; and (b) λ against C.

Furthermore, it is well known that thermal inertia is typically characterised using the thermal admittance μ=λC (Oke, 2017), which characterises the ability of a material to absorb and release heat. This suggests that there might also be a correlation between C and λ. Figure 4b shows the final values for C and λ plotted against each other, once more showing a clear relation between the two quantities. However, Fig. 3h shows that even though C and λ are clearly correlated, there is considerable uncertainty in μ, and it is not possible to replace C and λ for μ.

To conclude this section, we have learnt that a “naive” parameter estimation approach does not result in reliable parameter values. Firstly, the parameters can be correlated to each other, and we found that it is not possible to determine h and β individually. However, even though it is shown that C and λ are strongly correlated, the quantity μ does not capture this behaviour and we thus conclude that it is impossible to obtain reliable parameter values using temperature observations at a single soil depth.

3.2 Observations at two depths

In this section, parameter estimation using two sets of observations of temperatures is explored, once more using the synthetic data. In addition to T (t; x = 0.05 m), the surface temperature T_s will be used, since this is where the temperature amplitudes are largest. Instead of changing the parameter estimation strategy and using $h(\mathrm{1}+{\mathit{\beta }}^{-\mathrm{1}})$ as a model parameter, we will keep using h and β, but with the understanding that we should only look whether the quantity $h(\mathrm{1}+{\mathit{\beta }}^{-\mathrm{1}})$ has a value that is independent of its initial value.

Figure 3 shows that by using observations at two depths, it becomes possible to reliably estimate parameter values. The parameter values obtained from the inverse model are α = 0.199 ± 0.010, T_b = 292.99 ± 0.05 K, λ = 0.81 ± 0.11 $\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$, C = 2.24 ± 0.32  $\mathrm{MJ}{\mathrm{m}}^{-\mathrm{3}}{\mathrm{K}}^{-\mathrm{1}}$ and $h(\mathrm{1}+{\mathit{\beta }}^{-\mathrm{1}})$ = 24.90 ± 0.78 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}$. Obtained from 50 trials, the mean values of the parameters are within 2 % of their real value (see Table 1). Note that C and λ are the only quantities that have outliers up to 15 % larger than the mean value.

3.3 Application to Phoenix dataset

In order to apply the inverse model to real data, we use data from a flux tower managed by the CAP LTER programme. The tower is located in a residential suburban of Maryvale, in the city of Phoenix, United States. The forcing data (incoming shortwave radiation K^↓, re-emitted longwave radiation L^↓ and atmospheric temperature T_a) were collected over the calendar year of 2012 (Chow et al., 2014; Chow, 2017). The data was then bias-corrected and gap-filled by the Urban-PLUMBER project (Lipson et al., 2022b, 2024). Phoenix is known for its arid, desert-like climate. This minimises the effect of the moisture content of soil on heat flux exchange, which is not included in the equations of the current model. We extracted 200 h of flux data from 1 May 2012 00:00 (UTC−7 h) to 9 May 2012 08:00 (UTC−7 h). The first 150 h of soil temperature measurements are used for data assimilation to estimate the parameters, whilst the remaining 50 h of data are used for evaluation of the model and its parameters. A longer time series is used than the synthetic data for assimilation to account for the error in assuming a uniform initial temperature profile, which is unlikely to occur in real life. For detailed information about the measurement instruments in the West Phoenix flux tower, see Appendix A.

For this dataset, we aim to estimate four parameters, namely T_b, $h(\mathrm{1}+{\mathit{\beta }}^{-\mathrm{1}})$, λ and C. The value of α will not be determined using optimisation, as we found that the parameter values could not be determined reliably without specifying α, presumably since α controls the absorbed shortwave radiation which is the primary driver of the system. However, this value can be inferred directly from measured incoming and outgoing short-wave radiations. In this case, it is taken to be α=0.172, an average mid-day value estimated for the site by Lipson et al. (2022b, 2024). Once more, a total of 50 trials are conducted using initial parameter values that are randomly generated using the procedure discussed in Sect. 3.1. Table 2 shows the final parameter values from assimilating to temperature measurements at 5 and 15 cm, which are taken to be the mean of the optimised values from all 50 trials.

Table 2 Estimated parameters for the West Phoenix flux tower data set for the period of 1 May 2012 00:00 (UTC−7 h)–9 May 2012 08:00 (UTC−7 h) from assimilating soil temperatures at 0.05 and 0.15 m.

Comparing results of inverse analysis using real data to that of synthetic data when assimilated to measurements at two depths (Fig. 3), λ has a larger distribution in its optimised values, indicating higher uncertainties. The estimated thermal conductivity λ and volumetric heat capacity C are within the range of typical values for all three types of soils. The estimated boundary temperature T_b is also reasonable for an arid climate during summer (Table 1).

Figure 5 Comparison of modelled and recorded data for the West Phoenix flux tower dataset. (a) Temperatures at 0.05 and 0.15 m depth. (b) Ground heat flux at 0.02 m depth. (c) Outgoing long-wave radiation. (d) Combined sensible-latent heat flux.

In order to compare the modelled fluxes with the observations, the forward model is run using the mean of the obtained parameter values. Figure 5a shows good agreement of the model output for the two temperature time series with the final 50 h of observations, which are not used for the optimisation. More interesting are comparisons with other observed quantities, such as the outgoing longwave radiation L^↑, the combined sensible-latent heat flux Q_E+Q_H and conductive soil fluxes in the soil. Figure 5b shows the conductive soil flux measurements in gravel and sandy soil as observed together with the model outputs. For the model, the flux $G=-\mathit{\lambda }\partial T/\partial x$ is estimated using central differences. The modelled ground heat flux lies in between the ground heat flux measured in gravel and sandy soil, as it resembles the flux recorded in sandy soil more in the day, and flux recorded in gravel more at night. Figure 5c compares the observed upwelling longwave radiation L^↑ with the model output, once more showing good agreement between the two. It is not self-evident that this should be the case, since the long-wave radiation picks up information from buildings as well as the soil, and the surface temperature has a much larger amplitude than the temperatures deeper in the soil, thus allowing for extrapolation errors. Finally the combined sensible-latent heat fluxes Q_H+Q_E are calculated and compared with flux tower measurements in Fig. 5d. Here, the agreement is not as good as for the other surface energy balance components. We note that the differences between the observations and the model for the turbulent fluxes could be due to the difference in scale of the forcing data and target temperature measurements. Turbulent fluxes are measured at a height of 22.1 $\mathrm{m}\mathrm{a}.\mathrm{g}.\mathrm{l}.$ and represent an averaged value across the source area, which is highly variable as it changes with wind direction and atmospheric stability (Chow et al., 2014; Schmid et al., 1991), whilst soil temperature measurements are localised (taken only at one location beneath the flux tower).

4 Conclusions

We developed a land-surface model that enables end-to-end parameter estimation using gradient-based optimisation through Neural Physics, an approach which takes advantage of the tools from machine-learning libraries to provide differentiability whilst remaining entirely physics based. The approach also offers insights into parameter identifiability and confounding, and we provide a concrete, quantitative characterisation of equifinality for this application, identifying specific pathways through which it can be mitigated.

Results from a “naive” approach that assumes parameter independence show that temperature observations at a single depth are not sufficient to obtain reliable parameter values. In fact, the problem is shown to be ill-posed as multiple parameter combinations (e.g., h and β) collapse onto an identifiable product $h(\mathrm{1}+{\mathit{\beta }}^{-\mathrm{1}})$ and soil parameters (λ,C) exhibit strong correlations with a large spread in their individual estimates. We demonstrate that this non-identifiability is not intrinsic to the model alone, but depends on the observational configuration, as adding temperature observations at a second depth resolves this degeneracy and enables reliable recovery of thermal conductivity and heat capacity, independent of initial conditions. We also show that some confounding cannot be resolved without additional physics or data, for example, the sensible and latent heat fluxes remain inseparable without direct flux measurements. The model is then applied to observations of a 150 h period of the West Phoenix flux tower dataset and evaluated on a subsequent 50 h period of data. The estimated parameter values for the thermal conductivity, volumetric heat capacity and the combined sensible-latent heat transfer coefficient are physically plausible, and the agreement of the predictions for outgoing longwave radiation, conductive soil fluxes and combined sensible-latent heat fluxes with observations in the evaluation period is good.

In terms of parameter estimation, we have therefore learnt the following:

1. Data at two depths is required in order to obtain reliable estimates for thermal conductivity and heat capacity.

2. It is not possible to partition sensible and latent heat fluxes without specific measurements of either Q_H or Q_E.

3. It is necessary to specify the effective albedo of the site, as without it parameter values did not become independent of their initial estimate. We surmise that this is because the effective albedo is the primary forcing during day-time and all fluxes will be affected by this. Since this is one of the simplest quantities to measure, this is not considered a major limitation.

4. We recommend performing optimisation with at least 50 random initial parameter estimates to ensure the robustness of the obtained parameter values and to take the mean value over all the trials. However, our results suggest that it is possible to use a lower number.

Several challenges remain for future work. First, the LSM used to demonstrate the proposed approach is intentionally simple. The absence of a subsurface moisture model allows a detailed investigation of parameter identifiability and equifinality, but also motivates future work aimed at extending the framework to more realistic settings. In the present study, we therefore focus on a period from the West Phoenix flux tower dataset during which latent heat fluxes are negligible (i.e., no rainfall). Applying the approach to operational LSMs will require methods capable of estimating the larger number of parameters associated with more complex process representations. In addition, the inverse model is applied over a limited time period, avoiding the need to assume slowly varying parameter values. This assumption becomes particularly important when considering soil moisture, which is inherently time-dependent and typically requires multiple state variables for accurate representation (Grimmond et al., 2010). For longer time series, more sophisticated approaches may be required, potentially including the incorporation of a learned component. Although it is feasible to augment the current model with such methods, the parameter estimation results indicate that the present physics-based formulation is already capable of reproducing observations with reasonable accuracy.

The second challenge concerns practical adoption of the framework. LSMs are typically implemented in legacy languages such as Fortran, whereas machine-learning frameworks are commonly developed in Python, JAX or Julia. Therefore, applying the Neural Physics approach, or, indeed, adapting code to any framework with automatic differentiability, requires some degree of model refactoring. Although challenging, this process is becoming increasingly tractable with modern tools such as Github Copilot, Claude Code and ChatGPT (Zhou et al., 2024).

A third challenge is to deal with the non-smooth parameterisations arising from more complex LSMs. Many operational land-surface models include discontinuities arising from lookup tables, threshold-based parameterisations and conditional logic, which can complicate gradient-based optimisation in practice. Although modern automatic differentiation frameworks can handle piecewise-smooth operations, extending our approach to discontinuous or hybrid formulations is an important direction for future work and an active area of research.

Despite these challenges, the proposed framework offers a potential advantage when scaling to more complex systems. Dukes et al. (2026) identify differentiable hybrid physics–machine-learning models as a new and promising direction for parameter estimation in Earth System Models, but highlight the substantial computational burden associated with increasing model complexity. In particular, stability constraints often necessitate smaller time steps, leading to higher iteration counts and increased memory demands in inverse problems. In this context, the Neural Physics formulation provides a natural way of incorporating fully differentiable multigrid solvers through U-Net architectures (Chen et al., 2026, 2024; Phillips et al., 2023). Such approaches are well suited to reducing iteration counts, thereby mitigating a key computational bottleneck. Whilst not required for the simple case considered here, this capability represents a potential advantage of the proposed approach when applied to more complex LSMs.

Appendix A: West Phoenix flux tower data

Information about the West Phoenix flux data measurements is provided in Table A1.

(Lipson et al., 2022a, b, 2024)(Chow et al., 2014; Chow, 2017)

Table A1 Summary of West Phoenix flux data.

¹ Urban-PLUMBER gap-filled the forcing data using a combination of three methods: contemporaneous and nearby flux tower data where available, linear interpolation for small gaps (⩽2h), and bias-corrected ERA5 reanalysis meteorological data for larger gaps. ERA5 combines satellite, atmospheric and ground-based observations to produce globally consistent forcing data at 0.25° spatial and hourly temporal resolutions; ² Soil temperature measurements are only available from March 2012.