Deriving the Vertical Variations in the Diffuse Attenuation Coefﬁcient of Photosynthetically Available Radiation in the North Paciﬁc Ocean from Remote Sensing

: Diffuse attenuation coefﬁcient of photosynthetically available radiation (PAR), K PAR , is a key product of ocean color remote sensing. Current ocean color algorithms generally detect only the average K PAR within one optical depth, K RSPAR . Due to the marked vertical variations of K PAR , knowledge of K RSPAR is insufﬁcient to accurately evaluate the submarine light ﬁeld. By using ﬁeld in situ observations, a two-step approach, based on the development of an ocean color algorithm for K RSPAR and the relationships between K RSPAR and the average K PAR from the surface down to depth Z ( K ZPAR ), was developed to remotely estimate the vertical variations in K ZPAR in the North Paciﬁc from the MODerate-resolution Imaging Spectrometer (MODIS) imagery. The root mean square difference of log( K ZPAR ) in depths within the euphotic zone was around ± 0.059 (in unit of m − 1 for K ZPAR ), which corresponded to a deviation of ± 15% for the estimated K ZPAR and the penetration depths of PAR. Our study may provide a promising approach to detect the vertical variations of K ZPAR and underwater PAR distributions in the North Paciﬁc Ocean.


Introduction
Knowledge of the vertical distribution of photosynthetically available radiation (PAR), which is associated with approximately the visible part of the solar radiation spectrum from 400 to 700 nm, is essential to understanding the underwater biogeochemical processes and the response of marine ecosystem to environmental changes [1,2]. Playing a key role in driving water-column photosynthesis, the vertical distribution of PAR, especially in the euphotic layer, is needed for accurately quantifying and modeling the aquatic primary production, which, as a whole for the global oceans, accounts for about half of global carbon fixation and thus alters the global oxygen and carbon balances [3][4][5]. In addition to being a controlling factor in the water-column photosynthesis, the vertical distribution of PAR is fundamental to predictions of vertical subsurface heating rates, given that more than 98% of total incoming solar irradiance as short wave is converted to heat. Subsurface heating can lead to instabilities of the upper layer, which in turn affects the physical and K PAR , and thus K Z PAR , can be measured in situ relatively simply, but such data are severely limited in their spatial and temporal coverage. On the other hand, satellite remote sensing provides a synoptic, large-scale, and continuous view of the oceans to augment the direct field in situ observations [18]. Efforts were made in past decades to develop and validate ocean color algorithms for K PAR [10,12,[19][20][21]. Three types of algorithms are commonly developed for remotely estimating K PAR : chlorophyll-based approach [12,19,22], IOP-based approach [10,13], and empirical algorithm based on the band ratios of remote sensing reflectance (R rs ) [20]. The first two types require remote sensing products, e.g., chlorophyll a (Chl_a) concentration and IOPs of absorption and backscattering coefficients, respectively. An additional uncertainty is likely introduced when propagating these products to estimate K PAR . Furthermore, the vertical variations in Chl_a and IOPs, especially in stratified waters, may also cause difficulty in evaluating the vertical distributions in K PAR [10,22].
The third type, namely, the R rs -based empirical algorithm, has been routinely used to estimate the diffuse attenuation coefficients in the global oceans from space [20]. For historic reasons, algorithms for the diffuse attenuation coefficient of downwelling plane irradiance at 490 nm (K 490 ), rather than K PAR , are commonly developed [20,21]. For example, the remotely sensed K 490 (K RS 490 ) by the MODerate-resolution Imaging Spectrometer on Aqua (MODIS-Aqua) is estimated as log K RS 490 − 0.0166 = −0.8515 − 1.8263X + 1.8714X 2 − 2.4414X 3 − 1.0690X 4 (2) where X = log[R rs (488)/R rs (555)]. In order to convert K RS 490 to the remotely sensed K PAR (K RS PAR ), an empirical algorithm is commonly used such that [12] K RS PAR = 0.0864 + 0.884K RS Although the combination of Equations (2) and (3) provides an operational approach to estimate K RS PAR from the MODIS-Aqua imagery, two questions are raised when applying it to evaluate the vertical distributions in K PAR . First, this approach requires the estimation of K RS 490 , and an additional uncertainty may be introduced when propagating K RS 490 to estimate K RS PAR . Secondly, there exists a gap between K RS PAR and the vertical variations in K PAR  depth at which PAR is reduced to about 37% (e −1 ≈ 37%) of its surface value [13], since more than 90% of water-leaving radiance is reflected from the upper one optical depth [23]. K RS PAR is thus valid to evaluate the penetration of PAR in the layer above one optical depth only. Even with vertically uniform profiles of IOPs, K PAR can vary by a factor of two; thus, its value below one optical depth may be much lower than K RS PAR [10,13]. As such, the penetration depths of PAR in the lower layer may be substantially underestimation if using a constant K PAR equal to K RS PAR for the whole water column. Therefore, the relationship between K RS PAR and K Z PAR , in addition to the ocean color algorithms for K RS PAR , has to be elucidated in order to accurately quantify the submarine penetration of PAR from space.
Although previous works have reported the estimations of K Z PAR on few specific layers, e.g., from the surface down to the base of the euphotic zone [10,14,15,17] or two optical depths [12], an approach to evaluate the vertical profiling in K Z PAR continuously has seldom been documented [22].
In this study, we report an attempt to develop an operational approach for remotely estimating the vertical distributions of K Z PAR in the oceanic water by using field observations collected in the North Pacific ( Figure 1). The approach was designed in a two-step method: (1) to develop an algorithm for estimating K RS PAR directly from the R rs band ratio, and (2) to build the relationships between K RS PAR and K Z PAR . Section 2 describes the data and methods employed in this study. Section 3 shows the results and discusses them in terms of algorithm development and validation. Section 4 shows the application of the algorithms in evaluating the distributional patterns of K Z PAR in the North Pacific. Section 5 draws the conclusions.
Remote Sens. 2023, 15, x FOR PEER REVIEW 3 of 17 optical depth only. Even with vertically uniform profiles of IOPs, KPAR can vary by a factor of two; thus, its value below one optical depth may be much lower than [10,13]. As such, the penetration depths of PAR in the lower layer may be substantially underestimation if using a constant KPAR equal to for the whole water column. Therefore, the relationship between and , in addition to the ocean color algorithms for , has to be elucidated in order to accurately quantify the submarine penetration of PAR from space. Although previous works have reported the estimations of on few specific layers, e.g., from the surface down to the base of the euphotic zone [10,14,15,17] or two optical depths [12], an approach to evaluate the vertical profiling in continuously has seldom been documented [22].
In this study, we report an attempt to develop an operational approach for remotely estimating the vertical distributions of in the oceanic water by using field observations collected in the North Pacific ( Figure 1). The approach was designed in a two-step method: (1) to develop an algorithm for estimating directly from the Rrs band ratio, and (2) to build the relationships between and . Section 2 describes the data and methods employed in this study. Section 3 shows the results and discusses them in terms of algorithm development and validation. Section 4 shows the application of the algorithms in evaluating the distributional patterns of in the North Pacific. Section 5 draws the conclusions.

Study Area
The general distributional patterns of the biophysical properties, as shown in the climatologically seasonal distributions in spring in the MODIS-Aqua-derived sea surface temperature (SST) and Chl_a, in the North Pacific are shown in Figure 1. SST generally increases southwards from~4 • C at 50 • N to~29 • C at the equator, primarily reflecting the latitudinal dependence of heat loss/gain across the air-sea interface. The major circulations, namely, the Kuroshio, the North Pacific Current, the California Current, and the North Equatorial Current, define the approximate boundaries of the North Pacific Subtropical Gyre (NPSG), which is the most extensive gyre of the Earth. The SSTs at the north and south boundaries of the NPSG are about 15 • C and 28 • C, respectively (Figure 1a). Corresponding to the southward-increasing SST, which may indicate the strengthened vertical stratification and thus reduced nutrient supply for phytoplankton growth, Chl_a generally decreases southwards, varying from~0.4 mg m −3 at 50 • N to~0.1 mg m −3 at the equator (Figure 1b). High Chl_a, e.g., over 1 mg m −3 , can be found in the coastal waters, while the most oligotrophic waters are found in the NPSG, with Chl_a typically around 0.05 mg m −3 . In terms of PAR, as exampled in its climatologically seasonal distribution in spring in the North Pacific ( Figure A1 in Appendix A), it generally increases southwards from 30 E m −2 d −1 at 50 • N to~45 E m −2 d −1 at the equator. Latitudinal maximum values of PAR typically occur at about 10-17 • N. Seasonally, PAR generally increases from the winter to the summer, then following a decrease to the winter. Seasonal variations are typically more pronounced in the high-latitude region than in the low-latitude region. Such a distributional pattern may primarily reflect the spatial and temporal variability in the incoming total solar irradiance, in addition to modulation of the atmospheric transmission associated with cloud cover and aerosol distributions.  Table 1). At each station, vertical distributions of water temperature and salinity were recorded with a SeaBird SBE9/11 conductivity-temperature-depth (CTD) recorder. The mixed layer depth (MLD) was defined as the depth at which the temperature was 0.5 • C lower than the surface temperature [24]. Discrete water samples were collected with depth by using 20 L SBE bottles mounted onto a Rosette sampling assembly. Samples for the determination of the concentrations of Chl_a were collected through 47 mm GF/F Whatman glass fiber filters by filtering~2.5 L seawater onboard the ship. The concentrations of Chl_a were determined by fluorimetry [25], though such a historically widely used technique may produce a relatively high uncertainty in determining Chl_a, e.g., tending to underestimate Chl_a when chlorophyll b is unusually distributed in the water. Only the surface concentrations of Chl_a, whose samples were collected at 10 m depth, are reported here.

Field In Situ Observations in the Experiment of Pac_2017
At each optical station conducted in the daytime, vertical profiling of the spectral downwelling irradiance [E d (λ, Z)] and upwelling radiance [L u (λ, Z)], along with simultaneous measurements of the sea surface solar irradiance [E s (λ)], between 400 and 700 nm in 2 nm intervals was determined through a Satlantic HyperPro II system by following the manufacturer's manual (https://www.seabird.com/hyperpro-ii/product-downloads? id=54627923897 (accessed on 1 December 2021)). PAR(Z) (µE m −2 s −1 ) was calculated as follows: The above-water spectral remote sensing reflectance [R rs (λ)] was determined from E d (λ, Z), L u (λ, Z), and E s (λ) using a data processor in compliance with the Ocean Optics Protocols established for in situ observations in support of ocean color remote sensing measurements [26].   Table 1). The extracted data included PAR(Z), E d (λ, Z), L u (λ, Z), and E s (λ). Similarly, K Z PAR and R rs (λ) were determined from these observations by following the methods as used in Section 2.2. To reduce the effect of the variation of the incident solar light on accurately determining K Z PAR and R rs (λ), e.g., overpassing clouds might markedly reduce the incoming PAR(0 − ) and thus PAR(Z), casts whose coefficients of variation in E s (λ) during sampling were >1% were excluded. Then, 745 sets of observations in K Z PAR and R rs (λ) were available. The extracted data were pooled into our observations, and, finally, there were 770 sets of observations. About 80% (617) of stations, which were randomly selected by using Excel ® "Rand" function, were used for algorithm development, and the remaining 20% (153) of stations were used to validate the algorithms.

Statistics
The performance of the algorithm of the remotely sensed K Z PAR and the penetration depth of PAR was evaluated by the root mean square difference (RMSD) such that where n, F, and D are the number of samples, the observations, and the corresponding derived products, respectively. Unless noted, uncertainty used in this study means RMSD.

Hydrographic Properties Observed In Situ
The in situ observations followed similar patterns to the remotely sensed data (Figure 1), as exampled in the latitudinal variations of SST, surface salinity, MLD, EZD, and surface Chl_a along 143-149 • E observed during our cruise (Pac_2017) conducted in April-May 2017 ( Figure 2). According to reports in the literature [28][29][30][31][32][33], this area might be subdivided into three primary hydrographic regimes, the Kuroshio Extension (north of 34 • N), the NPSG (15-25 • N) and the western Tropical Pacific (south of 10 • N), and two transition zones (25)(26)(27)(28)(29)(30)(31)(32)(33)(34) • N and 10-15 • N). The SST, 14.2 ± 3.0, 27.0 ± 1.5, and 29.9 ± 0.4 • C in the Kuroshio Extension, the NPSG, and the western Tropical Pacific, respectively, increased southwards. Surface salinity generally increased from 34.6 ± 0.2 in the Kuroshio Extension to 34.9 ± 0.3 in the NPSG and then decreased to 34.0 ± 0.2 in the western Tropical Pacific, consistent with previous reports [34,35]. The MLD, 138 ± 72 m, was deep and was much deeper than the EZD, 53 ± 21 m, in the Kuroshio Extension, suggesting a replete supply of nutrients by vertical mixing and thus resulting in a large surface Chl_a, 1.12 ± 0.30 mg m −3 . In the NPSG, as reported in the literature [36], the MLD, 29 ± 10 m, became very shallow and was much shallower than the EZD, 115 ± 17 m. Because the permanent anticyclonic circulation leads to an Ekman downwelling, a deep pycnocline, and a limited supply of nutrients by vertical mixing [36,37], the surface Chl_a, 0.07 ± 0.02 mg m −3 , was very low in the NPSG. In the western Tropical Pacific, the MLD, 76 ± 20 m, was shallower than the EZD, 103 ± 16 m, and the surface Chl_a, 0.07 ± 0.02 mg m −3 , was also low. Nevertheless, the MLD, the EZD, and the surface Chl_a were all significantly (p < 0.01) correlated with the SST (r = −0.53, 0.88, and −0.84), suggesting that the temperature-associated status of water stratification is important, if not dominant, in controlling the distributions of the physical, hydrographic, and biological qualities in the North Pacific.
where n, F, and D are the number of samples, the observations, and the corresponding derived products, respectively. Unless noted, uncertainty used in this study means RMSD.

Hydrographic Properties Observed In Situ
The in situ observations followed similar patterns to the remotely sensed data ( Figure  1), as exampled in the latitudinal variations of SST, surface salinity, MLD, EZD, and sur face Chl_a along 143-149°E observed during our cruise (Pac_2017) conducted in April-May 2017 ( Figure 2). According to reports in the literature [28][29][30][31][32][33], this area might be sub divided into three primary hydrographic regimes, the Kuroshio Extension (north of 34°N) the NPSG (15-25°N) and the western Tropical Pacific (south of 10°N), and two transition zones (25- 34°N and 10-15°N). The SST, 14.2 ± 3.0, 27.0 ± 1.5, and 29.9 ± 0.4 °C in the Kuro shio Extension, the NPSG, and the western Tropical Pacific, respectively, increased south wards. Surface salinity generally increased from 34.6 ± 0.2 in the Kuroshio Extension to 34.9 ± 0.3 in the NPSG and then decreased to 34.0 ± 0.2 in the western Tropical Pacific consistent with previous reports [34,35]. The MLD, 138 ± 72 m, was deep and was much deeper than the EZD, 53 ± 21 m, in the Kuroshio Extension, suggesting a replete supply o nutrients by vertical mixing and thus resulting in a large surface Chl_a, 1.12 ± 0.30 mg m −3 In the NPSG, as reported in the literature [36], the MLD, 29 ± 10 m, became very shallow and was much shallower than the EZD, 115 ± 17 m. Because the permanent anticyclonic circulation leads to an Ekman downwelling, a deep pycnocline, and a limited supply o nutrients by vertical mixing [36,37], the surface Chl_a, 0.07 ± 0.02 mg m −3 , was very low in the NPSG. In the western Tropical Pacific, the MLD, 76 ± 20 m, was shallower than the EZD, 103 ± 16 m, and the surface Chl_a, 0.07 ± 0.02 mg m −3 , was also low. Nevertheless the MLD, the EZD, and the surface Chl_a were all significantly (p < 0.01) correlated with the SST (r = −0.53, 0.88, and −0.84), suggesting that the temperature-associated status o water stratification is important, if not dominant, in controlling the distributions of the physical, hydrographic, and biological qualities in the North Pacific.  actuality, a positive correlation (r = 0.77) was found between K RS PAR and the surface Chl_a, suggesting that this study area is well classed into the Case-1 waters in which the optical properties are determined primarily by the phytoplankton and their associated living and inanimate materials.
Depth-dependence in K Z PAR was pronounced in all water conditions. Vertical variations of PAR(Z)/PAR(0 − ) and K Z PAR , as exampled in a well-mixed station (149 • E, 38 • N) and a stratified station in the NPSG (143 • E, 22 • N) during the Pac_2017 cruise, are shown in Figure 3. More examples can be found in the Appendix A ( Figure A2). At both stations, although PAR(Z)/PAR(0 − ) followed a generally exponential decay function with depth, it was reduced faster in the upper layer than in the lower layer, resulting in a generally decreasing K Z PAR with depth, consistent with previous reports [10,12,13,17]. In the well- were 42 and 92 m, which were shallower than the observations of 49 and 117 m by 14% and 21% in these two stations ( Figure 3). These results indicated that knowing K RS PAR alone is insufficient to retrieve the submarine penetration of PAR. In order to accurately estimate the penetration depths of PAR at different light levels from space, the relationships between K Z PAR and K RS PAR have to be built.
Horizontal variations of KPAR, as exampled in the latitudinal variations of culated from the in situ PAR observations in the first optical depth during the P cruise in April-May 2017, are shown in Figure 2c. Ranging between 0.05 and 0 generally decreased southward, following a similar trend to the surface C actuality, a positive correlation (r = 0.77) was found between and the surfac suggesting that this study area is well classed into the Case-1 waters in which th properties are determined primarily by the phytoplankton and their associated liv inanimate materials.
Depth-dependence in was pronounced in all water conditions. Vertic tions of PAR(Z)/PAR(0 − ) and , as exampled in a well-mixed station (149°E, 38 a stratified station in the NPSG (143°E, 22°N) during the Pac_2017 cruise, are s Figure 3. More examples can be found in the Appendix (Figure A2). At both stati hough PAR(Z)/PAR(0 − ) followed a generally exponential decay function with d was reduced faster in the upper layer than in the lower layer, resulting in a gene creasing with depth, consistent with previous reports [10,12,13,17]. In t mixed station, , ranging between 0.115 and 0.094 m −1 , varied by a factor of 1.2 3a). In the stratified station, , ranging between 0.056 and 0.038 m −1 , varied by of 1.5 (Figure 3b). Without considering the vertical variations in , i.e., by vertically constant to , the penetration depths of PAR would be mark derestimated. For example, the estimated EZDs by using constant equal were 42 and 92 m, which were shallower than the observations of 49 and 117 m and 21% in these two stations ( Figure 3). These results indicated that knowing is insufficient to retrieve the submarine penetration of PAR. In order to accurately the penetration depths of PAR at different light levels from space, the relations tween and have to be built.

Algorithm Development for Remotely Sensing K Z PAR in the North Pacific
In this study, we adapted a two-step method for remotely estimating K Z PAR : developing an ocean color algorithm for K RS PAR , and building the relationships between K Z PAR and K RS PAR . For the first step, K RS PAR could be related to R rs band ratio such that (Figure 4) log K RS PAR = (−0.697 ± 0.006) − (0.951 ± 0.013)X n = 617, r 2 = 0.896 (6) In this study, we adapted a two-step method for remotely estimating ing an ocean color algorithm for , and building the relationships between . For the first step, could be related to Rrs band ratio such that (Figu log( ) = (−0.697 ± 0.006) − (0.951 ± 0.013) n = 617, r 2 = 0.896 Again, X = log[Rrs(488)/Rrs (555)]. The RMSD of the derived log( ) was unit of m −1 for ), which corresponded to a deviation of about ±16% in t . This uncertainty was reduced to 2/3 of that by applying the commonl proach, the combined Equations of (2) and (3), which introduced the RMSD o the derived log( ) and the corresponding deviation of about ±27% in . the use of the combined Equations (2) and (3) generally introduced an overest . is given in unit of m −1 . Black solid line-best fit line of Equation (6); da one RMSD (±0.065) from the best fit line; red solid line-relationship by using the com tions (2) and (3).
For the second step, our data indicated that for a given light level (f) wit photic zone whose corresponding depth is Zf so that f = PAR(Zf)/PAR(0 − ) × 100 erage KPAR from the surface down to Zf, f , could be linearly related to Figure 5), while the derived slope (A) was related to log(f) by a third-order p function such that (r 2 = 0.996; Figure 6 The RMSD of log( f ) derived from Equation (7) was around ±0.041 (r tween ±0.015 and ±0.055 for f varying from 1% to 70%; f in unit of m −1 ), r the corresponding deviation of the derived f to be around ±10%. By defin Again, X = log[R rs (488)/R rs (555)]. The RMSD of the derived log(K RS PAR ) was ±0.065 (in unit of m −1 for K RS PAR ), which corresponded to a deviation of about ±16% in the derived K RS PAR . This uncertainty was reduced to 2/3 of that by applying the commonly used approach, the combined Equations of (2) and (3), which introduced the RMSD of ±0.103 in the derived log(K RS PAR ) and the corresponding deviation of about ±27% in K RS PAR . Moreover, the use of the combined Equations (2) and (3) generally introduced an overestimation on K RS PAR . For the second step, our data indicated that for a given light level (f) within the euphotic zone whose corresponding depth is Z f so that f = PAR(Z f )/PAR(0 − ) × 100%, the average K PAR from the surface down to Z f , K Z f PAR , could be linearly related to K RS PAR (r 2 > 0.96; Figure 5), while the derived slope (A) was related to log(f) by a third-order polynomial function such that (r 2 = 0.996; Figure 6)  By combining Equations (6) and (7), given the light level of f, f could be determined from the Rrs band ratio, and the corresponding Zf could then be further determined from Equation (10). The RMSD of log( f ) derived from this approach was around ±0.059 ( f in unit of m −1 ), which corresponded to a deviation of ±15% in both the derived f ( Figure 7a) and Zf (Figure 7b). The above results suggest that, given substantially accurate Rrs spectra by satellites, the remote estimates in f and Zf at any light levels within the euphotic zone could be accurately estimated from space with a deviation of ±15% only.   By combining Equations (6) and (7), given the light level of f, f c mined from the Rrs band ratio, and the corresponding Zf could then be furth from Equation (10). The RMSD of log( f ) derived from this approach was ( f in unit of m −1 ), which corresponded to a deviation of ±15% in both the ( Figure 7a) and Zf (Figure 7b). The above results suggest that, given substan Rrs spectra by satellites, the remote estimates in f and Zf at any light le euphotic zone could be accurately estimated from space with a deviation o The RMSD of log(K Z f PAR ) derived from Equation (7) was around ±0.041 (ranging between ±0.015 and ±0.055 for f varying from 1% to 70%; K Z f PAR in unit of m −1 ), resulting in the corresponding deviation of the derived K Z f PAR to be around ±10%. By definition, Remote Sens. 2023, 15, 3023 10 of 17 Z f could thus be determined such that By combining Equations (6) and (7), given the light level of f, K Z f PAR could be determined from the R rs band ratio, and the corresponding Z f could then be further determined from Equation (10). The RMSD of log(K Z f PAR ) derived from this approach was around ±0.059 (K Z f PAR in unit of m −1 ), which corresponded to a deviation of ±15% in both the derived K Z f PAR (Figure 7a) and Z f (Figure 7b). The above results suggest that, given substantially accurate R rs spectra by satellites, the remote estimates in K Z f PAR and Z f at any light levels within the euphotic zone could be accurately estimated from space with a deviation of ±15% only. (m −1 ) and (b) the associated penetration depths (Zf, m) by using our approach with the field in situ observations for selected light levels of f = 50% (○), 10% (∆), and 1% (◊) for data used for algorithm development. Solid lines-1:1 plots; dashed lines-±15% deviated from the 1:1 plots.

Validation on the Derived f
The performance of the proposed approach was further validated by comparing f and Zf derived from Rrs(λ) to the observations in the dataset used for validation. The resulting RMSD in log( f ) for light levels within the euphotic zone was around ±0.057 ( f in unit of m −1 ), which corresponded to a deviation of ±14% in both the derived f ( Figure 8a) and Zf (Figure 8b). These values were even slightly smaller than those found in the algorithm development, suggesting that there was a general good agreement on the remotely sensed and the observed f and Zf. Comparisons of the MODIS-Aqua-derived f and Zf to the field in situ observations, by following the protocols [38], are shown in Figure 9. For each station, the main criteria included minimum 5 valid pixels in the 3 × 3 pixel array centered on the field station, coefficient of variation of valid satellite pixels lower than 0.15, and time difference within ±3 h between satellite overpass and field observation. The resulting RMSD in log( f ) for light levels within the euphotic zone was around ±0.073 ( f in unit of m −1 ), which corresponded to a deviation of ±18% in both the derived f ( Figure 9a) and Zf (Figure 9b). These values were slightly larger but still quite comparable to those found in

Validation on the Derived K Z f PAR
The performance of the proposed approach was further validated by comparing K Z f PAR and Z f derived from R rs (λ) to the observations in the dataset used for validation. The resulting RMSD in log(K Z f PAR ) for light levels within the euphotic zone was around ±0.057 (K Z f PAR in unit of m −1 ), which corresponded to a deviation of ±14% in both the derived K Z f PAR (Figure 8a) and Z f (Figure 8b). These values were even slightly smaller than those found in the algorithm development, suggesting that there was a general good agreement on the remotely sensed and the observed K PAR and Z f to the field in situ observations, by following the protocols [38], are shown in Figure 9. For each station, the main criteria included minimum 5 valid pixels in the 3 × 3 pixel array centered on the field station, coefficient of variation of valid satellite pixels lower than 0.15, and time difference within ±3 h between satellite overpass and field observation. The resulting RMSD in log(K  (Figure 9a) and Z f (Figure 9b). These values were slightly larger but still quite comparable to those found in Figures 7 and 8, suggesting that there was a reasonable agreement between the remotely sensed and the observed K Z f PAR and Z f . Nevertheless, since only 13 match-up data points were available for this evaluation, further validation when additional field observations became available would be desirable.

Application: Distributions of Remotely Sensed f and Zf in the North Pacific
The operational approach developed in this study can be applied to estimate the vertical variations of KPAR from space. Using MODIS-Aqua monthly Level-3 products (Reprocessing R2022.0) of Rrs(488) and Rrs(555) extracted from the NASA Ocean Color Web (http://oceancolor.gsfc.nasa.gov (accessed on 1 December 2022)), the climatologically (2002-2022) seasonal distributions in spring in the satellite-derived , 1% , and EZD  (Figure 10a). Such a distributional pattern generally follows the pattern of the latitudinal variations in SST (Figure 1a) such that high K RS PAR is associated with low SST and vice versa. Enhanced nutrient supply and, hence, primary production by vertical mixing to the surface layer in low SST (Figure 2) may account for high K RS PAR in low-SST conditions. High values of K RS PAR can also be found in the coastal waters, in which land-sourced inputs and/or southward currents are marked, and in the cyclonic circulation gyres, e.g., region surrounded by the Oyashio, the North Pacific Current, and the Alaska Current at about 160 • E-135 • W, 40-60 • N [32]. The upwelling may enhance vertical transport of nutrients to the surface layer to support phytoplankton growth, resulting in the increased turbidity and, thus, high K RS PAR . In contrast, the NPSG is characterized as oligotrophic waters with very low K RS PAR because of the occurrence of a permanent anticyclonic circulation gyre. K RS PAR is slightly enlarged south of the NPSG, possibly due to the slightly enhanced vertical mixing by the blowing of the trade winds [32]. K Z 1% PAR , whose value is about 81% of K RS PAR , follows a similar distributional pattern to K RS PAR (Figure 10b). The distributional pattern in the EZD follows a reverse trend to K Z 1% PAR (Figure 10c) according to Equation (9). As such, the EZD may vary from <50 m north of 40 • N to >110 m in the NPSG. the surface layer to support phytoplankton growth, resulting in the increased turbidity and, thus, high . In contrast, the NPSG is characterized as oligotrophic waters with very low because of the occurrence of a permanent anticyclonic circulation gyre. is slightly enlarged south of the NPSG, possibly due to the slightly enhanced vertical mixing by the blowing of the trade winds [32]. 1% , whose value is about 81% of , follows a similar distributional pattern to (Figure 10b). The distributional pattern in the EZD follows a reverse trend to 1% (Figure 10c) according to Equation (9). As such, the EZD may vary from <50 m north of 40°N to >110 m in the NPSG.
The application of the remotely sensed and the associated vertical distributions of the PAR may improve the accuracy in the assessment of the carbon fixation in the North Pacific. For example, the remotely sensed net primary production in the euphotic zone (PPeu) is thought to be linearly related to the euphotic zonal depth (EZD) [5]. The EZD is typically underestimated, e.g., by 14 ± 16% using the models integrated in the PPeu model [5] or by 16 ± 11% estimated by treating vertically constant to ( Figures  3 and A2). As a result, PPeu may be substantially underestimated in the North Pacific. A reassessment on the contribution of the North Pacific to the global carbon cycle may be needed.
Although an improved approach to remotely estimate in the North Pacific has been proposed in this study, it is necessary to note that it is designed to be applied to the Case-1 waters, given that most of the in situ observations were conducted in the offshore regions. For inshore locations and marginal seas at land-sea boundaries, which are typically characterized as the Case-2 waters, regionally tuned algorithms are generally needed. The application of the remotely sensed K Z PAR and the associated vertical distributions of the PAR may improve the accuracy in the assessment of the carbon fixation in the North Pacific. For example, the remotely sensed net primary production in the euphotic zone (PP eu ) is thought to be linearly related to the euphotic zonal depth (EZD) [5]. The EZD is typically underestimated, e.g., by 14 ± 16% using the models integrated in the PP eu model [5] or by 16 ± 11% estimated by treating K Z PAR vertically constant to K RS PAR (Figures 3 and A2). As a result, PP eu may be substantially underestimated in the North Pacific. A reassessment on the contribution of the North Pacific to the global carbon cycle may be needed.
Although an improved approach to remotely estimate K Z PAR in the North Pacific has been proposed in this study, it is necessary to note that it is designed to be applied to the Case-1 waters, given that most of the in situ observations were conducted in the offshore regions. For inshore locations and marginal seas at land-sea boundaries, which are typically characterized as the Case-2 waters, regionally tuned algorithms are generally needed.

Conclusions
This study demonstrated that, by adapting a two-step approach, that is, developing an ocean color algorithm for K RS PAR and building the relationships between K Z PAR and K RS PAR , the vertical variations in K Z PAR in the North Pacific may be assessed from space. The deviation in the remotely sensed K Z PAR in depths within the euphotic zone, as well as the corresponding penetration depths of PAR, is about ±15%. The significant underestimation in the remotely sensed euphotic zonal depth from methods adapted in current primary production models suggests a substantial underestimation of the primary production in the North Pacific. The variations in K Z PAR , together with those in other hydrographic properties (e.g., SST and Chl_a), suggest that our dataset covers various water conditions and, thus, may be a good representative to the Case-1 waters in the North Pacific Ocean. A similar approach may be developed for the global oceans when field observations in the global oceans become available.