Determination of global Earth outgoing radiation at high temporal resolution using a theoretical constellation of satellites

New, viable, and sustainable observation strategies from a constellation of satellites have attracted great attention across many scientific communities. Yet the potential for monitoring global Earth outgoing radiation using such a strategy has not been explored. To evaluate the potential of such a constellation concept and to investigate the configuration requirement for measuring radiation at a time resolution sufficient to resolve the diurnal cycle for weather and climate studies, we have developed a new recovery method and conducted a series of simulation experiments. Using idealized wide field‐of‐view broadband radiometers as an example, we find that a baseline constellation of 36 satellites can monitor global Earth outgoing radiation reliably to a spatial resolution of 1000 km at an hourly time scale. The error in recovered daily global mean irradiance is 0.16 W m−2 and −0.13 W m−2, and the estimated uncertainty in recovered hourly global mean irradiance from this day is 0.45 W m−2 and 0.15 W m−2, in the shortwave and longwave spectral regions, respectively. Sensitivity tests show that addressing instrument‐related issues that lead to systematic measurement error remains of central importance to achieving similar accuracies in reality. The presented error statistics therefore likely represent the lower bounds of what could currently be achieved with the constellation approach, but this study demonstrates the promise of an unprecedented sampling capability for better observing the Earth's radiation budget.


Introduction
The Earth reflects part of the incoming solar radiation back to space and responds to the remaining absorbed energy by adjusting its temperature and emitting thermal infrared radiation out to space. Observing these energy flows exiting the top-of-the-atmosphere (TOA), referred to as Earth outgoing radiation (EOR), has advanced our understanding of fundamental climate parameters such as the planetary brightness [Vonder Haar and Suomi, 1971], the greenhouse effect [Dickinson, 1985], and the zonal and meridional heat transports required by the atmosphere and ocean to redistribute regional energy imbalances [Rasool and Prabhakara, 1966;Charney, 1975]. EOR observations also play a crucial role in studying climate forcing and feedbacks [Futyan et al., 2005;Loeb et al., 2007;Brindley and Russell, 2009;Ansell et al., 2014], global energy imbalance and its implication in the hydrological cycle [Trenberth and Fasullo, 2011;Stephens et al., 2012;Allan et al., 2014;Hegerl et al., 2015], and in the development, improvement, and validation of weather and climate models [Forster and Gregory, 2006]. [Smith et al., 1977;Barkstrom, 1984;Kyle et al., 1993], the current Clouds and the Earth's Radiant Energy System (CERES) [Wielicki et al., 1996], and the Geostationary Earth Radiation Budget (GERB) experiment . CERES has proven invaluable for studying global energy balance; the uncertainty on the monthly net TOA irradiance determined from CERES ranges from À2.1 to 6.7 W m À2 [Loeb et al., 2009]. Improved absolute accuracy is attainted by incorporating ocean heat content observations [Lyman et al., 2010;Church et al., 2011]. Unlike CERES, operated mainly in Sun-synchronous orbits, the geostationary capability of GERB provides outgoing radiation products at a 15 min temporal resolution, but over the European and African continents and surrounding areas only. interpolated data set [Doelling et al., 2013[Doelling et al., , 2016] that has been used in studies of cloud and aerosol radiative forcing [Taylor, 2012;Su et al., 2013], explaining TOA diurnal cycle variability [Taylor, 2014;Dodson and Taylor, 2016], and to make comparisons with models and reanalyses [Itterly and Taylor, 2014;Hill et al., 2016]. However, incorporating different geostationary observations with unique sensor characteristics and varying degrees of quality can produce significant artifacts resulting in unnatural spatial patterns in radiation fields [Doelling et al., 2013]. It is clear that current missions are not designed for these high temporal resolution applications, despite the apparent need for global EOR at high temporal resolution, preferably sufficient for resolving the diurnal cycle.

Measurements of EOR have been made from dedicated missions since 1975; examples include the early Earth Radiation Budget (ERB) experiments
Such EOR measurements require a new observation strategy. Recently, because of a technology revolution in small satellites and sensor miniaturization [Sandau et al., 2010], a constellation approach has drawn considerable attention and has been applied in monitoring surface wind speeds over oceans and within tropical cyclones [Komjathy et al., 2000;Ruf et al., 2012], surface albedo [Nag et al., 2015], and providing surface imagery for natural disaster monitoring [Underwood et al., 2005]. Its potential for providing high temporal resolution EOR observations, however, has not been fully explored.
This paper aims to evaluate whether the new constellation concept can potentially provide an EOR observational data set with high temporal resolution and accuracy for weather and climate studies and to identify the key factors affecting the performance of the constellation. While various advanced sensors like scanning narrow field-of-view (FOV) radiometers can be potential candidates for a constellation, we focus on applications using wide FOV (WFOV) broadband radiometers. This type of sensor is low cost and, more importantly, has features preferable for any Earth radiation budget mission, such as mature instrumentation technology and no moving parts.
WFOV measurements have proven invaluable for monitoring planetary EOR [Smith et al., 1977;Barkstrom, 1984;Kyle et al., 1993], but their large footprint size makes it less straightforward to obtain EOR at a synoptic scale. To enhance the spatial resolution of WFOV observations, Hucek et al. [1987Hucek et al. [ , 1990 applied spherical harmonic analysis to 6 day measurements of a single WFOV radiometer and showed significant improvement not only in spatial resolution but also in the root-mean-square error (RMSE) of the global mean albedo. Their recovery analysis greatly enhanced the EOR product from monthly mean to weekly mean; however, the enhancement reached a limit, because a time interval of 6 days was necessary to accumulate sufficient instantaneous measurements from a single radiometer. Clearly, an increase of available satellites using a constellation can help push this time limit further and enhance the temporal resolution, but it is unclear what the required configuration of the constellation will be if we aim to observe EOR hourly to monitor fastevolving systems such as dust storms and cyclones in the tropics and extratropics. An enhanced ability to track these systems allows us to capture their evolution and understand the underlying radiative processes, further constraining future warming of our climate system. We will perform extensive simulation experiments to explore the capability of the constellation concept for observing EOR. Similarly to Hucek et al. [1987Hucek et al. [ , 1990, Han et al. [2008], and Salby [1988a, 1988b, we use spherical harmonic analysis for detailed EOR recovery from WFOV measurements, but necessarily, different constraints are applied to make the recovery work for a much higher temporal resolution. In this paper, section 2 details the recovery method and the design of simulation experiments, while section 3 presents the results from a baseline constellation and related sensitivity tests. Finally, section 4 summarizes our key findings and highlights the upcoming opportunities that this work can be directly applied to.

Simulation of Measurements
To develop and evaluate our recovery method, synthetic data sets were generated using output from the Met Office global numerical weather prediction (NWP) model [Walters et al., 2014]. The model was run from an operational 00 Z analysis using 1024 longitudinal and 769 latitudinal grid points, with a coarsest spatial resolution of approximately 39 km × 26 km along the equator. Normally, the time step used for global NWP with this configuration is 12 min and for reasons of computational expense, full radiation calculations are only done every hour [Manners et al., 2009]. In this simulation, however, the time step was reduced to 5 min and to represent evolution of the radiation fields better, the radiation scheme [Edwards and Slingo, Journal of Geophysical Research: Atmospheres 10.1002/2016JD025514 1996 was called on every one of those time steps. Arbitrarily, we chose a 1 day model output from 12 Z on 28 August 2010 (T + 12 to T + 36) to simulate satellite WFOV measurements for our experiment. Figure 1 shows an example of the modeled outgoing shortwave (SW) and longwave (LW) irradiance fields at the TOA of 80 km, for 00 Z on 29 August 2010. These irradiance fields contain clear structures associated with meteorological phenomena including large regions of tropical convection in the central and western Pacific and a midlatitude frontal system in the northern Pacific. Meteorological phenomena of this type are large in spatial scale (~1000 km), evolve quickly (approximately hourly) and appear to dominate the variability of the outgoing irradiance fields, and will therefore be of interest to observe when recovering outgoing irradiance fields.
Using these modeled irradiance fields, the instantaneous irradiance measured by each satellite, F sat , at colatitude θ and longitude λ, is then the integration of radiation received from all points within the WFOV, given as where I TOA (θ ', λ ') is the radiance reaching the satellite from TOA location (θ ', λ ') within the satellite FOV Ω, A(θ, λ, θ 0 , λ 0 ) is a factor to account for the laboratory-measured angular response of the instrument, β is the angle between the satellite nadir and the line of satellite to location (θ ', λ '), and e represents systematic and random instrument error, which is investigated in section 3.2.2. For the sake of simplicity, we assume that the instrument has a flat spectral response that is insensitive to the spectral composition of the observed scene and a flat angular response to omit A(θ, λ, θ 0 , λ 0 ) from hereon in. Since the model output provides irradiance rather than radiance, I TOA were generated using an isotropic assumption. In other words, equation (1) can now be rewritten as where F TOA is the instantaneous outgoing irradiance at the TOA. This assumption inevitably introduces uncertainty in the recovery, which will be evaluated in section 3.2.1 using angular distribution models from the CERES [Loeb et al., 2003].
Using equation (2), synthetic measurements are generated at the sampling rate of the satellite instruments. At each time step, a satellite will have progressed through space, giving a different value of θ and λ and thus a different FOV and F sat . This forms nonuniformly distributed yet dense satellite measurements, which rely on spherical harmonic analysis to recover the global distribution of the outgoing irradiance field, as explained next.

Recovery Method
Spherical harmonic analysis serves as a homogenization of multiplatform observations to provide the complete distribution of the irradiance field on the entire Earth. This distribution can then be used to interpolate the irradiance field to smaller scales by optimally exploiting the dense measurements. At a given location (θ, λ), F TOA can be represented exactly by a spherical harmonic series as or be approximated by a truncated series as where L is the truncation limit, leading to (L + 1) 2 terms on the right-hand side. The spherical harmonic functions Y C lm and Y S lm represent the spatial distribution, defined as where P lm is the associated Legendre function. The coefficients S l0 are always zero, while the first term containing the harmonic coefficient C 00 finds the global mean outgoing irradiance. In contrast, the higher degree harmonics, C lm and S lm , represent the amplitude of structures at finer resolutions. Denoting R E as the Earth radius, the subscript l represents structures at a wavelength of 2πR E /l in spherical harmonic analysis and thus a resolution (half of the wavelength) of approximately 20,000/l km. In other words, the spatial resolution of the recovered irradiance field is determined by the truncation limit L; a value of 20 leads to a resolution of 1000 km, recovered from 441 (i.e., (L + 1) 2 ) spherical harmonic coefficients.
Substituting equation (4) into equation (2) gives where the truncation inequality has been dropped; Y C lm and Y S lm represent the spatially integrated spherical harmonic functions within the WFOV and respectively relate to Y C lm and Y S lm in equation (5) by For M satellite measurements, equation (6) becomes a system of equations and can be written in matrix form as where F sat is a M × 1 matrix comprising all satellite measurements, Y is a M × N matrix containing N number of the spatially integrated functions (i.e., Y C lm and Y S lm ) for each simulated satellite measurement, c is a N × 1 matrix comprising the unknown coefficients C lm and S lm , and e is a M × 1 matrix containing the measurement errors. Using a constrained least squares approach, an estimate of c in equation (8), ĉ, is given by where ε 2 is empirically chosen to be 10 À4 and I is an N × N identity matrix. The ε 2 factor was chosen to be large enough to improve the condition number of the normal matrix and therefore stabilize the solution, while being 6 orders of magnitude smaller than the mean of the simulated satellite measurements and therefore not introduce significant bias in the recovery. Once equation (9) is solved, we place the coefficients, ĉ, back into equation (4) to recover the TOA outgoing SW or LW irradiance F TOA (θ, λ) for all global locations.

10.1002/2016JD025514
Assuming that the errors in F sat are uncorrelated and have constant variance σ 2 , the variance-covariance matrix of ĉ can be given as [Hastie et al., 2009, p. 47] var where Then, the variance in the recovered outgoing irradiance field at each location (θ, λ) can be estimated by where Y(θ, λ) ia a 1 × N matrix comprising the spherical harmonic functions Y C lm and Y S lm in equation (5), and the square root of it (i.e., the standard deviation representing about 68% of the normal distribution) is used as the recovery uncertainty estimate.
Note that the spherical harmonic analysis above uses one set of matrices (i.e., equation (8)) to recover one irradiance field (i.e., equation (4)), which implicitly assumes that the irradiance field does not change during the analysis. This assumption is a particular concern for the SW spectral region, because SW irradiance fields evolve fast for two reasons. One is that the illuminated portion of the Earth changes rapidly due to Earth's rotation and solar geometry; the other is that albedo varies with changes in the radiative properties of the atmosphere and surface. For a WFOV, the former is found to be the dominant factor for SW radiation evolution and needs to be incorporated into the recovery method. Therefore, instead of instantaneous satellite measurements, we used time-averaged values in equation (9) through the following process. For each WFOV satellite measurement, the instantaneous albedo, A INS , is calculated as the ratio of this measurement and the instantaneous incoming solar irradiance over the WFOV. The average incoming solar irradiance within the measurement collection period (e.g., 1 h), S, can be also calculated. The time-averaged satellite measurement in the SW, F ' sat;SW , can then be estimated by and the set of measurements generated in this way replaces F sat in equation (9) when performing the recovery. The LW measurements remain unchanged.

Experiment Design and Setup
The configuration of a satellite constellation is defined by several parameters including the inclination angle, altitude, number of satellites, number of orbit planes, phasing of the satellites between planes, instrument sampling rate, and collection time period used to recover the irradiance field. While possible configurations can be unlimited, they unavoidably depend on launch opportunities, available finances, and, importantly, the required accuracy of the recovered irradiance.
To provide a general guideline, we focus on a baseline constellation configuration that consists of 36 satellites spread evenly across six orbital planes ( Figure 2a), with no phase difference between satellites in adjacent planes. Each orbital plane resides in an 86.4°non-Sun-synchronous orbit, allowing every satellite to precess throughout 24 h and provide coverage all the way to the poles, while maintaining plenty of separation to avoid any possibility of collision. The altitude of the satellites is chosen to be 780 km, the upper end of typical small satellite altitudes [Maessen, 2007;Lücking et al., 2011]. This altitude enables a FOV as large as possible, corresponding to a horizon-horizon footprint of approximately 6000 km in diameter for an instrument FOV of 126°. For simplicity, we use circular orbits with all satellites at the same altitude, but the performance for varied satellite altitudes and other combinations of orbits is consistent given that the same sampling density shown in Figure 2b is provided. Note that this baseline configuration typically requires six launches, but it is possible to use fewer launches and let satellites drift and spread out, although it may take several months for satellites to reach the specified planes [Bingham et al., 2008].
Additionally, we use a 5 s sampling rate to be consistent with past WFOV radiometer response times [Hucek et al., 1987]. We also focus on measurements over a 1 h collection period because of the following trade-off. For scientific purposes, it is crucial to enhance our ability to observe the diurnal cycle of the global irradiance field; therefore, the collection time period must be as short as possible. However, in practice, it is necessary to accumulate sufficient measurements to provide dense global coverage that is required for the success of the recovery. To optimize this collection time period, we used the model output to examine how many times each grid point needs to fall within satellite's FOVs (i.e., to be seen by the satellites) to produce a stable recovery. As expected, the most challenging region is at lower latitudes. As shown in Figure 2b, while it is not difficult to sample high-latitude regions densely, the sample size at low latitudes is at least 5 times smaller than that at high latitudes during a 1 h time period. Interestingly, several regions along the equator with small sample sizes somehow coincide with satellite ground tracks (see Figure 2a), which is counterintuitive. The reason is that the "gaps" between satellite tracks have a better chance of being observed from both sides; in other words, the overlapping FOVs from satellites in adjacent orbital planes make the regions between tracks better sampled than those along the tracks. In general, we found that a minimum density of 500 samples is required for stable recovery, which can be achieved for the entire globe by a 1 h collection time period as shown in Figure 2b. A much longer duration is inappropriate because it would result in time aliasing, strongly violating the assumption of a static irradiance field as discussed in section 2.2, and introduce a source of recovery instability.
The performance of the recovery method is measured by three metrics. In addition to commonly used global mean bias and RMSE, we also use the signal-to-noise ratio (SNR) from power spectra to quantify the recovery errors dependent on spatial resolution, investigate how fast the recovery errors increase with spatial resolution, and when the recovery noise becomes dominant and thus the recovery outcome is no longer reliable. For each spatial resolution, the SNR is defined as the ratio of the recovery amplitudes and error amplitudes. The average recovery amplitude associated with a spatial resolution of 20,000/l km, A r,l , is given by where Ĉ lm and Ŝ lm are the recovered coefficients in ĉ lm and the corresponding error amplitude, A e,l , is calculated by We used a value of SNR of unity to identify the spatial scale at which fast-evolving meteorological phenomena can be reliably recovered.
The performance of the baseline configuration is also used to compare and contrast recovery errors introduced by other factors, such as the number of satellites, the configuration of the satellites, instrument performance/calibration, and the isotropic assumption. We conducted a series of sensitivity tests to identify which of these factors play the most crucial roles and to provide guidelines about what needs to be considered when designing an optimal configuration.

Recovered Outgoing Irradiance From the Baseline Constellation
The performance of the recovery method depends on the required spatial resolution. The SNR of the recovered irradiance field, calculated from the power spectra analysis in Figure 3, generally decreases as Journal of Geophysical Research: Atmospheres 10.1002/2016JD025514 the spatial resolution becomes finer. These variations are determined by the spatial structures in the radiation field and are therefore nonmonotonic in nature. For example, the error at the largest planetary scales in the SW is relatively large due to the need to recover the position of the terminator. At a 1000 km spatial resolution, the SNR remains greater than unity in both the SW and the LW, suggesting that features of the atmosphere and surface are still recovered well at this resolution.
Taking 1000 km resolution as an example, Figure 4 shows the truth irradiance field averaged from the finer-resolution model output (i.e., Figure 1) over the 1 h observation period, along with the synthetic WFOV satellite measurements and the recovery from the baseline constellation. As shown in Figures 4a  and 4b, although the truth irradiance field at a coarse resolution loses fine structures of the atmospheric systems, synoptic-scale features such as tropical cyclones in the central Pacific and midlatitude frontal systems that were previously identified as dominant features of the outgoing irradiance field are retained. Due to the WFOV of the satellites, these features are not directly revealed in satellite measurements at the native resolution (Figures 4c and 4d). However, once the spherical harmonic analysis is applied, the recovered irradiance fields (Figures 4e and 4f) over a 1 h collection time period show a remarkable resemblance with the truth fields, which reconfirms the added value of performing the spherical harmonic analysis.
The uncertainty associated with the recovery is shown in Figures 4g and 4h, which has a spatial distribution similar to that of the sample density in Figure 2b. Interestingly, the regions with large recovery uncertainty do not coincide with the gaps of the sample density but fall between the gaps. Recall that the gaps of the sample density coincide with satellite ground tracks. Compared to the samples close to the edge of the FOV, nadir samples along the satellite track have much larger radiative contribution imparted to the WFOV measurement due to the cosine factor as seen in equations (1) and (2). Therefore, even though the corresponding sample density along the satellite tracks is relatively lower than neighboring regions, nadir radiation is retained well in WFOV observations and can be recovered with small uncertainty. In contrast, although regions between ground tracks are seen more frequently in the WFOV than along the tracks, their radiative contribution is small and thus the corresponding irradiance has larger uncertainty. Additionally, the regions with large recovery uncertainty in the Eastern Hemisphere are generally to the south of the equator, while those in the Western Hemisphere are to the north of the equator. This is as a result of sampling. Note that the satellites are going north in the Eastern Hemisphere and going south in the Western Hemisphere. Since all satellites have an orbital period longer than our collection time of 1 h and have not completed a full orbit yet, some regions are therefore not sampled as well as others, leading to larger recovery uncertainty.
We evaluate the recovery performance further via scatterplots and histograms of the recovery errors, calculated by subtracting the recovered outgoing irradiance from the truth on a 1000 km grid. Consistent with Figure 4, the scatterplots in Figure 5 show that the truth and the recovered outgoing irradiance are correlated well in the both SW and LW regions, and the majority of the data points fall onto the 1:1 line. Compared to the SW, the LW recovery errors in Figure 5d are smaller (within 10%), because the LW irradiance field is smoother and evolves more slowly, leading to less truncation error and time aliasing. As a result, the absolute errors in the LW recovery (Figure 5f) are distributed more evenly around the globe but increase generally toward the equator as the magnitude of the outgoing LW irradiance increases. In contrast, while 94% of recovery in the SW agrees with the truth within 25% (Figure 5c), Figure 5e shows that large isolated errors can exist, particularly on the daylight side of the terminator and in the regions where tropical convection is common.
As shown in Figures 4 and 5, the recovery errors are significantly larger than the estimated uncertainty, which calls into question the sources of the errors. To investigate the issues involved, we first looked into two convection systems that are associated with similar SW outgoing irradiance of over 600 W m À2 but have very different recovery errors (one less than 25 W m À2 and the other~120 W m À2 ). The case with small errors in Figure 6a is mainly composed of an isolated convective system at the center of the domain. There are many surrounding cloud bands, but they are small and far from the main convection system, leading to a good recovery in both location and magnitude. Unlike the case in Figure 6a, the case with large errors in Figure 6b is composed of many complex cloud systems that surround each other closely, making it difficult to capture the fine gaps between systems and to recover the exact boundaries of convections. Consequently, although the recovery is able to capture the main features of the truth, the magnitude and location of the main system centered at the domain are off, leading to large recovery errors surrounding the main convection system. These case studies indicate that the recovery includes two types of errors: one is the omission error that is owing to a limited spherical harmonic degree (i.e., a spatial resolution issue) and the other is the commission error mainly due to producing a static field from measurements collected over some finite time interval (i.e., a temporal resolution issue). To diagnose their relative impact on the recovery, the following experiments were performed. We started with a low resolution of 1000 km, static irradiance field (averaged over 1 h as shown in Figures 4a and 4b) as input to simulate synthetic satellite measurements. As expected, this leads to an excellent recovery with negligible errors of less than 0.01 W m À2 for all grid points. We then put spatial complexity back into the input irradiance field by replacing the static, low-resolution field with a field averaged over 1 h at the original high resolution. The corresponding recovery (Figure 7a) reveals similar errors to those found in Figure 5e in the western and central Pacific region, suggesting that the spatial resolution is the primary issue for the large errors found in the daylight side of the terminator (except the edges). For the temporal resolution issue, we estimate the largest possible errors by taking the difference between 00 Z and 01 Z irradiance fields. As shown in Figure 7c, the rapid change of AE100 W m À2 in SW occurs not only on the edges of the terminator but also in the regions approximately 45°east and west extended from the edges, making the recovery prone to large commission errors. Note that we take an extreme estimate in Figures 7c and 7d and that the actual errors are much less, as shown in Figures 5e and 5f. Based on these analyses, adding low-inclination satellites to improve samples in the tropics and shortening the collection time period will be effective ways to reduce both the omission and commission errors.
The recovery errors over all grid points in Figures 4e and 4f are 0.5 AE 16.4 W m À2 and À0.02 AE 5.94 W m À2 (mean with one standard deviation) in the SW and LW, respectively. To examine whether the reported bias and errors for this particular time period are representative, Figure 8 shows a 24 h time series of the hourly recovered SW and LW global mean biases. Recall that the previous results and discussions are based on the 1 h observations from 00 Z to 01 Z, around the middle of the time series in Figure 8. Overall, on this day, the hourly global biases in the SW have a mean 0.16 W m À2 with a standard deviation of 0.45 W m À2 .   Stephens et al. [2012]), demonstrating that a constellation can potentially greatly enhance our capability to observe global EOR with much smaller uncertainty even at the hourly and daily time scale. However, the results from this experiment, which we will refer to as the control experiment, are clearly based on theoretical simulations containing some simplifying factors when compared with the real observational problem. To identify the factors that are most important in limiting the performance of the recovery method, including the assumption of isotropic radiation (section 3.2.1), instrument performance (section 3.2.2), and the number of satellites (section 3.2.3), we next present results from a series of sensitivity tests.

Sensitivity Tests 3.2.1. Assumption of Isotropic Radiation
As described in section 2.1, an isotropic assumption was made to convert modeled irradiance output to radiance during the generation of the WFOV satellite measurements. To evaluate the magnitude of the recovery errors introduced by the isotropic assumption, we capitalize on the Angular Distribution Models (ADMs) outlined by Loeb et al. [2003] to account for anisotropy in our synthetic satellite measurements. These ADMs consist of a set of anisotropic factors that empirically relate the isotropic radiance to that observed under certain conditions. In the SW, the anisotropic factors depend on the solar-viewing geometry, the surface type, and the meteorology of the scene (i.e., near-surface wind speed, cloud phase, cloud optical depth, and cloud fraction) which we also obtain from the Met Office NWP model output. In the LW, anisotropic factors span a smaller range and are typically much closer to 1 [Gupta et al., 1985], so we focus our attention on the SW here where the most serious consequences of making an isotropic assumption are likely to exist.
To include a more realistic radiance distribution, we simulate the 1 h set of measurements from the baseline constellation again but this time apply an appropriate anisotropic factor to each and every grid point within each WFOV measurement (i.e., insert this factor within the integral in equation (2)). We found that producing this new set of simulated measurements introduces a systematic bias to the measurements of 1.61 W m À2 which will of course hamper the recovery accuracy. To determine the origin of this bias, we simulated another set of measurements by randomly including the meteorology in the anisotropic factors to isolate the influence of meteorology and the viewing geometry. Those two sets of measurements are only offset by 0.10 W m À2 (Figure 9a), suggesting the bias is largely related to the viewing geometry. We call the set of anisotropic factors that fully account for the scene meteorology "Full-ADM," and the set that randomly include the meteorology as "Random-ADM." To alleviate the anisotropic issue, we perform the recovery on the simulated measurements using the Full-ADM (i.e., the most realistic estimate of what the measurements would be) but use the Random-ADM in the recovery process. This is a reasonable step to take since at any instant in time we will always know the solar-viewing geometry and surface type, even if simultaneous observations of the scene meteorology are not available. As shown in Figures 9b and 9c, the SW anisotropic factors used in this experiment can change the radiance by more than a factor of 10 in extreme cases. However, the resulting bias in the hourly recovered global mean outgoing SW irradiance is just 0.56 W m À2 larger than the control experiment and is 0.93 W m À2 smaller than the experiment that uses Full-ADM measurements but isotropic radiation in the recovery process. Interestingly, this suggests that the instrument's WFOV of 6000 km in diameter encompasses such Journal of Geophysical Research: Atmospheres 10.1002/2016JD025514 a large area that contains a range of meteorological regimes and makes the meteorological dependence of the ADM less critical for obtaining the hourly global mean. If one wishes to consider regional information, however, the errors increase significantly using the Random-ADM in the recovery process. In this case we recommend performing the recovery assuming that the radiation is isotropic, which only increases the RMSE by 6% compared to the control experiment. In practice, the regional error could be reduced by using climatological distributions of the meteorological variables or better still by incorporating the simultaneous meteorology using geostationary imagery or observations similar to the Moderate Resolution Imaging Spectroradiometer.
Note that Su et al. [2015] have recently provided a new generation of ADMs. The new ADMs incorporate additional anisotropic effects, such as those of aerosols and sastrugi. Although the change in regional monthly mean instantaneous irradiances over some regions can be large (up to 5 W m À2 ), the global mean irradiances do not change dramatically (less than 0.5 W m À2 ). Again, since our 5 s instantaneous measurements are generated from a WFOV that samples a large area and a wide range of viewing angles, regional differences would have a limited effect on our results and, therefore, the difference between the old and new ADMs is unlikely to affect our estimate of the anisotropic impact significantly.

Instrument Performance and Calibration
Because there is no existing constellation for us to properly estimate potential instrument performance, and because development and demonstration of miniaturized broadband radiometers are ongoing research activities [Swartz et al., 2015[Swartz et al., , 2016, we test a number of possible scenarios to investigate how fast recovery errors grow with varying instrument performance characteristics to provide a level of uncertainty that engineering development should aim for. These scenarios include adding random instrument noise and systematic calibration errors to synthetic satellite measurements for the case used in Figure 4. If these additional uncertainties are applied based on certain probability density functions, we repeat each test 10 times to ensure that the mean result is robust. A summary of these sensitivity results is listed in Table 1, which applies to both the SW and LW.
First, a random instrument noise is included by adding 0.1 W m À2 of white noise to the simulated measurements, more than an order of magnitude larger than that for similar early instruments [Kyle et al., 1985]. This could account for unknown variations in the angular and spectral response across the WFOV, wobbles in instrument pointing, blurring of scenes due to insufficient instrument response time, or other imperfections in the instrument or satellite performance. Since we applied white noise, the resulting difference in average global mean irradiance from the control experiment is essentially zero with a standard deviation of approximately 0.002 W m À2 (i.e.,~0.002% in SW and less than 0.001% in LW). The reason that this difference is almost negligible is that there are a huge number of individual simulated measurements that are fed to the recovery inversion such that any random noise quickly averages to zero, one of the key advantages of a large constellation approach. When this white noise is increased to 0.25 W m À2 and 0.5 W m À2 , the average global mean difference remains essentially zero, but the associated standard deviation increases to approximately 0.01 W m À2 . This change is still an order of magnitude smaller than the absolute bias derived from the recovery performance (0.16 W m À2 in SW and 0.13 W m À2 in LW; see Figure 8) and therefore is tolerable.
Second, we consider the influence of a systematic calibration bias on the recovered global mean outgoing irradiance. A flat systematic bias of 0.5 W m À2 across the constellation is found to result in a change in the recovered global mean of 0.63 W m À2 . Including 0.1 W m À2 of instrument noise on top of this systematic bias results in the same mean with a standard deviation of 0.003 W m À2 , consistent with that from instrument noise alone. Finally, instead of a systematic calibration bias, we included a standard deviation of 0.1 W m À2 between individual satellites, with a mean bias of 0.5 W m À2 and 0.1 W m À2 of instrument noise; this results in an average change of 0.657 W m À2 in the recovered global mean outgoing irradiance with a standard deviation of 0.043 W m À2 . Overall, the error growths due to compounding sources of uncertainty generally follow a linear behavior. Table 1 suggest that maintaining any systematic bias to a minimum would be crucial for a satellite constellation to measure the EOR accurately, consistent with most climate monitoring missions. To achieve it, in addition to prelaunch and on-flight calibrations, several steps that are unique to the constellation approach can potentially help facilitate further calibration. The footprints of satellites in adjacent orbital planes have considerable simultaneous overlap at the poles, and the footprints of adjacent satellites in the same orbital plane are almost identical just minutes apart, providing the perfect opportunity to carry out frequent cross calibrations as suggested by Wiscombe and Chiu [2013]. In addition, the slowly evolving, bright, and largely homogeneous polar regions present an ideal natural calibration target. Absolute calibration, required to track and correct for long-term drift, would require a small subset of the constellation to have a highly accurate onboard calibration system or have absolute calibration transferred to the constellation from an external source, both of which are active research areas [Wielicki et al., 2008;Swartz et al., 2015Swartz et al., , 2016.

Number of Satellites
Since the number of satellites available in reality may differ from the 36-satellite baseline constellation, we also test how sensitive the recovery performance is to the number of satellites using two distinct configurations: satellite limited and orbit plane limited. The satellite-limited configuration is defined as having six Journal of Geophysical Research: Atmospheres 10.1002/2016JD025514 equally spaced orbit planes while varying the number of satellites in each plane, useful for when launch opportunities are plentiful but fewer satellites are desirable. Conversely, the orbit plane-limited configuration has six satellites in each plane while the number of orbit planes is varied, useful for when mass production of satellites is possible (i.e., CubeSats) but launch opportunities are restricted. Since these two different configurations directly lead to sampling differences, we performed recoveries every hour throughout a 24 h time period to ensure that we capture any potential worse-case scenarios within the diurnal cycle.
As the number of satellites increases and the Earth becomes better sampled, Figure 10 shows that the absolute global mean bias in the recovered irradiance decreases. The improvement in absolute bias from 36 satellites to 42 satellites is small in both satellite limited and orbit plane limited, and both in the SW and LW, because the recovered irradiance from 36 satellites already achieves a satisfactory SNR for a recovery to a 1000 km spatial resolution (see Figure 3). For fewer than 36 satellites, the density criterion of 500 samples per hour is no longer met globally; mathematical instability begins to influence the recovery process, contaminating the global mean bias and resulting in a sharp increase in this bias. The instability becomes more widespread as the number of satellites is further reduced below 30, increasing both the average bias and the RMSE. To overcome this problem with fewer satellites would require relaxation of the 1000 km spatial resolution requirement. Additionally, while the bias and RMSE are comparable in the SW and LW for 36 and 42 satellites, the SW bias generally becomes much larger for fewer satellites. This is due to the nature of the SW irradiance field evolving more quickly in space and time, and spanning a larger range of magnitudes, resulting in more pronounced instabilities when they do occur.
When comparing the global mean bias and RMSE between the satellite-limited and orbit plane-limited configurations, the orbit plane-limited configuration generally yields improved results. This improvement is most apparent for fewer satellites and suggests that if one needs to reduce the total number of satellites from the baseline constellation, limiting the total number of orbit planes would provide the most scientific value with the same number of satellites, recovering the global mean outgoing irradiance more accurately and better representing spatial patterns.

Conclusions and Summary
We investigated the capability of a new constellation concept to determine global EOR with sufficient temporal resolution and accuracy, crucial for understanding fundamental atmospheric processes and feedbacks  (Figures 10a and 10c) and an orbit plane-limited configuration (Figures 10b and 10d). Dots and error bars represent the means and the 25th and 75th percentiles, respectively, from 24 recovered irradiance fields.
Journal of Geophysical Research: Atmospheres 10.1002/2016JD025514 at both weather and climate scales. Capitalizing on technology revolutions in small satellites and sensor miniaturization, the proposed baseline constellation comprises 36 identical WFOV radiometers, evenly distributed in six non-Sun-synchronous orbit planes. The WFOV feature is desirable to provide a global coverage over a short time scale, allowing us to track fast-evolving synoptic phenomena such as cyclones and dust storms, for understanding their interactions with radiation and for model evaluations. To investigate the errors associated with the baseline constellation, and the error growth with spatial resolution and time scale, we developed a recovery method and performed extensive simulation experiments.
The baseline constellation provided sufficient 5 s instantaneous WFOV measurements for a stable recovery of hourly outgoing irradiance fields at both planetary and synoptic scales. At a spatial resolution of 1000 km, assuming isotropic radiation and perfect calibrations, the mean hourly global errors averaged for 1 day are respectively 0.16 AE 0.45 W m À2 and À0.13 AE 0.15 W m À2 for SW and LW. These results are unprecedented when compared with current observational products but are based on simplified radiance fields and idealized instrument characteristics. To identify the influence of these simplifications and changes in the constellation setup, we performed a series of sensitivity tests.
First, the influence of the anisotropic nature of the radiance field was investigated by generating a new set of simulated measurements using angular distribution models that take into account the directional effects of the solar-viewing geometry, the surface type, and the meteorology of the scene. We found that the global mean outgoing shortwave irradiance could be recovered to within 0.56 W m À2 of that in the control experiment by randomly accounting for the directional effects of the meteorology. This method significantly increases regional errors, but low regional errors can be obtained by assuming that the radiation is isotropic during the recovery process, which kept the RMSE to within 6% of the control experiment. It is also possible to reduce the error by incorporating angular distribution models directly into the core of the recovery method, but in that case, ancillary observations will be needed to help classify the scene and determine its associated angular distribution model. Second, we tested the impact of a systematic calibration error of 0.5 W m À2 , with and without random instrument noise up to 0.5 W m À2 . Overall, the recovery errors in hourly global means due to compounding sources of uncertainty grow linearly, and thus, one can adjust the total error estimate linearly if the actual calibration errors are smaller or larger than the tested value. A systematic calibration error of 0.5 W m À2 was found to have the largest individual impact of 0.63 W m À2 on the recovered SW and LW irradiance, while random instrument noises tend to have a negligible impact.
Third, the recovery errors in hourly global mean irradiance increase with decreasing number of satellites. Since irradiance fields in the SW tend to be more inhomogeneous and evolve faster than that in the LW, the recovery errors grow faster in the SW when the number of satellites is reduced. When a reduction in the number of satellites is necessary in reality, our results suggest that keeping the same number of satellites per orbital plane but reducing the total number of orbit planes would be more effective to retain the observation capability in measuring global mean outgoing irradiance.
Finally, we note that other constellation configurations beyond those that were tested in this study are possible, including hosted payloads on commercial satellites and a single dedicated launch with differing satellite drifts. Similarly, global mean irradiance fields using various configurations can be obtained using our recovery method, although we envisage that some modification or prior constraints may be necessary to ensure a stable recovery. Nevertheless, this study demonstrates the great potential of a new constellation concept for monitoring global Earth energy flows in the climate system at a time scale shorter than monthly, which is challenging to achieve using existing radiation budget measurements.