Assessments of Solar, Thermal and Net Irradiance from Simple Solar Geometry and Routine Meteorological Measurements in the Pannonian Basin

: In this paper, we discussed several different procedures for calculating irradiation from routine weather measurements and observations. There are between four and eight frequently used parameterizations of radiation balance components in meteorological preprocessors, and we investigated them. First of all, the estimated and measured solar and net irradiance were compared. Af-terwards, the estimated and measured longwave irradiance were investigated. Then, we recalculated the net irradiance from the sum of global solar irradiance, longwave downwelling irradiance, reflect solar irradiance and upwelling longwave irradiance. Statistical estimates of the described methods were also recalculated compared with each shortwave and longwave radiation budget component measured separately with WMO first-class radiation instruments (Kipp&Zonen CMP6 and CMP11 and CGR3 and CGR4) in the Agrometeorological Observatory Debrecen, Hungary during one-year time period. Finally, we compared the calculated and measured values for longer periods (2008–2010 and 2008–2017) through statistical errors. The suggested parameterizations of the net radiation based on the separately parameterized all radiation balance components were: Foken’s calculation for clear sky solar global irradiance, Beljaars and Bosveld parameterization for albedo, Dilley and O’Brien methodology for the clear sky incoming longwave (LW) irradiance and Holstlag and Van Ulden cloudiness correction for all sky incoming LW and for the LW outgoing irradiance.


Introduction
In meteorology and climatology, a routine meteorological observation program generally selects only state variables of the atmosphere for measurements. At the same time, the growing research on the problem of climate change, air pollution dispersion or wind energy resource has increased demand for the best possible information about surface layer parameters (SLP): turbulence fluxes of momentum ( ), sensible ( ) and latent heat ( ), trace gases ( ) and aerosol particles ( ) near the ground [1][2][3][4]. It has been known for some time already that surface fluxes of momentum, heat and moisture ( ) are essential for determining atmospheric steady states [5,6]. These fluxes can be both directly measured and calculated by, e.g., eddy covariance, profile, gradient and Bowen ratio methods or calculated from routine measurements [7][8][9]. Therefore, the development and checking models for these calculations are an important area of research where the Monin-Obukhov similarity theory (MOST) is a generally accepted framework for describing interactions between the Earth surface and atmosphere [10].
When dealing with problems inside the surface layer using MOST, Monin-Obukhov (MO) scaling variables are the key variables in a dimensional analysis. There are: length scale ( ), scales of velocity or friction velocity ( * ), velocity scale for the convective planetary boundary layer (PBL) dependent on the sensible heat flux ( * ), roughness length ( ) and temperature scale ( * ). The most important scaling variable in MOST is the turbulent length scale-, which dictates stratification of the atmosphere. In Western Europe, many authors have investigated processes inside the surface and boundary layer, using the Cabauw experiment and data, suggesting ten classes of atmospheric stability, depending on the turbulent length scale, . These classes are: extremely stable, very stable, moderate stable, light stable, neutral, light unstable, moderate unstable, very unstable and free local convection. It has been discovered that the universal function of MOST does not work properly in the case of extremely stable stratification of the atmosphere. Comparing the number of hours with extremely stable stratification of atmosphere on the Cabauw location and on the locations of Novi Sad and Debrecen, both in the Pannonian region (Figure 1d), enormous differences were noticed. While the Cabauw location has about 3% percent of these situations, both stations in the Pannonian Basin have more than 30%. Due to the significances of the , the authors who offered calculations of this parameter from routine meteorological measurements gave huge support to the investigation of the PBL. The modeling of these scales from standard weather observation FM12-SYNOP was developed and described by Holtslag, Van Ulden and De Bruin from 1982 to 1988 and partly modified by Foken and Göckede recently [1, 5,[11][12][13][14][15]. Their model needs only astronomical conditions and routine SYNOP measurements for determining SLP, so it can be widely accepted where there are no special flux or profile measurements. The first step in these models, which still determine the development of the parameterization methods [16][17][18], is assessing a global solar radiation ( ↓) in the cases when this parameter is not measured on a routine basis. The next steps are the assessments of albedo ( ) and, finally, crucial parameter net radiation ( * ), which is usually not included in standard meteorological measurement programs. However, these models/parameterizations can be used only after careful consideration of all the relevant facts. Namely, the parameters that are used were empirically determined by experiments, and they are a function of the place, time and weather situation in the period where the experiment was done. Therefore, these models may have different representativeness in different regions. Since the calculation of from routine SYNOP measurements essentially depends on a sensible heat flux, and a sensible heat flux essentially depends on upwelling and downwelling shortwave and longwave radiation, the basic question is to find out which are the most suitable radiation models in the Pannonian Basin. It is the first step in the classification of stability. The stability class is crucially important in the region that has a specific distribution of atmospheric stability in comparison with other regions. There are already several papers for checking the parameterization procedures of the individual radiation balance components, especially for extreme geographical locations in the recent years. For example, longwave parameterization procedures were investigated over the Tibetan Plateau by Zhu et al. [19] and Liu et al. [20]. Short-and longwave radiation parameterizations were tested with satellite measurements over an alpine glacier in Italy (Senese et al. [21]). Stettz et al. [22] analyzed the shortwave radiation parameterizations in a tropical area. Lindauer et al. [23] developed a new general parameterization for incoming solar radiation dependent only on screen-level relative humidity with site-specific astronomical information. In addition to the surface measurements, satellite observations are increasingly used in radiation balance modeling. Our work is different, as we take into account all of the components of the radiation balance, focusing on the modeling methods used in meteorology for a lowland area of Europe where a high-density surface measurement network is available over eight countries of the Carpathian Basin.
The results of previous research of radiation and surface energy budget investigations in Southeast Europe and the Carpathian region were also illustrated by a few examples [24][25][26][27][28]. Two of the SVAT models are highlighted. The land-surface flux model (PROGSURF) was designed jointly at the Universities of Vienna and Budapest [24]. This model solves the surface energy budget equation using the Penman-Monteith approach. Net radiation is modeled separately for bare and vegetation-covered soil through the radiation balance components. The model comprises one vegetation layer and three soil layers. Surface temperature prediction is made by the heat conduction equation in conjunction with the force-restore method. The ground-air parameterization system (LAPS) also relies on data from synoptic stations. It was developed at the University of Novi Sad [25]. It is structured according to similar principles as the previous one [24]. The model is suitable for calculations over both heterogeneous and nonheterogeneous surfaces. It follows the methodology of De Bruin and Holtslag [10] in the modeling of radiation balance components.
The parameterization of the radiation balance components plays an important role in the meteorological preprocessor prepared for Hungary, which determines the daily course of the PBL thickness and the surface energy budget's components for air quality purposes [26]. The development is based on the methodology of Holtslag and Van Ulden [5] and COST Action 710 [1].
For many tasks, e.g., energy production forecasting of solar collectors, a detailed shortwave radiation estimate is required. The next two articles on this topic are significant in the Carpathian Basin. Based on radiation measurements in Romania, computing global and diffuse solar hourly irradiation in a clear sky was reviewed and tested using 54 parameterizations on daily and hourly scales [27]. The empirical models for estimating solar insolation were developed in Serbia by using meteorological data on cloudiness [28]. Solar energy utilization has been a priority in the countries of the Carpathian Basin in recent years. This requires reliable global radiation estimates at least in hourly time resolutions. Here, the primary data sources are standard meteorological measurements. It is important to select the optimal parameterization method to be used in the Carpathian Basin to quantify the application uncertainties. Modeling of the radiation balance components in daily and longer time scales, as well as applications of GIS-based modelling for solar energy estimation in different time scales and different climate regions are beyond the scope of this paper (see references [29][30][31][32][33] for an overview).
This paper aims to contribute to the present understanding of how successful the methods for the assessment of shortwave and longwave radiations are in the Pannonian Region.
Two solar radiation models, often used in meteorology, are described. Further, we investigate the relationship between downwelling shortwave radiation and the meteorological parameters that modify it. Besides that, many methods that calculate the upwelling shortwave radiation hourly, as well as a component of longwave radiation as a function of a state of the ground and a state of the atmosphere, are studies. This model of the approximation has become common in SVAT (Soil-Vegetation-Atmosphere-Transport) modeling, and these types of radiation parameterization procedures are also often used in one-dimensional PBL models.
We believe that the results of our study are significant for data quality control. Further, we believe they may help in understanding and bridging differences between the point measurements and atmospheric models for all radiation components, SLP, as well as the state variables, especially temperature and wind, in our region.

Terminology
Solar radiation is electromagnetic energy originating from the Sun [34,35]. Of the light that reaches Earth's surface, infrared radiation makes up around 50% of it, while visible light provides about 42% and ultraviolet radiation just over 8%, respectively.
Thermal or terrestrial radiation emitted by the Earth's surface and atmosphere is in the range of 4-100 μm.
Irradiation is the energy received per area. It is the process by which an object is exposed to radiation; in this case, coming from the Sun, from the atmosphere or from the Earth surface. Its SI unit is (J m -2 ), including time-averaged values (for example, hourly or daily irradiation).
Irradiance is the radiant flux received by a surface per unit area (W m -2 ).

Description of Key Variables
Although solar radiation has a leading role in almost all processes in the atmosphere it was not part of standard reports. During the last century, the instruments that measure irradiance were rare and usually placed on high mountains or in deserts, so as to minimize the atmospheric influence. The growing interest for renewable energy sources is making this measurement more and more popular. In this paper, the following irradiance measurements were used on a horizontal plane (without surrounding obstacles and buildings) at the Earth's surface: Global solar irradiance is radiant flux emitted from the Sun and received at the Earth's surface separated in two basic components: direct and diffuse. Global solar irradiance is a measure of the rate of total incoming solar energy, both direct and diffuse, on a horizontal plane at the Earth's surface. It depends on position of the Sun in the sky, season, time of the day and turbidity of the atmosphere. Turbidity mostly depends on the cloudiness, humidity, content of aerosol particles and, of course, from the pressure (amount of the air column). Reflected solar irradiance is part of global solar irradiance that is reflected from the Earth's surface. It depends on the global solar irradiance and the surface albedo (function of the angle of solar elevation and characteristics of the ground surface). Incoming longwave (LW) irradiance is a downward flux of thermal radiation emitted from the atmospheric molecules (such as H2O, CO2 and O3); aerosol particles and clouds per unit horizontal area in a given time period. It depends, first of all, on cloudiness, temperature precipitable water and turbidity of the atmosphere. Outgoing (upwelling) longwave irradiance represents a redistribution of the absorbed global solar irradiance. The power of this energy emitted from the Earth's surfaces per unit area and in the given time is called thermal (terrestrial) irradiance. Besides the global solar irradiance, it depends on the temperature of the Earth's surface (or atmospheric temperature). Net irradiance is often convenient to split into four components: , , and . Therefore, net irradiance is a sum of these compo- (1) In our case, the signs of all the surface radiation balance components were positive or zero.
Net irradiance is necessary fuel for all motion processes in the atmosphere. Energy budget on the Earth surface in the case of the quasi-steady state of PBL, when there is no heating or cooling, sees the net irradiance as a sum of the sensible heat flux, , latent heat flux, , and flux into or from the ground, .
Accordingly, is key connection between radiation and stability of the atmosphere.
In order to estimate the upper irradiances, the next variables are measured or derived from routine SYNOP reports and the metadata about the location where measurements are conducted: ϕ The angle of solar elevation (rad) Hour angle (rad) ℎ Hour angle used in [5] (rad) ℎ Hour angle used in [11] (rad)

Description of the Key Method
The measurement mentioned in Section 1 helps to better understand how solar energy changes the state of the atmosphere and how the state of the atmosphere influences the global solar irradiance and thermal irradiance on the Earth's surface. There are some investigations that try to track the interaction between the Sun and Earth at some locations hourly through the relationship between: global solar irradiance, ↓, reflected solar irradiance, ↑ , incoming longwave, LW, irradiance, ↓, and outgoing LW irradiance, ↑ , in the surface and meteorological parameters. The methodologies that estimate the net radiation and all its components used in this paper are listed in Table 1 below. Each methodology is represented through its source, i.e., reference (ref.), abbreviation (abbr.) and equations number (eq.), which describe the components of irradiance as a function of variables written in parentheses in the last column.  Table 11 row (1) = f (T, e) [44] ↓ Table 11 (2) = f (T, e) [45] ↓ Table 11 (3) = f (T) [46] ↓ Table 11 (4) = f (T, e) [47] ↓ Table 11 (5) = f (T) [48] ↓ Table 11 (6) = f (T, e) [43] ↓ Table 11 (7) = f (T, e) [46] ↓  Cloudiness parameterization for downwelling LW radiation based on Iziomon et al. [51] Cloudiness parameterization for downwelling LW radiation based on Swinbank [47] and Dilley and O'Brien [48] Cloudiness parameterization for downwelling LW radiation based on Niemelä et al. [43] Each an estimated radiation component is compared with the measurements and discussed through statistical errors: BIAS, MAE, RMSE and correlation coefficient ( ). Statistical errors for the estimated and measured values are calculated by: where and are the mean values of the modeled and measured hourly values, respectively, and n is the number of cases (measured hours).

The Description of the Datasets
The Agrometeorological Observatory in Debrecen is on a flat, mostly homogenous terrain in agricultural surroundings with geographic coordinates: = 47.5291° (latitude) and = 21.6397° (longitude).
There is a meteorological 10-m mast with profile measurements of ten-minute average values of: • temperature, relative humidity and wind speed, at 1, 2, 4 and 10 m; • infrared ground surface temperature; • two levels of soil temperature and humidity (5 and 10 cm) and • all radiation balance components: global solar radiation or , reflected solar radiation , incoming longwave (LW) radiation and outgoing LW radiation .
Together with the station's pressure, precipitation and snow cover obtained from routine meteorological measurements and observations, we used a 2-m temperature and relative humidity and 10-m wind speed from the mast. The momentum, CO2 and sensible and latent heat fluxes were also measured directly using the CSAT-3 and LI-7500 instruments with the Campbell data acquisition program.
Standard meteorological measurements, including cloudiness and visibility, were provided from the SYNOP station in Debrecen Airport (12882) 8.5 km from the Agrometeorological Observatory [54,55].
The new generation of micrometeorological measurements started in 2008. The instruments were initially factory-calibrated, and a year-long dataset was investigated in the year 2009. That is why the first full year has been chosen for model development. For the long-term comparative study of the developed calculation methods, two test periods were selected. The first was 3 years between 2008 and 2010 when cloud and present weather observations were made hourly at Debrecen Airport (SYNOP station No. 12882). The second, longer test period was 7 years between 2011 and 2017. Namely, after that, from 2011, only partial cloud observations took place at night, so we checked the parameterization of the shortwave radiation balance components and, on occasion, the longwave radiation balance components. The entity responsible for the operation of Agrometeorological Observatory Debrecen was changed in 2017. Data with quality control has been available to us until this year. Since 2018, the measuring site has been operating as a "backup" site with limited data access because there has been a new measuring site installed.
Long-term agroclimatological investigations and reorganization of the high-precision radiation and micrometeorological measurements are in progress at Agrometeorological Observatory Debrecen. The location of the observation in the Pannonian Basin, main measurement items and the instrumentation are presented in Figure 1 and Table 2.
The measurement of solar radiation balance components is part of the measuring program of the surface observing system of the Hungarian Meteorological Service. The operation of observing the system and calibration of sensors is regulated by the ISO9001:2015 quality control system. The calibration of solar radiation sensors is done according to the relevant working instructions; therefore, the calibration of the solar radiation sensors of Agrometeorological Observatory Debrecen is done the same way. Regarding the shortwave sensors, the reference instrument is an HF-type absolute cavity pyrheliometer (s.n. 19746) that has been calibrated at International Pyrheliometer Comparisons (IPC-s) at the World Radiation Centre (WRC) in Davos every five years since 1980. In the case of longwave sensors, the reference sensors are an Eppley PIR (s.n. 29582F3) that was calibrated two times at WRC Davos between 2006 and 2014 and a Kipp&Zonen CGR4 (s.n. 080066) that was calibrated during IPC-XII in 2015. The calibration of the radiation sensors was done with parallel measurements between reference and site instruments, which happened in every two years between 2008 and 2013 and every three years after 2013. The history of the calibrations shows that, in the case of shortwave down-and upward sensors, the averaged differences among the values of reference and site instruments did not exceed ±1.5% and ±1%, respectively, if the level of irradiance was over 200 W m -2 . Regarding the downward and upward longwave sensors, the differences, because the downward sensor was not shaded and ventilated, were between ±3.5% and 2%, respectively. The addition of the personal checking of measuring sites was done every 3 months when the leveling of the sensors was checked and cleaning was done. The detailed calculations of uncertainty of the calibration factors of the solar sensors is not part of the relevant working instructions.

Assessments of Solar Irradiance
The Earth's surface, the basic energy transfer area for atmospheric processes, is heated by the absorption of Sun shortwave radiation, ↓. A small part of this radiation, depending on the complexity of the surface, the state of the ground surface, the Sun's position on the sky, visibility and cloudiness, is reflected back to the sky, ↑.
The mean annual daily variation of calculated ↓ and ↑, as well as measured and irradiance, for the clear and cloudy sky are represented by taking the different procedures for their calculations. The statistical errors (See Section 2.3) of their calculations represent the reliability of these calculations.

Assessments of Downwelling Shortwave Solar Radiation for Clear Sky
The astronomical quantities and downwelling solar radiation are calculated using Foken [11] and Holtslag and Van Ulden [5] methodology. Both well-known methods estimate downwelling solar radiation as a function of , the angle of solar elevation. This angle is a function of: solar declination, , latitude, , and hour angle, ℎ (all angles are in radians): where , the solar latitude, is and is the Day of the year. Holtslag and Van Ulden [5] used the hour angle as a function of , east longitude in radians and time as UTC time. The equation of hour angle, ℎ , is calculated as follows: Foken [11] took the hour angle, ℎ , as a function of Δ , the duration of full rotation of the Earth (between two midnights = 86,400 s), and , the time distance from to culmination of the Sun in seconds.
is the difference between the true and averaged local time and is the Central European Time. From tables of time equations [56] for 15° E and 50° N, Göckede [57] calculated an approximation equation for , and the coefficients , , … , were listed in Foken [11]. After checking the agreement between Equations (10) and (11), one notices that the resultant angles of solar elevation (see, also, Equation (7)) are very similar but not same. They have a maximum at the same hour but, obviosly, not at the same moment. At the hour with maximum solar elevation, the Foken's value of the solar elevation angle is about 0.6° greater than Holstlag's during the day of the sumer solsticium and about 0.2° greater than Holstlag's during the winter solsticium. Concerning the daily circling, it looked that Foken's solar elevation lasted about 20-25 min after Holstlag's. This deviation was established arbitrary based on Figure 2a, on which the x-axis represents "hours in the year" and y-axis represents solar elevation. This deviation was not the same every day, and 20-25 min was the annual maximum. It happened on August 2, when the differences between the angles of solar elevation were a little smaller than 6.5°. In 285 out of 8760 h, these differences were greater than 6°; in 1104 h, they were greater than 5° and, in 4145, they were greater than 3°. Using the angle of solar elevation , which is calculated with the hour angle described by Equations (10) and (11), downwelling the clear sky solar radiation ↓ following Holstlag's calculations based on reference [5] is: where and are the mean and actual distances from the Earth to the Sun. The ratio of both can be determined according to Hartmann (1994) (see reference [11]) as a sum of the first three components of the Fourier series. The argument of the Fourier series is 2 265 ⁄ (denominator is 266 in the case of a leap year), and the coefficients of the Fourier series are listened in reference [11]. The angle of solar elevation is calculated with the hour angle described by Equation (11).
is the solar constant ( = 1368 W m ). Although during the year 2009, there were only 201 h both with > 0 and with clear skies, based on the SYNOP station at the airport in Debrecen (12882) (Figure 2), statistical errors for this set of data have also been investigated (Table 3.) Based on the 201 clear sky hours, the Foken's calculation made an error over 100 W m -2 in 11-h data, while Holtslag's calculation had nine such situations. The errors over 50 W m -2 took up about one-quarter of the hourly data, mostly during the winter period, and were fewer for Holtslag's than Foken's calculations.

Assessments of Downwelling Shortwave Solar Irradiance for All Sky Condition
To include cloudiness in the calculation of shortwave incoming solar irradiance, we applied two widely used methodologies based on Kasten and Czeplak [39]: and Burridge and Gadd [40]: where is total cloud cover that goes from 0 for clear skies to 8 octa for total cloudy skies, and , and are cloud covers of high, middle or low clouds, respectively.
The first method was described by Holtslag and Van Ulden [5] and Nyren and Gryning [37]. This method was also popular in climatological investigations and applications requiring rough estimations of solar radiation on horizontal surfaces [58][59][60]. The second method was described by Stull [61] and used in different micrometeorological modeling works for radiation parameterization, such as in surface hydrological modelling or snow evolution system modeling [62,63]. The number of SYNOP stations with complex cloud observations, especially the cloudiness and cloud types, decreased with the installation of automatic surface weather stations from the mid-1990s [64]. Observations of cloud types are still generally performed by human observers and naturally have a few limitations [65]. Although it seems clear that the second method includes cloudiness in a more sophisticated way than the first, the lack of information about cloudiness for different levels: , and , makes the first method more practical. Namely, in the case when the total cloud cover was 8 octa and the upper cloud levels were invisible, we had to approximate and and took = = 5. Besides that, when = 8 and < 8, we supposed that = 8 or = 8, which depended on the type of clouds. On the other hand, instead of a total cloud cover in the first method, effective cloudiness = ( + )/2 is used as an alternative solution in the calculation of ↓ whenever is available from SYNOP.
Statistical errors for global solar irradiance for the four investigated models, comparing 4022 h during the year 2009, when all data necessary for the calculations was available ( > 0 W m ), are shown in Table 4. The errors were absolutely the same using either total or effective cloud cover. About 3000 of 4022 h of data had absolute errors in estimation fewer than 100 W m -2 , but about 500 of them had errors greater than 200 W m -2 . There is no comment here whether these errors were the result of the crude method for estimation of the input solar irradiance or even sporadic errors in the measurements. Namely, the radiation balance components are measured with first-class radiation sensors separately using the Campbell CR1000 datalogger, so the measurements are quality-controlled and mostly reliable. Therefore, without analyzing the reason, all situations when the differences between the measured and estimated solar radiations were fewer than: 50, 100, 200,..., 700 W m -2 were countered, and the results are shown in Table 5. After this counting, we recalculated statistical errors for 3500 h of data when all the methods had absolute errors fewer than 200 W m -2 . In that case, the MAE became smaller around 25% and RMSE around 50%, and the coefficient of correlation grew to 0.96, while the BIAS stayed similar. The few greatest absolute errors (>200 W m -2 ) seemed to be connect to the Kasten and Czeplak [39] method for describing cloudiness, because after excluding these situations, the RMSE and MAE became 10% smaller for this method than for the Burridge and Gadd [40]-type cloudiness parameterizations for both clear sky calculations.
Relative errors were fewer than 10% for only 18-25% of the hourly data. Both methods look very similar through statistical errors (Table 4). Looking at the absolute errors, one notices that the Kasten and Czeplak [39] cloud parameterization has a greater number of hours with errors smaller than 50 W m -2 , although the few greatest AE relate to this parameterization. All methods have relative errors below 25% in about 45% of data hours. The usual errors are connected with the short periods when the sky fast becomes either overcast or clear, or when solar elevation is low, especially during the winter period.
Although it is not very obvious which clear sky model and model for cloudiness is better (Figure 3 and Tables 4 and 5), we clearly notice that: (1) Foken's clear sky model, together with Kasten and Czeplak [39] methodology, were the best during unstable stratification when the solar elevation was high, although a few of the greatest absolute errors happened then (Figure 3b).
(2) Holtslag's clear sky model and Burridge and Gadd [40] cloudiness parameterization were the best for low solar elevation, when, especially during sunset, Foken's calculation permanently overestimated the measurements.
One can conclude that both clear sky methods and both cloudiness methods are very useful tools for data quality control through intercomparison tests that check the relationship between cloudiness and global solar radiation at midday and separately at sunset or sunrise.

Assessments of Upwelling Shortwave Solar Irradiance
The upwelling solar radiation is part of downwelling solar radiation, which is reflected from the Earth's surface and the atmosphere. The ratio of reflected to incoming shortwave radiation is called albedo. It depends on the angle of solar elevation, cloudiness and composition and state of the Earth's surface. It is not measured on a routine basis. There are hourly values of the albedo as the fraction of the measured incoming shortwave radiation that is reflected in a Debrecen dataset.

Assessments of Albedo
In 1995, Geiger et al. [36] offered a classification where high solar elevation albedo depended on the surface type and varied from 0.024 for rough seas to 0.98 for clean snow.
Since there was no data in Debrecen for "state of the ground with/without snow or measurable ice cover"-WMO SYNOP code tables 0901 E and 0975 E', the hourly data about precipitation and present and past weather, as well as the temperature, must be used in order to classify the "state of the ground". Namely, the nearest stations that included the "state of the ground" in the SYNOP reports at 06 UTC were: the station WMO 13174 Kikinda (in Serbia) and WMO 11968 Kosice (in Slovakia). In this paper, the albedo A was calculated for three different values of A ∶ = 0.11 for wet soil (precipitation in present or past weather), = 0.25 for dry soil (nice weather) and A = 0.85 for snow (snow cover).
The equation that was adopted after reference [66] and used by Nyren and Gryning [37] expresses albedo, , as a function of the solar elevation angle (in degrees): The other similar empirical equation was used by references [67,68]. This type of methodology for the calculation of albedo that depends on the solar elevation and the soil wetness has been applied up to the present day [41].
Statistical errors for proposal albedos are the same for the first three decimal places, using both Equations (16) and (17). Due to that, only one set of results is presented. The calculated albedo compared with the measurements is illustrated in Table 6.
Without hours with precipitations and snow cover, a total of 3473 out of 4022, the statistical errors were smaller, while the 126 h with precipitation and 423 with snow cover without precipitation had expectedly greater errors.  [37] model of albedo for all data and separately for nice weather without snow cover ( = 0.25), for weather with snow cover and without precipitation ( = 0.85) and, finally, for precipitation weather ( = 0.11). Besides this calculation, we tried to find albedo as a function of solar elevation in Debrecen in a similar way as Duynkerke [69], as well as Beljaars and Bosveld [38]. They suggest that albedo depends on the angle of solar elevation and on the relation between diffuse and total incoming solar radiation, i.e., cloud cover.

Statistical Error for Albedo
Following their idea about the dependences of the albedo of the cloud cover, firstly, a linear function = ( ) for clear skies is calculated, using observations of albedo and Foken's calculation for solar elevation during clear sky conditions ( = 0). For the 178 h without snow cover, this linear dependency is: Following idea described in reference [38] alongside this equation, cloudiness is included to describe albedo in nonclear sky conditions using the least squares method for the different forms of the functions without and with snow cover surface. Beljaars and Bosveld [38] introduced the relation of downwelling diffuse incoming solar radiation, , and totally downwelling solar radiation, , as a factor with leading influence on the albedo. Since there is no solar radiation instrument equipped with a shadow band at Agrometeorological Observatory Debrecen, taking that diffuse radiation as 10% of the total downward solar radiation for clear skies ( = 0) and 100% for cloudy skies ( = 8), we supposed that their ratio rises linearly between 0.1 and 1.0 with the rising cloudiness. Therefore, it used (0.125 • − 0.1), instead of ( / − 0.1), given by reference [38]. For 3473 situations without snow cover and precipitation, when ≥ 0, we obtained: For 423 data with snow cover, we obtained: For hourly data both with precipitation and snow, we took = 0.2. The statistical errors obtained from these calculations were significantly smaller than in previous calculations, except in the case of precipitation (Table 7). Useful, statistical significance correlation was formed in the case of all data ( = 0.79), but in all three independent cases, the absolute values of the correlation coefficients were below 0.3, so there were no useful connections for the practices among the measured and calculated albedo. The absolute values of the BIAS, except for the few precipitation cases, were below 0.025, which is an acceptable low value. The shortwave reflected solar irradiance is computed as a multiplication of global solar radiation and albedo. The precision for evaluated albedo is very important, and we represented this fact by comparing the mean annual daily variation of calculated and measured upwelling shortwave radiations and through the statistical errors for two mentioned methods: first-Nyren and Gryning [37], Equation (16), and second-Beljaars and Bosveld [38], Equations (20) and (21), for the assessments of albedo (  Through all the statistical errors, the evaluation of albedo that includes cloudiness is much better than other methods. The advantage of the Beljaars and Bosveld [38] method is also obviously visible in Figure 4, which represents the mean annual daily variation of measured and calculated reflected solar irradiance.

Assessments of Net Irradiance
The net radiation of the ground surface is given by the sum of the shortwave and longwave radiation components. However, contrary to its importance, it is very rarely included in the standard meteorological measurements. In this section, a mean annual daily variation and statistical errors for measured and estimated * are discussed, where the estimated value depends on the global solar radiation, albedo, cloudiness and, in one parameterization, on the temperature. We analyze the frequently used classical parameterizations in the SVAT (Soil-vegetation-atmosphere transfer) models in the daytime ( ↓> 0) and in the nighttime ( ↓= 0) separately for a better understanding of the accuracy of the parameterizations in the Pannonian Basin.

Daytime
To estimate the net irradiance ( * , (W m -2 )) from the routine meteorological measurements, an empirical equation was taken, such as Holtslag and Van Ulden [5]. * = 0.9 (1 − ) where : atmospheric temperature (K) from a 2-m measuring height, : the Stefan-Boltzmann constant ( = 5.670 • 10 W m K ), = 5.3 • 10 W m K , = 60 W m and = 0.12. In this equation, we used calculations for Holtslag's clear sky calculations and changed both the cloudiness corrections and both albedos to obtain the statistical errors. Besides that, the empirical equation following the idea of Göckede and Foken [14] was used. They used the global solar irradiance as an input parameter [11].
In this equation is the measured global solar irradiance and can be substituted with ↓, is albedo, as earlier, ↓ is described in Sections 2 and 4.1 as the downwelling solar radiation for clear sky situations based on Foken [11] and is the air density (kg m -3 ). Based on the ideal gas low, where measures the station pressure in (hPa), = 287.06 J kg K is the specific gas constant for the dry air and is the virtual temperature in (K).
where is the specific humidity (kg kg -1 ) that is calculated as = 0.622 is the water vapor pressure in the same unit as the pressure (hPa), which is calculated based on the relative humidity ( ℎ) and the saturated water vapor pressure ( ) in the temperature, , is calculated as follows:  (27) for ice.
The water vapor pressure is calculated as = • ℎ/100, and finally, we give the formula of the specific heat of moist air at a constant pressure: Furthermore, we compared Holtslag and Van Ulden [5] and Foken [11] equations for * , using measured global solar radiation and Beljaars and Bosveld [38] methodology for calculation of the albedo.
The following statistical errors are calculated and shown in Table 9: BIAS, MAE and RMSE, together with the percentage of cases in which the relative errors were fewer than 25% (Rerr < 25%) for the calculation for clear sky global solar radiation, cloudiness calculations and albedo, where albedo ′ does not depend on the cloudiness, whereas does. The second numbers in the column represent statistical errors of the net radiation but only for periods with high solar elevation when the sensible heat flux ( ) was positive. The moment when the sensible heat flux changed its direction was calculated, according to De Bruin and Holtslag [10], from equation: where = 20 W m is the constant, is the ground heat flux transported into the soil, = / and the psychometrics constant is = / , where λ is the latent heat of vaporization. The rough estimation of the ground heat flux transported into the soil is = 0.1 • * [70]. There are a few approximated relationships between and [11,12,37]. We checked all the approximations in our dataset and chose the second method [12]. Table 9. Statistical errors and percentages of cases in which the relative errors were smaller than 25% (Rerr < 25%) in the calculations of the net irradiance, * . We used the next symbols in the table: Foken [11] net radiation model * and Holtslag and Van Ulden [5] net radiation model * . The low index used for the measured global solar radiation and , , and for its assessment depending on Foken's or Holtslag's clear sky calculation and cloudiness correction.
′ represent the Nyren and Gryning [37] and Beljaars and Bosveld [38]-type parameterizations of albedo. The second column inside the main columns represents statistical errors of the net radiation but only for period with high solar elevation when the estimated sensible heat flux based on De Bruin and Holtslag [10] was positive. Gray-methods with higher correlation coefficients. Although a significant difference is not noticed between Holtslag and Van Ulden [5] and Foken [11] calculations for net radiations, or between Kasten and Czeplak [39] and Burridge and Gadd [40] mathematics for cloudiness, it is clear that the inclusion of cloudiness in the calculation of albedo, as recommended in Beljaars and Bosveld [38], improves the estimation of net radiation. It is obviously confirmed that the inclusion of measured global solar irradiance, as suggested in Göckede and Foken [14], significantly improved the previous estimations of net radiation when these measurements were not available, as is expected.

Net
Furthermore, describing a successful estimation of net radiation can help the estimation of daytime sensible heat flux from routine meteorological measurements. It is especially useful at locations where special profile or flux measurements do not exist. In addition, these methods can be very useful tools for data quality control for special micrometeorological measurements.

Nighttime
The nighttime net thermal radiation depends on the temperature and cloudiness. We used an equation that was adapted from Burridge and Gadd [40]: where is the measured temperature at a 2-m height (K), = 285 K , * = −91 W m and = ( + )/16 is the effective cloud cover, where is in octa (0, 1,…, 8).
Although nighttime net thermal radiation is not usually used for the calculation of sensible heat flux from routine measurements, we offered statistical errors for its estimation: = −24 W m , = 26 W m , = 31 W m and the coefficient of correlation was 0.74. The small difference between the absolute value of and represents a systematic underestimation, as we can see in Figure 5. We represented only * net radiation models with Foken's clear sky calculations, * net radiation models with Holtslag's clear sky calculations and both models with measured , emphasizing the differences between ′ and albedo based on Nyren and Gryning [37] and Beljaars and Bosveld [38] methodology separately. Figure 5 represents the mean annual daily variation of net radiation, which clearly confirms the facts that were discussed through statistical errors and suggests the following conclusions: (1) During the midday period, when sensitive heat flux is directed upwards, Foken's calculation is clearly better than Holtslag's calculation.
(2) Holtslag's method is slightly better for low solar elevation when a sensible heat flux is directed downwards.
(3) The difference between the measured ( ) and calculated ( * ) net irradiance is very small when we compared both methods for mathematically described cloudiness.
(4) An albedo that includes cloudiness ( ) is much better than one ( ′) that does not.
(5) Measured global solar radiation makes an estimation of the net irradiance significantly better in comparison to when this value is estimated.
(6) Using the same value for global solar radiation, Foken's estimation for net irradiance is slightly better than Holtslag's. (7) The estimation of nighttime net thermal radiation is crude and has systematic errors.

Assessment of Longwave Irradiance
The absorbed solar energy on the Earth's surface is emitted into the atmosphere as upwelling LW radiation ( ↑). The part of this radiation and part of solar radiation absorbed by cloudiness and atmospheric particles and gases go back and/or forward to the Earth's surface as downwelling LW radiation ( ↓). The fluxes of longwave radiation >4 μm are key terms in the surface energy budget and vitally important in applied meteorology. Its knowledge is essential for the forecast of frosts, fogs, temperature variations, etc. [71,72].
Although parameterization schemes for longwave radiation have their background in the irradiative transfer theory, most of them, as in the case of shortwave radiation, are based on empirical relationships derived from observed fluxes. Empirical relations usually describe LW radiation well, especially in climatic conditions similar to those from which they were derived and in clear sky conditions. Cloudiness correction, again as in the case of SW radiation, has to be included in the calculations.

Assessment of Upwelling LW Irradiance
Outgoing longwave radiation from the Earth's surface is driven by the effective radiative, infrared surface temperature, : where ~0.98 is the Earth emissivity, and is Stefan-Boltzmann's constant. is very rarely available from routine observations. Due to this, it is necessary to substitute this value with the screen level atmospheric temperature, , which is always smaller during the day and greater during the night. Sellers [52] and Holtslag and van Ulden [5] described upwelling LW radiation using a Taylor series centered at : Since we consider the above equation over short grass, we assume = 1; hence, the effect of a gray body is neglected. Since feels the change of with a delay, especially during the day period, we used atmospheric temperature with a one-hour delay when the solar elevation was positive. For the daytime case, the last term can be approximated using either ↓ or .
Since more and more routine meteorological measurements are including global solar radiation, the last approximation is used more often. According to Foken [11] and Offerle et al. [53], who modified the first form of Holtslag and van Ulden [5] approximation by including albedo ( ): We also used the same mathematical expression approximated by * according to Holtslag and van Ulden [5]: Instead of a constant 0.12 in this equation, the same authors used values that depended on atmospheric and surface conditions, especially soil moisture: where and were mentioned in previous sections, and , is the modified resistance for sensible heat flux. The default value : , = 80 s m , as noted in reference [5], is a psychrometric constant, and is the Pristly-Taylor coefficient, which is a function of the soil moisture, as described in reference [73].
= 0.1 and = • * represents the part of the net radiation that goes into the soil.
Although the soil moisture, as well as the surface soil temperature, is not part of the routine weather observations, in this paper, it is compared by measuring the outgoing LW radiation with the estimated ↑ in a few different cases: There are 7933 data hours calculated and compared. Contrary to expectations, we found that the calculations with measured did not have the best evaluation of ↑. Unexpectedly, it is evident that including soil moisture measurements through did not improve our calculations (Table 10). In that case, three numerical experiments were done: the first experiment = 1.0, the second experiment: = 0.95 from April to September and = 0.85 from October to March and the third experiment = 0.65. The last one was a little better than the first and the second, but it was worse than when we used the constant. According to the assessment of upwelling LW from the routine measurements, it was compared to the Offerle et al. [53] parameterization, which used (Equation (33)) with measured ( ↑ ) and calculated ↓ ( ↑  ) with Holtslag's parameterization (Equation (34)) with * , where the net radiation was calculated again with the measured ( ↑ * ) and calculated ↓ ( ↑ * ). It is evident that Holtslag and van Ulden [5] parameterization with net radiation is the best in both cases: with and without measured ( Figure 6 and Table 10). Figure 6. Mean annual daily variation of ↑ for four models and measurements. Models that used routine meteorological measurements included, separately: ↓ or . Abbreviations of the models are same as in Table 9. In this paper, eight methods for modeling clear sky downwelling LW radiation as a function of screen level temperature, , expressed in (°K), and water vapor pressure, , expressed in (hPa), are compared with measurements through statistical errors. The applied parameterizations are presented in Table 11.  [43] parameterization, ↓, was the best (Table 12). Similar to ↓, similar errors have parameterizations from Idso [45], ↓, although their errors were the smallest during the day period and almost the worst during the night, when clear sky conditions were more frequent. Firstly, in order to include cloudiness in our calculations, the eight parameterizations for clear sky , ,… , ↓ are combined with three different methods that describe the effects of cloudiness. The first group of parameterizations was constructed on the basis of cloudiness correction, following Niemelä et al. [43]. There are three parameterizations, where: Based on Jacobs [49] and Maykut and Church [50], separately. The third cloudiness correction method was for a lowland site, following Iziomon et al. [51]: All clear sky downwelling LW radiation parameterizations ( , ,…, ↓), (see Section 6.1.1) were applied in cases ( ↓) , ( ↓) & and ( ↓) . There were (3 • 8 = 24) different cases for hourly downwelling LW irradiance calculations for all the sky conditions.
The downwelling LW radiation, as described by Holtslag and Van Ulden [5], was calculated using clear sky parameterizations, where the downwelling irradiance was a function of from Swinbank [47] and Dilley and O'Brien [48] separately: Finally, following Niemelä et al. [43], the downwelling LW irradiance was calculated for all sky conditions, with a parameterization that was developed using clear sky downwelling irradiance, cloudiness and upwelling LW irradiance: There were eight cases for each hour. Based on the intention to work with routine measurement data and considering the results from previous sections, the upwelling LW irradiance in Equation (41) was calculated as a function of and * , while the net irradiance was calculated by Foken's clear sky model, Kasten and Czeplak [39]-type cloudiness parameterizations, Beljaars and Bosveld [38] albedo and Göckede and Foken [14] calculations for net irradiance, ↑ * .
In general, comparing by the statistical errors only for cloudiness correction, we concluded that the Maykut and Church [50] method was slightly better than Jacobs [49]. At the same time, the cloudiness correction described by Iziomon et al. [51] was between them with: BIAS between 0 W m -2 for ( ↓) and 17 W m -2 for ( ↓) , MAE between 18 W m -2 for ( ) and 26 W m -2 for ( ↓) , RMSE between 22 W m -2 for ( ↓) and 32 W m -2 for ( ) and a correlation coefficient 0.91 for almost all parameterizations except ( ↓) and ( ↓) , which were less. Slightly better statistical errors than all three above-mentioned assessments of downwelling LW irradiance were described by Niemela et al. [43] and shown in the grey column of Table 13. The number of hours with absolute errors greater than 100 W m -2 for Niemelä et al. [43] parameterization was very stable and varied between 24 for Dilley and O'Brien [48] to 29 for Swinbank [47] and Iziomon et al. [51] clear sky downwelling LW parameterization. According to the relative errors, it was noticed that, for almost all parameterizations, the relative errors were smaller than 25%. However, in order to distinguish slight differences between them, percentages for the relative errors smaller than 10% are offered in Table 14. The columns represent different clear sky parameterizations, while the rows represent additional cloudiness corrections. ). It is noticeable for the downwelling irradiance, according to Niemelä et al. [43], that all sky cloudiness conditions, together with Dilley and O'Brien [48] ( ↓), Niemelä et al. [43] ( ↓), Iziomon et al. [51] ( ↓) or Pratta [5] ( ↓) parameterization of the clear sky LW irradiance, gave an excellent estimate (Table 12 and Figure 7).  Table 8. Two parameterizations: ( ( ↓) and ( ↓) ) based on Holtslag and Van Ulden [5] are represented in colored lines. See, also, Tables 11,12 and 14.

The Statistical Errors and Validation of the Results
The stability validation of our previous results was done by comparing the statistical errors calculated for one year with the same values calculated for longer periods. Namely, we compared BIAS, MAE, RMSE and the coefficient of correlation ( ) calculated for the year 2009 with the same errors for the two longer periods. The first period was from 2008 to 2010 and included 17,964 data hours for net radiation and 9814 data hours when > 1 W m −2 . The second period was from 2008 to 2017 and included 61,880 data hours for LW radiation and 40,869 data hours with > 1 W m −2 (Table 15). We used the parameterization with the proven best results.  The results of the interpretation were separated, because there was the transition from manual to automatic measurements after the year 2010 at Debrecen Airport. It had the consequence of missing cloudiness measurements during the night. Namely, hours with a lack of cloudiness had to be excluded from the calculations, as well as the hours with accidentally missed data. During the period 2008-2017, there were night measurements of cloudiness only in the first half of the period-in the beginning completely and, later, only partly.
Since 2018, the measuring site Debrecen Airport WMO 12882 has been operating as a "backup" site with limited data access. Finally, the comparison between the estimated * and measured irradiance was conducted in four different cases, and the results are shown in Tables 15 and 16 and   Comparing statistical errors for a 10-year time period, 2008-2017, the stability of the presented calculations was clearly confirmed (Table 17). The same statistical errors for the same variables for the time period from 2008 to 2010 are shown in Table 16 and Figure 8. It is obvious that the correlation coefficient is better for a longer period for all the variables, while the statistical errors MAE and RMSE looked clearly better only for * , * and * . Again, one emphasizes the significance of the measurements of global solar irradiance for the estimation of net irradiance. Namely, Q * and Q * have smaller statistical errors and better coefficients of correlation than * and Q * , although we believe that the estimation of global solar radiation can become better if we include astronomical parameters for the EQT (Equations of Time) for our geographical location instead of = 15° and = 50°. We recalculated the solar elevation angles using EQT (in minutes) from reference [74] and compared the original Foken's [11] and Holtslag's [5] methodology with modified Foken's methodology through the application EQT from reference [74] in Equation (11) during the year with an hourly time step. The differences between the angles of solar elevation for this new calculation and original Foken's methodology were much smaller compare with Holtslag's methodology. The absolute differences were 0.84 ± 0.72° and 3.0 ± 1.6° C, respectively. We underlined the importance of the precise calculations of solar geometry.
Besides that, it was shown that net radiation ( * ) is little better if it is calculated through the sum of its components than if it is calculated from the first step of parameterization-directly, Tables 16 and 17 and Figure 9.

Conclusions
Comparisons of widely used parameterizations of shortwave and longwave radiation balance components in micrometeorology were provided based on the long-term measurements in the Northeastern part of the Pannonian region in Agrometeorological Observatory Debrecen. Quality-controlled datasets from WMO first-class radiation sensors were used during the time period 2008-2017.
The method for estimating the global solar irradiance (K ↓) can be used for data quality control of the measured parameters and for intercomparison tests between the global solar irradiance ( K ↓) and cloudiness ( ). All the methods had relative errors below 25% in around 45% of the data hours. The usual errors were connected with the short periods when the sky fast became either overcast or clear or when solar elevation was low, especially during the winter period. Individual overview of the data, synoptic situation, weather and cloudiness measurements, which were provided by Debrecen Airport (12882) 8.5 km from the Agrometeorological Observatory, are necessary in situations when the errors between the measured and calculated K ↓ are significant. In these cases, decisions can be made whether the measured or modeled data are acceptable.
Comparing nearly 62,000 h of observations with all sky conditions, we concluded that:

•
Foken's clear sky calculations [11], together with Kasten and Czeplak's [39] methodology, are the best during unstable stratification when solar elevation is high.

•
Measured global solar irradiance makes the estimation of the net irradiance significantly better in comparison to when this value is estimated.

•
The estimation of nighttime net thermal irradiance is crude and has systematic errors.

•
Dilley and O'Brien's [48] clear sky parameterization with Holtslag and Van Ulden [5] cloudiness corrections had the smallest statistical errors for the estimation of downwelling longwave irradiance.

•
Holtslag and Van Ulden [5] modified by Offerle [53] had the smallest statistical errors for the assessment of upwelling LW irradiance.
The method for estimating net irradiance is significant, because it is a variable that dictates the sensible heat flux and flux of the momentum when the measured surface wind is taken into consideration in the surface layer of PBL. Processes inside this layer are typically described by the laws of similarity theory, which lean on the Buckingham π-theorem. The most famous of these theories of Monin-Obukhov use the mentioned fluxes to determine the characteristic turbulent length scale on which classes of atmospheric stability are typically based. It would be valuable to know this scale and the components of the radiation and other surface layer parameters (SLP) by the same values from the numerical models, as well as to discuss the spotted differences.
We believe that the observed differences emphasized the importance of special flux and radiation measurements in the Pannonian Plain, where the wind is weak and the duration of extremely unstable and extremely stable stratifications of PBL, which is still a discussion topic, are far more common than in Western Europe. The chosen methods can be used for data quality control and filling the gaps in these measurements.