A Numerical Study of Effects of Radiation on Deep Convective Warm Based Cumulus Cloud Development with a 3-D Radiative Transfer Model

The effects of radiation heating and cooling on cumulus cloud development have been the focus of considerable attention for many years. However, it is still not clear how radiation impacts cloud droplet growth. Since cloud inhomogeneity has a great influence on radiation transmission, we coupled the 3D atmospheric radiative transfer model using the spherical harmonic discrete ordinate method with WRF-LES, which can improve the simulation accuracy of the inhomogeneous effect of clouds on radiation compared with that of the 1D radiation method. The shortwave and longwave radiation fluxes for upward and downward directions were simulated with different solar zenith angles. The comparison of 1D and 3D radiative solvers for deep convective cloud cases shows that the 3D radiative solver provides an accurate structure of solar and thermal radiation characteristics and the spatial distribution field. The solar radiation heating is likely to increase perpendicular to the solar incidence direction. For longwave radiation, the cooling effect on the cloud top and the heating effect on the cloud base are both more intense in the 3D radiation model. This study focuses on 3D cloud-radiative interactions in an inhomogeneous cloud field in a large eddy simulation, and the results suggest that compared with the widely used 1D radiative solver in WRF, the 3D radiation model can provide a precise description of the radiation field in an inhomogeneous atmosphere.


Introduction
Radiation and cloud interactions have been widely recognized as one of the key factors influencing the microphysical processes in atmospheric science [1][2][3][4][5][6][7][8][9][10][11][12][13]. The clouds absorb and scatter solar shortwave and longwave radiation, and emit longwave radiation simultaneously, which affect the short-term radiation flux change and long-term radiation balance. In return, a number of numerical and observational studies have suggested that the impacts of radiation on clouds are based on the gradient of the heating rates which includes the gradient inside the cloud, between two clouds or between the clouds edge and clear sky, resulting in the changing of cloud formation [14][15][16][17][18][19][20]. Therefore, it is increasingly important to improve the ability to describe cloud microphysical processes by accurately describing cloud radiation patterns. processes and the time-consuming calculations, the effect of coupling 3D and 1D radiation models on the calculation accuracy of mesoscale models has rarely been verified.
Although the 3D radiation model computes radiation characteristics that can achieve high precision and considers the inhomogeneous cloud medium with a higher angular resolution, it also has several flaws such as a large amount of calculation and a complicated radiation transfer equation [44]. The research on the 3D radiation model has developed in recent decades. Kattawar preliminarily attempted to calculate 3D radiance using the Monte Carlo method and discussed the dependence on the parameters for cloud models and radiation parameters [45]. Liu et al established a 3D Monte Carlo transport model of polarized radiation in the microwave spectrum to study the influence of rain clouds on the 3D microwave brightness temperature and pointed out that brightness temperatures are reduced in the 3D case due to the net leakage of radiation from the sidewalls of the cloud [46]. Gautier et al used a Monte Carlo based 3D radiative transfer model to investigate the solar radiation within a cloudy atmosphere [4]. The analysis of a tropical cloud field showed that the cloud morphology influences the absorption of ultraviolet radiation and the absorption of water vapor in the overhead sun. Klinger used large eddy simulations to investigate the effect of 1D and 3D longwave radiation on cloud droplet growth in cumulus clouds over six hours [23]. The results showed that the effect of 3D thermal radiation on rain production is stronger than that of 1D thermal radiation.
At present, it is still not clear how radiation impacts cloud droplet growth. The Weather Research and Forecasting (WRF) model is a numerical weather prediction (NWP) system which is developed by NCAR [47]. In regard to the detailed analysis of the cloud structure and cloud radiation characteristics in the 3D model, there are few studies on the coupling of the 3D radiation model and (Weather Research and Forecasting) WRF mesoscale model for large eddy simulations. The reasons for the lack of adoption of the 3D radiation model in most large eddy simulation models are the computation speed and calculation amount considerations. In contrast to 1D atmospheric radiation transfer, solving the radiative transfer equation in a 3D medium is a five-dimensional (three spaces, the zenith angle and azimuth angle) boundary value problem [44]. Due to the high spatial and temporal resolution of large eddy simulations, the spatial resolution is less than one kilometer and the temporal resolution less than one minute. Coupling a full 3D radiation transfer solver to an LES model is extraordinarily computationally expensive. Moreover, the radiation transmission must be integrated within the solar and thermal spectra. Only computationally efficient, accurate, and flexible 3D radiative models are appropriate for large eddy simulations. The spherical harmonic discrete ordinate method (SHDOM) radiative solver is a widely used 3D radiative model that can compute unpolarized spectral band radiative transfer in 1D, 2D, and 3D media for either collimated solar or thermal emission sources of radiation with high efficiency [48][49][50][51]. The correlated-k distribution is used for the integration, and scattering tables are adopted in the model, which improve the computational efficiency. The variable grid is available, and the calculation accuracy is adjustable, which allow exploring the trade-off between calculation speed and accuracy.
The (Advanced Research Weather Research and Forecasting) ARW-WRF is the ARW dynamics solver and other components of the WRF system compatible with the solver used to produce simulations. In this paper, we coupled a 3D radiative solver (SHDOM) in ARW-WRF with a mixed-phase microphysical scheme. The high-resolution large eddy simulations are used to study the radiation characteristics of cumulus clouds. The current study extends previous work by responding to sufficiently complex physical realities. This paper is arranged as follows. Section 2 briefly introduces the WRF model and the SHDOM radiative solver. Section 3 gives a description of the experimental design. Section 4 then proceeds to compare the 1D and 3D simulation results with the newly coupled radiation model. Section 5 is discussion part and Section 6 presents the conclusions.

The Model Introduction
The Weather Research and Forecasting (WRF) model is a numerical weather prediction (NWP) and atmospheric simulation system which is a non-hydrostatic compressible atmospheric model with a capability of grid nesting. The advanced research WRF (ARW) modelling system used in our study is the ARW dynamics solver together with other components of the WRF which was developed at the National Centre for Atmospheric Research (NCAR) and the National Oceanic and Atmospheric Administration (NOAA) [52][53][54][55] and has been widely used for both meteorological research and numerical weather prediction for the past two decades. The model has several initialization programs for both idealized and real-data simulations cases and provides several parameterization options such as land surface, surface layer, planetary boundary layer (PBL), microphysics, cumulus parameterization. WRF solves the fully compressible, nonhydrostatic Euler equations formulated on the terrain following hydrostatic-pressure vertical coordinates. The equations for the momentum, potential temperature, subgrid-scale (SGS) turbulence kinetic energy (TKE) and other scalars has a form coupled with the dry air mass. The basic predicted variables include velocity, pressure, temperature, humidity, and SGS TKE. In microphysical parameterization scheme, variables such as mass mixing ratio and cloud droplet concentration are also predicted. WRF uses the Runge-Kutta second and third-order time integration schemes with smaller time step for acoustic and gravity-wave modes. During each Runge-Kutta step, the high frequency mode of horizontal propagation is combined with the forward and backward scheme of the acoustic time step, which is usually one order of magnitude smaller than the physical time step, the high frequency mode of vertical propagation adopts implicit mode. For horizontal grid, Arakawa C-grid staggering is adopted, geopotential is located at the vertical velocity point, while pressure and density are located at the scalar point. The fifth-order advection scheme is used in advection calculation, and the third-order advection scheme is used in vertical advection calculation. [56]. The diffusion due to turbulent mixing is represented by the eddy diffusivity (K-theory) in the height coordinates of the coupling with the dry air mass, which involves a coordinate transformation from the mass coordinates to the height coordinates of the horizontal derivative. A 1.5-order TKE closure is used as a SGS turbulence parameterization.
The Version 4.0 is used in this study and has been available for download since June 2018. The monotonic technique is employed for advection of scalar and moist variables. WRF is suitable for use in a broad range of applications across scales ranging from kilometers to meters, including ideal case and real case simulations, and a number of options for parametric radiation, the planetary boundary layer, microphysics, dynamics, and land surface interactions is available as the need arises. The large eddy simulation of WRF (WRF-LES) used in our study is a technique in which large scales of atmospheric turbulence are resolved, while the smaller-scale turbulence from the flow field is removed using a spatial filter [57]. WRF-LES is able to provide a sophisticated model framework, including multiple types of physical schemes to parameterize microphysics, dynamics, radiation, and cumulus processes, and reproduce high-resolution numerical simulation results to the greatest extent while saving computing resources for the simulations. The detailed description of WRF-LES can be found in Stevens et al [58]. The microphysical scheme and the newly coupled radiation scheme used in this study are described below.

Microphysics
The spectral bin microphysics (SBM) scheme adopted in this study was developed by Khain at the Hebrew University of Jerusalem [59,60]. This scheme is a kinetic equation system based on solving the size distribution function of several types of aqueous particles, including water drops, three types of ice crystals, aggregates, graupel, and hail, as well as the distribution function of condensation nucleus (CCN) aerosol particles. The microphysics model is an equation system based on solving a size (number) distribution function f i , the size distribution function f ik of a kth mass category is described as follows: where V tik is the terminal velocity of hydrometeors of type i belonging to the kth mass bin, the right terms represent the rates of the size distribution function change due to nucleation, condensation or evaporation of cloud drops, deposition or sublimation of ice particles, collisions, freezing or melting, and breakup, respectively [59]. The CCN activation spectrum of the initial time step is given by a semi-empirical formula: where S is supersaturation (%), and N 0 and k are measured constants, the different magnitudes of No and k is used within different ranges of supersaturations by empirical dependence. Supersaturation is the work of Tzivion and other predecessors, and it is analyzed and extended by using an accurate calculation method [61]. Using supersaturated values, the critical size of the aerosol particles was calculated as activated droplets. Unlike the standard parameterization bulk scheme, the cloud size distribution of hydrometeors and aerosols is not predetermined but calculated in the process of model integration. Each size distribution function contains 33 bins. The fast version of SBM was used in our study, and the structures of microphysics and dynamics are similar to the full SBM, requiring approximately 40% of the computing power [62]. The outstanding advantage of this scheme is that each aqueous particle is described by a defined size distribution function on a grid containing 33 bins. The provided number concentration of cloud droplets is crucial to the calculation of cloud droplets' effective radius. A description and details of the parameterization were noted by Khain [63,64]. The effective radius is defined as the area weighted mean radius, as shown in Equation (3) below: where r is the radius, and n(r) is the gamma distribution of cloud droplets with respect to r. α and β specify the gamma distribution in which the typical value for water clouds is α = 7 and β is related to the concentration and content of cloud droplets. For water clouds, the cloud water content is determined by Equation (4) where ρ w is the density of the cloud water. In this case, the effective radius of cloud particles is a function of the cloud water mixing ratio and cloud number concentration. The cloud particle distribution will affect the amount of solar radiation reflected back to free space. In particular, with a constant cloud water mixing ratio, smaller cloud drops reflect more solar radiation to space and cause the cooling of the atmosphere, whereas larger droplets reflect less solar radiation to space and allow more solar radiation to penetrate to the ground, resulting in more radiation being absorbed by the ground [65]. A reasonable microphysical process can provide a more accurate background field for the radiation process, which is critical to the accuracy of radiation calculations.

Radiation Scheme
Observational and simulation results indicate that cumulus clouds play an important role in the radiation budget [66][67][68]. In addition, radiation feedback influences the formation and evolution of cumulus clouds [31]. However, due to complexity and computational efficiency, the radiation scheme generally adopts the assumption of local plane-parallel atmospheric radiation in WRF, the radiation intensity and atmospheric parameters only change in the vertical direction and the horizontal interaction is ignored, which is problematic and not suitable for radiation transmission, particularly in an inhomogeneous atmosphere [69]. To that end, one such strategy is to couple the 3D radiative solver to WRF online to achieve high precision, considering the inhomogeneous cloud medium and performing computations in more angular directions. The spherical harmonic discrete ordinate method for 3D atmospheric radiative transfer (SHDOM) is a widely used 3D radiative transfer model in atmospheric science that has been coupled to WRF in our study. SHDOM uses the extinction coefficient, the single scattering albedo, the temperature, and the expansion coefficients of the phase function in terms of Legendre polynomials as input parameters. The integral form of the radiative transfer equation shows that The extinction k and source function (Jk) are assumed to vary linearly with distance s. The radiative transfer equation of radiation quantity I(s) in the direction of the ray at the distance s along the ray is expressed in the differential form of dI(S) where k is the volume extinction coefficient. The extinction varies linearly with distance across the grid cell: SHDOM can compute unpolarized monochromatic or spectral band radiative transfer in a one-, two-or three-dimensional medium for either collimated solar or thermal emission sources of radiation. The radiative solver uses a correlated-k method to calculate the absorption and scattering of molecules and cloud hydrometeors in broadband spectra. Furthermore, the spherical harmonic and discrete ordinate is used to represent the source function and radiation field during the computation, which is efficient to save calculation time.
The source function in spherical harmonics space J lm is related to the source function in angular coordinates: where Y lm are orthonormal spherical harmonic functions, µ is the cosine of the zenith angle, and φ is the azimuth angle. Because the scattering is mainly determined by the scattering angle, the scattering integral of the source function in the spherical harmonic function coordinate system is transformed into the following simple multiplication: where l is the radiation coefficient of the spherical harmonic function, χ 1 is the coefficient of the Legendre phase function, and ω is the reflectance of single scattering. The properties of the medium can be specified completely generally, i.e., the extinction, single scattering albedo, Legendre coefficients of the scattering phase function, and temperature for the particular wavelength or spectral band may be specified at each input grid point. The stability of the model is controlled by a Picard iteration, and for each iteration.
(1) The source function is transformed to discrete ordinates from spherical harmonics at each grid point.
(2) The source function is integrated using the radiative transfer equation in discrete ordinates to obtain the radiation field at each grid point. (3) The radiation field is transformed to spherical ordinates. (4) The source function is computed from the radiation field in spherical harmonics. To balance the accuracy and efficiency of the calculation, the radiance is calculated in a certain order in SHDOM, and the specific content and algorithm will be introduced.
The longwave and shortwave rapid radiative transfer models (RRTM) were developed at Atmospheric and Environmental Research and have been proven to be highly accurate compared with the line-by-line radiative transfer model [70]. The RRTM method is applied to SHDOM by the k-distribution for atmospheric gases treatment, which produces a file that includes absorption molecules integrated with each broadband spectrum. The shortwave and longwave broadband spectral wavelength ranges follow the RRTM k-distributions and are displayed in Table 1. The longwave broadband spectra are divided into 16 bands from 10 microns to 3000 microns, the shortwave broadband spectra are divided into 14 bands from 0.2 micron to 12 microns. For the main trace gases, H 2 O, CO 2 , O 3 , N 2 O, and CH 4 are considered. To reduce the amount of calculations, the single scattering properties of the gamma distribution of spherical particles, including aerosol, cloud water, and cloud ice, were calculated by referring to the Mie scattering table prepared in advance. If the particle type is water or ice, the integral can be performed across the specified wavelength band. If the aerosol particle type is selected, the refractive index of the aerosol is specified. Using the Mie scattering table can save a substantial amount of calculation time and is effective in large eddy simulations. The three-dimensional grid containing the liquid water content and the effective radius of the droplet are stored in the cloud field. The effective radius was calculated by the cloud water content and the number concentration of each grid point provided by the double-moment bulk mixed-phase microphysical scheme in WRF. SHDOM combines the atmospheric gas information with the cloud particle information to generate a medium property, including extinction, single scattering albedo, Legendre coefficients of the scattering phase function, the specific wavelength or the spectral band temperature, which can be stored in the property file calculated at each grid point. Finally, the radiation transmission results of each broadband spectrum are obtained by using the radiation transmission equation and Picard iteration to calculate the radiation intensity and heating rate. The radiation transmission results of 14 shortwave broadband and 16 longwave broadband spectra are separately output or accumulatively output, as required.
SHDOM maintains a balance between computational speed and accuracy and is proven to be superior to Monte Carlo radiative transfer methods when many radiative quantities are considered, which is suitable for cloudy sky conditions in WRF. The boundary conditions can be open or periodic. Selecting adaptive grids can increase the number of grids in important areas. It also provides a good tool for studying the influence of 1D radiation in comparison with 3D radiation by calculating with the independent pixel mode. The amount of calculation of the 3D radiative solver in SHDOM is mainly affected by the various angular and spatial resolutions. For example, in the same spatial resolution, the calculation speed is reasonably fast under a low angular resolution (such as four-stream), whereas an increase in the angular resolution will exponentially increase the amount of computation. Therefore, it is very important to choose appropriate angular and spatial resolutions to balance the calculation amount and precision. The details of SHDOM are described by Evans [44].

Experimental Design
On the basis of a detailed analysis of the cumulus structure, the 3D radiation scheme (SHDOM) was coupled to WRF online to study the radiation characteristics of cumulus clouds, including the developmental life history of cumulus clouds. For radiation settings, the concentration of CO 2 , CH 4 and N 2 O are set to 370 ppmv, 1.7 ppmv and 0.3 ppmv, respectively, at the surface, and the ozone density profile follows the American standard atmosphere. Because of the complex shape and radiation characteristics of ice crystal particles, the mixed phase states of water droplets and ice crystals will affect the analysis of simulation results to some extent. We focus our study on the cloud droplets in clouds, which do not contain ice crystals. The Mie scattering table of spherical water droplets is set as the gamma size distribution within the effective radius. The gamma distribution shape parameter is set to 7, and the typical values for the effective radius range from 0.5 microns to 25 microns, which is suitable for the most effective radius since the maximum radius in the distribution is 75 microns, which is large enough to fit large aqueous particles in the situation [71][72][73][74]. The horizontal boundary condition in SHDOM is open in x-coordinates and y-coordinates, and the radiation in the open horizontal boundary condition that leaves the sides does not wrap around. Since the accuracy of the model calculation is driven by the angular resolution, the number of zenith angles is set to 8 in both hemispheres, and the number of azimuth angles is set to 16 to ensure the accuracy of the simulation. To compare the 1D model with the 3D model and avoid introducing other errors from the model initialization, the 1D model is chosen in SHDOM and all other variables are set the same as the 3D model. The spectral bin microphysics (SBM) scheme adopted in this study was developed by Khain at the Hebrew University of Jerusalem [63,64]. A more complete microphysical process can provide a better background field and accurate feedback for the radiation process.
In this section, we initiate the large eddy simulation with a thermal bubble perturbation and employed the radiation setting in WRF. The sounding data are from Yau [75], which is a typical deep convective atmospheric junction, and Figure 1 shows the temperature and dew point temperature profile of the initial sounding. The temperature is very close to the dew point temperature from the surface to the model top, especially at the height of 910 hPa which is conducive to condensation. Moist air and a high dew point temperature in the whole model layer are favorable for condensation of water vapor. To better verify the coupling effect, a time period with strong solar shortwave radiation accompanied by longwave radiation is selected. The large eddy simulation begins at 12:00 UTC and ends at 12:30 UTC, such that the model runs for 30 min and includes the lifetime of the thermal bubble from formation to collapse. The simulation time step was one second, the number of horizontal grids was 81 with a resolution of 150 m × 150 m, number of vertical layers was 81 with the model top by 12 km. Moreover, we interpolate three extra levels, including 15 km, 20 km, and 30 km, to ensure the calculation of atmospheric optical thickness starts from the stratosphere in the radiation process. In the 1D radiation model, solar radiation is vertically downward and upward and is independent of the zenith angle and azimuth angle, whereas in the 3D radiation model, the solar incident direction is highly related to the zenith angle and azimuth angle. As a consequence, in order to avoid the difference caused by the solar incident direction, we conducted two experiments. In case one, the zenith angle is set as zero in both the 1D model and 3D model, and the scattering structure of solar radiation is highly symmetrical. In this case, we can better compare the difference between the two models in a symmetric structure. In case two, the solar cosine zenith angle is outputted by WRF in real time from 12:00 UTC to 12:30 UTC, and the solar zenith angle varies from 22.99 • to 23.89 • , which can reflect a more realistic and complex situation in 3D radiation and better reflect the difference between the 1D radiation and 3D radiation models from the solar incident direction. The azimuth angle is 90 degrees in both cases, and the radiation solver is called every one min.

Simulation Results
There are two cases described in the following section. To better evaluate the difference between the solar zenith angle and the radiation transmission process, the solar zenith angle was set as zero, and the real zenith angle was obtained from 12:00 UTC to 12:30 UTC. The former describes in detail the scattering process of solar radiation under ideal conditions. The latter reflects the effect of the actual solar incident angle on the radiation process. Figure 2 shows the cloud water mixing ratio and cloud water number concentration in the 3D radiation case. The development of cumulus clouds was divided into three main parts: the generation period from 0 to 5 min, the development period from 5 to 15 min and the decline period from 15 to 30 min. At 5 min, the clouds began to form, and the cloud top was approximately 1.6 km high with a low cloud water mixing ratio and cloud water number concentration. At 15 min, the cloud top reached up to 5 km high with the strong cloud water content and cloud water number concentration, and the maximum cloud water content reached 7.0 g/kg. At 20 min, the clouds grew higher and reached 7.5 km high, and the cloud water content decreased, while the cloud water number concentration continued to increase. The completely symmetry of the initial conditions and perturbation conditions for the large eddy simulations resulted in the symmetric distribution of the cloud water content generated by WRF. Transmittance and absorptivity are strongly dependent on the geometry of the cumulus cloud and the solar zenith angle. The vertical solar incidence is adopted by setting the solar zenith angles as zero. When the solar energy is vertically downward, the transmittance and absorption occur at the edge of the cloud top because of the spilling of the solar energy over the top of the cloud area, and the distribution characteristics of the radiation simulation results should be highly symmetrical. In this case, the influence of the solar incident angle on the simulation results obtained by the 1D and 3D radiation models will be eliminated. Compared with the clear sky atmosphere, the existence of the cloud particles has a great impact on the incoming shortwave radiation by absorbing and scattering, such that the shortwave shadowing is strongest when the zenith angle is zero and results in a strong cooling effect. Figure 3 presents the downward and upward shortwave radiation flux at 20 min. In the 1D radiation model, the downward shortwave radiation decreases from 1500 W/m 2 to 300 W/m 2 from the cloud top to the cloud base, and the downward shortwave radiation is less than 100 W/m 2 below the cloud base. In the cloudless region, there is a large horizontal gradient between the downward radiation of ambient air and the cloud region, and the large value region of the downward shortwave radiation is located directly above the cloud top. Since the assumption of the two-stream approximation is adopted in the 1D model, each column is independent and scatters only in the vertical direction; therefore, neither optical variability nor horizontal cloud particles are taken into account. The interactions between two cloud regions or between clouds and ambient air were not considered. Previous studies have shown that the multiple scattering of shortwave radiation is the main source of increased downward shortwave radiation [76]. In this case, the increased value of shortwave radiation exists directly above where the cloud droplet scattering occurs. Ignoring the horizontal transport will accumulate errors over time, and the horizontal interactions of photons will be significant particularly when a large amount of downward transmission is carried from the top of the model. In the 3D radiative model, the radiation solver calculates the radiation transmission of upward and downward fluxes, and additionally for sideward streams, the photons in the 3D model can travel in the air until encountering the new cloud particles that once deviated to the horizontal direction. The secondary scattering of the scattered radiation from the cloud top and the multiple scattering of the shortwave radiation transmitted horizontally from other directions contributed to large-value regions of the cloud top and cloud side that reached 1500 W/m 2 , and the horizontal transmission and multiple scattering also caused a radiation disturbance in the noncloud region near the surface. The upward shortwave radiation mainly comes from the scattering of downward solar radiation by cloud particles in the cloudy area, and the scattering intensity is closely related to the number concentration and mass mixing ratio of cloud particles. Due to the symmetrical distribution of structures inside the cloud, the maximum cloud water content exists in the cloud center with the strongest scattering ability and decreases from cloud center to cloud edge. In the 1D radiation model, there is no horizontal scattering except for the vertical direction in the two-stream approximation. As a result, the horizontal transportation of scattered solar radiation between two grid points is inhibited, the maximum upward shortwave radiation is located in the center of the cloud top and extends straight up to the model top, and the upward shortwave radiation gradually decreases from the cloud top to the cloud base inside the cloud. Furthermore, the defect of the 1D radiation scheme will lead to large values in the vertical direction above the cloud top and no horizontal radiation scattering on the cloud sides. In the extreme case, the upward shortwave radiation reaches 1400 W/m 2 at the model top as shown in Figure 3c which is not consistent with the reality. This faulty radiative transmission will lead to the accumulation of errors, especially in the presence of an adjacent cloud. In the 3D radiation model, the scattering intensity and area are enhanced by the growth of cloud droplets, the internal structure of the cloud is clear, and the reflected radiation shows a symmetrical structure. The large-value region is located at the cloud top as high as 1400 W/m 2 , and the upward shortwave radiation decreases layer by layer from the cloud top to 300 w/m 2 . The main reason for this phenomenon is that the scattering intensity of cloud particles depends on the phase function as a function of direction in the 3D model. When the solar radiation reaches the cloud top, cloud particles scatter solar radiation in both horizontal and vertical directions, which maximizes the upward shortwave radiation from the cloud top. Moreover, the scattered shortwave radiation containing a horizontal component travels horizontally and is scattered multiple times in the ambient air, which results in a large scattering region of shortwave radiation within a certain range above the cloud top.

The Effects of Vertical Incidence on Shortwave Radiation
The gradient of the net radiation flux affects the atmospheric heating rate and further influences the cloud microphysical process and dynamic process. The net shortwave radiation flux is shown in Figure 4, where the downward direction is positive. In the 1D radiation model, due to the attenuation of the downward shortwave radiation by the atmosphere, the net shortwave radiation vertically decreases from the model top to the surface in the cloudless region and there is no horizontal gradient in the net shortwave radiation flux in ambient air. In addition, in the cloud region, the spurious upward shortwave radiation above the cloud top results in less net shortwave radiation extending towards the top of the model. In the 3D radiation model, the cloud field is surrounded by clear and well-structured net shortwave radiation scattering, and the net shortwave radiation gradually spreads out from the cloud center. The downward shortwave radiation reaching the bottom of the cloud is reduced by intense absorption and scattering within the cloud field, and the net shortwave radiation is less than 100 W/m 2 . It can also be found that the horizontal scattering of cloud water on solar radiation causes the cloud to have an influence on the noncloudy area and there is a horizontal gradient between the surrounding air and the cloud region around the cloud area. Radiation solvers change ambient temperatures through atmospheric heating rates, which in turn affect cloud microphysical and dynamic processes [38,40,42]. The influence of the atmospheric heating rate on clouds is a complex process. The rising of temperature in clouds increases buoyancy and accelerates the development of clouds. Figure 4c,d presents the shortwave heating rate in the 1D and 3D radiation model. The high value region of the shortwave heating rate is located at the cloud top and indicates that the absorption of shortwave radiation is the strongest. In the 1D radiation case, the maximum shortwave heating rate reaches 1.6 k/h at the cloud top. Cloud fields affect the atmospheric heating rate in the vertical direction but not in the horizontal direction. The heating rate above the cloud top is higher than that of the 3D radiation scheme, while the heating rate on the side of the cloud is lower than that of the 3D radiation scheme. In the 3D radiation scheme, the heating rate of ambient air produces a disturbance due to the horizontal transmission beside the cloud fields.

Longwave Radiation Simulation Results
Many studies have shown that the effects of 3D thermal radiation have a large impact on the development of cumulus clouds, and the cooling rate at the cloud top and slight heating rate at the cloud base change the cloud circulation by creating an unstable motion [23,26,35]. Figure 5a,b presents the downward longwave flux in the 1D and 3D radiation case, respectively. In the 1D radiation case, the downward longwave flux decreases with height from 425 W/m 2 to 50 W/m 2 in the cloud-free region, which is mainly from the scattering of the upward longwave radiation emitted by the ground and the downward longwave flux emitted by the atmospheric environment gases. In the cloudy area, the downward longwave flux is larger than that of the environmental atmosphere horizontally since the cloud water particles are the main source of downward radiation contributing to the increased value. The presence of clouds, however, has no effect on the surrounding air horizontally and there is a clear distinction between the cloudy area and the area without clouds given that the downward longwave flux in the cumulus cloud region is independent of the cloud-free region. The longwave radiation effects of the clouds on the environment or adjacent clouds are ignored in the 1D radiation model. In the 3D radiation case, it is not constant in horizontal direction in the cloud-free area, the downward longwave flux close to the cloud area has a large horizontal gradient, and the horizontal scattering from the cloud field makes the cloud and free regions linked to each other and excessive, allowing for the interaction between two clouds. Figure 5c,d presents the upward longwave radiation flux in the 1D radiation model and 3D radiation model. Unlike shortwave radiation, atmospheric air also strongly attenuates longwave radiation. Cloud droplet particles and surrounding air are important media for longwave radiation scattering and absorption. The upward longwave radiation flux is mainly contributed by the radiation flux emitted by the surface. In the 1D radiation case, the following is observed. In the cloud-free region, the upward longwave flux decreases linearly from 440 W/m 2 to 240 W/m 2 from the ground to the model top. In the cloudy area, the upward longwave flux is lower and is due to the attenuation of longwave radiation by cloud droplets, due to the strong absorption of the cloud field, the upward longwave radiation above the cloud top is less than 240 W/m 2 and extends to the top of the model. In the 3D radiation case, the following is observed. The downward longwave flux decreases with the height from the surface ground to the model top. In the cloudy area, the upward longwave flux is lower than the environment atmosphere caused by cloud water absorption, the upward longwave radiation is less than 240 W/m 2 at the cloud top, and in the cloud-free area, the upward longwave flux has a large horizontal gradient near the cloud area. In the horizontal direction, the upward longwave radiation from the grid point near the cloud is less than that away from the cloud, which may be caused by single scattering and multiple scattering shadows in a variety of directions. The thermal source containing the horizontal component cannot reach the cloud sides due to the shadow effect. There are two major differences between 1D radiation and 3D radiation, as shown in 20 min. In the cloud-free area, there is no horizontal gradient in the upward longwave flux in the 1D radiation model, and cloud particles only affect the vertical grid points. In the 3D radiation model, there exists a large horizontal gradient near the cloud area, and the longwave radiation decreases when it is near the cloud area. This may be because the shadow effect of the cloud prevents some longwave radiation from reaching the boundary of the cloud. In cloudy areas, the minimum value of longwave radiation in the 1D radiation model is less than 240 W/m 2 and extends to the top of the model; however, the minimum value of longwave radiation is located in the cloud top in the 3D radiation model. This occurs because in the 1D radiation model, there is only vertical transmission, and given the absorption and scattering of particles inside the cloud, the longwave radiation above the cloud top is relatively small, and the longwave radiation cannot omit the cloud properties in the vertical direction. In the 3D radiation model, given the contribution of the horizontal transmission, the radiance that contains the horizontal component can reach every corner of the model through multiple scattering. The longwave flux can be transmitted from other grid points and avert the cloud top to the model top, given that the value is smallest at the cloud top because of the strongest absorption by cloud water. The thermal radiation drives the generation of clouds and affects the cloud formation, therefore ignoring the influence of clouds on the ambient air in the 1D radiation model will inevitably affect the accuracy of the simulation results. Figure 6a,b presents the net longwave radiation flux in the 1D radiation model and 3D radiation model. The upward direction is positive. The net longwave radiation flux is less than 20 W/m 2 inside the cloud field in both cases, so our focus is on the cloud-free region. In the 1D radiation case, the net longwave flux clearly presents two distinct parts; in the cloud-free region in the vertical direction, the net longwave flux of each layer is constant horizontally and there exists a large gradient between the cloud field and cloud-free field. Above the cloud area, the net longwave flux is small, which is mainly due to absorption and scattering of upward longwave radiation by cloud particles. Since there is only vertical transmission, the upward longwave radiation cannot reach the model top. As a result, the net radiation flux is relatively small at the model top. In the 3D radiation case, the cloud-free region shows continuous changes related to the cloud region and there is a large horizontal gradient in cloudless regions, especially near the cloud top and cloud edge boundaries. In the horizontal direction, the net radiation flux is small near the cloud area and occurs because the cloud region suppresses some of the longwave radiation transmitted upwards. The presence of clouds not only affects the change of radiation flux inside the cloud but also affects the change of ambient air radiation flux through horizontal scattering. Figure 6c,d presents the longwave heating rate in the 1D radiation model and 3D radiation model. The longwave radiation heating rate causes strong local cooling at the cloud top accompanied by a modest warming at the cloud base in two cases, which leads to a destabilization of the cloud layer. The constant longwave heating rate will modify the development of the cloud process. The decrease of the cloud top temperature and the increase of the cloud bottom temperature accelerate the convection. The cooling effect on the cloud top and the heating effect at the cloud bottom are both more intense in the 3D radiation model.

The Radiation Simulation Results by the True Zenith Angle
The Effects of the True Zenith Angle Incidence on Shortwave Radiation Figure 7 shows the cloud water mixing ratio and cloud number concentration in 1D and 3D radiation cases. The energy balance of shortwave radiation is often affected by illumination and shadowing, and the sun position is an important aspect. In the 1D radiation model, the direction of solar radiation is vertically downward regardless of the zenith angle and azimuth angle, and the shadow of the cloud is by definition always directly underneath it. It can be seen from Figure 8a that even though the real zenith angle is adopted, the solar incidence direction is the same as that of the vertical incidence case, only with a different value. This treatment seems to be numerically reasonable. The influence of the solar incident direction on clouds is ignored. It was found that the solar incidence angle provides the possibility of organizing the cloud through solar radiation heating. Furthermore, the effect of the solar zenith angle on solar radiation scattering and absorption will accumulate gradually with the increase of simulation time. In contrast, in Figure 8b the downward shortwave radiation is sloping by the incident solar radiation in the 3D radiation case. The flux facing the solar incident direction is greater than that deviating from the incident direction at the cloud top. The direction of solar incidence not only affects the energy budget inside the cloud but also influences the distribution of ground heat flux through the cloud shadow. The cloud shadow is tilted according to the solar incidence angle. Moreover, 3D radiation transmission allows the energy to move horizontally and to correctly transmit to the cloud's shadow based on the position of the sun.  Figure 8c,d presents the upward shortwave radiation in the true zenith angle incidence case. It can be seen that the solar incident direction has little influence on the reflected solar radiation in 1D and 3D radiation cases, and the structure of upward shortwave radiation is highly symmetrical.
The influence of radiation feedback on clouds is eventually reflected in the temperature change in the ground and ambient air as well as in the occurrence and development of clouds. Figure 9 presents the net shortwave radiation and shortwave heating rate in the true zenith angle incidence case. In 3D radiation model, the following is observed. Above the cloud top, the structure is similar to the vertical incidence case, while below the cloud base, the net shortwave radiation is inclined by the influence of the incident direction of the sun, and the position of the cloud shadow is obviously perpendicular to the sun's incidence. However, the shadows always fall directly below the clouds by definition in the 1D radiation model. The changing rate of net radiation flux will affect the atmospheric heating rate, and thus affect the development of the cloud. In the 1D radiation model, the shortwave heating rate is similar to the case of vertical incidence, with a maximum value of 1.3 k/h. In the 3D radiation case, the heating rates on both sides of the cloud are asymmetric, and the shortwave heating rate facing the incident radiation direction is greater than that deviating from the incident direction at the cloud top. The results indicate that solar radiation heating is likely to increase the heating rate of the part perpendicular to the solar incidence angle. Moreover, 3D radiation also causes the horizontal variation of the heating rate in cloud-free regions, and the heating rate in the cloud shadow is lower than that of nonshaded regions. The horizontal displacement of the cloud shadow affects the structural distribution of the atmospheric heating rate. Figure 10 shows the heating rate in 1D and 3D radiation models in which the solar radiation that is very intense at 12:00 is superimposed with the cooling rate by longwave radiation and indicates there is cooling at the cloud top and heating below the cloud top. The cloud top cooling rate is stronger in the 3D radiation scheme, especially in terms of the deviation from the solar incident direction. Furthermore, the difference between the two radiation cases is asymmetrically distributed after subtraction. This may be due to the uneven distribution of energy in the 3D radiation model of the real zenith angle. This complex relationship between the asymmetric distribution structure of the atmospheric heating rate and radiation propagation may facilitate or prevent the development of clouds.

Discussion
The 3D radiation effect of deep convective cumulus clouds in ARW-WRF model is discussed above; however, there are also other atmospheric models to choose from such as LES (large eddy simulation) model, GCM (general circulation model), cloud tomographic reconstruction, etc. It is important to consider the 3D radiation effects by parameterizing the 3D radiation model with them. The applications for SHDOM is portable and tolerable. It can be used as an independent radiation parameterization scheme coupled with other models and compared with other radiation schemes. The cloud model can provide input parameters for the radiation scheme, and the radiation scheme will give radiation results. In this case, the radiation scheme can be coupled online or offline with cloud model as needed. In this process, it is important to select the input parameters, output parameters and the radiation parameters should be set according to the actual case. One of the difficulties is determining the desired trade-off between accuracy and computation cost in choosing parameters. For input parameter, the cloud model should provide the information including temperature, pressure, atmospheric composition, and cloud property information at each grid point. The cloud property should contain the content of cloud water or cloud ice, effective radius, and cloud number concentration. The output format includes radiance, heating rates, hemispheric flux, source function, and medium properties. These variables can be output as needed. The radiation parameters include the setting of boundary conditions, solution accuracy, angular resolution, and spatial resolution. On this basis, the radiation model can be coupled with other cloud models for further research.

Summary and Conclusions
This work aims to study the 3D radiation effect of deep convective cumulus clouds. To that end, using a 3D radiation transmission model, the spherical harmonic discrete ordinate method (SHDOM) is applied online with radiation feedback to ARW-WRF. The overall description of the WRF model and the technical challenges encountered in migrating the 3D radiative solver to WRF are discussed. To conduct the verification tests and explore the 3D radiation effect of deep convective cumulus clouds, we initiated a large eddy simulation with a thermal bubble perturbation under high humidity. The differences between 1D and 3D radiation schemes at different solar zenith angles are compared. The most important findings can be summarized as follows: (1) In the 1D radiation model, grid points are independent of each other in the calculation, the absorption and scattering only occur in the vertical direction, and neither the optical variability nor the horizontal cloud particles are taken into account. In this case, the interactions between two cloudy regions or between clouds and ambient air are not considered in the simulation. As a result, there is a clear distinction between the cloudy grid point and the cloudless grid point.
In the vertical direction, the radiation cannot pass through the cloudy area due to the shadow of the cloud, which results in spurious high-value or low-value regions that exist at the model top. This deviation can accumulate over time and affect cloud microphysics and dynamic processes. (2) In the 3D radiation model, the solver computes the radiation transmission, including upward and downward streams, and additionally for sideward streams, the scattering intensity of cloud particles depends on the phase function of the direction, such that photons scattered in multiple directions by cloud particles can travel through the air until encountering new cloud particles. The scattered radiation structure is reasonable, and the distribution of scattered radiation corresponds to the cloud field. In addition, the horizontal transmission and multiple scattering can also cause radiation disturbance in noncloudy regions. As a result, the radiation presents a clear structure around the cloud field, and the cloud area affects the cloud-free area through horizontal transmission, which makes the cloud and cloud-free regions linked to each other and excessive, allowing for the interaction between two clouds. (3) In the true zenith angle case, the most obvious difference between 1D and 3D radiation transfer is the distribution of radiation to the cloud top and the displacement of the sun's shadow on the surface ground. In 1D radiation transmission, the distribution of radiant energy is completely symmetric regardless of the zenith angle, and the shadow of the cloud is directly underneath the cloud. In 3D radiation, the solar flux facing the incident radiation direction is greater than that deviating from the incident direction at the cloud top, and the cloud shadow is based on the position of the sun. (4) The influence of radiation feedback on the cloud is eventually reflected in the atmospheric heating rate and ultimately affects the development of the clouds. For shortwave radiation, in the 1D radiation case, the cloud fields affect the atmospheric heating rate in the vertical direction but not in the horizontal direction. The heating rate above the cloud top is higher than that of the 3D radiation scheme, while the heating rate on the side of the cloud is lower than that of the 3D radiation scheme. In the 3D radiation case, the solar radiation heating is likely to increase perpendicular to the solar incidence direction in the cloud field, and the heating rate of the ambient air produces a disturbance due to horizontal transmission beside the cloud fields. The horizontal displacement of the cloud shadow also affects the surface heat flux through the land surface model. This complex relationship between the asymmetric distribution structure of the atmospheric heating rate and radiation propagation may promote or hinder cloud formation and development. For longwave radiation, the cooling effect at the cloud top and the heating effect at the cloud bottom are both more intense in the 3D radiation model.
Through the description of the model details and the analysis of the experimental results, the coupling of the model attains a certain reliability, and the distribution of the radiation field is significantly different between the 1D radiation scheme and 3D radiation scheme in the large eddy simulation. The effects of radiation on clouds are multifaceted and complex and there are also some inadequacies in this study. Thirty minutes is not long enough for the difference between the radiation fields to have a significant impact on the clouds. Furthermore, the ideal cases in WRF are also inadequate to represent real atmospheric radiation conditions. The influence of the 3D radiation feedback on the cloud field needs to be further studied. Subsequent work will allow for a long simulation time with the WRF real case in the presence of adjacent clouds to better compare the effects of 3D radiation patterns on the cloud's macrostructure and evolution. The coupled model is fully open to all researchers, and email requests will be accommodated.

Conflicts of Interest:
The authors declare no conflict of interest.