Towards a simple representation of chalk hydrology in land surface modelling

. Modelling and monitoring of hydrological processes in the unsaturated zone of chalk, a porous medium with fractures, is important to optimize water resource assessment and management practices in the United Kingdom (UK). However, incorporating the processes governing water movement through a chalk unsaturated zone in a numerical model is complicated mainly due to the fractured nature of chalk that creates high-velocity preferential ﬂow paths in the subsurface. In general, ﬂow through a chalk un-saturated zone is simulated using the dual-porosity concept, which often involves calibration of a relatively large number of model parameters, potentially undermining applications to large regions. In this study, a simpliﬁed parameterization, namely the Bulk Conductivity (BC) model, is proposed for simulating hydrology in a chalk unsaturated zone. This new parameterization introduces only two additional parameters (namely the macroporosity factor and the soil wetness threshold parameter for fracture ﬂow activation) and uses the saturated hydraulic conductivity from the chalk matrix. The BC model is implemented in the Joint UK Land Environment Simulator (JULES) and applied to a study area encompassing the Kennet catchment in the southern UK. This parameterization is further calibrated at the point scale using soil moisture proﬁle observations. The performance of the calibrated BC model in JULES is assessed and compared against the performance of both the default JULES parameterization and the uncalibrated version of the BC model implemented in JULES. Finally, the model performance at the catchment scale is evaluated against independent data sets (e.g. runoff and latent heat ﬂux). The results demonstrate that the inclusion of the BC model in JULES improves simulated land surface mass and energy ﬂuxes over the chalk-dominated Kennet catchment. Therefore, the simple approach described in this study may be used to incorporate the ﬂow processes through a chalk unsaturated zone in large-scale land surface modelling applications.


Introduction
Chalk can be described as a fine-grained porous medium traversed by fractures (Price et al., 1993).Previous studies showed that the unsaturated zone of the chalk aquifers plays an important role in groundwater recharge in the UK (e.g. Lee et al., 2006;Ireson et al., 2009).Therefore, both monitoring (e.g.Bloomfield, 1997;Ireson et al., 2006) and modelling (e.g.Bakopoulou, 2015;Brouyère, 2006;Ireson andButler, 2011, 2013;Sorensen et al., 2014) strategies have been adapted previously to understand the governing hydrological processes in the chalk unsaturated zone.
In chalk, the matrix provides porosity and storage capacity, while the fractures greatly enhance permeability (Van den Daele et al., 2007).Water movement through the chalk matrix is slow due to its relatively high porosity (0.3-0.4) and low permeability (10 −9 -10 −8 m s −1 ).A fractured chalk system, in contrast, conducts water at a considerably higher velocity because of the relatively high permeability (10 −5 -10 −3 m s −1 ) and low porosity (of the order 10 −4 ) of fractures (Price et al., 1993).
Simulating water flow through the matrix-fracture system of chalk has been the subject of research for some time.Both conceptual (e.g.Price et al., 2000;Haria et al., 2003) and physics-based (e.g.Mathias et al., 2006;Ireson et al., 2009) models have been proposed previously to describe water Published by Copernicus Publications on behalf of the European Geosciences Union.
flow through the chalk unsaturated zone.The physics-based models mentioned above were developed based on a dualcontinua approach and required relatively large numbers of parameters (i.e. of the order of 20-30 parameters) that were calibrated via inverse modelling using observed soil moisture and matric potential data (e.g.Ireson et al., 2009;Mathias et al., 2006).
In recent years, representation of chalk has gained attention in land surface modelling.For example, Gascoin et al. (2009) applied the Catchment Land Surface Model (CLSM) over the Somme River basin in northern France.A linear reservoir was included in the TOPMODELbased runoff formulation of the CLSM to account for the contribution of chalk aquifers to river discharge.Le Vine et al. (2016) applied the Joint UK Land Environment Simulator (JULES; Best et al., 2011) over the Kennet catchment in southern England to evaluate the hydrological limitations of land surface models.In that study, two intersecting Brooks and Corey curves were proposed, which allowed a dual-curve soil moisture retention representation for the two distinct flow domains of chalk (i.e.matrix and fracture) in the model.Considering this dual Brooks and Corey curve, a three-dimensional groundwater flow model (ZOOMQ3D; Jackson and Spink, 2004) was coupled to JULES to demonstrate the strong influence of representing chalk hydrology and groundwater dynamics on simulated soil moisture and runoff.
The above-mentioned studies illustrate the importance of representing chalk in land surface modelling.However, including chalk hydrology in large-scale land surface modelling using the contemporary dual-porosity concept can be complicated due to the large number of additional parameters.In this context, we propose a new parameterization, namely the Bulk Conductivity (BC) model, as a first step towards a simple chalk representation suitable for land surface modelling.In order to test the proposed parameterization, the BC model is included in JULES (version 4.2), which, by default (i.e.uniform soil column representation using a general soil database as typically applied in land surface models), does not represent any chalk feature.In this study, the BC model (included in JULES) is applied at two distinct spatial scales (i.e. point and catchment).At the point scale, the proposed parameterization is calibrated using observed soil moisture profile data.This is achieved by randomly sampling the parameter space and extensively running the model in order to minimize the differences between observed and simulated soil moisture variability at different depths.Finally, the proposed model is applied to the Kennet catchment in southern England and the fluxes and states of the hydrological cycle are simulated for multiple years.The simulation results are evaluated using observed latent heat flux (LE) and runoff data to assess the performance of the BC model in simulating land surface processes at the catchment scale.

A model of flow through a chalk unsaturated zone
In this study, the BC model based on the work by Zehe et al. (2001) is incorporated into JULES to represent the flow of water through the fractured chalk unsaturated zone.According to this approach, if the relative saturation (S) exceeds a certain threshold (S 0 ) at a soil grid, the saturated hydraulic conductivity of the chalk matrix (K s ) is increased to a bulk saturated hydraulic conductivity (K sb ) as follows: with where f m is a macroporosity factor (-), θ is soil moisture (m 3 m −3 ), θ s is soil moisture at saturation (m 3 m −3 ) and θ r is the residual soil moisture (m 3 m −3 ).Note that S ranges from 0 in the case of completely dry soils to 1 for fully wet soils.
At the first step of evaluation, the K s , S 0 and f m parameters are estimated based on the existing literature to assess the performance of the uncalibrated BC model.For the matrix saturated hydraulic conductivity (K s ), we use K s = 1.0 mm day −1 following Mathias et al. (2006).In addition, Eq. ( 1) indicates that the onset of water flow through the fracture system of chalk is controlled by the threshold S 0 .According to Wellings and Bell (1980), water flow through fractures dominates over matrix flow in chalk when the pressure head in soil becomes higher than −0.50 m H 2 O.We consider a value of S 0 = 0.80 for the uncalibrated BC model, which is based on the observed soil moisture-matric potential relationship in the study area.
Finally, in Zehe et al. (2001), f m was defined as the ratio of the saturated water flow rate in all macropores in a model element to the corresponding value in the soil matrix, which can be determined based on the density and length of fractures on small scales.In addition, f m has also been considered as a calibration parameter previously (e.g.Blume, 2008;Zehe et al., 2013).In this study, we define f m as a characteristic soil property reflecting the influence of fractures on soil water movement (Zehe and Blöschl, 2004) and estimate it from the relative difference of permeability between chalk matrix and fractured chalk system that can be of the order 10 4 -10 6 according to Price et al. (1993).Consequently, we consider a macroporosity factor of f m = 10 5 for the uncalibrated BC model.In the following step, the BC model is calibrated to minimize the differences between the variability of observed and simulated soil moisture at individual depths.The calibration strategy will be discussed elaborately in Sect. 3.5. Hydrol. Earth Syst. Sci., 21, 459-471, 2017 www.hydrol-earth-syst-sci.net/21/459/2017/   3 Methods

Study area
The study area encompasses the Kennet catchment located in southern England, with an area of about 1033 km 2 (Fig. 1a).Generally, Kennet is rural in nature with scattered settlements and has a maximum altitude of approximately 297 m (above ordnance level).The River Kennet discharges into the North Sea through London.The major tributaries of this river are Lambourn, Dun, Enborne, and Foudry Brook.An average annual rainfall of approximately 760 mm was recorded in the catchment over a 40-year period from 1961 to 1990.The solid geology of the Kennet catchment is dominated by chalk, which is overlain by a thin soil layer.While lower chalk outcrops along the northern catchment boundary, progressively younger rocks are found in the southern part.In general, surface runoff production is very limited over the regions of the catchment where chalk outcrops.The flow regime shows a distinct characteristic of slow response to groundwater held within the chalk aquifer (Le Vine et al., 2016).According to Ireson and Butler (2013), the unsaturated zone of chalk shows slow drainage over summer and bypass flow during wet periods in this catchment.

Field measurements and remotely sensed data
Table 1 summarizes the field measurements and remote sensing data used in this study.We use in situ soil moisture and runoff measurements along with remotely sensed LE data to assess model performance in simulating the mass and energy balance components of the hydrological cycle.Point-scale soil moisture measurements at two adjacent sites (∼ 20 m apart) at Warren Farm (Fig. 1) were provided by the Centre for Ecology and Hydrology (CEH).A Didcot neutron probe was used at these locations to measure fortnightly soil moisture at different depths below the land surface (10 cm apart down to 0.8 m, 20 cm apart between 0.8 and 2.2 m, and 30 cm apart between 2.2 and 4.0 m) (Hewitt et al., 2010).
The National River Flow Archive (NRFA) coordinates discharge measurements from the gauging station networks across the UK.These networks are operated by the Environmental Agency (England), Natural Resources Wales, the Scottish Environment Protection Agency, and the Rivers Agency (Northern Ireland).We use discharge measurements provided by the NRFA to assess model performance in simulating runoff over the Kennet catchment in this study.
The MOD16 product of the Moderate Resolution Imaging Spectroradiometer (MODIS) is a part of the NASA/EOS project that provides estimation of global terrestrial LE.The LE estimation from MOD16 is based on remotely sensed land surface data (e.g.Mu et al., 2007).In this study, the 8-day and monthly LE data products from MODIS are used to evaluate the model performance in simulating land surface energy fluxes.

Land surface model
In this study, we use the Joint UK Land Environment Simulator (JULES; e.g.Best et al., 2011;Clark et al., 2011) version 4.2.JULES is a flexible modelling platform with a modular structure aligned to various physical processes developed based on the Met Office Surface Exchange Scheme (MOSES; e.g.Cox et al., 1999;Essery et al., 2001).Meteorological data including precipitation, incoming shortand long-wave radiation, temperature, specific humidity, surface pressure, and wind speed are required to drive JULES.Each grid box in JULES can comprise nine surface types (broadleaf trees, needleleaf trees, C3 grass, C4 grass, shrubs, inland water, bare soil, and ice) represented by respective fractional coverage.Each surface type is represented by a tile and a separate energy balance is calculated for each tile.
Subsurface heat and water transport equations are solved based on finite difference approximation in JULES as described in Cox et al. (1999).Moisture transport in the subsurface is described by the finite difference form of Richards' equation.The vertical soil moisture flux is calculated using Darcy's law.While the top boundary condition to solve the Richards' equation is infiltration at the soil surface, the bottom boundary condition in JULES is free drainage that contributes to subsurface runoff.
Surface runoff is calculated by combining the equations of throughfall and grid box average infiltration in JULES.In order to direct the generated runoff to a channel network, river routing is implemented based on the discrete approximation of a one-dimensional kinematic wave equation (e.g.Bell et al., 2007).In this approach, a river network is derived from the digital elevation model (DEM) of the study area and different wave speeds are applied to surface and subsurface runoff components and channel flows (e.g.Bell and Moore, 1998).A return flow term accounts for the transfer of water between subsurface and land surface (e.g.Dadson et al., 2010Dadson et al., , 2011)).

Model configurations and input data
In this study, simulations are performed at two distinct spatial scales, namely point and catchment.At the point scale, JULES is configured to simulate the mass and energy fluxes at the Warren Farm site (Fig. 1a).A total subsurface depth of 5 m is considered in the model with a vertical discretization ranging from 10 cm at the land surface to 50 cm at the bottom of the model domain.Note that this discretization is consistent with the soil moisture measurement depths mentioned in Sect.3.2.The vegetation type is implemented as C3 grass using the default parameters in JULES.Point-scale simulations were performed over 2 consecutive years from 2003 to 2005 at an hourly time step.Except for precipitation, hourly atmospheric forcing data to drive JULES were obtained from an automatic weather station operated by the CEH at Warren Farm.In order to estimate hourly precipitation data to run JULES, rain gauge measurements from the Met Office (Met Office, 2006) were used.The inverse distance interpolation technique (e.g.Garcia et al., 2008;Ly et al., 2013) was applied to rainfall measurements from the 13 gauges closest to Warren Farm (the distance varies from 25 to 60 km) to obtain hourly precipitation for the point-scale simulations.
At the catchment scale, JULES is configured over a study area encompassing the Kennet catchment (Fig. 1a) considering a uniform lateral grid resolution of 1 km with 70 × 40 cells in the x and y dimensions, respectively.The total subsurface depth and vertical discretization are identical to those of the point-scale simulations.Spatially distributed vegetation-type information for the study area (Fig. 1b) is obtained from the Land Cover Map 2007 (LCM2007) data set (Morton et al., 2011).Simulations were performed over The hydraulic properties for chalk used in this study are summarized in Table 3.These properties are obtained based on the existing literature as a first step when evaluating the uncalibrated BC model.The BC model parameters are subsequently calibrated to minimize the differences between observed and simulated θ (Sect.3.5) at various soil depths.
In this study, we consider two different model configurations, namely default and macro (Fig. 2).The default configuration corresponds to the standard parameterizations of JULES that does not represent chalk hydrology in the model.In this configuration, each soil column in JULES is considered to be vertically homogeneous with the soil properties defined in Table 2, which is motivated by the Met Office JULES Global Land 4.0 configuration described in Walters et al. (2014).The macro configuration, in contrast, explicitly represents chalk by applying the BC model starting at 30 cm below the land surface to the bottom of the model domain (i.e.500 cm).Therefore, the soil column in the macro configuration can be divided into topsoil (0-30 cm) and chalk (30-500 cm).The topsoil depth of 30 cm is defined based on several augured soil samples collected during a field campaign at Warren Farm in 2015 (Fig. 2).This depth is corroborated by additional information from the British Geological Survey (BGS) operated borehole records (http://www.ukso.org/pmm/soil_depth_samples_points.html), which show that topsoil depths vary from 10 to 40 cm over the study area.We apply the macro configuration assuming a spatially homogeneous topsoil depth of 30 cm for both point-and catchmentscale simulations.Note that except for this inclusion of chalk, default and macro configurations are identical in terms of model set up and input data.It should also be emphasized that default represents a "naïve" configuration deprived of model calibration.Moreover, this configuration does not represent chalk, which, according to previous studies (e.g.Le Vine et  1c).Saturated hydraulic conductivity (K s ) and porosity data are obtained from Rawls et al. (1982).The Van Genuchten parameters are acquired from Schaap and Leij (1998) al., 2016), substantially affects the hydrology of the study area considered here.

Calibration of the BC model
We calibrate the BC model at the point scale to minimize the differences between observed and simulated soil moisture variability ( θ ) at different depths.The root mean squared error (RMSE) is used as the objective function to optimize the BC model parameters (e.g.Ireson et al., 2009): where nd is the number of soil layers, nt is the number of soil moisture observations available for a layer d, θ obs is the observed variability of soil moisture and θ sim is the simulated variability of soil moisture.Note that we consider θ for this optimization because of its relevance to the water flux and recharge through a chalk unsaturated zone (e.g.Ireson and Butler, 2011).Equation (1) reveals that the calibration of the BC model involves optimizing three parameters, namely the saturated hydraulic conductivity of the chalk matrix (K s ), the saturation threshold (S 0 ), and the macroporosity factor (f m ).Ireson et al. (2009) suggested a range of 0.2-2.0mm day −1 for K s .On the other hand, Price et al. (1993) argued that, in general, K s is around 3-5 mm day −1 for most chalk soils.Therefore, we consider a range of 0.2-5.0mm day −1 in optimizing K s .We consider the S 0 range 0-1, representing the entire physical domain for soil wetness from fully dry to fully wet, respectively.For f m , a range of 10 4 -10 6 is considered, which, as discussed earlier, is consistent with the relative difference between the permeability of fractured chalk and chalk matrix according to Price et al. (1993).The Latin hypercube sampling technique (e.g.Mckay et al., 2000) is used to generate 2000 random samples for each BC model parameter within the ranges discussed above.Note that for the K s parameter, the random sampling was performed from a logarithmic distribution (Ireson et al., 2009).We perform simulations using these random samples and calculate model performance (Eq. 3) to select the optimum parameter values for the BC   from the default and macro configurations at 2 m below the land surface.Note that the macro configuration uses the chalk hydraulic parameters collected from the existing literature (Table 3).This figure shows that the default configuration considerably underestimates θ throughout the simulation period, which is improved remarkably in the case of macro.Figure 3b plots observed and simulated soil moisture variability ( θ) from the default and macro configurations ( θ default and θ macro , respectively) at the Warren Farm site.
In general, both configurations show discrepancies with observed θ with macro showing relatively better model performance.
The results show that despite the macro configuration improves simulated θ , it shows considerable discrepancies with observed θ , which is consistent throughout the whole chalk profile (results from other model layers are not shown).In order to minimize the differences between observed and modelled θ from the macro configuration, we calibrate the BC model following the methodology described in Sect.3.5.The optimization results are summarized in Fig. 4. Note that for each combination considered in the optimization, 2000 model runs were performed using randomly sampled parameters as discussed in Sect.3.5.In addition to the default and macro cases, the calibrated cases in Fig. 4 correspond to the results from the model runs yielding the lowest RMSE for each parameter combination evaluated.
The RMSE between observed and simulated θ for the model configurations considered in the optimization is shown in Fig. 4a.This figure illustrates that the RMSE of the default configuration is larger than that of macro, indicating better model performance in reproducing θ for the latter (corresponding to a reduction of 15 % in RMSE compared to the default case).Therefore, the uncalibrated BC model (i.e.macro configuration) better reproduces the soil moisture variability compared to the default.Concerning the calibration of single BC model parameters, Fig. 4a shows that S 0 results in a 46 % reduction of RMSE compared to the macro configuration.Calibrating K s or f m individually yields only about a 25 % reduction of RMSE compared to macro.
Optimizing both K s and S 0 simultaneously shows the largest reduction (50 %) of RMSE compared to macro, which coincides with the total RMSE reduction found when all parameters are calibrated.Arguably, the BC model can be implemented in other chalk regions by constraining only the S 0 parameter.Such a result could potentially be advantageous for transferability to other regions in the UK in order to assess chalk hydrology on a large scale.However, since this is the first time the BC model is introduced, we decide to take a conservative approach and select the macro configuration with optimized K s and S 0 (macro opt hereafter) to simulate chalk hydrology over the study area that ensures the best overall model performance.
The lower three panels in Fig. 4 presents the BC model parameter values for the default and uncalibrated macro cases as well as for different combinations of parameters cali-Hydrol.Earth Syst.Sci., 21, 459-471, 2017 www.hydrol-earth-syst-sci.net/21/459/2017/ brated.The red bars in Fig. 4b-d highlight the cases in which a given parameter is constrained by optimization.In those cases, the calibrated parameter values are obtained from model runs producing the lowest RMSE.An interesting feature in Fig. 4b (calibrating K s individually) is that the optimization suggests a compensation mechanism in which K s is increased remarkably in order to physically represent the "effective" flow through the chalk fractures in the BC model.This is not surprising and arguably the simplest way to attempt to improve model performance.For macro opt , the value used for K s is relatively lower than that of the uncalibrated macro case nevertheless consistent with previous estimates (e.g.Ireson et al., 2009).Figure 4c clearly shows the dominance of S 0 in the BC model as all the relatively low RMSE bars in Fig. 4a are associated with S 0 calibration (see red bars in Fig. 4c).In addition, calibrated S 0 values for all cases show a consistent constraint around 0.50.Finally, Fig. 4d indicates the lack of influence for f m parameter on model performance.
Figure 5 compares θ default , θ macro , and θ from the macro opt configuration ( θ opt ) with observed soil moisture variability ( θ obs ).As mentioned earlier, θ default and θ macro show considerable discrepancies with θ obs , while the macro configuration exhibits relatively better performance (Fig. 3). Figure 5 illustrates that the overall agreement between observed and simulated θ improves substantially in the case of macro opt compared to the default and macro, which is pronounced, especially in the deeper chalk layers.Therefore, this figure indicates that the performance of the BC model in simulating θ is further improved by optimizing the K s and S 0 parameters simultaneously at the Warren Farm site.
As mentioned earlier, efficiently reproducing soil moisture variability over the profile is important due to the fact that θ significantly affects water flux and recharge through the chalk unsaturated zone.The drainage flux through the bottom of the soil column (d b ) of a land surface model can be considered the potential recharge flux to groundwater (e.g.Sorensen et al., 2014).Figure 6 compares the daily sum of d b from the default and macro opt configurations at the Warren Farm site.Daily rainfall at this site over the simulation period is shown in Fig. 6a.In Fig. 6b, the macro opt configuration shows considerable d b during the colder months, while relatively slow drainage is observed in summer.In contrast, the default configuration shows relatively high d b in summer compared to the colder months.In general, the recharge rate through the chalk unsaturated zone during the warmer periods of the year is lower than that in the winter months (Wellings and Bell, 1980;Ireson et al., 2009).Therefore, the macro opt configuration appears to be more consistent with the recharge mechanism in chalk compared to the default.
In this section, the BC model was evaluated at the point scale.The results showed that, in general, the macro configuration outperforms the default case in simulating θ .In order to improve the model performance even further, model parameter calibration was performed to minimize the differences between observed and simulated θ at the point scale.In the next sections, the optimized model (macro opt ) is evaluated at the catchment scale.

Catchment-scale simulations
At the catchment scale, simulation results from the default and macro opt configurations are compared with the observations over the Kennet catchment.In order to assess the differences between LE from the default and macro opt configurations at the catchment scale, Fig. 7 plots spatially averaged 8-day composites of LE from MODIS (LE MOD ) against the LE from these configurations (LE default and LE opt , respectively) over Kennet.The agreement between simulated LE and LE MOD is evaluated using the coefficient of determination (R 2 ; see the Appendix) and the mean bias.Comparison between LE default and LE MOD shows a coefficient of determination of R 2 default = 0.78 and a mean bias of  bias default = 10.5 W m −2 .The agreement between simulated LE and LE MOD improves in the case of the macro opt configuration, which is reflected by an increased coefficient of determination of R 2 opt = 0.80 and a reduced mean bias of bias opt = 7.1 W m −2 .
Figure 7 shows considerable differences between LE default and LE opt for relatively high LE, indicating discrepancies, especially during the warmer months of the year.Spatially averaged time series of monthly LE MOD , LE default , and LE opt are presented in Fig. 8a.This figure shows that the differences between LE default and LE opt increase in summer compared to the colder months of the year, which is consistent with Fig. 7. Consequently, the default configuration under- estimates LE in summer compared to LE MOD , which is improved in the case of the macro opt configuration.In contrast, the differences between LE default and LE opt are negligible during the colder months of the year.In addition, Fig. 8b compares the observed and simulated monthly average discharge from the two model configurations at the Kennet at Theale gauging station (Fig. 1a).This figure shows that the default configuration generally overestimates discharge at this gauging station, which is improved considerably in the case of macro opt .We use the Kling-Gupta efficiency criterion (KGE; Gupta et al., 2009) to compare the performance of the two model configurations in reproducing observed discharge variability.As mentioned above, the default configuration overestimates discharge with KGE default = −0.17.On the other hand, the macro opt configuration improves the agreement between observed and simulated discharge, which is reflected by KGE opt = 0.51.
In order to summarize the results at the catchment scale, Table 4 compares observed and simulated runoff from the two model configurations over the Kennet catchment from 2006 to 2011.The runoff ratio (RR; see the Appendix), which is equal to the mean volume of flow divided by the volume of precipitation (e.g.Kelleher et al., 2015), assesses the partitioning of precipitation into runoff over the catchment.The default configuration (RR = 0.82) shows considerably higher RR compared to observation (RR = 0.40), indicating overestimation of runoff by the model that is consistent with Fig. 8b.Including chalk hydrology in the model remarkably improves the agreement between observed and simulated mean runoff over the Kennet catchment, which is assessed from a runoff ratio of RR = 0.46 for the macro opt configuration, which is much closer to the observed RR value than the default.
In Table 4, the relative bias ( µ) of 1.04 between observed and simulated runoff from the default configuration again indicates the overestimation by the model.In comparison, macro opt shows a smaller relative bias of µ = 0.12, indicating improved agreement between observed and simulated mean runoff volume compared to the default.The relative difference in standard deviation ( σ ; see the Ap-  4 relating directly to the seasonal change in runoff.This comparison shows that the default configuration overestimates the variability of runoff over the Kennet catchment ( σ = 2.04), which is improved in the case of macro ( σ = 0.65).This improvement in reproducing flow variability is also clearly observed in Fig. 8b.In this section, the BC model is evaluated using observed mass and energy fluxes over the Kennet catchment.The default configuration suggested relatively low summertime LE over the catchment.The agreement between observed and simulated LE was improved in the case of the macro opt configuration compared to the default.It was also observed that the overall runoff prediction was considerably improved by macro opt compared to the default.Given its simplicity, our results indicate that the proposed parameterization is suitable for use in land surface modelling applications.

Summary and conclusions
In this study, we proposed a simple parameterization, namely the Bulk Conductivity (BC) model, to simulate water flow through the matrix-fracture system of chalk in large-scale land surface modelling applications.This parameterization was implemented in the Joint UK Land Environment Simulator (JULES) and applied to the Kennet catchment located in the southern UK to simulate the mass and energy fluxes of the hydrological cycle for multiple years.Two model configurations, namely default and macro, were considered, with the latter using the BC model to simulate chalk hydrology.
The proposed BC model is a single-continuum approach to modelling preferential flow (e.g.Beven and Germann, 2013) that involves only three parameters, namely the saturated hydraulic conductivity of the chalk matrix (K s ), the macroporosity factor (f m ), and the relative saturation threshold (S 0 ).Initially, these parameters were estimated from the existing literature to assess the performance of the uncalibrated BC model.Finally, the BC model parameters were optimized to minimize the differences between observed and simulated soil moisture variability.Our results indicated that S 0 is by far the most influential parameter in the model when repre-senting water movement through a soil-chalk column.This highlights the simplicity of the proposed BC model for largescale studies and potential ease in transferability.In comparison, K s and f m showed secondary (low) sensitivity in the model performance.Since this study introduces the BC model, we decided however to take a conservative approach.We optimized K s and S 0 simultaneously for our catchmentscale simulations since this combination resulted in the best overall model performance.
At the catchment scale, the proposed BC parameterization improved simulated latent heat flux (especially in summer) and the overall runoff compared to the default.Note that the complexity (i.e. the number of parameters) of the BC model for simulating water flow through a chalk unsaturated zone is substantially lower compared to more commonly used models for this purpose (e.g.dual-porosity models).Despite its simplicity, the proposed parameterization considerably improves the key hydrological fluxes simulated by JULES at the catchment scale.Therefore, the BC model can potentially be useful for land surface modelling applications over largescale chalk-dominated areas.

Data availability
The soil moisture profiles and atmospheric information used for the point-scale simulations (Warren Farm site) are available from the Centre for Ecology and Hydrology (CEH) upon request.Rain gauge measurements to obtain the precipitation estimates at the Warren Farm site are available upon request from the NCAS British Atmospheric Data Centre (http://catalogue.ceda.ac.uk/uuid/bbd6916225e7475514e17fdbf11141c1).For the catchment-scale simulations, the soil texture information is freely available at http://www.fao.org/soils-portal/soil-survey/soil-maps-and-databases/ harmonized-world-soil-database-v12/en/.The Moderate Resolution Imaging Spectroradiometer (MODIS) data are freely available at https://modis.gsfc.nasa.gov/data/dataprod/mod16.php.Land cover information for the catchment-scale simulations can be acquired upon request from the CEH (https://www.ceh.ac.uk/services/land-cover-map-2007).The daily atmospheric variables over the Kennet catchment can also be obtained upon request from the CEH (https://catalogue.ceh.ac.uk/documents/ 8baf805d-39ce-4dac-b224-c926ada353b7).Daily river discharge data (from the National River Flow Archive, NRFA) at the Kennet at Theale station are freely available at http://nrfa.ceh.ac.uk/data/station/info/39016.
Warren Farm; b Locations are shown in Fig. 1a.

Figure 1 .
Figure 1.(a) Location, (b) vegetation cover and (c) soil texture over the study area.The red line in (a) outlines the Kennet catchment boundary, while the river network is shown in blue.The black triangle in (a) shows the location of the discharge gauging station at the catchment outlet and the black square corresponds to the Warren Farm location where point-scale simulations are carried out.The black line in (c) encloses the area of the catchment where chalk is present.

Figure 2 .
Figure 2. (a) Example of soil profiles collected at Warren Farm during a field campaign in 2015 and (b) the two model configurations.

Figure 3 .
Figure 3.Comparison between observed and simulated (a) soil moisture (θ ) and (b) change in soil moisture ( θ ) from the default and macro configurations at a depth of 2 m below the land surface at the Warren Farm site.The shaded areas constructed from two soil moisture probes at the Warren Farm site denote the range of observed data in these plots.

Figure 4 .
Figure 4. (a) Model performance in reproducing observed and simulated θ, (b) K s , (c) S 0 and (d) f m for the different parameter combinations considered in the optimization.For each parameter (i.e.b, c, and d), the red bars show cases in which the relevant parameter is calibrated (either individually or in combination with others), while the blue bars correspond to cases in which the selected parameter is not calibrated (i.e. a fixed value according to the literature as in the macro case).Note that except for the default and macro, the simulation yielding the lowest RMSE (out of 2000 model runs) is presented in this plot.

Figure 5 .
Figure 5.Comparison between observed and simulated θ from default, macro and macro opt configurations at various depths below the land surface.The shaded areas, which are constructed from two soil moisture probes at the Warren Farm site, denote the range of θ .

Figure 6 .
Figure 6.(a) Daily precipitation and (b) daily drainage through the bottom of the soil column at Warren Farm over the 2 simulated years (2003-2005).

Figure 7 .
Figure 7. Catchment average 8-day composites of MODIS estimated LE (LE MOD ) against simulated LE from the default and macro configurations (LE default and LE macro , respectively) along with the linear models fitted for LE default (black line) and LE macro (red line).The 1 : 1 line is shown in grey, which represents the perfect fit between LE MOD and simulated LE.

Figure 8 .
Figure 8.(a) Spatially averaged monthly latent heat flux (LE) from MODIS, default and macro opt configurations over the Kennet catchment and (b) monthly average observed and simulated discharge from the default and macro opt configurations at the Kennet at Theale gauging station.

Table 1 .
Field measurements and remote sensing data.

Table 2 .
Hydraulic properties for different soil types (refer to Fig. .

Table 3 .
Hydraulic properties of chalk.

Table 4 .
Comparison between observed and simulated daily average runoff from the two configurations over the Kennet catchment.Metrics include the runoff ratio (RR), relative bias ( µ), and relative difference in standard deviation ( σ ) (refer to the Appendix for further information).